一元配置分散分析 - 統計学備忘録一元配置分散分析

例題
歩行が自立している高齢者のバランストレーニングの効果を確認するために、30名を対象に床反力計を使用して足圧中心(COP)の総軌跡長を30秒間評価した。年齢別に差があるか有意水準5％で検定してみましょう。

僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。

データの準備と要約
回帰分析と分散分析
多重比較
- 分散分析表の理解
- 平方和の分解
分散分析の検定方法

データの準備と要約

使用するパッケージ

library(psych)
library(multcomp)

以下のファイルを読み込んでください（ファイルの読み込み方）

anova-01 データのダウンロード

dat <- read.csv("anova-01.csv", header=T, fileEncoding = "UTF-8")

縦のデータセットです

head(dat)

グラフ(赤い線は全体の平均値です)

グラフの描き方はここを参照ください

各グループの要約

psych::describeBy(dat$data, dat$age)

情報が多すぎるので・・・

dat2 <- psych::describeBy(dat$data, dat$age)
lapply(
    dat2,
    "[",
    c("n", "mean", "sd", "median", "min", "max", "se")
)

これくらいでいいでしょう

IDとageはカテゴリー変数に変更しましょう（重要）

dat$ID <- as.factor(dat$ID)
dat$age <- as.factor(dat$age)

回帰分析と分散分析

回帰分析

分散分析表の前に単回帰分析をやります（ただしageは0 or 1の2値変数となります）

fit <- lm(data~age, data=dat)
summary(fit)

lm(data~age, data=dat)

Rで単回帰分析を実行するときのプログラムです。y＝data, x=ageです。ただしここではageが名義変数ですので解釈には注意が必要です。

この結果から以下のことが分かります

age65-70の平均=1092
age70-75の平均=1092+83.9=1175.9
age75-80の平均=1092+234=1326

また一番下の行が分散分析の結果です

F検定のｐ値が<0.05なので３群のどこかに有意差があることが分かります

回帰分析でもう少し解説する予定です

分散分析

回帰分析の結果をanova関数に渡すことで分散分析表が出力されます（自由度、F値、p値は同じです）

anova(fit)

分散分析表

F値は回帰分析の結果と同じです

群間に有意な差があることが分かりました

多重比較

多重比較の方法のみ紹介します

Dunnet法

この検定はコントロール群と他の群との比較になります

ano <- aov(data ~ age, data = dat)
summary(multcomp::glht(ano, linfct=mcp(age = "Dunnett")))

65-70 vs 75-80 に有意差があります

これは、65-70群がコントロール群になっています

75-80群をコントロールに置きたいときは

dat$age <- relevel(dat$age, ref="75-80") 
ano <- aov(data ~ age, data = dat)
summary(multcomp::glht(ano, linfct=mcp(age = "Dunnett")))

65-70 vs 75-80 と 70-75 vs 75-80 に有意差があります

Tukey‒Kramer法

TukeyHSD(aov(data ~ age, data = dat))

65-70 vs 75-80, 70-75 vs 75-80 に有意差があります

Bonferroni法

pairwise.t.test(dat$data, dat$age, p.adj = "bonf")

同じく65-70 vs 75-80, 70-75 vs 75-80 に有意差があります

結果のみでよければ、ここまででOKです

ここからはもう少し詳しく勉強してみましょう

分散分析表の理解

データサイエンス入門ですので、分散分析表をもう少し詳しく見てみましょう

平方和 (square sum, sum of squares)について

平方和（偏差$^2$の和）はテキストによって様々な書き方がありますが、ここでは以下の名称を使用します

全体平方和
群平方和
群内平方和（残差平方和）

下記の図のなかのそれぞれの偏差$^2$を足し合わせた合計になります

65-70歳をA群、70-75歳をB群、75-80歳をC群とします

A群、B群、C群は独立したサンプルです（ここ重要です）

$i=1,\, 2,\, 3$

$i=1=A$
$i=2=B$
$i=3=C$

なんだかややこしいのですが数式を組み立てるための準備です…最初から１群、２群、３群とすればこの手間は除けますがやや意味合いが異なってきます

$j=1,\cdots,n_i$ （$n_1=10, \, n_2=10, \, n_3=10$）

データ：$y_{ij}$ （例えば $y_{ij}=y_{27}$ はB群の7番目となります）

各群の平均：$\bar{y_i}$

全体の平均：$\bar{y}=1197.967$

mean(dat$data)

群平方和の偏差は各群に1つです。例えば、A群の群平方和＝（A群の偏差$^2$ ✕ n1）となります。数式になったら分かりずらいのですが、群平方和＝（A群の偏差$^2$ ✕ n1）＋（B群の偏差$^2$ ✕ n2）+ （C群の偏差$^2$ ✕ n3）となります。

全体平方和$\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{n_i} (y_{ij}-\bar{y})^2$

（自由度：$n-1=29$）

群平方和$\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{n_i} (\bar{y_{i}}-\bar{y})^2=\sum_{i = 1}^{3} n_i(\bar{y_{i}}-\bar{y})^2$

（自由度：$群の数-1=2$）

群内平方和（残差平方和）$\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{n_i} (y_{ij}-\bar{y_i})^2$

（自由度：$n-群の数=27$）

分散分析では群平方和と群内平方和を使います。分散分析表には、全体平方和の出番はありません。でも、全体平方和＝群平方和＋群内平方和が成り立つことから、それぞれの平均平方和の比を検定して群間差を検証することになるので、全体平方和を改めて記載する場合もあります。

平方和の分解

全体平方和＝群平方和＋群内平方和を証明しておきましょう

ここで $\bar{y_i}-\bar{y}=\alpha_i$ を群の効果、$y_{ij}-\bar{y_i}=e_{ij}$ を残りの効果（残差）として考えます

$\blacksquare$ 一元配置分散分析の統計モデル

${y_{ij} = \bar{y}+\alpha_i+e_{ij}}$

グラフから 全体の偏差 = 群の偏差 + 群内の偏差 ということが理解できます

数式に直すと以下のようになります

$(y_{ij}-\bar{y}) = (\bar{y_i}-\bar{y}) + (y_{ij}-\bar{y_i})$

統計モデルを書き直すと、こうなります

$(y_{ij}-\bar{y}) = \alpha_i+e_{ij}$

両辺の二乗和を求めてみましょう(例題のi=1,2,3を適用してます)

$\displaystyle \sum_{i=1}^{3}\sum_{j=1}^{n_i}(y_{ij}-\bar{y})^2$

$\displaystyle =\sum_{i=1}^{3}\sum_{j=1}^{n_i}(\alpha_i+e_{ij})^2$

$\displaystyle =\sum_{i=1}^{3}\sum_{j=1}^{n_i}(\alpha_i^2+2\alpha_ie_{ij}+e_{ij}^2)$

ここで残差の和は0になることから

（$\displaystyle \sum_{i=1}^{3}\sum_{j=1}^{n_i}e_{ij}=0)$

$\displaystyle \sum_{i=1}^{3}\sum_{j=1}^{n_i}(y_{ij}-\bar{y})^2=\sum_{i=1}^{3}\sum_{j=1}^{n_i}(\alpha_i^2+e_{ij}^2)$

つまり

$\displaystyle \sum_{i=1}^{3}\sum_{j=1}^{n_i}(y_{ij}-\bar{y})^2= \sum_{i = 1}^{3} n_i(\bar{y_{i}}-\bar{y})^2+\sum_{i = 1}^{3} \sum_{j = 1}^{n_i} (y_{ij}-\bar{y_i})^2$

証明終わり

分散分析の検定方法

帰無仮説：$\alpha_1=\alpha_2=\alpha_3$

平方和から平均平方和を算出し、群平均平方和と残差平均平方和の比がF分布に従うことを利用して検定します

統計量：$F＝\dfrac{群の平均平方和}{群内の平均平方和（残差平均平方和）}$

$\displaystyle F \sim F_{群の平方和の自由度、群内の平方和の自由度}$

例題のF値は $F_{2、27}$ に従います

$F > F_{2, 27}(0.95)$ となる場合には帰無仮説が棄却される検定です

群の平均平方和=$\dfrac{群の平方和}{自由度}=\dfrac{群の平方和}{2}$

群内の平均平方和（残差平均平方和）=$\dfrac{群内の平方和}{自由度}=\dfrac{群内の平方和}{27}$

群の平均平方和(gun)

gun <- sum(
    (mean(dat[dat$age=="65-70", "data"]) - mean(dat$data))^2*10 +
    (mean(dat[dat$age=="70-75", "data"]) - mean(dat$data))^2*10 +
    (mean(dat[dat$age=="75-80", "data"]) - mean(dat$data))^2*10
)/2

dat[dat$age==”65-70″, “data”]

これはdatというデータフレームのage列が65-70に該当する対象者のdata列の値を抜き出したいときに使います。便利ですよ！

残差平均平方和(zan)

zan <- sum(
    (dat[dat$age=="65-70", "data"] - mean(dat[dat$age=="65-70", "data"]))^2 +
    (dat[dat$age=="70-75", "data"] - mean(dat[dat$age=="70-75", "data"]))^2 +
    (dat[dat$age=="75-80", "data"] - mean(dat[dat$age=="75-80", "data"]))^2
)/27