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

統計学備忘録 since2016

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

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

更新日2017.3.13
更新日2017.6.22


        f:id:yoshida931:20170721165122p:plain

一元配置分散分析とは、一つの因子に基づく“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 )
xx<- data.frame ( A , B , C , D )
xx<- t ( xx )
colnames ( xx )<- c ( 1:5 )

f:id:yoshida931:20170622170844p: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

f:id:yoshida931:20170622172124p:plain

 

水準間変動  ( A ~ D で4水準なので df = 4 - 1 = 3  )
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


水準内変動   (A~Dで4水準, 繰り返い数は5回なので df=4×(5-1)=16 )

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
    水準間平均平方 ÷ 水準内平均平方
= 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

正規分布のグラフ

グラフ初歩の初歩

手探り状態でやってます

赤文字のみRで実行

 

正規分布
#1正規分布の乱数を100個生成
x<-rnorm(100)
#ヒストグラム,freq=Fで確率密度,ylimでy軸範囲
hist(x,freq=F,ylim=c(0,0.6))
#枠の処理、lty線種は実践、btyで左と下のみ
box(lty=1,bty = "l")
#確率密度のグラフ,add=Tで重ねる、x軸範囲-3~3
curve(dnorm(x),add=T,-3,3)


これらをオリジナルの関数として記憶させます
mynorm<-function(n,g)
{
x<-rnorm(n)
hist(x,freq=F,ylim=c(0,0.6),xlab=g)
box(lty=1,bty = "l")
curve(dnorm(x),add=T,-3,3)
}


mynorm(100,"n=100")

 

f:id:yoshida931:20170212153358p:plain

 

#四つのグラフを2行2列のまとめてみます
par(mfrow = c(2,2))
mynorm(30,"n=30")
mynorm(100,"n=100")
mynorm(300,"n=300")
mynorm(1000,"n=1000")

f:id:yoshida931:20170212153327p:plain

二項分布のグラフ

こんな説明でどうでしょうか?
  f:id:yoshida931:20170524161843p:plain

コイン投げ(いかさまなし)
表の出る確率確率1/2
3回コイン投げを実施してみた場合・・・
表が出る回数(x)は0回、1回、2回、3回
xの確率は𝐵𝑖(3 , 0.5) に従う.
その確率分布を二項分布といいます.

    f:id:yoshida931:20170524161933p:plain

f:id:yoshida931:20170524161906p:plain

 

 

村上 正康,安田 正実; 統計学演習, 培風館,1989,p51のグラフをRで作ってみます
B(10,0.2),B(10,0.5),B(10,0.8)を重ねたグラフです
赤文字のみRで実行してください

 (foo.bar@baz.com様よりコメントいただき修正しております。簡単で理解しやすい方法です。)

t0.2 <- 0:6
P0.2 <- dbinom(t0.2, 10, 0.2)
t0.5 <- 1:9
P0.5 <- dbinom(t0.5, 10, 0.5)
t0.8 <- 4:10
P0.8 <- dbinom(t0.8, 10, 0.8)
plot(t0.2, P0.2, type="b", lty=2, xlab="", ylab="", xlim=c(0,10), ylim=c(0,0.4))
lines(t0.5, P0.5, type="b", lty=2)
lines(t0.8, P0.8, type="b", lty=2)

f:id:yoshida931:20170121213325p:plain

 

 

 

途中経過を忘れないために残しておきます

 
n=10 、確率0.2 
t0.2 <- 0:6;p<-0.2;n<-10
y<-factorial(n)*(p^t0.2)*(1-p)^(n-t0.2)
r<-factorial(t0.2)*factorial(n-t0.2)
y/r

n=10、確率0.5
t0.5 <- 1:9;p<-0.5;n<-10 
c<-factorial(n)*(p^t0.5)*(1-p)^(n-t0.5)
d<-factorial(t0.5)*factorial(n-t0.5)

n=10、確率0.8
t0.8 <- 4:10;p<-0.8;n<-10 
e<-factorial(n)*(p^t0.8)*(1-p)^(n-t0.8)
f<-factorial(t0.8)*factorial(n-t0.8)
e/f


p=0.2,0.5,0.8の結果をそれぞれ算出
間違わないように名前付け
P0.2<-(y/r)
P0.5<-(c/d)
P0.8<-(e/f)

linesでグラフを重ねる

plot(t0.2,P0.2,type="b",lty=2,xlab="",ylab="",xlim=c(0,10),ylim=c(0,0.4))
lines(t0.5,P0.5,type="b",lty=2)
lines(t0.8,P0.8,type="b",lty=2)

 

グラフに色をつける

library(RColorBrewer)
#RColorBrewerパッケージのサンプル
display.brewer.all()
 
ヒストグラムを塗ってみる
どの色セットを使用するかを指定する
cols <- brewer.pal(8,"Pastel1")   # brewer.pal(何色、パレット名)

y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mfrow = c(3,3),mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 
pos.x <- barplot(q,ylim=c(0,6),col=cols[1])
pos.x <- barplot(q,ylim=c(0,6),col=cols[2])
pos.x <- barplot(q,ylim=c(0,6),col=cols[3])
pos.x <- barplot(q,ylim=c(0,6),col=cols[4])
pos.x <- barplot(q,ylim=c(0,6),col=cols[5])
pos.x <- barplot(q,ylim=c(0,6),col=cols[6])
pos.x <- barplot(q,ylim=c(0,6),col=cols[7])
pos.x <- barplot(q,ylim=c(0,6),col=cols[8])
pos.x <- barplot(q,ylim=c(0,6),col=cols[9])

 

f:id:yoshida931:20170114232836p:plain

 

棒グラフの中央に散布図をプロットする方法 (pos.x)

y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 
pos.x <- barplot(q,ylim=c(0,6))
points(pos.x, p)

 

f:id:yoshida931:20170113120942p:plain

Rことはじめ

1、Rのインストール方法と使い方

2、以下のデータの母平均の95%信頼区間を求めてみます

73.0 75.8 74.0 72.3 76.9 82.5 69.1 80.9 74.9 71.6 73.5 84.6 81.3 82.5
78.1 71.3 76.3 74.9 74.8 72.5 78.5 71.5 70.0 81.1 76.6 77.9 75.5 70.8
81.9 75.1 85.4 82.0 73.4 76.9 78.5 77.0

 以下のファイルをご覧ください

Rの使い方.pdf - Google ドライブ

 

 

サンプル

対馬栄輝研究室 - 弘前大学医学部保健学科・大学院保健学研究科

http://www.hs.hirosaki-u.ac.jp/~pteiki/

統計解析→Rコマンダー

 

級内相関係数 ICC

修正日 2017.3.2
修正日 2017.6.21


ICC (Intraclass correlation coefficients)

ICCは信頼性の指標の一つとして利用されています.
ICCには3つの分類があり、6Caseに分類されています.

f:id:yoshida931:20170621175544p:plain

注)ICCの記載方法は統一されていない

ICCは分散分析の結果を活用しています
一元配置分散分析の重要な仮定
水準Aiに対してXij = μ+ Ai + Eijというモデルを考えた場合
μ+ Ai=水準Aiの母平均、Eij=誤差
条件1:標本は無作為抽出
条件2:誤差は母平均0の正規分布に従う(標本は正規分布に従う)
条件3:水準間の母分散はすべて等しい(等分散性の検定が必要)
ICCの場合、条件3は必ず保証されなければいけないものではない(対馬

一般可能性理論
一般可能性理論は一般化可能性研究(generalizability study:G研究)と決定研究(decision study:D研究)からなります.G研究は、検査の測定誤差に着目しており、その変動要因の成分とバラツキの大きさの推定値を求めます.そして各変動要因や交互作用について検討していく研究方法です.D研究では、検査の信頼性を検討します.G研究で得られた分散分析の推定値を利用し、一般化可能性係数を算出します.そして、どの程度の評価項目や評価者が必要なのかを検討する研究です(山西2005).分散分析では、水準の分散と誤差分散の推定値の比を見ることで、信頼性を推定します。これは一般化可能性理論と考えられます.

f:id:yoshida931:20170621195330p:plain

 

ICCは分散分析から求められる平均平方和の期待値を利用します.
全ては下記の検査結果の一覧表から始まります.

f:id:yoshida931:20170302184239p:plain


各自由度
被験者間:n—1
被験者内:n(k-1)
検者間:k-1
残差:(n-1)(k-1)

Case1の例
1名の検者で複数の被験者にそれぞれ数回検査を実施します.検者の誤差は考慮しません.ICC ( 1,1 ) は一回の測定に関する検者内信頼性、ICC ( 1,k ) はk回繰り返し測定した平均値に対する検者内信頼性. 

f:id:yoshida931:20170626153155p:plain

Case2の例
複数の検者で複数の被験者にそれぞれ数回検査を実施します.検者の誤差が発生します.ICC ( 2,1 ) はランダム効果による、つまり任意の検者に対する検者間信頼性.ICC ( 1,k ) の1回目、2回目、3回目を検者に置き換えると、k人の検者がn人の患者を検査した2元配置 ( 繰り返しなし ) の結果となる.

f:id:yoshida931:20170626153022p:plain
ICC ( 2,k )  は、複数(k人:3人)の検者が、複数(5人)の被験者を複数回(4回)繰り返した平均値に対する検者間信頼性.

f:id:yoshida931:20170626153428p:plain

Case3の例
複数の検者で複数の被験者にそれぞれ数回検査を実施します.検者が混合モデルとなります.ICC ( 3,1 ) は、固定効果による、つまり測定に参加した検者群そのものが対象となる.測定の精度(一致性:agreement)というよりも整合性(一貫性:consistency)を確かめるもの.
f:id:yoshida931:20170626161532p:plain

 

Rによるケース2(Shrout&Fleiss,1979)の算出方法
赤文字のみRにペーストして実行するのみ
二要因変量モデル

ICC(2,1)のRによる出力は以下のような分類になります McGraw & Wong (1996)
ICC(A,1) 一致性を示す級内相関係数(検者間の信頼性)
ICC(C,1) 一貫性を示す級内相関係数(検者内の信頼性)

準備
行列の配置を確認
行:被験者(測定の対象)
列:検者(評定者)
irrを読み込む
library(irr)

村山航;級内相関係数についての覚書p4のサンプル
例1
x<-c(1,2,3,4,5)
y<-c(1,2,3,4,5)
例2
x<-c(1,2,3,4,5)
y<-c(5,6,7,8,9)
例3
x<-c(1,2,3,4,5)
y<-c(3,5,7,9,11)

行列形式に変換cbind()してiccを算出

 

ICC(A,1) 一致性を示す級内相関係数
icc(cbind(x,y), model ="twoway",type = "agreement")
ICC(C,1) 一貫性を示す級内相関係数
icc(cbind(x,y), model ="twoway",type = "consistency")

  

ICCが0.7以上であれば信頼性は良好と考える
ICCは点推定であるので区間推定を提示することが重要

0.0-0.20    slight
0.21-0.40   fair
0.41-0.60   moderate
0.61-0.80   substantial
0.81-1.00   almost perfect

 

 

例題
対馬栄輝;頼性指標としての級内相関係数p.23の膝屈曲のサンプル

A B C D
126 122 131 125
137 143 141 141
113 119 115 105
153 143 135 144
146 157 150 149
161 157 160 160
110 109 105 113
145 151 152 156
126 141 132 122
114 126 130 125

x<-read.table("clipboard",header=T)

 

Case1 :一要因モデル (A=1回目、B=2回目、C=3回目、D=4回目とする)
検者内信頼性 ( Intra-rater reliability )

膝角度計測の検者信頼性を求める
一人の理学療法士が10人の患者に対してそれぞれ4回膝の角度を計測し信頼性を求める

ICC ( A,1 ) = ICC ( 1 , 1 ) 一致性を示す級内相関係数

icc(x, model ="oneway",type = "agreement")

Single Score Intraclass Correlation

Model: oneway
Type : agreement

Subjects = 10
Raters = 4
ICC(1) = 0.909

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,30) = 40.9 , p = 2.06e-14

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


参考)

繰り返しのある一元配置分散分析
       Df Sum Sq Mean Sq F value Pr(>F)

被験者内 treat 3 76 25.4 0.894 0.457
被験者間 pat 9 10319 1146.6 40.421 2.25e-13 ***

ICC(A,1) =
    ( 群間平均平方ー 郡内平均平方 ) ÷ ( 群間平均平方 + ( 繰り返し数 - 1 )×郡間平均平方 )
     = ( 1146.6 - 25.4 ) ÷ ( 1146.6 + ( 4 - 1 ) × 25.4  ) = 0.916912

 

4回づつ測定した10名分の平均値から検者内信頼性 ( Intra-rater reliability )を求める
ICC ( C,1 )ICC ( 1,k ) 一貫性を示す級内相関係数
icc(x, model ="oneway",type = "consistency")
結果

Single Score Intraclass Correlation

Model: oneway
Type : consistency

Subjects = 10
Raters = 4
ICC(1) = 0.909

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,30) = 40.9 , p = 2.06e-14

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

 

 Case2 :二要因モデル (4人の理学療法士をA,B,C,Dとする)
検者間信頼性 ( Inter-rater reliablity )

膝角度計測の検者信頼性(ばらつきの程度)を求める
4人の理学療法士が10人の患者に対してそれぞれ1回膝の角度を計測し信頼性を求める
ICC ( 2 , 1 ) 一致性を示す級内相関係数
icc(x, model ="twoway", type = "agreement")

Single Score Intraclass Correlation

Model: twoway
Type : agreement

Subjects = 10
Raters = 4
ICC(A,1) = 0.909       #一致性を示す級内相関係数

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,30) = 40.4 , p = 2.46e-14

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

 

膝角度計測の検者信頼性(ばらつきの程度)を求める
4人の理学療法士が測定した10人の患者の平均値から信頼性を求める
ICC ( C,1 ) = ICC ( 2 , 1 ) 一致性を示す級内相関係数
icc(x, model ="twoway", type = "consistency")

Single Score Intraclass Correlation

Model: twoway
Type : consistency

Subjects = 10
Raters = 4
ICC(C,1) = 0.908         #一貫性を示す級内相関係数

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,27) = 40.4 , p = 2.25e-13

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


ICC全てを一気に計算する
psychをインストールする
library ( psych )

ICC ( x )                        #大文字

Call: ICC(x = kneef)

Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute  ICC1   0.91  41  9  30  2.1e-14  0.79  0.97
Single_random_raters ICC2   0.91  40   9  27   2.3e-13   0.79   0.97
Single_fixed_raters ICC3   0.91   40   9   27  2.3e-13   0.78   0.97
Average_raters_absolute ICC1k   0.98   41   9   30   2.1e-14   0.94   0.99
Average_random_raters ICC2k   0.98   40   9   27   2.3e-13   0.94   0.99
Average_fixed_raters ICC3k   0.98   40   9   27   2.3e-13   0.93   0.99

Number of subjects = 10 Number of Judges = 4

 

反省
分散分析を復習してから再更新した方がよさそうです