glm 関数と分散分析 | グループ情報をデザイン行列に記述することで、分散分析を一般化線形モデルとして解くことができる

分散分析は、一般化線形モデルに帰着することができ、一般化線形モデルの問題として解くことができる。分散分析を一般化線形モデルの問題として解く際に、グループ情報を独立変数としてデザイン行列に記述し、各グループにかかる重みをパラメーター（係数）とする。また、観測値を従属変数とする。

このページでは分散分析を一般化線形モデルの問題としてとく例を示す。ここで使用するデータは、ctrl、trt1、trt2 の 3 種類の処理で処理した植物の乾燥重量を記録したデータを使用する。解析目的として、3 つの処理で、植物の乾燥重量に変化があったかどうかを調べることとする。

head(PlantGrowth)
##   weight group
## 1   4.17  ctrl
## 2   5.58  ctrl
## 3   5.18  ctrl
## 4   6.11  ctrl
## 5   4.50  ctrl
## 6   4.61  ctrl

この PlantGrowth データセットに対して、一般化線形モデルとして解くためのモデルを構築する。一般化線形モデルは次のような式で書き表すことができる。

\[ E[\mathbf{Y}] = g(\mathbf{\mu}) = \mathbf{X}\mathbf{\beta} \]

このモデル式を PlantGrowth データセットに当てはめるためには、パラメーター（回帰係数）として 2 つ用意する。まず、ctrl 群の平均値を表すパラメーターとして β₀ を用意する。次に、trt1 群の平均値と ctrl 群の平均値の差を β₁、trt2 群の平均値と ctrl 群の平均値の差を β₂ とおく。こうすることで、PlantGrowth データセットを解析するためのモデルは以下のように構築できる。

\[ E \begin{bmatrix} y_{ctrl} \\ y_{ctrl} \\ \vdots \\ y_{trt1} \\ \vdots \\ y_{trt2} \\ \vdots \end{bmatrix} = \begin{bmatrix} 1&0& 0 \\ 1&0&0 \\ \vdots & \vdots & \vdots \\ 1&1&0 \\ \vdots & \vdots & \vdots \\ 1&0&1 \\ \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} \beta_{0} \\ \beta_{1} \\ \beta_{2} \end{bmatrix} \]

この式をみると、E[y_ctrl] = β₀、E[y_trt1] = β₀ + β₁、E[y_trt2] = β₀ + β₂ の関係が得られる。これらの y に実験で得られた値を代入することでパラメーター β₀、β₁、β₂ を求めることができるようになる。

まず、ctrl、trt1、trt2 グループ情報をデザイン行列に変換する。

x <- model.matrix(~ PlantGrowth$group)
head(x)
##   (Intercept) PlantGrowth$grouptrt1 PlantGrowth$grouptrt2
## 1           1                     0                     0
## 2           1                     0                     0
## 3           1                     0                     0
## 4           1                     0                     0
## 5           1                     0                     0
## 6           1                     0                     0

次に、glm 関数でモデル解析を行ってみる。

y <- PlantGrowth$weight
m <- glm(y ~ x)
summary(m)
## Call:
## glm(formula = y ~ x)
## 
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -1.0710  -0.4180  -0.0060   0.2627   1.3690
## 
## Coefficients: (1 not defined because of singularities)
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)              5.0320     0.1971  25.527   <2e-16 ***
## x(Intercept)                 NA         NA      NA       NA
## xPlantGrowth$grouptrt1  -0.3710     0.2788  -1.331   0.1944
## xPlantGrowth$grouptrt2   0.4940     0.2788   1.772   0.0877 .
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for gaussian family taken to be 0.3885959)
## 
##     Null deviance: 14.258  on 29  degrees of freedom
## Residual deviance: 10.492  on 27  degrees of freedom
## AIC: 61.619
## 
## Number of Fisher Scoring iterations: 2

Coefficinets 項の Pr 列に示されている p 値は各係数（β₁, β₂）がゼロかどうかの検定結果を示している。この場合、危険率 5% のもとで trt1 の係数がゼロであるから、ctrl と trt1 の間に有意差が見られない。また、trt2 に関しても同様に、ctrl と trt2 の間に有意差が見られないことがわかる。ただし、このモデルの立て方として、trt1 と trt2 の差を検出することが不可能で、分散分析とまったく同じ効果が得られないことに注意。