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

読者です 読者をやめる 読者になる 読者になる

統計学備忘録 since2016

理学療法、医療、福祉、小児

繰り返し for

いつも忘れるので書いておこう
#x=3に2を3回足す式
#3+2+2+2=9

myf<-function(){
x<-3
for(i in 1:3){
x<-x+2}
return(x)
}
myf()

 

#Xに2を3回加える式
#x+2+2+2
#xに5を代入すると5+2+2+2=11


myf<-function(x){
a<-x
for(i in 1:3){
a<-a+2}
return(a)
}
myf(5)


#Xに3を掛けることを3回繰り返す
#X*3*3*3
#xに5を代入すると
5*3*3*3

myf<-function(x){
for(i in 1:3){
x<-x*3}
return(x)
}
myf(5)

 

#10にxを3回加える式
#10+x+x+x
#x=2のとき10+2+2+2=16

myf<-function(x){
a<-10
for(i in 1:3){
a<-a+x}
return(a)
}
myf(2)

 

級内相関係数 ICC(1,1), ICC(1,k)

1人の検者により検査結果のばらつきを求める

ICC(1,1)

= ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数 ) ÷

 ( ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数) - 被験者内の平均平方)

 

被験者=水準
被験者間の平均平方
 =∑(被験者平均―総平均) ^2 ÷ (被験者数-1)

被験者内の平均平方

     =∑∑(X―被験者平均) ^2 ÷被験者数×(繰り返し数-1)

 

「一元配置分散分析 (対応なし) F値の算出方法」より
http://yoshida931.hatenablog.com/entry/2017/02/19/144239

 

仮想データ

一人の理学療法士が被験者A~Dに対してそれぞれ5回づつ検査を実施した結果

f:id:yoshida931:20170313152805p:plain

Rで算出してみます
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)
x<-data.frame(A,B,C,D)

パッケージirrのダウンロードが必要です
library(irr)
icc(t(x))   #行に各被験者のデータが並ぶように設定します

Single Score Intraclass Correlation
Model: oneway
Type : consistency
Subjects = 4
Raters = 5
ICC(1) = 0.24
F-Test, H0: r0 = 0 ; H1: r0 > 0
F(3,16) = 2.58 , p = 0.0898
95%-Confidence Interval for ICC Population Values:
-0.079 < ICC < 0.877

 

ICCの値と95%信頼区間の求め方
被験者数4、繰り返し数5
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)
ALL<-c(A,B,C,D)  #全データ 
mean(ALL)     #総平均


被験者間の平方和
y<-5*{(mean(A)-mean(ALL))^2+
(mean(B)-mean(ALL))^2+
(mean(C)-mean(ALL))^2+
(mean(D)-mean(ALL))^2}
=66.15

被験者内の平均平方
水準間平均平方=∑(水準平均―総平均) ^2 ÷ (水準数-1)
66.15÷(4-1) = 22.05


被験者内の平方和
sum( (A-mean(A))^2)+
sum( (B-mean(B))^2)+
sum( (C-mean(C))^2)+
sum( (D-mean(D))^2)
=136.8

被験者内
平均平方=∑∑(X―水準平均) ^2 ÷水準数×(n-1)
136.8÷{4☓(5-1)} = 8.55

 

ICC(1,1) 
= ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数 ) ÷
( ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数) - 被験者内の平均平方)
= ( (22.05-8.55)/5)/( ( (22.05-8.55)/5)+ 8.55)
=0.24

 

ICC(1,1)の95%信頼区間の求め方
被験者間の平均平方÷被験者内の平均平方
= 2.578947
qf(0.975,16,3)
= 14.23152
qf(0.975,3,16)
= 4.076823
( (2.578947/4.076823)-1)/ {(2.578947/4.076823)+(5-1)}
= -0.07931044
( (2.578947*14.23152)-1)/( (2.578947*14.23152)+(5-1))
=0.877157
95CI = [ -0.07931044 0.877157 ]

 


1人の検者により検査の信頼性を調べる
ICC(1,k)
= (被験者間の平均平方-被験者内の平均平方)÷ 被験者間の平均平方

 

Rで算出してみます
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)
x<-data.frame(A,B,C,D)
icc(t(x),"oneway","consistency", "average")

Average Score Intraclass Correlation
Model: oneway
Type : consistency (整合性、一貫性)
Subjects = 4
Raters = 5
ICC(5) = 0.612

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(3,16) = 2.58 , p = 0.0898

95%-Confidence Interval for ICC Population Values:
-0.581 < ICC < 0.973

 

ICCの値と95%信頼区間の求め方
ICC(1,5)
= (被験者間の平均平方-被験者内の平均平方)÷ 被験者間の平均平方
= (22.05-8.55) / 22.05
= 0.6122449

ICC(1,5)の95%信頼区間の求め方
1-(1/(2.578947/4.076823))= -0.5808091
1-(1/(2.578947*14.23152))= 0.9727538
95CI = [-0.5808091 0.9727538]

 


次のような関数で全てのICCが算出可能になります
パッケージpsychを使用します
library(psych)
ICC(t(x))

Call: ICC(x = t(x))
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.24 2.6 3 16 0.09 -0.079 0.88
Single_random_raters ICC2 0.24 2.6 3 12 0.10 -0.084 0.88
Single_fixed_raters ICC3 0.24 2.6 3 12 0.10 -0.094 0.88
Average_raters_absolute ICC1k 0.61 2.6 3 16 0.09 -0.581 0.97
Average_random_raters ICC2k 0.61 2.6 3 12 0.10 -0.631 0.97
Average_fixed_raters ICC3k 0.61 2.6 3 12 0.10 -0.752 0.97

Number of subjects = 4 Number of Judges = 5

 

f:id:yoshida931:20170313152828p:plain

今回使用した結果は
Single_raters_absolute、Average_raters_absoluteの部分になります

 

対馬栄輝;頼性指標としての級内相関係数

http://www.hs.hirosaki-u.ac.jp/~pteiki/research/stat/icc.pdf

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

データは仮想です

繰り返し数がすべて等しい場合 (=5)

   

因子B

 

 

B1

B2

B3

因子A

A1

7

9

21

11

11

22

8

10

19

10

9

18

9

11

22

A2

23

12

9

21

16

8

26

16

11

27

17

10

22

15

9

エクセルにペーストして、青文字のみをコピー


Data<-c(x$B1,x$B2,x$B3)x<-read.table("clipboard",header=T)


以下赤文字のみRで実行

inshiA1<-c(rep("A1",5),rep("A2",5))

inshiA<-c(rep(inshiA1,3))

inshiB<-c(rep("B1",10),rep("B2",10),rep("B3",10))

summary ( aov ( Data ~inshiA* inshiB) )

              Df Sum Sq Mean Sq F value   Pr(>F)   

inshiA         1   67.5    67.5   21.89 9.40e-05 ***

inshiB         2   73.3    36.6   11.88 0.000259 ***

inshiA:inshiB  2  850.2   425.1  137.87 6.94e-14 ***

Residuals     24   74.0     3.1                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

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

 

Df

SumSq

MeanSq

Fvalue

inshiA

1

67.5

67.5

21.89

inshiB

2

73.3

36.6

11.88

inshiA:inshiB

2

850.2

425.1

137.87

Residuals

24

74

3.1

 

 

Df:自由度、SumSq:平方和、MeanSq:平均平方(平方和÷自由度)、

Fvalue:F値(…の平均平方÷水準内変動の平均平方)

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

 

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

 

inshiA:因子Aの水準間の変動を調べる

inshiB:因子Bの水準間の変動を調べる

inshiA:inshiB:因子Aと因子Bの交互作用の変動を調べる

Residuals:水準内の変動を調べる

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

 

Df(自由度)

inshiA:a-1

inshiB:b-1

inshiA:inshiB:(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

 

因子Aの水準A1・A2間の変動(平方和)

inshiASumSq=67.5         
A1<-c(Data[1:5],Data[11:15],Data[21:25]) #(A1,B1),(A1,B2),(A1,B3)全てのデータ
mean(A1)
A2<-c(Data[6:10],Data[16:20],Data[26:30])  #(A2,B1),(A2,B2),(A2,B3)全てのデータ
mean(A2)
n*b*( (mean(A1)-mean(Data))^2+(mean(A2)-mean(Data))^2)

 

f:id:yoshida931:20170310132759p:plain

 

因子Bの水準B1・B2・B3間の変動(平方和)
inshiB SumSq73.3
B1<-c(Data[1:10])  #(A1,B1),(A2,B1)全てのデータ
mean(B1)
B2<-c(Data[11:20])  #(A1,B2),(A2,B2)全てのデータ
mean(B2)
B3<-c(Data[21:30])  #(A1,B3),(A2,B3)全てのデータ
mean(B3)
n*a*( (mean(B1)-mean(Data))^2+(mean(B2)-mean(Data))^2+(mean(B3)-mean(Data))^2)

 

交互作用の変動
inshiA:inshiBのSumSq=850.2
水準内平均‐水準A平均-水準B平均+総平均
n*(mean(Data[1:5])-mean(A1)-mean(B1)+mean(Data))^2+
 n*(mean(Data[6:10])-mean(A2)-mean(B1)+mean(Data))^2+
 n*(mean(Data[11:15])-mean(A1)-mean(B2)+mean(Data))^2+
 n*(mean(Data[16:20])-mean(A2)-mean(B2)+mean(Data))^2+
 n*(mean(Data[21:25])-mean(A1)-mean(B3)+mean(Data))^2+
 n*(mean(Data[26:30])-mean(A2)-mean(B3)+mean(Data))^2

f:id:yoshida931:20170310142941p:plain

 

各水準の平均値をプロットしたグラフ
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

 

水準内の変動
ResidualsSumSq74
sum( (Data[1:5]-mean(Data[1:5]))^2)+
 sum( (Data[6:10]-mean(Data[6:10]))^2)+
 sum( (Data[11:15]-mean(Data[11:15]))^2)+
 sum( (Data[16:20]-mean(Data[16:20]))^2)+
 sum( (Data[21:25]-mean(Data[21:25]))^2)+
 sum( (Data[26:30]-mean(Data[26:30]))^2)

f:id:yoshida931:20170310140533p:plain


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

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

 

相関係数の区間推定

サンプル:日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012,p125

 

身長のデータ
大学生
x<-c(172,167,184,175,176,175,170,180,170,179,167,175,174,162,165,163,170,169,165,175)
父親
y<-c(165,165,178,176,150,171,172,175,170,156,163,170,165,160,163,170,163,165,160,172)

 

cor.test関数を使用します
同じ数の要素をもつベクトル2つの相関関係、また相関係数の 95% 信頼区間,p値を求めます.Pearsonの相関係数をもとめるのであれば次の式を実行するのみでOKです.


  cor.test(x,y)

 

Pearson's product-moment correlation

data: x and y
t = 1.4587, df = 18, p-value = 0.1619
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1371017 0.6711054
sample estimates:
cor
0.3251458

cor.test(x,y)の引数
x, y:同じ大きさの数字ベクトル
alternative(hypothesis): "two.sided", "greater", "less".
method : "pearson(相関係数cor)", "kendall"(相関係数tau), "spearman"(相関係数rho).

   (r-de様のコメントを参考に修正しました)

 

次に上記関数結果の95%信頼区間[-0.1371017 0.6711054]の求め方を忘れないように残しておきます

大学生の標準偏差
Vx<-(sum( (x-mean(x) )^2) )/20
Sx<-sqrt(Vx)
Sx

父親標準偏差
Vy<-(sum( (y-mean(y) )^2) )/20
Sy<-sqrt(Vy)
Sy

大学生と父親の共分散
Sxy<-(sum((x-mean(x) )*(y-mean(y) ) ) )/20

標本の相関係数r
Sxy/(Sx*Sy)
r=0.325

 

相関係数ρの区間推定

 

標本の相関係数rの分布を考えてみます.
r=0の場合
t=r*(sqrt(n-2)/sqrt(1-r^2) )は自由度n-2のt分布に従う

r≠0の場合(FisherのZ変換
rの分布はρが0から離れると非対称性が強くなります.
Z={log((1+r)/(1-r) ) }/2と変換した場合、Zの分布は対称になります.
またZは近似的にN({log( (1+ρ)/(1-ρ) ) }/2,1/(n-3) )に従うことが知られています.

 

r=0.325をZ={log((1+r)/(1-r) ) }/2に代入してZ=0.337を求めます.

{log((1+ρ)/(1-ρ) )=Rとおくと

(0.337-R)=±1.96×標準偏差
0.337-1.96/√(17)<R<0.337+1.96/√(17)
2×(0.337-1.96/√(17))<log((1+ρ)/(1-ρ)<2×(0.337+1.96/√(17))
‐0.138×2< log((1+ρ)/(1-ρ) )<0.671×2
したがってρの95%信頼区間
(exp(-0.276)-1)/(exp(-0.276)+1)
-0.137
(exp(1.626)-1)/(exp(1.626)+1)
0.671

 

また勉強して更新していきますのでコメントお願いします.

なぜt検定を繰り返してはいけないのか?

sample は「石村卓夫;分散分析のはなし,東京図書,1992,p137」

A<-c(12,14,16)
B<-c(13,15,17)
C<-c(15,17,19)
D<-c(17,19,21)

x<-c(rep("A",3),rep("B",3),rep("C",3),rep("D",3))
y<-c(A,B,C,D)
stripchart(y~x,vertical=T)

f:id:yoshida931:20170224165734p:plain

視覚的にはAD間、BD間には差があるように見えます

 

重要な前提です
Aの母平均をμa
Bの母平均をμb
Cの母平均をμc
Dの母平均をμd とします
「μa=μb=μc=μd」を検定する場合にt検定を繰り返してはいけないのです.
ある特別なペアの違いなどを確認したい場合にはt検定を使用します.研究内容により異なるので間違いないように.

 

参考 http://www.gen-info.osaka-u.ac.jp/MEPHAS/ave.html
2種類の既存薬AとBを組み合わせた配合薬Cの配合効果を評価する場合.この場合、既存薬A,Bのそれぞれの効果と配合薬Cの効果を比較します.ここで いいたいのはCが既存薬A、Bの両方よりも効果があるということです.「CがAよ りも優れている、かつCがBよりも優れている」ということを 示します.つまり帰無仮説はC=AかつC=Bとなり、A=B=Cとは異なります.このような場合には2標本t検定を繰り返して用いるほうが適切だと考えられます.

 

一元配置分散分析
oneway.test(y~x,var.equal = T)
#F = 3.6875, num df = 3, denom df = 8, p-value = 0.06214
#有意水準5%で差が「ない」という結果になりました

 

#全ての群の母分散は等しいと仮定してt検定を繰り返してみます
AB<-t.test(A,B,paired=FALSE, var.equal=T, conf.level=0.95)
AC<-t.test(A,C,paired=FALSE, var.equal=T, conf.level=0.95)
AD<-t.test(A,D,paired=FALSE, var.equal=T, conf.level=0.95)
BC<-t.test(B,C,paired=FALSE, var.equal=T, conf.level=0.95)
BD<-t.test(B,D,paired=FALSE, var.equal=T, conf.level=0.95)
CD<-t.test(C,D,paired=FALSE, var.equal=T, conf.level=0.95)

 

#それぞれの検定結果のP値のみを取り出してみます
str(AB3);str(AC3);str(AD3);str(BC3);str(BD3);str(CD3)
AB 0.573
AC 0.14
AD 0.0376
BC 0.288
BD 0.0705
CD 0.288

 やはりAD間には有意水準5%で「差がある」という結果になりました.しかし分散分析では差がないという結果でした.

 

分散分析の帰無仮説
有意水準0.5において「μa=μb=μc=μd」と仮定する

t検定の帰無仮説
有意水準0.5において「μa=μb」、
かつ有意水準0.5において「μa =μc」、
かつ有意水準0.5において「μa =μd」、
かつ有意水準0.5において「μb=μc」、
かつ有意水準0.5において「μb=μd」、
かつ有意水準0.5において「μc¬=μd」と仮定した検定になります.

したがってt検定を使った「μa=μb=μc=μd」の確率は(1-0.05)^6となります.少なくとも一組には差がある確率は1-(1-0.05)^6=0.2649081となり、有意水準26.5%の検定になってしまいます。つまりαエラーの危険性が高くなるのです.このような矛盾を解消するための手法が多重比較です(しかし、本当にt検定は繰り返して使ってはいけないのだろうか? 本当に有意水準の高い検定になっているのだろうか? 2群間比較にはt検定しかないのだが…疑問が残ります)

多重比較(Tukey)を行ってみます

TukeyHSD(aov(y ~ x),ordered = TRUE)
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered

Fit: aov(formula = y ~ x)

$x
  diff  lwr    upr   p adj
B-A 1 -4.2294189 6.229419 0.9252929
C-A 3 -2.2294189 8.229419 0.3245304
D-A 5 -0.2294189 10.229419 0.0609499
C-B 2 -3.2294189 7.229419 0.6297636
D-B 4 -1.2294189 9.229419 0.1441838
D-C 2 -3.2294189 7.229419 0.6297636

有意水準0.05ではどのペアにも差はありませんでした.

等分散性の検定

準備としてカイ二乗分布とF分布が必要になります

カイ二乗分布
確率変数Z1,Z2,…Znが互いに独立で標準正規分布にしたがうとき
W=Z1^2+Z2^2+…+Zn^2
の従う分布を自由度nのカイ二乗分布とよび、χ^2(n)で表します.
期待値E[W]=n, 分散V[E]=2n

カイ二乗分布の性質
① 二つの確率変数W1,W2が独立にχ^2(n1), χ^2(n2)に従うとき
W1+ W2は、χ^2(n1+n2)に従います.(カイ二乗分布の再生性)
② 母集団N(μ,σ^2)からの無作為標本X(X1,X2,…X3)について
Xiを標準すると Zi= (Xi-μ)/ σとなる
Ziは互いに独立で標準正規分布に従うので
∑{(Xi-μ)/ σ}^2 はカイ二乗分布χ^2(n)に従うことになります
③ 標本不偏分散
s^2={∑{(Xi-X平均) ^2 }/(n-1) より
(n-1) ×s^2=∑(Xi-X平均) ^2
母集団N(μ,σ^2)からの無作為標本Xi=(X1,X2,…X3)は互いに独立でN(μ,σ^2)に従うので
W=∑{(Xi-X平均) ^2}/ σ^2 = {(n-1) ×s^2}/σ^2
母平均μをX平均で置き換えると自由度がn-1になります
したがって{(n-1) ×s^2}/σ^2は自由度n-1のカイ二乗分布χ^2(n-1)に従います.

f:id:yoshida931:20170222184024p:plain

http://tips-r.blogspot.jp/2015/01/r.htmlより)


2標本の母分散の比の検定
散らばりを表す分散は非負の値なので、その大きさを考える場合には差ではなく比で考える方が適切です.

例)単位はmm
A群の身長の標準偏差はσa=5 (分散25)
B群の身長の標準偏差はσb=10 (分散100)
σa/σb=0.5 A群の方がB群に比べると2倍散らばりの小さい集団である
現実的には標準偏差が便利で理解しやすいのですが、推定での取り扱いは上述のような分散を使います.

 

F分布
条件
・確率変数U:自由度mのカイ二乗分布に従う
・確率変数V:自由度nのカイ二乗分布に従う
・UとVは独立である

フィッシャーの分散比 F= (U/m)/(V/n)

Fが従う確率分布を自由度(m,n)のF分布F(m,n)という

2標本の標本分散s1^2, s2^2を考えます
標本不偏分散s^2={∑{(Xi-X平均) ^2 }/(n-1) 
標本s1^2:{(m-1) ×s1^2}/σ1^2は自由度m-1のカイ二乗分布χ^2(m-1)に従います
標本s2^2:{(n-1) ×s2^2}/σ2^2は自由度n-1のカイ二乗分布χ^2(n-1)に従います
s1^2とs2^2は独立です

F={(m-1) ×s1^2}/(σ1^2)*(m-1)÷{(n-1) ×s2^2}/(σ2^2)*(n-1)
= (s1^2/σ1^2)÷(s2^2/σ2^2)
= (s1^2/s2^2)*(σ2^2/σ1^2)
の従う分布は、自由度(m-1,n-1)のF分布F(m-1,n-1)となります

 

 

それでは始めます
赤文字のみRで実行

2標本の場合

母分散の比の検定になります
二つの母集団からの標本をXi~N(μx,σx^2),Yi~N(μy,σy^2)とします.
帰無仮説:σx^2=σy^2、(σx^2)/(σy^2)=1
対立仮設:σx^2≠σy^2、(σx^2)/(σy^2)≠1

サンプルは「石村卓夫;分散分析のはなし,東京図書,1998,p69」
Xi<-c(12.2,18.8,18.2)
Yi<-c(26.4,32.6,31.3)
σx二乗<-sum((Xi-mean(Xi) )^2)/2;σx二乗
σy二乗<-sum((Yi-mean(Yi) )^2)/2;σy二乗

T<-σx二乗/σy二乗
1.25
F分布を使用します

qf(0.025,2,2,lower=F) # F 分布の上側5%点
39
qf(0.975,2,2,lower=F) # F 分布の下側5%点
0.026

T=1.25は棄却域には入らないので帰無仮説は棄却できない
標本をXiとYiの母分散は等しいと仮定できることになります

 

SPSSではlenene検定が使用されていますので以下に算出しておきます
Levene's Test もF分布を使用していますが、上記方法とは異なります
leveneTest(x,group)

Levene's Test for Homogeneity of Variance (center = median)
    Df  F value  Pr(>F)
group  1  0.0031  0.9585
    4

 

 

 

 

3標本以上の場合

 

繰り返し数が等しいとき:Hartleyの検定
各水準の分散(Si^2 )について
Si^2の最大分散÷最小分散が1に近づくことで等分散性を仮定します
最大分散/最小分散=max{Si^2}/min{Si^2}

 

f:id:yoshida931:20170222184102p:plain

http://www2.vmas.kitasato-u.ac.jp/lecture0/statistics/hartley.pdf


検定表の見方:Fmax(k , f)(k) = Fmax(水準数,繰り返し数-1)(水準数)
帰無仮説:各水準の分散は等しい

サンプルは「石村卓夫;分散分析のはなし,東京図書,1998,p108」
A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

x<-c(A,B,C,D,E)

SuppDistsをインストトールします
library(SuppDists)
factor()カテゴリーを要素としたベクトル
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))
levels()でグループ化をstr()で構成を確認できます
fmax=max( by(x, group, FUN=var) ) / min( by(x, group, FUN=var) )  
5.70

by関数の引数(data,INDICES,FUN=,simplify=)
data:data.frameやmatrixもOK
INDICES:factorをsubsetに指定
FUN :引数dataに適用する関数を指定。ここでは分散varを指定。
simplify :結果のデータ形式。TRUE:scalar、FALSE:list or array。

 

積分布pmaxFratio
pmaxFratio( fmax, 繰り返し数の自由度, 水準数)
1-pmaxFratio(fmax, 2, 5)  #p値
0.8154272

帰無仮説は棄却されず、等分散性を仮定することができる

 

 

繰り返し数が異なるとき:Bartlettの検定
Rヘルプより
繰り返し数は異なりませんが同じサンプルで試行してみます
以下実施方法のみ記載
A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

x<-c(A,B,C,D,E)
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))

bartlett.test(x,group)


Bartlett test of homogeneity of variances
data: x and group
Bartlett's K-squared = 1.2291, df = 4, p-value = 0.8733

 


Leveneの検定
最もよくみかける等分散性の検定だと思われます
水準の繰り返し数に縛りはありません

A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

統計量(Wikipediaより)
W={(N-k)/(k-1)}×{∑Ni(Zi.-Z..)^2/∑∑(Zij-Zi.)^2}

数値はsample
k:比較する群の数水準
水準数i  1,2,3,4,5
繰り返しj  1,2,3
N:総サンプル数
Ni:第i群のサンプル数
自由度N-k=15-5=10
自由度k-1=5-1=4
Z.. 全データの平均値
Zi. i番目の水準内の平均値 
Yij:第i群のj番目の変数
Zij:| Yij – Yiの平均|、| Yij – Yiの中央値|


carパッケージをインストールします
x<-c(A,B,C,D,E)
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))

leveneTest(x,group)

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 0.1256 0.9698
10

p値=0.9698で帰無仮説は棄却されず、等分散性を仮定した検定が可能になります

 

引用・参考
バイオインフォマティクス入門 http://bio-info.biz/
石村卓夫;分散分析のはなし,東京図書,1998

一元配置分散分析 (対応なし) F値の算出方法

(2017.3.13更新)

赤文字の部分だけRにペーストして実行

一元配置分散分析とは、一つの因子に基づく“3つ以上の母平均の差の検定”です
因子Factor:実験結果に影響を与えると考えられる要因
水準Lebel:因子をいくつかに分ける基準 (水準=群)

例題)
治療因子に基づく、患者20名の満足度データを比較する(繰り返し数=5)
因子Factor=治療、水準Level=治療A, 治療B, 治療C, 治療D
因子<c(水準)としてベクトルを作成
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)

f:id:yoshida931:20170313105833p:plain

繰り返し数(n)=5(順序が入れ替わっても問題はない)

 

ALL<-c(A,B,C,D)  #全データ 
mean(ALL)     #総平均
Levels<-c(rep("A",5),rep("B",5),rep("C",5),rep("D",5))    # 水準oneway.test(ALL~Levels,var.equal = T)         #一元配置分散

One-way analysis of means
data: AllD and AllL
F = 2.5789, num df = 3, denom df = 16, p-value = 0.08979

それではF = 2.5789を求めてみます

 

分散分析の理解
分散分析は、観測データにおける変動(誤差変動、水準内変動)と各要因による変動(水準間変動)に分解し効果を判定する方法です。(2元配置では交互作用による変動も考えます)。
1、 データの分散成分の平方和(全変動)を分解し、誤差による変動と要因の効果による変動を分離します。
2、 平方和を自由度で割った平均平方を算出します。
3、 統計量はF値になります。「要因の効果によって説明される平均平方(水準間変動)」を分子、「誤差によって説明される平均平方(水準内変動)」を分母にとった値が分散分析の統計量です。
4、 F検定。F分布により各効果の有意性を判定します。

以下、水準間変動、水準内変動という表現を使用します

 

まず各水準の平均間に差がない場合を考えてみます
赤線は全平均
F = 0.30185, num df = 3, denom df = 16, p-value = 0.8236

f:id:yoshida931:20170219143818p:plain


残っている変動は水準内の変動になります。つまり、「全変動=水準内の変動」ということになります。偏差の2乗和(平方和)を全2乗和、または全変動といいます。

逆に考えると水準間変動>水準内変動の場合には、「各水準間に差がある」と考えることができます。したがって「水準間変動÷水準内変動」が分散分析の統計量になります

ここで一元配置分散分析の数学モデル(X=x1,x2,…)を考えてみます
X = 総平均 + (水準平均―総平均) + (X―水準平均)
(X-総平均) = (水準平均―総平均) + (X―水準平均)

f:id:yoshida931:20170313111752p:plain


ここで差の平方をとって考えていきます
(X-総平均)^2 、 (水準平均―総平均) ^2 、(X―水準平均) ^2
例) 20個       4個         20個 
平方和∑にすると次の式が成り立ちます
n=各水準に含まれている繰り返し数(例題ではn=5)

∑(X-総平均)^2=n∑(水準平均―総平均) ^2 + ∑(X―水準平均) ^2

f:id:yoshida931:20170313112921p:plain

 

全変動=∑∑(X-総平均)^2
水準間変動=n∑(水準平均―総平均) ^2
水準内変動=∑∑(X―水準平均) ^2

 

水準間変動
y<-5*{(mean(A)-mean(ALL))^2+
(mean(B)-mean(ALL))^2+
(mean(C)-mean(ALL))^2+
(mean(D)-mean(ALL))^2}
=66.15

平均平方和
水準間平方和:n∑(水準平均―総平均) ^2
各水準平均は正規分布N( 総平均+(水準平均―総平均)、σ^2/n)に従います。
したがって∑(水準平均―総平均) ^2÷(σ^2/n)は自由度(水準数-1)のχ二乗分布に従います。
水準間平均平方=∑(水準平均―総平均) ^2 ÷ (水準数-1)
66.15÷(4-1) = 22.05


水準内変動

sum( (A-mean(A))^2)+
sum( (B-mean(B))^2)+
sum( (C-mean(C))^2)+
sum( (D-mean(D))^2)
=136.8

平均平方和
水準内平方和:∑(X―水準平均) ^2
それぞれの水準における∑(X―水準平均) ^2 ÷ σ^2 はn-1のχ二乗分布に従う統計量になります。したがって∑∑(X―水準平均) ^2 ÷ σ^2は自由度、水準数×(n-1)のχ二乗分布に従います。
水準内平均平方=∑∑(X―水準平均) ^2 ÷水準数×(n-1)
136.8÷{4☓(5-1)} = 8.55


ポイント
誤差を正規分布に従う確率変数と考えます。したがって例題の水準A標本(10,6,10,9,10)は、正規母集団(6.55+9.0、σ^2)からの大きさ5の標本ということになります。∑(A-6.55)^2/σ^2 は自由度4のχ二乗分布に従う統計量となります。

分散分析の統計量F=水準間平均平方÷水準内平均平方
{∑(水準平均―総平均) ^2 ÷(水準数-1)} ÷ {∑∑(X―水準平均) ^2 ÷ n(水準数-1)}
22.05÷8.55=2.5789


Rで分散分析表を求めてみます
anova(aov(ALL~Levels))

Analysis of Variance Table
Response: ALL
       Df Sum Sq Mean Sq F value Pr(>F)
Levels 3 66.15 22.05 2.5789 0.08979 .
Residuals 16 136.80 8.55
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

f:id:yoshida931:20170313123012p:plain