一般線形モデルを踏まえて、ここからは変量効果が導入される。まず、変量効果がどのようなものかを示すために、単純な分散分析の固定効果の部分を変量効果に置き換えて、両者を対比している。このステップが理解をすすめるためのポイントのように思われる。すなわち、変量効果は何も特別なことをしているわけではなく、固定効果ならば変数の各水準の効果を検討する代わりに、変量効果ではそれぞれの水準が目的変数に及ぼす効果を正規分布からランダムに取り出されたものとして推定しているに過ぎないということがわかってくる。そうすると、変量効果の使いどころもおのずと見えてくるのではないか。

具体的には、階層化されたデータにおいて、上位レベルのグループごとに下位レベルの変数の誤差項が独立ではない時とか、上位グループの違いを統制して分析をしなければならない時、一般線形モデルでは統制変数として所属する上位グループをダミー変数として投入する。このときグループによる違い、すなわちグループという水準間の差に関心があれば、そのままダミー変数を固定効果として検討すればよい。しかし、グループごとに誤差項が独立ではないことが想定されるが、グループが異なることによる違いそのものには関心がなく、かつグループ数が大量にある時には、これらを固定効果として推定するのは（自由度の消費など、、）問題が大きい。そのような場合、上位グループを変量効果としてモデルに投入し、各グループの効果を水準として推定するのではなく、それぞれの効果は平均が0の正規分布に従ってランダムに生じたとみなして、その分散だけを推定することで倹約的にモデルを推定する。これが混合モデル、ないしはマルチレベルモデルの1つの捉え方だといえるのではないか。

1元配置分散分析の固定効果を、そのまま変量効果に置き換えてみる。一般線形モデルにおける「m1.2.1」の1元配置分散分析と比較していく。

> m2.1.1m <- lmer(ACHIEVE ~ 1 + (1|DEPART), data=dat, REML=F)
> summary(m2.1.1m)
Linear mixed model fit by maximum likelihood t-tests use Satterthwaite
  approximations to degrees of freedom [lmerMod]
Formula: ACHIEVE ~ 1 + (1 | DEPART)
   Data: dat

     AIC      BIC   logLik deviance df.resid
   222.6    227.3   -108.3    216.6       33

Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.55900 -0.91239  0.03589  0.78482  1.90087

Random effects:
 Groups   Name        Variance Std.Dev.
 DEPART   (Intercept)  6.477   2.545
 Residual             20.050   4.478
Number of obs: 36, groups:  DEPART, 6

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)
(Intercept)   82.528      1.279  6.000   64.51 9.33e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> AIC(m2.1.1m)
[1] 222.5667

まずRのlmer()関数について。モデル式において、変量効果はカッコで表現する。この場合、固定効果は切片のみ（モデル式では「1」が切片を示す）で、学科（DEPART）が変量効果（学科といっても、本書での想定はランダムにクラス分けされたグループのようなものを想定するが）となっている。カッコ内を見ると、縦棒「|」の前に「1（＝切片）」があり、後ろにDEPARTがある。したがって、DEPAERTを変量効果とし、DEPAERTの各水準ごとに切片を推定するということを示す。すなわち、ランダム切片モデルである。ただし、分散分析はそもそも切片しかないが。

モデル式の後はデータの指定と推定方法の指定。推定方法は、完全最尤法と制限最尤法があり、標準では制限最尤法となる。ここでは「REML=F」としているので、完全最尤法で推定している。完全最尤法では、全てのパラメータを最尤法で推定するが、制限最尤法では固定効果のパラメータに一般線形モデルの値を用いる。完全最尤法では、「尤法は、固定効果(fixed effects、母数因子)を推定することによる自由度の減少を考慮しないので変量効 果(random effects、変量因子)の分散推定量には偏りが生じる」ということらしい。ただし、制限最尤法では、モデルを比較する際に用いるAICなどを求める際に固定効果の変数の数が考慮されないため、固定効果の変数のみが異なるモデルの比較ができない。AICを比較する際は、完全最尤法で推定する必要がある。本書では、両方で推定しAICは完全最尤法の結果を用いて検討している。

> m2.1.1r <- lmer(ACHIEVE ~ 1 + (1|DEPART), data=dat, REML=T)
> summary(m2.1.1r)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to
  degrees of freedom [lmerMod]
Formula: ACHIEVE ~ 1 + (1 | DEPART)
   Data: dat

REML criterion at convergence: 214.1

Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.51008 -0.87509  0.04156  0.73409  1.86322

Random effects:
 Groups   Name        Variance Std.Dev.
 DEPART   (Intercept)  8.441   2.905
 Residual             20.050   4.478
Number of obs: 36, groups:  DEPART, 6

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)
(Intercept)   82.528      1.401  5.000   58.89 2.67e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> AIC(m2.1.1r)
[1] 220.1479

同モデルを制限最尤法で推定したもの。

書籍に示されている結果と比較すると、主要な数値は一致しているが検定に関する出力は異なっている。これは、混合モデルの検定について標準的な方法が定まっていないため、実装者の考え方の違いを反映しているようだ。推定される変量効果の分散や、固定効果の値については一致した値が得られているが、Rでは変量効果に関する検定が行なわれていない。変量効果の分散が「0」かどうかの検定は、変量効果によって説明される成分が「0」ではないことを検定するものであるため、やはり必要だと思う。Rでどのように検討すれば良いのかは今後の課題として調べたい。lmerTestパッケージにはrand()関数があり、lmer()関数の結果を入力すると変量効果の検定らしきものが出力されるがSPSSの結果とは異なっている（p値も異なるので、違う見方で等価な検定をしているということでもなさそう）。

この他にも、書籍に添付されているRのサンプルスクリプトへのコメントでも指摘されているが、AICも異なっている。AIC()関数に結果を投入するとAICの値が出力されるが、完全最尤法の場合には一致するが制限最尤法ではSPSSと一致しない。本書では、制限最尤法でも固定効果の数までカウントすることが原因ではないかと推測している。AICの値はモデル比較の際にしか利用しないので、制限最尤法での値に問題があっても特に不便はないだろう。

次に、本書では「ネストした混合分散分析」としているが、ちょっと違うような気もする。ここでは、学部が固定効果で学科が変量効果となっている。つまり、学部による違いは注目するが、学部内のサブグループはランダムに振り分けらえたクラスであって、クラス間の違いには特に関心はなく統制だけしたいというモデルのようだ。ネストした混合分散分析というならば、変量効果がネストしているような場合ではないだろうか。マルチレベルモデルでは、上位のグループレベルの変数を変量効果とすることが多いが、ここではこれが逆転しているように見える。こうした場合、推定された固定効果や変量効果は何を示しているのだろうか？

回帰式の性質から考えれば、変量効果となっている学科の分散の推定値は、学部が同一だとした場合の分散ということになるのか。しかし、学部によって学科間の分散が大きく異なる場合、このモデルでの分析は適当ではないようにも思える。学科の効果が学部ごとに異なるような場合は、学部と学科の交互作用項を変量効果にするとか、学部ごとに分けて分析するとか、そういったことが必要かもしれない。この辺は、考え出すとよくわからなくなるので今後ゆっくり考えることにして、とりあえずは上位のグループ変数を変量効果とするモデルをまず理解しよう。

> m2.2.1 <- lmer(ACHIEVE ~ FACULTY + (1|DEPART), data=dat, REML=T)
> summary(m2.2.1, ddf="Kenward-Roger")
Linear mixed model fit by REML t-tests use Kenward-Roger approximations to
  degrees of freedom [lmerMod]
Formula: ACHIEVE ~ FACULTY + (1 | DEPART)
   Data: dat

REML criterion at convergence: 209

Scaled residuals:
     Min       1Q   Median       3Q      Max
-1.42052 -0.87863 -0.02588  0.69950  1.98494

Random effects:
 Groups   Name        Variance Std.Dev.
 DEPART   (Intercept)  7.626   2.762
 Residual             20.050   4.478
Number of obs: 36, groups:  DEPART, 6

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)
(Intercept)   80.944      1.912  4.000  42.334 1.86e-06 ***
FACULTYS       3.167      2.704  4.000   1.171    0.307
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr)
FACULTYS -0.707

とりあえず、学科を変量効果として学部を固定効果としたモデル。推定値は、固定効果の検定まで書籍と一致。ここでは、固定効果の検定のためlmerTestパッケージを使っている。この時、summaryを出力する際に「ddf=“Kenward-Roger”」オプションを与えて、Kenward-Roger法で固定効果の自由度の推定を行っている。ddfを指定しなければ、Satterthwaiteの近似が使われるようだ。ただし、ここでのモデルではいずれを用いても結果は同一。SPSSではどうなっているのだろう（今度調べよう）。