glm 関数と分散分析

分散分析は、観測値を従属変数とし、各グループの平均をパラメーターとし、各グループのグループ名をデザイン行列にすることによって、一般化線形モデルとしてあらわせる。

ここで、サンプルデータとして、5 つのグループ(At, S5, S4, S3, S2)の比較データを利用する。

dat <- read.table("https://stat.biopapyrus.net/data/tseed.csv", header = T, sep = ",")
head(dat)
##     At   S5   S4   S3   S2
## 1 1.29 2.25 3.91 4.59 5.58
## 2 0.82 2.31 4.24 4.95 4.26
## 3 0.87 2.40 2.81 4.78 4.19
## 4 1.12 2.49 4.25 4.29 5.78
## 5 1.09 2.58 3.75 5.11 4.10
## 6 0.99 2.28 3.49 3.51 5.20

一般化線形モデルとして解くためには、次のようなモデルを構築する必要がある。

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

上で見たデータをこのモデルにフィットさせるためには、以下のように変形する必要がある。すなわち、At グループの真の平均値を β0 とおくと E[YAt] = β0 になる。次に、S5 の真の平均値を At の真の平均値の差を β1 と置けば、E[YS5] = β0 + β1 となる。以下同様に β2, β3, β4 を考える。

\[ E \begin{bmatrix} y_{At} \\ y_{At} \\ \vdots \\ y_{S5} \\ \vdots \\ \vdots \\ y_{S2} \end{bmatrix} = \begin{bmatrix} 1&0& 0 & 0 & 0 \\ 1&0&0&0&0 \\ 1&\vdots & \vdots & \vdots & \vdots \\ 1&1&0&0&0 \\ 1& \vdots & \vdots & \vdots & \vdots \\ 1& \vdots & \vdots & \vdots & \vdots \\ 1&0&0&0&1\end{bmatrix} \begin{bmatrix} \beta_{0} \\ \beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{bmatrix} \]

R で読み込んだ行列型のデータを、このモデルと同じ配置となるように変形する。その際に reshape2 ライブラリーを利用すると便利。デザイン行列を作成する際に、モデル式では 5 列となっているが、At グループの平均を基準(定数)として考えれば、デザイン行列は 4 列で十分である。そこで、R ではモデル式中にある全部 1 からなる列を取り除いたものをデザイン行列として作成する。

library(reshape2)
y <- melt(dat)
y <- y[!is.na(y[, 2]), ]         # 欠損値 NA を除く
head(y)
##   variable value
## 1       At  1.29
## 2       At  0.82
## 3       At  0.87
## 4       At  1.12
## 5       At  1.09
## 6       At  0.99

X <- as.data.frame(matrix(0, ncol = 4, nrow = nrow(y)))
X[y[, 1] == "S5", 1] <- 1
X[y[, 1] == "S4", 2] <- 1
X[y[, 1] == "S3", 3] <- 1
X[y[, 1] == "S2", 4] <- 1

# デザイン行列中の 0 と 1 を因子型に変換する
# そうしないと 0 と 1 がグループラベルと識別されずに、回帰分析と同じ効果になってしまう
for (i in 1:4) X[, i] <- as.factor(X[, i])

あとで解析しやすいように観測値とデザイン行列を束ねる。

dat.df <- data.frame(y[, 2], X)
colnames(dat.df) <- c("y", "S5", "S4", "S3", "S2")
head(dat.df)
##         y S5 S4 S3 S2
## [1,] 1.29  0  0  0  0
## [2,] 0.82  0  0  0  0
## [3,] 0.87  0  0  0  0
## [4,] 1.12  0  0  0  0
## [5,] 1.09  0  0  0  0
## [6,] 0.99  0  0  0  0

glm 関数を利用してモデルの解析を行う。Coefficinets 項の Pr 列に示されている p 値は各係数(β1, β2, β3, β4)がゼロかどうかの検定結果を示している。。この場合、4 つとも p 値が非常に小さいため、各係数がゼロでないことになる。すなわち、S5, S4, S3, S2 のグループの平均値が基準(すなわち At グループの平均値)との間に有意差が認められることになる。

v <- glm(y ~ S5 + S4 + S3 + S2, data = dat.df)
summary(v)
## Call:
## glm(formula = y ~ S5 + S4 + S3 + S2, data = dat.df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.32333  -0.23231  -0.01077   0.31795   1.45308  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.0438     0.1714   6.092 9.62e-08 ***
## S51           1.4769     0.2423   6.095 9.51e-08 ***
## S41           2.3295     0.2473   9.418 2.72e-13 ***
## S31           3.5131     0.2423  14.497  < 2e-16 ***
## S21           3.5670     0.2473  14.422  < 2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for gaussian family taken to be 0.3817213)
## 
##     Null deviance: 136.76  on 62  degrees of freedom
## Residual deviance:  22.14  on 58  degrees of freedom
## AIC: 124.9
## 
## Number of Fisher Scoring iterations: 2

aov 関数を利用した分散分析

aov(value ~ factor(variable), data = y)
## Call:
##    aov(formula = value ~ factor(variable), data = y)
## 
## Terms:
##                 factor(variable) Residuals
## Sum of Squares         114.61596  22.13984
## Deg. of Freedom                4        58
## 
## Residual standard error: 0.617836
## Estimated effects may be unbalanced