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

統計学備忘録 since2016

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

二元配置分散分析 分散分析表の理解

投稿日2017.3.8
更新日2017.6.26

繰り返し数がすべて等しい場合 (=5)
水準A1 × 5回 = (A11 ,A12 ,A13 ,A14 ,A15 ) 
水準A2 × 5回 = (A21 ,A22 ,A23 ,A24 ,A25 )
つまり水準B (1,2,3) はA1、A2をそれぞれ5回繰り返したことになります

f:id:yoshida931:20170626180739p:plain
(仮想データ使用)

B1<-c(7,11,8,10,9,23,21,26,27,22)
B2<-c(9,11,10,9,11,12,16,16,17,15)
B3<-c(21,22,19,18,22,9,8,11,10,9)  
dat<-c(B1,B2,B3)
rown<-rep(c("A1","A2"),times=c(5,5))       #行のラベル名
A<-rep(rown,3)
B<-rep(paste0("B",1:3),each=10)
summary ( aov ( dat ~ A * B) )

     Df Sum Sq Mean Sq F value Pr(>F)
A  1   67.5  67.5    21.89    9.40e-05 ***
B  2   73.3  36.6       11.88    0.000259 ***
A:B 2   850.2  425.1    137.87   6.94e-14 ***
Residuals 24 74.0 3.1

分散分析表の計算方法を忘れないように記載しておきます!

f:id:yoshida931:20170626185319p:plain

Fvalue:F値(…の平均平方÷水準内変動の平均平方)
Df:自由度、SumSq:平方和、MeanSq:平均平方(平方和÷自由度)、

(平均平方、F値は平方和から算出可能.よって、平方和を調べるべし!)

 

http://yoshida931.hatenablog.com/entry/2017/02/19/144239

 

A:因子Aの水準間の変動を調べる
B:因子Bの水準間の変動を調べる
A:B:因子Aと因子Bの交互作用の変動を調べる
Residuals:水準内の変動を調べる

(変動:各測定値から平均値へのずれ(バラつき)が、「0」にならないように平均値からの変位の二乗和をとったもの.分散分析では平方和.)

 

Df(自由度)
A:a-1
B:b-1
A:B:(a-1)×(b-1)
Residuals:a×b×(n-1)

因子Aの水準数a、因子Bの水準数b、繰り返し数n、データ総数a*b*n
総数= length(data)、総平均=mean(data)とします
a=2; b=3; n=5; mea<-mean(dat)


因子A
の水準A1・A2間の変動(平方和)
ASumSq      
A1<-c( B1[1:5] , B2[1:5] B3[1:5] ) 
A2<-cB1[6:10] , B2[6:10] B3[6:10] ) 
n*b*( (mean(A1)-mea)^2+(mean(A2)-mea)^2)
=67.5

 

f:id:yoshida931:20170310132759p:plain

 

因子Bの水準B1・B2・B3間の変動(平方和)
B SumSq
n*a*( (mean(B1)-mea )^2+ ( mean (B2) - mea )^2 + (mean ( B3)-mea )^2 )
=73.26667

 

交互作用の変動
inshiA:inshiBのSumSq
水準内平均‐水準A平均-水準B平均+総平均
n*(mean ( B1[1:5] ) - mean(A1) - mean(B1) + mea )^2+
 n*( mean ( B1[6:10] ) - mean(A2) - mean(B1) + mea )^2+
 n*( mean ( B2[1:5] ) - mean(A1) - mean(B2) + mea )^2+
 n*( mean ( B2[6:10] ) - mean(A2) - mean(B2) + mea )^2+
 n*( mean ( B3[1:5] ) - mean(A1) - mean(B3) + mea )^2+
 n*( mean ( B3[6:10] ) - mean(A2) - mean(B3) + mea)^2

=850.2

 

f:id:yoshida931:20170310142941p:plain

 

各水準の平均値をプロットしたグラフ
A1<- c ( mean(B1[1:5] ) , mean ( B2 [1:5] ), mean ( B3[1:5] ) )
A2<- c ( mean(B1[6:10] ) , mean ( B2 [6:10] ), mean ( B3[6:10] ) )
x<- data.frame ( A1 = c(9,10,20.4) , A2 = c (23.8,15.2,9.4) )
plot ( x$A1 , type="o", lty=3, xaxt="n", xlim = c(0.8,3.2), ylim=c(7,25), ylab="", xlab = "", bty = "l" ) 
lines (x$A2, pch=19 , type="o")
axis (side = 1, at =1:3 , labels=c("B1","B2","B3") )
legend ( 2.7 , 25 , legend = colnames(x) , pch = pchs, lty = ltys )

f:id:yoshida931:20170310172329p:plain

 

水準内の変動
ResidualsSumSq
sum( ( B1[1:5] - mean ( B1[1:5] ) ) ^2 )+
 sum( ( B1[6:10] - mean (B1[6:10] ) )^2)+
 sum( ( B2[1:5] - mean (B2[1:5] )^2)+
 sum( ( B2[6:10] - mean ( B2[6:10] ) ) ^2)+
 sum( ( B3[1:5] - mean (B3[1:5] ) ) ^2 )+
 sum( ( B3[6:10] - mean (B3[6:10] ) ) ^2 )
=74

 

f:id:yoshida931:20170310140533p:plain


注意)水準内変動はいわゆる各水準内のバラつきです.
F値は水準内平均平方を分母にとるため、バラつきが
小さいと交互作用や因子内の水準間差が出やすくなります

例)
交互作用の変動:F=交互作用の変動の平均平方 ÷ 水準内変動の平均平方

 

ベクトルとデータフレームの関係

更新日2017.6.26
ここのページは、ちょこちょこ更新する予定

物忘れが激しいので何でも書いときます
ブログ名通りの使い方です

ベクトルをフレームにしたら列になる!
data.frame (列名1 = ベクトル1, 列名2 = ベクトル2, ... )

まずはここを忘れないように…


A<-c(10,6,10,9,10)
B<-c(10,5,5,12,4)
C<-c(5,4,11,4,6)
D<-c(9,5,2,3,1)
xx<-data.frame(A,B,C,D)

見慣れているデータフレームで確認してみます
xx<-data.frame(A,B,C,D)

f:id:yoshida931:20170626144828p:plain
                                              ( 画面はRstudio )

行ラベルを変更
rownames(xx)<-paste0("P", 1:5)     # 列ラベル変更colnames()

f:id:yoshida931:20170626145118p:plain


全て水準と要素をデータフレームに取り込みます

x<-c(A,B,C,D)
pat <- paste0("P", 1:5)
treat <- gl(4,5,labels = LETTERS[1:4])      #4つのラベルを5回繰り返し
xy <- data.frame(treat, pat, dat=x)

   

   f:id:yoshida931:20170626145636p:plain


r-de-r様よりいただきましたコメントをもとに修正しました
かなりスッキリしました!

 

Rstudioでのデータ保存、読み込み

よく使うので書き留めておきます

事前準備
Rstudio 四つの画面から編成されていますので、勝手に名前をつけます.
①書き込み画面、②Consol画面、③環境画面、④ファイル画面

f:id:yoshida931:20170620165251p:plain


まずRスクリプト上に書いているベクトルを保存してみます
①にデータを書き込みます
t3.4<-c(150,132,144,139,118,135,123,133,152,136)

f:id:yoshida931:20170620164100p:plain

  出典)出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010, 表3.12

 

保存するために①に
save ( t3.12 , file = "表3.12.RData" )

f:id:yoshida931:20170620164424p:plain

 

④にデータファイルができてます

f:id:yoshida931:20170620164601p:plain

 

完成したファイルをクリックすると③にデータを読み込みます.

f:id:yoshida931:20170620164731p:plain

 最後に①にt3.12と入力するだけでベクトルになります.
データファイルの保管場所は、Rprojectのフォルダになります
もちろんデータフレームも同じ要領で保存できます.

x<-c(LETTERS[1:10])
y<-t3.12
z<-data.frame(x,y)
save(z,file = "表3.12名前入り.RData")

f:id:yoshida931:20170620170456p:plain

Zの部分をクリックすると①に次の表が出力されます

f:id:yoshida931:20170620170607p:plain

Zという一つのアルファベットの中に上記のデータセットが格納されているのです.
便利です!

 

 

 

 

比率の信頼区間

比率の信頼区間

例)有権者から2,400人を無作為抽出した結果、1,250人は支持していたことがわかった.有権者の支持率の95%信頼区間を求めよ.
出典 日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012, p149

二項分布の問題です.
考え方)試行回数n=2,400回、成功回数X=1,250回、成功確率p=1,250÷2,400 

成功回数Xは二項分布 B ( n , p ) にしたがう.
母集団の成功確率をpとした場合の推定値は

   f:id:yoshida931:20170619182944p:plain

二項分布の正規近似よりXを標準化した

   f:id:yoshida931:20170619183856p:plain

nが大きい場合、次のZは近似的に標準正規分布 N ( 0 ,1 ) にしたがう.

   f:id:yoshida931:20170619184137p:plain
したがって95%信頼区間

         f:id:yoshida931:20170619184527p:plain

n=2400, X=1250 
p1 <-  ( X/n )  -  qnorm(0.975)  * sqrt ( ( ( X/n )*( 1 - X/n ) ) / n )
p2 <-  ( X/n ) +  qnorm(0.975)  * sqrt ( ( ( X/n )*( 1 - X/n ) ) / n )
       p1 < 95%信頼区間 < p2

Rの関数では
binom.test(1250, 2400)

Exact binomial test

data: 1250 and 2400
number of successes = 1250, number of trials = 2400, p-value = 0.04328(両側確率)
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.5006234 0.5409923
sample estimates:
probability of success
0.5208333

 

 

比率の差の信頼区間
A群はn1人の中のn1'人が支持している、B群はn0人の中のn0'人が支持している.両群の支持率に差の95%信頼区間を求めよ.

f:id:yoshida931:20170619194012p:plain

f:id:yoshida931:20170619191059p:plain

もとめるものは… f:id:yoshida931:20170619194047p:plain

 

Rの関数では
prop.test(c(n1', n0'), c(n1,n0), correct=F)

例)
prop.test(c(50, 38), c(76,79), correct=F)

2-sample test for equality of proportions without continuity correction

data: c(50, 38) out of c(76, 79)
X-squared = 4.9384, df = 1, p-value = 0.02627
alternative hypothesis: two.sided
95 percent confidence interval:
0.02353533 0.33022883
sample estimates:
prop 1 prop 2
0.6578947 0.4810127

 
出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010, p142

 

 

並行箱ひげ図

エクセルでデータのみコピー

t<-read.table("clipboard")
x<-t$V1
y<-t$V2
boxplot(x,y,xaxt="n",main="タイトル")   #x軸のラベルを消去
name<-c("ラベル1","ラベル2")     #x軸のラベルを定義
axis(side=1,at=c(1,2),labels=name)   #x軸のラベルを更新

 

f:id:yoshida931:20170615115707p:plain

このような小技が意外と役立つ!

二項分布からの最尤推定

二項分布から最尤推定に挑戦

二項分布のグラフ - 統計学備忘録 since2016

例)試行回数3回・確率1/2に二項分布
       𝐵𝑖(3 , 0.5)

f:id:yoshida931:20170613163424p:plain



表が出る確率をp=1/2とすると、次のような分布が考えられます.
確率pのときに表が〇回出る確率(同時確率

    f:id:yoshida931:20170613163444p:plain

これらの標本から得られる同時確率は、右に記載している確率です.
L(P) = 1回目表の確率×2回目表の確率×3回目表の確率

例)表、表、ウラの場合を考えると
L(P) = p×p×(1-p) = p^2×(1-p) となります

ここでpが未知の場合のp=0.2とp=0.8を仮定してみます.

L(P) = p×p×(1-p) = p^2×(1-p)
L(0.2)=0.032
L(0.8)=0.128
p=0.8のL(P)が大きく、p=0.8が尤もらしい値と言えます.

 

L(P) 尤度関数 ( maximum function )
pのいろいろな値における尤もらしさを表す関数である.
尤度関数L(P)を最大にする値が最尤推定値、関数としては最尤推定量である.

求めるためには微分が必要になります・・・
L(P) = p×p×(1-p) = p^2×(1-p) を微分して=0としたときのpが最尤推定値となります.
dL(P)/dp = (p^2 – p^3)’ = 2p-3p^2 = 0
最尤推定値p=2/3
    ↓
確率密度関数fp(y)をもつ母集団から
大きさ3の確率標本を取り出す{1回目, 2回目, 3回目}
同時確率密度:fp(y1,y2,y3) = p×p×(1-p)
これをpの関数とみなしたものが尤度関数です
L(p)= fp(y1)* fp(y2)*,…,* fp(yn)


つまり次ようになります
θを未知の母数とした確率密度関数fθ(y)の母集団から
大きさnの確率標本 {Y1,Y2,…,Yn} を取り出した場合 
同時確率分布:fθ(yn)
同時確率密度:fθ(y1,y2,…,yn) = fθ(y1)* fθ(y2)*,…,* fθ(yn)
これをθの関数とみなしたものが尤度関数です
L(θ) = Πfθ(y) = fθ(y1)* fθ(y2)*,…,* fθ(yn)


{Y1,Y2,…,Yn}を観測して{a1,a2,…,an}が得られた場合、
尤度関数を最大にするθのことをY1= a1、Y2= a2、…、Yn= anが観測されたときのθの最尤推定値( maximum likelihood estimate )といいます.


例)
表の出る確率がpで{Y1=表、Y2=表、Y3=ウラ}が観測された場合
標本から得られる確率は L(p) = p^2*(1-p)^3-2
これを微分して=0としたpの値が最尤推定値となります
dL(p)/dp= (p^2 – p^3)’ = 2p-3p^2 = 0
p=2/3
表が2回観察されたときのpの最尤推定値は2/3


対数尤度 ( log likelihood function )

尤度関数を考えるとき、積の形なので数学的に扱うのが不便であるため対数をとる
logL(θ) = ∑log{ fθ(y1)* fθ(y2)*,…,* fθ(yn)}
二項分布の例) 
L(p) =Π(p^ai×(1-p)^n-ai ) = p^T×(1-p)^n-T    (T=∑ai)
logL(p)= log { p^T×(1-p)^n-T)} = Tlogp+(n-T)log(1-p)
最尤推定値を求めるためにやはり微分して=0とします
dlogL(p)/dp= { Tlogp+(n-T)log(1-p) }’ = T/p – (n-T)/(1-p) = 0
p=T/n=(1/n)*∑ai (相対頻度)
例)p = (1/3)*(1+1) = 2/3
pの最尤推定値は相対頻度

 

参考文献
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991

t分布からの信頼区間

自由度5のt分布より
確率点:qt(0.025,5,lower.tail = F)
積分布:pt(2.570582,5,lower.tail = F)
確率密度:dt(2.570582,5)


拡張期血圧を6回測定して95%信頼区間を求める
DBP<-c(86,92,88,94,89,88)
vDBP<-var(DBP)

t値
 { mean(DBP)-X } ÷ {sqrt(vDBP/6)}

 

95%信頼区間の求め方
mean(DBP)-qt(0.025,5,lower.tail = F)*sqrt(vDBP/6)
mean(DBP)+qt(0.025,5,lower.tail = F)*sqrt(vDBP/6)

関数で確認
t.test(DBP)

One Sample t-test

data: DBP
t = 74.326, df = 5, p-value = 8.352e-09
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
86.40461 92.59539
sample estimates:
mean of x
89.5

出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010,p132

 


t分布のグラフ
x <- seq(-4, 4, 0.1)
plot(x, dt(x, 20), type="l")

for(i in c(2,3,4)){
for(i2 in 2:4){
curve(dt(x, 6-i), col=i, type="l",add=T, )
}
}
labels<-c("自由度","黒:20","赤:4","緑:3","青:2")
legend("topleft", legend = labels)

f:id:yoshida931:20170613195843p:plain