2009.1.3
2009.1.22 最終修正
ランダム効果モデルで(たぶん)もっとも簡単な場合、一元配置分散分析でサンプル数が揃っている場合の計算例とその説明です。たとえば、a個の”もの”をn回ずつ測定した場合で、ランダム効果しかないので、混合モデルとは言えないでしょう。推定するのは、全体の平均、”もの”間の分散、同じ”もの”の測定間の分散となります。
分散分析表は以下のようになります:
−−−−−−−−−−−−−−−−−−−−−−−−−−
変動源 自由度 平方和 平均平方
”もの” a-1 SSamong σ^2+n・(σA)^2
測定 a(n-1) SSwithin σ^2
−−−−−−−−−−−−−−−−−−−−−−−−−−
σ^2が同じ”もの”の測定間の分散、(σA)^2が”もの”間の分散です。平方和の部分は固定効果の場合と同じです。
例は、スネデカーとコクランの『統計的方法』10.13節に登場する、表10.13.1のかぶら菜のデータです。計算はRで行なっています。分散分析の仮定(等分散の正規分布)は満たしているとします。以下のようなデータフレームにデータを入れておきます。4枚の葉について4回ずつ測定したというデータで、このページでの表記で言うとa=4,n=4となります。xf11は順序のないfactorで葉の番号を表し、yy1が目的変数です。
> mydata13
xf11 yy1
1 1 3.28
2 1 3.09
3 1 3.03
4 1 3.03
5 2 3.52
6 2 3.48
7 2 3.38
8 2 3.38
9 3 2.88
10 3 2.80
11 3 2.81
12 3 2.76
13 4 3.34
14 4 3.38
15 4 3.23
16 4 3.26
固定型の分散分析をして見ます。この分散分析は自由度、平方和、平均平方までランダム効果モデルの場合と同じで、Fも同じです(F検定の意味は同じでないので注意)(広津千尋『分散分析』など参照)。
> summary(aov(yy1~xf11,data=mydata13))
Df Sum Sq Mean Sq F value Pr(>F)
xf11 3 0.88837 0.29612 44.853 8.52e-07 ***
Residuals 12 0.07923 0.00660
---
(以下略)
ここから、先の分散分析表にあるように、分散を推定すると、
同じ葉の測定間の分散は0.00660で、
(0.88837/3-0.07923/12)/4= 0.0724
が葉のあいだの分散の推定値となります。
Rの関数lme()を使えば(nlmeパッケージにあります)、以下のようになります(上記の計算よりもわかりやすいでしょう)。
> gmydata13 <-groupedData(yy1 ~ 1 | xf11,data =mydata13)
> gmydata13
Grouped Data: yy1 ~ 1 | xf11
xf11 yy1
1 1 3.28
2 1 3.09
3 1 3.03
4 1 3.03
5 2 3.52
6 2 3.48
7 2 3.38
8 2 3.38
9 3 2.88
10 3 2.80
11 3 2.81
12 3 2.76
13 4 3.34
14 4 3.38
15 4 3.23
16 4 3.26
> resm13.lme<-lme(yy1~1,data=gmydata13,random=~1|xf11)
> resm13.lme
Linear mixed-effects model fit by REML
Data: gmydata13
Log-restricted-likelihood: 9.277319
Fixed: yy1 ~ 1
(Intercept)
3.165625
Random effects:
Formula: ~1 | xf11
(Intercept) Residual
StdDev: 0.2690357 0.0812532
Number of Observations: 16
Number of Groups: 4
分散の平方根(標準偏差)が計算されているので、
0.2690357と0.0812532をそれぞれ二乗してみると、上で求めた分散と同じであることがわかります。
Rの関数としてはlmer()もこういったときによく使われます(パッケージlme4にあります)。この例でlmer()を使うと、以下のようになります、
> resm13.lmer<-lmer(yy1~1+(1|xf11),data=mydata13)
> resm13.lmer
Linear mixed model fit by REML
Formula: yy1 ~ 1 + (1 | xf11)
Data: mydata13
AIC BIC logLik deviance REMLdev
-12.55 -10.24 9.277 -20.74 -18.55
Random effects:
Groups Name Variance Std.Dev.
xf11 (Intercept) 0.0723802 0.269036
Residual 0.0066021 0.081253
Number of obs: 16, groups: xf11, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.166 0.136 23.27
結果は先に示したものと同じです。Fixed effectsの検定結果は、固定効果として1つまり定数項を指定しているので、平均が0かどうかを検定したものです。
この例は、パッケージnlmeの付属データセットRailでの計算例(Pinheiro&Batesの本の1.1節に出てくるもの)とよく似ています。