最終更新日 ： 2007.11.25
変量効果
いわゆる対応ありのデータに対する解析方法です。ある個体（心理学では被験者）から反復測定によって得られたデータをそのままglm()で今までと同じように解析するのは間違いです。なぜなら、いままでやってきたglm()による解析は「データセットがすべて独立した別の個体から得られたものである」と仮定できたものであったのに対して、今回のような例は独立性を仮定できないわけです。したがって、独立でないデータを独立であると仮定して解析することは間違いということです。
では、結局のところ、glm()で解析して何が悪いかというと、
ということです。要するにp値の計算の仕方が異なるということです。
Rではlme4というパッケージに含まれているlmer()という関数を用いることができます。lme4というパッケージは最初から実装されているわけではないので、使用するためにはこれをインストールする必要があります。
> install.packages("lme4") #lme4のインストール
> library(lme4) #パッケージlme4の呼び出し
p229のBOX12.1と同じ解析を行ってみましょう。まずはglm()で分散分析表を作成しましょう（ただし、これによって得られるp値をそのまま採用してはいけない）。
# データフレームを作成しておく
> mydata
LEAF CACONC
1 1 0.4567
2 1 0.6033
3 1 0.3918
4 1 0.3016
5 2 0.7291
6 2 0.6527
7 2 0.5165
8 2 0.5704
9 3 0.7727
10 3 0.6171
11 3 0.7330
12 3 0.8495
13 4 0.3526
14 4 0.4497
15 4 0.3459
16 4 0.5446
# 分散分析表の作成
> attach(mydata) #データフレームにパスを通す
> result <- glm(CACONC ~ LEAF)
> summary(aov(result))
Df Sum Sq Mean Sq F value Pr(>F)
LEAF 3 0.280862 0.093621 8.7021 0.002443 **
Residuals 12 0.129100 0.010758
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# LEAFを変量要因として解析
> model1 <- lmer(CACONC ~ (1|LEAF), family = gaussian, data = mydata)
> summary(model1)
Linear mixed-effects model fit by REML
Formula: CACONC ~ (1 | LEAF)
Data: mydata
AIC BIC logLik MLdeviance REMLdeviance
-12.15 -10.60 8.075 -19.49 -16.15
Random effects:
Groups Name Variance Std.Dev.
LEAF (Intercept) 0.020715 0.14393
Residual 0.010758 0.10372
number of obs: 16, groups: LEAF, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.55545 0.07649 7.262
> detach(mydata) #最後に通したパスを解除するのを忘れずに
ここで出力結果の見方を解説しておきます。Random effects:の部分にある値がBOX12.1でいうところの「調整平方和を使う分散成分」の値です。
つづいて正しいp値を算出するための手順を紹介します。とりあえず、F値を求めなければなりませんが、簡単にF値を算出するために以下のような関数を定義しました。
KHH <- function(var.comp, group.num){
vec.n <- length(var.comp)
dat <- c()
for(i in 1:vec.n){
dat[i] <- sum(var.comp[i:vec.n] * group.num[i:vec.n])
}
f.value <- c()
for(j in 1:vec.n){
f.value[j] <- dat[j] / dat[j+1]
}
f.value
}
関数KHHの引数であるvar.compは先ほど求めた分散成分の値です。これをベクトルとして指定してください。また、引数group.numはグループの数というか、繰り返し観測の回数です。教科書の例題では1枚の葉から4回の反復観測が行われたということですから、LEAFに対する反復回数は4です。そして、ベクトルの最後に必ず1を入れておいてください。
> vc <- c(0.020715, 0.010758) > gn <- c(4, 1) > KHH(vc, gn) [1] 8.702175 NA
これで"正しい"F値が求められました（今回はたまたまglm()での解析結果と一致しましたがね）。これよりp値を求めるためには次のようにします。
> 1 - pf(8.702175, df1=3, df2=12) [1] 0.002442767
以上で"正しい"p値を求めることができました。
BOX12.3の解析例
> mydata2
CACONC2 LEAF2 PLANT2
1 1.6647 1 1
2 2.0276 1 1
3 2.2036 2 1
4 2.2238 2 1
5 1.9945 3 1
6 1.9045 3 1
7 1.9579 4 2
8 1.7441 4 2
9 2.4959 5 2
10 2.5570 5 2
11 1.8732 6 2
12 1.9046 6 2
13 2.2992 7 3
14 2.6597 7 3
15 2.7374 8 3
16 2.9178 8 3
17 2.8447 9 3
18 2.4489 9 3
19 2.6918 10 4
20 2.7402 10 4
21 3.3845 11 4
22 2.8510 11 4
23 2.2435 12 4
24 2.7079 12 4
> attach(mydata2)
#glmによる解析
> result2 <- glm(CACONC2 ~ PLANT2 + LEAF2)
> summary(aov(result2))
Df Sum Sq Mean Sq F value Pr(>F)
PLANT2 3 2.71433 0.90478 21.4480 4.114e-05 ***
LEAF2 8 1.26200 0.15775 3.7395 0.02001 *
Residuals 12 0.50621 0.04218
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#lmerによる解析
> model2 <- lmer(CACONC2 ~ (1|PLANT2) + (1|LEAF2),
+ family = gaussian, data = mydata2)
> summary(model2)
Linear mixed-effects model fit by REML
Formula: CACONC2 ~ (1 | PLANT2) + (1 | LEAF2)
Data: mydata2
AIC BIC logLik MLdeviance REMLdeviance
21.39 24.92 -7.693 13.92 15.39
Random effects:
Groups Name Variance Std.Dev.
LEAF2 (Intercept) 0.057827 0.24047
PLANT2 (Intercept) 0.124305 0.35257
Residual 0.042182 0.20538
number of obs: 24, groups: LEAF2, 12; PLANT2, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.378 0.194 12.26
> detach(mydata2)
ちなみに、次のように指定しても同じ結果が得られる。こちらの方がより教科書に忠実な形で表示される。
> mod <- lmer(CACONC2 ~ (1|PLANT2/LEAF2),
+ family = gaussian, data = mydata2)
> summary(mod)
Linear mixed-effects model fit by REML
Formula: CACONC2 ~ (1 | PLANT2/LEAF2)
Data: mydata2
AIC BIC logLik MLdeviance REMLdeviance
21.39 24.92 -7.693 13.92 15.39
Random effects:
Groups Name Variance Std.Dev.
LEAF2:PLANT2 (Intercept) 0.057827 0.24047
PLANT2 (Intercept) 0.124305 0.35257
Residual 0.042182 0.20538
number of obs: 24, groups: LEAF2:PLANT2, 12; PLANT2, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.378 0.194 12.26
これより"正しい"F値とp値を算出します。引数に指定するベクトルvc2とgn2は以下のように必ず階層順に値を入れてください。
> vc2 <- c(0.124305, 0.057827, 0.042182) > gn2 <- c(6, 2, 1) > KHH(vc2, gn2) [1] 5.725348 3.741786 NA > 1 - pf(5.725348, df1=3, df2=8) #PLANT2に対するp値 [1] 0.02165530 > 1 - pf(3.741786, df1=8, df2=12) #LEAF2(PLANT2)に対するp値 [1] 0.01996275
以上の結果より、glm()によって得られたp値とlmer()を用いて算出したp値を比較してみましょう。今回は結果としてはどちらの場合においても有意差が認められますが、glm()を用いた場合（すなわち独立性を仮定した場合）は随分とp値が小さいようです。つまり、「実際には有意ではないのに有意であると結論づけてしまうことがあるかもしれない」ということです。
| glm() | lmer() | |
| PLANT2 | 0.00004114 | 0.02165530 |
| LEAF2(PLANT2) | 0.02001 | 0.01996275 |
BOX12.4の解析例
> glm.res <- glm(WGAIN ~ SIRE + DAM)
> lmer.res <- lmer(WGAIN ~ SIRE + (1|DAM))
> summary(aov(glm.res))
Df Sum Sq Mean Sq F value Pr(>F)
SIRE 4 1.03439 0.25860 3.1285 0.06533 .
DAM 5 1.28898 0.25780 3.1188 0.05930 .
Residuals 10 0.82658 0.08266
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> summary(lmer.res)
Linear mixed-effects model fit by REML
Formula: WGAIN ~ SIRE + (1 | DAM)
AIC BIC logLik MLdeviance REMLdeviance
29.79 35.77 -8.896 12.50 17.79
Random effects:
Groups Name Variance Std.Dev.
DAM (Intercept) 0.087170 0.29525
Residual 0.082769 0.28770
number of obs: 20, groups: DAM, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.2905 0.2535 5.090
SIRE2 -0.3465 0.3585 -0.966
SIRE3 0.2055 0.3585 0.573
SIRE4 -0.4130 0.3585 -1.152
SIRE5 -0.0830 0.3585 -0.231
Correlation of Fixed Effects:
(Intr) SIRE2 SIRE3 SIRE4
SIRE2 -0.707
SIRE3 -0.707 0.500
SIRE4 -0.707 0.500 0.500
SIRE5 -0.707 0.500 0.500 0.500
> vc3 <- c(0.087170, 0.082769)
> gn3 <- c(2, 1)
> KHH(vc3, gn3)
[1] 3.106344 NA
> 1 - pf(3.106344, df1=5, df2=10) #DAMに対するp値
[1] 0.05992412
BOX12.6以降の解析はよく分からない（というかよく読んでいない）。ここで扱っている内容は「一般化線形混合モデル」の範疇に含まれる内容なので、より詳しい解説は以下の参考文献にゆずる。他にも「Multilevel Modeling」などのキーワードでgoogleしてみるとよい。
