これはR言語 Advent Calendar2025 の23日目の記事です。
はじめに
今回はRで分散分析を実装するという内容を書きます。実装といってもANOVA表の平方和$(y_{ij}-y_{i.})^2$や自由度などをベタに計算して検定にかけるのではなく、もう少し一般化された枠組み(二次形式)を使って平方和や自由度を計算していきます。ベタに計算していく方法では、平方和の計算の仕方や自由度が一元配置か二元配置なのか、また交互作用を考えるか考えないのかによって大きく変わってくるので(かつあまり一般化できない)、モデルごとにスクリプトを書かないといけませんが、二次形式を使う方法ではすべてのモデルを線形モデルで統一的に扱うため、どんなモデルにも対応したスクリプトを書くことができるというメリットがあります。
この方法については(分散分析を線形モデルを使って一般化する)で理論的なことをまとめましたが、実際にその理論に基づいて計算した結果がANOVAと一致するかの確認はまだしていませんでした。そこで今回は、二次形式に基づいて計算をし、anova関数の結果と比較しながら、理論が正しことを確認していこうと思います(アドカレの締切まで意外と時間が無いので今回は一元配置だけでやります)。
理論
はじめに二次形式で分散分析を行うとはどういうことか軽く説明します。
分散分析のモデルは線形モデルであり、どんなモデルを考えていても次のように表すことができます
y=X\beta+e
ここで$\beta$は切片や各処理の効果などを含む固定効果です。この$X$と$\beta$の構造を変えることにより様々なモデルを表すことができる仕組みです。
さて、分散分析は検定の一種ですから、もちろん帰無仮説と対立仮説を考えます。特に分散分析の場合は、ある特定の処理の効果があるかないかを調べるため、「処理の効果が無い」が帰無仮説、「処理の効果がある」が対立仮説になります。ここではそれぞれに対応するモデルを次のように書くこととします。
帰無仮説: y=X_{null}\beta_{null}+e
対立仮説: y=X_{alt}\beta_{alt}+e
次に各モデルに対する残差平方和(R)というものを次のように定義します。
R=y^T(I-X(X^TX)^-X^t)y
もちろんこの残差平方和にも帰無仮説に対応するもの、対立仮説に対応するものがあり、それぞれ次のように表されます。
R_{null}=y^T(I-X_{null}(X_{null}^TX_{null})^-X_{null}^T)y
R_{alt}=y^T(I-X_{alt}(X_{alt}^TX_{alt})^-X_{alt}^T)y
また、当然ですが、帰無仮説と対立仮説に対応する残差平方和は値が異なりますので、その差分を次のよう定義します。
R_{null}-R_{alt}=R_{diff}=y^T(X_{alt}(X_{alt}^TX_{alt})^-X_{alt}^T)y-y^T(X_{null}(X_{null}^TX_{null})^-X_{null}^T)y
さて、このように定義した残差平方和、およびその差分は$y$の二次形式になっており、$y$が正規分布に従う時、残差平方和はいわゆるカイ二乗分布と呼ばれる分布に従います。そしてそのカイ二乗分布には自由度という概念が存在します。これは基底をうまく取り直した時に、どれだけの自由度が残るかという概念であり、残差平方和の真ん中にある行列$(I-X(X^TX)^-X^t)$のランクに対応しています。従って
df_{null}=rank\big(I-X_{null}(X_{null}^TX_{null})^-X_{null}^T\big)
df_{alt}=rank\big(I-X_{alt}(X_{alt}^TX_{alt})^-X_{alt}^T\big)
df_{diff}=rank\big((X_{alt}(X_{alt}^TX_{alt})^-X_{alt}^T-X_{null}(X_{null}^TX_{null})^-X_{null}^T\big)
この時、ANOVA表で言うモデル平方和は$R_{diff}$、残差平方和は$R_{alt}$、それぞれの自由度は$df_{diff},df_{alt}$に対応するため、それに従って計算をすればANOVA表と値が一致するはずです。以下ではそれを確かめていきます。
実装(PlantGrowth)
ここではRで提供されているPlantGrowthというデータセットを使って実装をしていきます。
> data( PlantGrowth )
> 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
> summary( as.factor( PlantGrowth$group ) )
ctrl trt1 trt2
10 10 10
PlantGrowthはこんな感じで一元配置のデータになっています。データ数は30あり、それぞれコントロール、処理1、処理2のいずれかに属するようです。
さて、二次形式による方法を試す前に、よく使われているanova関数で検定をしてみます。lm関数で回帰してそれの結果をanova関数に入力すれば結果が出ます。
> result <- lm( weight~group, PlantGrowth )
> anova( result )
Analysis of Variance Table
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.7663 1.8832 4.8461 0.01591 *
Residuals 27 10.4921 0.3886
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ここで、groupの行に対応するのが処理に対する効果、Regisualsは誤差に対応しています。私の考えが間違っていなければ、groupのDfは$df_{diff}$、ResidualのDfは$df_{alt}$、groupのSumSqは$R_{diff}$、RegidualのSumSqは$R_{alt}$になるはずです。これを実際に計算をしながら検証していきます。
まず、対立仮説についてですが、モデルは$y_{ij}=\mu+\alpha_i+e_{ij} $ですので、$\beta_{alt}^T=(\mu,\alpha_1,\alpha_2,\alpha_3)^T$であり、$X_{alt}$は次のように表されます。
X_{alt}= \begin{pmatrix} 1 & 1 & & \\ \vdots & \vdots & O & \\ 1 & 1 & & \\ \vdots & & \ddots & \\ 1 & & & 1\\ \vdots & & O & \vdots \\ 1 & & & 1 \end{pmatrix}
次に帰無仮説ですが、「処理の効果がない」を仮定するのでモデルは$y_{ij}=\mu+e_{ij} $で、従って$X_{null}=(1,1,...1)^T$となります。
この2つの$X$が定義できてしまえばあとは単純で、上記に記したような方法で$R_{diff},R_{alt},df_{diff},df_{alt}$を求めればよいだけです。ここにそのコードを乗せておきます。
# number of observation
nObs <- nrow( PlantGrowth )
# define y
y <- PlantGrowth$weight
# define design matrices for null and alternative hyphothesis
X_alt <-cbind( rep( 1, nObse ), model.matrix( ~ group - 1, PlantGrowth ) )
X_null <- matrix( rep(1, nObs )
# define the matrices between y^T and y
Z_alt = diag( nObs ) - X_alt %*% ginv( t( X_alt ) %*% X_alt ) %*% t( X_alt )
Z_null = diag( nObs ) - X_null %*% ginv( t( X_null ) %*% X_null ) %*% t( X_null )
Z_diff = X_alt %*% ginv( t( X_alt ) %*% X_alt ) %*% t( X_alt ) - X_null %*% ginv( t( X_null ) %*% X_null ) %*% t( X_null )
# calculate rank of Z
df_alt <- qr( Z_alt )$rank
df_null <- qr( Z_null )$rank
df_diff <- qr( Z_diff )$rank
# calculate sum of squre of residual
R_alt <- t( y ) %*% Z_alt %*% y
R_null <- t( y ) %*% Z_null %*% y
R_diff <- t( y ) %*% Z_diff %*% y
# divide by degree of free dom
MSS_diff <- R_diff / df_diff
MSS_alt <- R_alt / df_alt
# calculte p-value
p <- pf( MSS_diff / MSS_alt, df1 = df_diff, df2 = df_alt, lower.tail = FALSE )
こんな感じですかね。では実際にこのコードを回して出た結果を表示してみましょう
> df_diff
[1] 2
> df_alt
[1] 27
> MSS_diff
[,1]
[1,] 1.88317
> MSS_alt
[,1]
[1,] 0.3885959
> p
[,1]
[1,] 0.01590996
anova関数で求めたANOVA表が以下のような感じなのでこのように計算しても正しく分散分析を行うことができると分かりました。
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.7663 1.8832 4.8461 0.01591 *
Residuals 27 10.4921 0.3886
まとめ
今回は二次形式に基づいて分散分析の計算を行いました。この二次形式による計算方法は日本の分散分析教育においてあまりメジャーではありませんが、一般化する上では非常に重要です。これを機に二次形式流派が増えてくれれば幸いです。