理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

分散共分散行列 相関のあるサンプル作成

データセットから分散共分散行列を求めてみます

x1 <- c(151, 164, 146, 158)
x2 <- c(48, 53, 45, 61)
x3 <- c(8, 11, 8, 9)
data <- data.frame(x1,x2,x3)

#分散共分散行列
var(data) 
         x1        x2        x3
x1 62.25000 38.250000 10.333333
x2 38.25000 48.916667  4.333333
x3 10.33333  4.333333  2.000000

#相関行列
cor(data) 
          x1        x2        x3
x1 1.0000000 0.6931597 0.9260955
x2 0.6931597 1.0000000 0.4381055
x3 0.9260955 0.4381055 1.0000000

定義に従って相関係数を求めてみます. 例)cor(data)[2]を算出してみます

cor(data)[2]
var(data)[2]/sqrt(var(data)[1]*var(data)[5])

シミュレーション

平均や分散を指定した乱数を作成するのは簡単です.例えば正規分布であればrnorm(10,5,3)で、平均5分散9の乱数が10個生成されます.
しかし繰り返しデータの場合に、それぞれの回数の相関を設定した乱数を作成するのはちょっと難しくなります.
5回繰り返しデータ(X1~X5)を想定します.また、それぞれの平均値が(0,3,5,-2,-1)、分散8、共分散3となるように乱数を発生させます.

#まず、次のような共分散行列を作成します
   x1 x2 x3 x4 x5
x1  8  3  3  3  3
x2  3  8  3  3  3
x3  3  3  8  3  3
x4  3  3  3  8  3
x5  3  3  3  3  8
  
sigma <- matrix(rep(6,25),ncol = 5)   # 3を25個作成
sigma  

     [,1] [,2] [,3] [,4] [,5]
[1,]    6    6    6    6    6
[2,]    6    6    6    6    6
[3,]    6    6    6    6    6
[4,]    6    6    6    6    6
[5,]    6    6    6    6    6

#次に分散8を挿入します
diag(sigma)=rep(8,5)       #斜めに挿入
sigma

     [,1] [,2] [,3] [,4] [,5]
[1,]    8    6    6    6    6
[2,]    6    8    6    6    6
[3,]    6    6    8    6    6
[4,]    6    6    6    8    6
[5,]    6    6    6    6    8

#X1~X5の平均値を設定します
m <- c(0, 3, 5, -2, -1) 

#平均m, 分散8, 共分散3になるようなX1~X5を10組の乱数を発生させます.
y <- rmvnorm(10, m, sigma, method = "chol")      #乱数
colnames(y) <- name
#イメージのために書き出しておきます.乱数なので算出毎に数値は変化します.
 y
               x1         x2         x3         x4         x5
 [1,]  3.52189310  6.3816131  8.5906099  0.3448154  2.4660227
 [2,]  0.92752113  4.9547725  7.9595828  0.3717669  1.3875238
 [3,] -5.19449966 -0.3019742 -1.4509597 -4.3760477 -5.0233539
 [4,] -1.01591408  4.2919216  5.4850598 -1.0672082 -1.4707111
 [5,] -3.80629633 -0.6012603  4.7040327 -6.4885712 -3.7804199
 [6,] -3.10433704  5.7128433  2.2531806 -1.1011307  0.8055750
 [7,]  2.05904069  4.8775650  6.2192196  0.3335801  0.4234229
 [8,] -0.62623689  4.2532797  3.7670864 -1.5564523 -1.8985524
 [9,] -0.09027824  3.3205490  5.0517354 -1.2558948 -1.0637946
[10,] -1.30173038 -3.0114726  0.8776159 -8.0027438 -3.9780232

#各列に適当に色をつけて、プロットしてくれるmatplotを使って、作図します.
#デフォルトでは col の順番なのでx1=黒 ,x2=赤, x3=緑, x4=青, x5=水色
matplot(y,type="l") 

f:id:yoshida931:20180423183232p:plain:w500

互いに相関のある群(相関係数0.75)を作成

#分散共分散行列よりそれぞれの相関係数は、6/√8*√8=0.75 が読み取れます
   x1 x2 x3 x4 x5
x1  8  3  3  3  3
x2  3  8  3  3  3
x3  3  3  8  3  3
x4  3  3  3  8  3
x5  3  3  3  3  8

しかしグラフからは、X1~X5の相関が理解できません.
ここからは、上記yをコピーしたy2を使用して進めていきます.
相関行列を確認してみたら6/√8*√8=0.75に近い相関になっています.

cor(y2)
          x1        x2        x3        x4        x5
x1 1.0000000 0.6014473 0.8170968 0.6413491 0.7693525
x2 0.6014473 1.0000000 0.6684236 0.9623296 0.8923660
x3 0.8170968 0.6684236 1.0000000 0.6635989 0.7835190
x4 0.6413491 0.9623296 0.6635989 1.0000000 0.8511389
x5 0.7693525 0.8923660 0.7835190 0.8511389 1.0000000

nのサイズを大きくすれば、相関係数はさらに6/√8*√8=0.75に近づきます

y3 <- rmvnorm(50, m, sigma, method = "chol")
cor(y3)
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.8068844 0.7566464 0.7541858 0.7390697
[2,] 0.8068844 1.0000000 0.7345947 0.6725027 0.6586652
[3,] 0.7566464 0.7345947 1.0000000 0.7810799 0.7455624
[4,] 0.7541858 0.6725027 0.7810799 1.0000000 0.7226452
[5,] 0.7390697 0.6586652 0.7455624 0.7226452 1.0000000

y3をグラフにしてみます

matplot(y3,type="l") 

f:id:yoshida931:20180423223431p:plain:w700

x軸をX1~X5に移行することで相関が可視化できます

matplot(t(y3),type="l")    

f:id:yoshida931:20180423223335p:plain:w700
x1~x5まで、ほぼ同じ範囲に入っており相関の強さが視覚的に確認できました

復習のためにy3を残しておきます

             [,1]         [,2]       [,3]       [,4]       [,5]
 [1,]  0.10151443  2.499704243  5.9057675 -3.4943595 -1.3900978
 [2,]  1.50663613  5.280020390  6.0871375 -1.6427891  4.1792270
 [3,]  5.15807547  9.240978804 10.1524971  1.6396755  4.2384627
 [4,]  0.52636740  7.362142992  8.2131410 -2.7272272 -0.1657287
 [5,] -1.82768487  2.078275786  4.8978490 -1.3363880 -2.6741024
 [6,]  4.43524347  8.060044313  8.8144178  2.4425139  2.0623879
 [7,]  0.37328244  5.317456919  3.7818383 -1.0436226 -0.7600355
 [8,] -4.91121977 -0.627691232  3.4766397 -4.2966815 -1.9693309
 [9,]  3.08701302  4.490715808  8.0062255  1.3102477  3.4494656
[10,] -1.07464046  4.124389139  4.1597881 -3.2331504 -1.5478108
[11,] -4.98187707 -2.960175387 -1.1852294 -8.4042762 -6.7175592
[12,]  0.55455335  4.811030651  4.4954359 -2.1924783 -1.6559525
[13,] -3.25343072  0.216598109  2.5139412 -4.0381462 -2.4717920
[14,]  0.71190379  0.664704194  7.4396890  1.2862593 -0.7421359
[15,]  2.05860252  2.432888164  3.5347403 -1.8087834 -1.5556714
[16,] -4.35657559  0.316841111  1.6688628 -5.3447686 -4.5653558
[17,]  0.52655249  0.326658121  3.3274462 -5.0033717 -1.4797982
[18,] -0.51997740  1.522057422  3.9754390 -2.2568758 -3.4186446
[19,] -2.48896458  0.609983346  4.4815961 -4.4651604  0.5507060
[20,] -1.27166463  0.860673117  4.5316817 -1.6695152 -1.5116476
[21,]  2.36579566  8.115692097 10.1842162  1.5469158  1.2745295
[22,] -2.61885567  2.264994539  1.8073273 -4.6059818 -5.2668700
[23,] -2.87195778 -0.354779514  4.4273497 -2.1824915 -2.3847656
[24,] -0.31689104  4.490830114  7.2127744  0.8737950  3.2108026
[25,]  0.86603870  3.574589801  4.4276743 -1.6189355 -0.6536883
[26,] -1.08641447  1.650209214  0.2495224 -3.6113148 -3.8762204
[27,] -1.06718456  1.333869467  2.7389823 -6.0151694 -1.0285024
[28,]  2.92878284  3.572365454  7.5983883  4.4834977  1.1370509
[29,]  1.25628266  3.077046228  4.4049753 -3.9389680  0.8752170
[30,] -2.25266309  1.031792708  2.0312006 -5.2208742 -1.1417857
[31,] -1.83386252  1.464612026  5.4984764 -2.5442451 -2.5926546
[32,] -3.74237017  0.948552374  1.6345401 -5.3518601 -4.3701054
[33,]  1.57641371  2.949825751  7.6793848  2.8385187  2.5930688
[34,] -1.70707918  4.328429267  4.8674954 -2.3957196 -0.6474211
[35,] -0.68729957  0.186676655  6.3561554 -4.9844379 -2.6635800
[36,] -5.70176744  0.002329385  3.3561491 -3.6218192 -2.5791208
[37,]  0.73014712  2.200125852  3.0212866 -1.4821494 -3.4750934
[38,]  2.51743524  4.256512695  3.2907951  1.0304869  3.3052088
[39,] -3.71975203 -0.761618136  3.5496253 -2.0532424 -1.8430078
[40,]  4.00292393  6.456935876  9.6564018  1.9699281  0.9976009
[41,] -0.76466675  5.343766540  4.4333739 -0.7760534 -1.2543057
[42,] -5.23139631 -2.144231142  0.3253171 -7.3690128 -5.2907999
[43,]  0.22047353  0.598735545  5.6964360 -3.8599177  0.1664410
[44,]  0.25324554  4.186434455  3.3374806 -2.2359334 -2.5759674
[45,] -0.33142225  3.259191383  4.9409355 -4.1641440 -3.3400309
[46,]  4.28859451  7.946808493  9.4197735  1.7908695  1.2385851
[47,]  0.91987986  4.440995369  4.2562333 -4.3308892 -1.1173378
[48,] -3.44169133  0.376280288  2.3650166 -3.1425327 -3.2670228
[49,]  0.07415149  2.969230453  5.1548999 -1.9221235  2.5954580
[50,]  2.51582740  4.186680476  6.5627895  1.7226084  3.3340985

正規分布の重ね描き

text関数でグラフに文字の挿入

curve(dnorm(x, -2, 4), from=-10, to=10, ylim=c(0,0.4),ylab ="")
text(-5, 0.1, "N(-2,4)")
par(new=T)
curve(dnorm(x,3, 1), from=-10, to=10, ylim=c(0,0.4),ylab ="")
text(1.5, 0.3, "N(3,1)")
par(new=T)
curve(dnorm(x, 3, 4), from=-10, to=10, ylim=c(0,0.4),ylab ="")
text(6, 0.1, "N(3,4)")
dev.off()

f:id:yoshida931:20180419170138p:plain:w600

一部分を塗りつぶし

curve(dnorm(x,3, 1), from=-10, to=10, ylim=c(0,0.4),ylab ="")
# X4~10を100個の多角形(台形)に等分割
xvals <- seq(4, 10, length=100)  
# 対応するグラフの高さ
dvals <- dnorm(xvals,3, 1)      
# 塗りつぶす
polygon(c(xvals,rev(xvals)),c(rep(0,100),rev(dvals)),col="2")     
#x軸に4を追加
axis(side=1,at=c(4))
dev.off()

f:id:yoshida931:20180419171823p:plain:w400

データ取り込みと保存

クリップボードから一覧表の形式に取り込み

エクセルなどの一覧表から必要な部分をコピーします.

f:id:yoshida931:20180408100451p:plain:w300

その後、いかのような操作でRに取り込むことができます.

#x のなかに一覧表として取り込みます  
x <- read.table("clipboard",header = T)
#xを確認すると
  ID    A    B    C
1 x1 10.5  8.9  9.9
2 x2 11.2  9.9 10.6
3 x3  9.5 10.2 11.5
4 x4 10.6  8.5 12.3
5 x5 12.3 13.2 11.3

#x2のなかに1列目のIDを行名として取り込む場合
x2 <- read.table("clipboard",header = T, row.names=1)
#x2を確認するとIDというラベルが消えて、IDは行名として扱われます
      A    B    C
x1 10.5  8.9  9.9
x2 11.2  9.9 10.6
x3  9.5 10.2 11.5
x4 10.6  8.5 12.3
x5 12.3 13.2 11.3

列の取り出し

x2のなかのAとBをそれぞれ取り出します

x2$A
[1] 10.5 11.2  9.5 10.6 12.3
x2$B
[1]  8.9  9.9 10.2  8.5 13.2

取り出したAとBをt検定してみます

t.test(x2$A,x2$B)

    Welch Two Sample t-test

data:  x2$A and x2$B
t = 0.71918, df = 6.2608,
p-value = 0.498
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.610436  2.970436
sample estimates:
mean of x mean of y 
    10.82     10.14 

cvsファイルからの取り込み

あらかじめ作業ディレクトリにファイルを保管しておく必要があります.

y <- read.csv("sample01.csv",row.names=1)

データの保存

取り込んだyの一覧表を、テキストとして作業ディレクトリに保存してみます

write.table(y,file = "data.txt")
#dataというテキストファイルが作成されます

ベータ分布

完全独習 ベイズ統計学入門

完全独習 ベイズ統計学入門


この本を参考にベータ分布を勉強します.
ベータ分布:ベータ関数により導かれる分布. ベイズ推定の際に事前分布として利用します.

f(x)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} \ \ \  \ \ (0<x<1)

ベータ関数\ \ \ \ B(\alpha,\beta)=\int _0 ^1 x^{\alpha-1}(1-x)^{\beta-1} dt

期待値 \ \ \ \frac{\alpha}{\alpha+\beta}

分散\ \ \ \ \frac{\alpha \beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}

グラフを描いてみます
ベータを1に固定してαのみ大きくすると、1より大きくなると指数関数となります.

curve(dbeta(x,1,1),from = 0,to=1,ylim=c(0,5))     #黒
par(new=T)
curve(dbeta(x,2,1),from = 0,to=1,ylim=c(0,5),col=2)    #赤
par(new=T)
curve(dbeta(x,3,1),from = 0,to=1,ylim=c(0,5),col=3)    #緑
par(new=T)
curve(dbeta(x,4,1),from = 0,to=1,ylim=c(0,5),col=4)    #青
f:id:yoshida931:20180326140755p:plain:w400

αとβをランダムに

curve(dbeta(x,1,1),from = 0,to=1,ylim=c(0,3))    #黒
par(new=T)
curve(dbeta(x,2,2),from = 0,to=1,ylim=c(0,3),col=2)   #赤
par(new=T)
curve(dbeta(x,3,3),from = 0,to=1,ylim=c(0,3),col=3)    #緑
par(new=T)
curve(dbeta(x,5,3),from = 0,to=1,ylim=c(0,3),col=4)    #青
f:id:yoshida931:20180326141950p:plain:w400

ベータ分布の使用例1

以下は自作(架空)の問題です.
問)ある疾患の患者群で10中で3人が転倒した場合の、転倒確率をベータ分布によるベイズ統計で推測せよ.
10人中3人なので\frac{3}{10}といのが直観的な確率になるのですが、ベイズ推定では次のように考えることができます. まず転倒する確率をx、転倒しない確率を1-xとします.転倒する確率(x)+転倒しない確率(1-x)=1.xは連続無限に存在するので、確率密度となります(0≦x≦1).ここで、xの事前分布設定のために、xがどの値でも「同様に確からしい」と仮定します.つまりxの事前分布を下図のような一様分布と仮定します.

y=1\ \ \ \ (0 \leqq x \leqq1)

f:id:yoshida931:20180326145458p:plain:w300

「転倒した患者3名、転倒しなかった患者7名」という結果になる確率はP=x^{3} (1-x)^{7} となり、この数式は\alpha=4, \beta=8のベータ分布です.
したがって転倒確率は \ \ \ \frac{\alpha}{\alpha+\beta}=\frac{4}{12}=\frac{1}{3}となります.

ベータ分布の使用例2

A君は身体に障害をもっており、ただいま装具装着の練習中です.30秒以内で装着できる成功確率をxとして考えてみます. 今日は8回練習して、5回成功しました.A君の装具装着の成功率はどれくらいでしょうか? ベータ分布の使用例1と同じように、xの事前分布設定のために、xがどの値でも「同様に確からしい」と考えて一様分布を仮定します.
y=1\ \ \ \ (0 \leqq x \leqq1)
確率は \ \ \ \frac{\alpha}{\alpha+\beta}=\frac{6}{10}=\frac{3}{5}となります.現在のA君の装具装着確率は以下のようになります.

curve(dbeta(x,6,4),from = 0,to=1,ylim=c(0,3))
f:id:yoshida931:20180329150808p:plain:w450

事前分布が同様に確からしいということもありえないので、A君の今日の成功確率を事前分布として確率分布を推定してみます.ここでは正規化定数となるベータ関数を無視して確率分布のみを求めてみます.
P( 現状の成功確率 , 成功確率 ) = P(x) × 事前分布 = x (x^{5}) (1-x)^{3}=x^{6}(1-x)^{3}
\alpha=7, \beta=4となり以下のような確率分布(赤)となります.

curve(dbeta(x,6,4),from = 0,to=1,ylim=c(0,3))
par(new=T)
curve(dbeta(x,7,4),from = 0,to=1,ylim=c(0,3),col=2)
f:id:yoshida931:20180329154806p:plain:w450

感想)理学療法士が適切な事前分布を設定して、子どもたちの可能性を探るって重要なことだと思いました.

カッパ係数

個人的によく利用させていただいております以下のHPをもとに、今回はカッパ係数について少し勉強してみます
統計学入門−第5章
まずはHPに掲載してある次のサンプルデータを使用して、Rを使って処理してみます

分類数が2つの場合

rater01<-c(rep(1,40),rep(2,40),rep(1,10),rep(2,10)) 
rater02<-c(rep(1,40),rep(2,40),rep(2,10),rep(1,10))
tk <- table(rater01,rater02)
addmargins(tk)

       rater02
rater01   1   2 Sum
    1    40  10  50
    2    10  40  50
    Sum  50  50 100

次にRの関数を使用して求めてみます

library(irr) # irrの読み込み
dk <- cbind(rater01,rater02)
kappa2(dk)

 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 100 
   Raters = 2 
    Kappa = 0.6 

        z = 6 
  p-value = 1.97e-09 

z値から95%信頼区間を算出してみます

n <- length(rater01)
po <- (40+40)/100
pc <- (50/100)*(50/100)+(50/100)*(50/100)
vk <- (po*(1-po))/(n*(1-pc)^2)

as.numeric(kappa2(dk)[5])-1.96*sqrt(vk) 
as.numeric(kappa2(dk)[5])+1.96*sqrt(vk) 

0.44320.7568

参考HPと比較してください

分類数が3つの場合

rater02<-c(rep(1,19),rep(2,26),rep(3,4),rep(1,17),rep(1,7),rep(2,7),rep(2,5),rep(3,3),rep(3,12)) 
rater03<-c(rep(1,19),rep(2,26),rep(3,4),rep(2,17),rep(3,7),rep(1,7),rep(3,5),rep(1,3),rep(2,12))
tk2 <- table(rater02,rater03)
addmargins(tk2)

       rater03
rater02   1   2   3 Sum
    1    19  17   7  43
    2     7  26   5  38
    3     3  12   4  19
    Sum  29  55  16 100

d2 <-cbind(rater02, rater03)
kappa2(d2)

 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 100 
   Raters = 2 
    Kappa = 0.198 

        z = 2.8 
  p-value = 0.00508 

#95%信頼区間
n2 <- length(rater03)
po2 <- (19+26+4)/100
pc2 <- (43/100)*(29/100)+(38/100)*(55/100)+(19/100)*(16/100)
vk2 <- (po2*(1-po2))/(n2*(1-pc2)^2)

as.numeric(kappa2(d2)[5])-1.96*sqrt(vk2) 
as.numeric(kappa2(d2)[5])+1.96*sqrt(vk2)  

0.043905650.3520686

重み付けカッパ係数

上記の例で、解答1、2、3を順序尺度として考えているので、rankA=一致とした場合に、次のような順序が考えられます
rankA > rankB > rankC
1と3のペアよりも1と2のペアや2と3のペアの方が一致に近いという考え方です

     [,1]    [,2]    [,3]   
[1,] "rankA" "rankB" "rankC"
[2,] "rankB" "rankA" "rankB"
[3,] "rankC" "rankB" "rankA"

この重み付けでカッパ係数をRで算出してみます

kappa2(d2, "squared")   #自乗距離に応じた重み付け

 Cohen's Kappa for 2 Raters (Weights: squared)

 Subjects = 100 
   Raters = 2 
    Kappa = 0.196 

        z = 2 
  p-value = 0.0453 


kappa2(d2, "equal")   #均等に重み付け

 Cohen's Kappa for 2 Raters (Weights: equal)   #評価者間の不一致を均等に重み付け

 Subjects = 100 
   Raters = 2 
    Kappa = 0.197 

        z = 2.67 
  p-value = 0.0076 

データフレームからの抽出 2

準備

下のデータをコピーして、Rでフレームにします

実験A   10   6  10   9  10
実験B   10   5   5  12   4
実験C    5   4  11   4   6
実験D    9   5   2   3   1

コピーして、データフレームに取り込み

(x <- read.table("clipboard",row.names = 1))

      V2 V3 V4 V5 V6
実験A 10  6 10  9 10
実験B 10  5  5 12  4
実験C  5  4 11  4  6
実験D  9  5  2  3  1

列名を記載

cn <- LETTERS[1:5]
colnames(x) <- cn
x
       A B  C  D  E
実験A 10 6 10  9 10
実験B 10 5  5 12  4
実験C  5 4 11  4  6
実験D  9 5  2  3  1

データフレームをベクトルに

準備したフレームxから2列目のみを取り出してみます

#列の取り出しは簡単です
x[,2]
[1] 6 5 4 5

準備したフレームxから2行目のみを取り出してベクトルにします

#2行目
x[2,]
       A B C  D E
実験B 10 5 5 12 4

#2行目をベクトルに変換
(xv <- as.integer(x[2,]))
[1] 10  5  5 12  4

マトリックスにしてベクトルとして一括抽出

(xm <- as.matrix(x))
       A B  C  D  E
実験A 10 6 10  9 10
実験B 10 5  5 12  4
実験C  5 4 11  4  6
実験D  9 5  2  3  1

(as.vector(xm))
 [1] 10 10  5  9  6  5  4  5 10  5 11  2  9 12  4  3 10  4  6  1

因子ベクトル(カテゴリカルデータ) 

実験A 10  6 10  9 10  a
実験B 10  5  5 12  4  a
実験C  5  4 11  4  6  b
実験D  9  5  2  3  1  b

コピーしてRのフレームに取り込みます

(xc <- read.table("clipboard",row.names = 1))
      V2 V3 V4 V5 V6 V7
実験A 10  6 10  9 10  a
実験B 10  5  5 12  4  a
実験C  5  4 11  4  6  b
実験D  9  5  2  3  1  b

cn2 <- LETTERS[1:6]
colnames(xc) <- cn2
xc
       A B  C  D  E F
実験A 10 6 10  9 10 a
実験B 10 5  5 12  4 a
実験C  5 4 11  4  6 b
実験D  9 5  2  3  1 b

上記のような場合には文字が因子に自動的に変換されています

#5列目は整数、6列目は因子
class(xc[,5])
[1] "integer"
class(xc[,6])
[1] "factor"

A列、B列のF列bのみ取り出す

subset(xc,select = c(A,B),subset = F=="b")
      A B
実験C 5 4
実験D 9 5

A列、B列 の F列a VS b
matrixからベクトルに変換してa VS bの配列にします

(xxx <- c(as.vector(as.matrix(xca)),as.vector(as.matrix(xcb))))
xxx
[1] 10 10  6  5  5  9  4  5

(dat <- array(xxx, dim = c(2,2,2)))
, , 1

     [,1] [,2]
[1,]   10    6
[2,]   10    5

, , 2

     [,1] [,2]
[1,]    5    4
[2,]    9    5

項目名などを記入

name <- list("実験"=c("実験A","実験B"),"結果"=c("A","B"),F=c("a","b"))
dimnames(dat)<-name
dat

, , F = a

       結果
実験     A B
  実験A 10 6
  実験B 10 5

, , F = b

       結果
実験    A B
  実験A 5 4
  実験B 9 5

分散分析の基本

最終更新日:2018.3.5 まだ理解できていない.なので書き直し・・・ 一元配置分散分析

言葉の整理

要因(factor), 因子(factor):実験結果に影響を与える要素.それぞれの分野で使い分ける場合もあるので注意.このブログでは要因と因子の区別をせず「要因」で統一します.
〇元配置分散分析:〇は要因の数
水準:要因を分ける条件( 要因に含まれる項目 )
:=水準
平方和:偏差の二乗の合計
変動:=平方和
水準内変動:= 群内変動 = 誤差変動 = 残差変動
平均平方和:平方和を自由度で割った値(分散)
効果:目的変数への影響
主効果:水準によって効果が異なる場合、その要因に主効果があるという
交互作用:2つの要因が目的変数に対して互いに影響を及ぼしている場合
   相殺効果:一方の要因の効果が高いと他方の要因の効果がそれをうち消すこと
   相乗効果:2つの要因により効果が大きくなること
要因効果:主効果と交互作用の総称

分散分析とは

  • 分散分析は、3群以上のそれぞれの群の母平均が等しいという帰無仮説のもとに検定を行う分析方法です.t検定を繰り返して使用することはできません.例えば3群の比較を5%有意水準でそれぞれt検定を行った場合、帰無仮説は「3通り全ての検定において差がみられない」となり、対立仮説は「少なくとも一組には差がみられる」となります.つまり、有意水準は0.05ではなく1- 0.953 = 0.142 で検定したことになります.つまり、第一種の過誤が生じる危険性が14.2%となります.

  • 分散分析では、各群の母平均の推定を行うのではなく、各群が及ぼす効果(要因の効果)の分散を推定します.そのために群間の平均の変動、各群の郡内変動から平均平方和 (分散) を求め、要因の効果の有無をF分布から推測します.


一元配置分散分析(対応なし)

対応のない一要因で分類される多群の検定(要因分散分析).
実験A,B,Cのそれぞれにn人が含まれており、全部で(3×n)人の被験者からなる実験です.
群間に差があるかどうかを知りたい

        1   2   3  ・・・  n
実験A  x11 x12 x12 ・・・  x1n
実験B  x21 x22 x23 ・・・  x2n
実験C  x31 x32 x33 ・・・  x3n
実験D  x41 x42 x43 ・・・  x4n  
  
サンプル(samp)  繰り返し数 5
実験A   10   6  10   9  10
実験B   10   5   5  12   4
実験C    5   4  11   4   6
実験D    9   5   2   3   1

上記サンプル(samp)を使用して考えていきます.
実験結果のバラつきは分散として得られ、分散の原因として次の2点が考えられます.

  • 原因1)実験A、B、Cの効果によるもの(要因の効果
  • 原因2)被験者の問題、測定の誤差、偶然などによるもの(誤差項

重要 データの構造を表すモデル
誤差項(測定値 - 各群の平均 )は互いに独立で、正規分布に従うことを仮定しています.

f:id:yoshida931:20180305174523p:plain:w400

測定値 = 全平均 + (各水準の平均 - 総平均) + (測定値 - 各水準の平均 )より

全変動 = 水準間変動 + 水準内変動  
∑∑(測定値-総平均)^2 = n*∑(各水準の平均-総平均)^2+∑∑(測定値-各水準の平均 )^2 
f:id:yoshida931:20180305180648p:plain:w450

  ポイント:∑∑2(総平均-各群の平均)(測定値-各群の平均) = 0

分散分析の帰無仮説
実験A,B,C,Dの母平均は等しい.実験A,B,C,Dの母分散は0である.
したがいまして群間の差に着目して検定することになります

F分布
全ての群の母平均が等しい場合、統計量FはF分布に従います.

F=\frac{群間平方和}{群間自由度}\div\frac{群内平方和}{群内自由度}

サンプル(samp)の場合、群間平方和の自由度は3、群内平方和の自由度は5*4-4=16となります.したがって自由度(3,16)のF分布に従うことになります

curve(df(x, 3, 16), from = 0, to = 5, type = "l", ylim = c(0, 0.8))
f:id:yoshida931:20171031142612j:plain:w300

検定の結果として帰無仮説が棄却された場合には、要因の効果が認められる(群間に差がある)と推測します.

Rを使って一元配置分散分析(対応なし)

data <- c(10,10,5,9,6,5,4,5,10,5,11,2,9,12,4,3,10,4,6,1)
exp <- factor(c(rep(c("A","B","C","D"),5)))
summary(aov(data~exp))

            Df Sum Sq Mean Sq F value Pr(>F)  
exp          3  66.15   22.05   2.579 0.0898 .
Residuals   16 136.80    8.55                 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

結果、有意水準0.05で群間の差は認められませんでした.

F値の算出方法
一元配置分散分析 (対応なし) F値の算出方法 - 統計学備忘録 since2016


一元配置分散分析(対応あり)

対応のある1要因で分類される多群の検定(反復測定分散分析)
正確には繰り返しのない二元配置分散分析になります
サンプル:被験者4名(S1,S2,S3,S4)が5種類(testA,B,C,D,E)の検査を受けた結果
同じ被験者が繰り返して試験を受けることになる(繰り返し数は5回)
群間に点数差があるかどうかを知りたい
  ・被験者間に差があるか
  ・検査間に差があるか

サンプル(samp)  
data <- c(10,10,8,9,6,5,4,5,10,12,11,9,9,8,7,3,8,4,6,2)
data <- array(data, dim = c(4,5))
name <- list(Sub=c("s1","s2","s3","s4"), TEST=c("A","B","C","D","E"))
dimnames(data)<-name
data

    TEST
Sub   A B  C D E
  s1 10 6 10 9 8
  s2 10 5 12 8 4
  s3  8 4 11 7 6
  s4  9 5  9 3 2

Rを使って一元配置分散分析(対応あり)

data <- c(10,10,8,9,6,5,4,5,10,12,11,9,9,8,7,3,8,4,6,2)
sub <- factor(c(rep(c("s1","s2","s3","s4"),5)))
test <- factor(c(rep("A",4),rep("B",4),rep("C",4),rep("D",4),rep("E",4)))
summary(aov(data~sub+test))
  
            Df Sum Sq Mean Sq F value   Pr(>F)    
sub          3   24.2   8.067   3.681 0.043466 *  
test         4   99.7  24.925  11.373 0.000475 ***
Residuals   12   26.3   2.192                     
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

対応ありを無視したら…

summary(aov(data~sub))

#結果
            Df Sum Sq Mean Sq F value Pr(>F)
sub          3   24.2   8.067   1.024  0.408
Residuals   16  126.0   7.875 


summary(aov(data~test))

            Df Sum Sq Mean Sq F value  Pr(>F)   
test         4   99.7  24.925   7.403 0.00168 **
Residuals   15   50.5   3.367                   
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

自由度、平方和は同じ値で、もちろん平均平方も同じ値になっています.
しかし対応ありを無視したら、残差が大きなり、F値が小さくなっていることが確認できました.
対応ありを考慮すると有意差が出やすくなることが理解できました.


二元配置分散分析(対応なし)

対応のない2要因で分類される多群の検定 (要因分散分析)
2つの要因の組み合わせによる目的変数への影響を分析します

#サンプル:テストA、Bで各ポイント(p1,p2,p3,p4,p5)を獲得した人数
data <- c(10,10,8,9,6,5,4,5,10,12,11,9,9,8,7,3,8,4,6,2,4,10,8,8,4,6,2,10,12,13)
data <- array(data, dim = c(5,3,2))
name <- list(point=c("p1","p2","p3","p4","p5"), ","=c("","",""),TEST=c("A","B"))
dimnames(data)<-name
data

#二元配置分散分析(対応なし)のデータ配列のイメージ
#繰り返し数 3
TEST = A
point         
   p1 10  5 11
   p2 10  4  9
   p3  8  5  9
   p4  9 10  8
   p5  6 12  7

TEST = B
point        
   p1 3  4  6
   p2 8 10  2
   p3 4  8 10
   p4 6  8 12
   p5 2  4 13

分析はこのページで確認
二元配置分散分析 (1) - 統計学備忘録 since2016


二元配置分散分析(対応あり)

対応のある2要因で分類される多群の検定 (反復測定分散分析)
2つの要因の組み合わせによる目的変数への影響を分析します

#例)5人の被験者がAMとPMに試験(t1,t2,t3)を受けた時の点数
# 時刻(AM,PM)と試験(t1,t2,t3)
data <- c(10,10,13,9,12,9,11,9,10,12,7,9,9,8,7,3,8,4,6,2,4,10,8,8,6,9,8,10,12,10)
data2 <- array(data, dim = c(5,3,2))
name <- list(sub=c("s1","s2","s3","s4","s5"), "test"=c("t1","t2","t3"),"time"=c("AM","PM"))
dimnames(data2)<-name
data2  

#二元配置分散分析(対応あり)のデータ配列のイメージ
, , time = AM

    test
sub  t1 t2 t3
  s1 10  9  7
  s2 10 11  9
  s3 13  9  9
  s4  9 10  8
  s5 12 12  7

, , time = PM

    test
sub   t1 t2 t3
  s1  3  4  9
  s2  8 10  8
  s3  4  8 10
  s4  6  8 12
  s5  2  6 10

分析方法はちょっと特別な方法になります

#準備
sub2 <- factor(c(rep(c("s1","s2","s3","s4","s5"),6)))
tes2 <- factor(rep(c(rep("t1",5),rep("t2",5),rep("t3",5)),2))
tim2 <- factor(c(rep("AM",15),rep("PM",15)))

#被験者による効果、交互作用による効果を調べます
sub2
sub2:poi2
sub2:tes2
sub2:poi2:tes2

summary(aov(data~tes2*tim2+Error(sub2+sub2:tes2+sub2:tim2+sub2:tes2:tim2)))

#結果
Error: sub2
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4  19.53   4.883               

Error: sub2:tes2
          Df Sum Sq Mean Sq F value Pr(>F)
tes2       2  8.267   4.133   2.988  0.107
Residuals  8 11.067   1.383               

Error: sub2:tim2
          Df Sum Sq Mean Sq F value Pr(>F)  
tim2       1  45.63   45.63   11.75 0.0266 *
Residuals  4  15.53    3.88                 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Error: sub2:tes2:tim2
          Df Sum Sq Mean Sq F value  Pr(>F)   
tes2:tim2  2  81.07   40.53   11.47 0.00447 **
Residuals  8  28.27    3.53                   
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

交互作用の効果が有意になっています.したがって、単純効果の検定を行います.さらに単純効果が有意な場合には多重比較を実施します.次回、投稿する予定です.



参考)
山田 剛史, 杉澤 武俊, 村井 潤一郎; Rによるやさしい統計学, オーム社 , 2008
小野 滋;読めば必ずわかる分散分析の基礎  http://elsur.jpn.org/resource/anova.pdf