以前の続き。 - 知行合一日記

leaf<-rep(1:4,each=4)
disc<-rep(1:4,4)
conc<-runif(16)*runif(16)*sqrt(leaf)
data<-data.frame(leaf=factor(leaf),disc=factor(disc),conc)

ここで，４つの葉っぱからそれぞれ４つの穴を開けて，カルシウム濃度を調べる。で葉っぱごとにどれだけばらつくか（分散）を調べたい場合を考えよう。どのようなときにこの解析を使うかといわれれば，調べる際に非常にバラツキが大きくなるような場合は，サンプリング回数を多くしなければならないし，慎重に行う必要がある。このデータ自体は自分で作ったが，大元はこちらの本を参考にしました。
さて，実際に解析する。

> model1<-lm(conc~leaf,data)
> coef(model1)
(Intercept)       leaf2       leaf3       leaf4 
 0.19492586 -0.02008269  0.54902054 -0.03600447 
> model2<-aov(conc~leaf,data)
> coef(model2)
(Intercept)       leaf2       leaf3       leaf4 
 0.19492586 -0.02008269  0.54902054 -0.03600447

どちらの関数でも同様の結果が得られる。でもより単純にコンセプトを学習するために対比を用いる。

> op<-options(contrasts=c("contr.sum","contr.treatment"))

contrasts
aov または lm といったモデル当てはめで用いられる既定のcontrasts。
長さ 2 の文字列で、最初は順序無し因子とともに使われる関数を、
第二は順序付き因子とともに使われる関数を与える。

http://www.is.titech.ac.jp/~mase/mase/html.jp/temp/options.jp.htmlより引用。
つまり，対比行列をどのように扱うかという問題。
ここでは，ANOVAにおける切片を全体平均にするために，デフォルトの"contr.treatment"から"contr.sum"にしている。polyが何を意味するか全然わからんかったけど，順序尺度に対して用いるのね。やっとわかったよ。で本題に戻りまして，

> model3<-aov(conc~leaf,data)
> coef(model3)
(Intercept)       leaf1       leaf2       leaf3 
  0.3181592  -0.1232333  -0.1433160   0.4257872 
> model4<-lm(conc~leaf,data)
> coef(model4)
(Intercept)       leaf1       leaf2       leaf3 
  0.3181592  -0.1232333  -0.1433160   0.4257872 
> mean(data$conc)
[1] 0.3181592

と全体平均が切片になっている。
さらに，ここでは葉ごとにどれだけばらつくかを調べたいので，

> summary(model3)
            Df  Sum Sq Mean Sq F value  Pr(>F)  
leaf         3 0.96951 0.32317  3.5809 0.04676 *
Residuals   12 1.08299 0.09025

前回も書いたけど， $\sigma^2_e$ は残差の平均平方が不偏推定量。
で，leafの $\sigma^2_S$ は｛（leafの平均平方）−（残差の平均平方）｝／leaf内の計測数で
推定される。この場合は，

> (0.32317-0.09025)/4
[1] 0.05823

と推定される。この0.058というleafの分散成分（ $\sigma^2_S$ ）の不偏推定量。なぜこの式で導かれるのかなかなか理解できなかった。がやっとわかった。通常の固定効果の場合 $\sigma^2_A$ の推定値は，帰無仮説が正しいなら $\sigma^2_e$ 。帰無仮説が棄却されるなら $r*\sigma^2_A$ になる。rは処理内の計測数。ちなみにこれはバランスデータを仮定。固定効果の場合は処理群間の平均値 $a$ を比較するのが目的なので， $\sigma^2_e$ と処理の期待平均平方 $\sigma^2_A$ の比を比較するだけでよい。

変量効果の場合，平均平方の期待値はr* $\sigma^2_S+\sigma^2_e$ となる。つまり $E(SSA)$ ；期待平方和 $=(a-1)(n\sigma^2_S+\sigma^2_e)$ から式を単純に移項するだけで， $\sigma^2_S=(MSA-MSE)/n$ となる。
要するに，固定効果における一元配置分散分析には分散成分は必要とされない。効果を検定する際にも，平均平方の比からF統計量を用いるためである。
一方，変量効果の一元配置分散分析では $\sigma^2_S=0$ を帰無仮説とするので，この $\sigma^2_S=0$ を検定するために， $\sigma^2_S$ を推定する必要がある。

ただし，この分散成分を推定する方法には，推定値が負になってしまうことがある場合やアンバランス型データでは推定量にバイアスが生じる場合もあるので，問題となっている。
（詳細はまた）

ちなみに，ランダム切片モデルとしてlmerで解析すると

> modellmer<-lmer(conc~1+(1|leaf),data)
> summary(modellmer)
Linear mixed model fit by REML 
Formula: conc ~ 1 + (1 | leaf) 
   Data: data 
   AIC   BIC logLik deviance REMLdev
 19.09 21.41 -6.545    10.99   13.09
Random effects:
 Groups   Name        Variance Std.Dev.
 leaf     (Intercept) 0.058228 0.24131 
 Residual             0.090250 0.30042 
Number of obs: 16, groups: leaf, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)   0.3182     0.1421   2.239

となり，変量効果であるleafの分散成分が先ほど求めた値と等しい。
つまり，制限付き最尤法と先のANOVA推定量は等しくなる（バランスデータの場合ね）。