例題 地域在住の60代と80代の女性を対象として、自主的に毎日体操しているA群、運動指導を受けているB群、運動指導を受けてかつ毎日自主的に体操しているAB群の肩関節屈曲可動域を計測します。肩関節屈曲可動域に年齢や運動により可動域の差はあるか、また年齢と運動の交互作用の効果について有意水準5％で検証します。
僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。平方和の分解などは、以下の二元配置分散分析（繰り返し測定なし）のページをご参照ください。
データの準備と要約
使用するパッケージ（パッケージのインストール）
library(psych)
library(tidyr)
library(ggplot2)
library(dplyr)
データセット anova03 をdatに格納します（ファイルの読み込み）
要因A：method(A, B, AB),
要因B：age(60代, 80代)
帰無仮説1：method（要因A）によってrangeの平均値は変動しない
帰無仮説2：age（要因B）によってrangeの平均値は変動しない
帰無仮説3：methodとageの間に交互作用はない
サンプルは以下のような横持ちのデータセットになっています
縦のデータセットに変換（縦持ち）
dat2 <- tidyr::gather(
dat,
key=method,
value=range,
-age)
head(dat2)
重なりが分かるように少しずらしてます
#グラフの描き方
g1 <- ggplot2::ggplot(
dat2,
aes(x = method, y = range)
)
g2 <- g1 + geom_jitter(
height=0, width =0.15, size = 3, aes(colour = age)
) + theme_test()
g2 + theme(axis.text.x = element_text(size = 12))
見にくいので、A → B → AB の順番に並び変えましょう
#methodをファクター変数に変更します
dat2$method <- as.factor(dat2$method)
#これでレベルが付与されました
levels(dat2$method)
#transform関数でレベルを入れ替えましょう
dat3 <- transform(
dat2,
method = factor(
method, levels = c("A", "B", "AB")
)
)
#これでレベルの順序が変更されているか確認
levels(dat3$method)
#グラフの描き方、dat2をdat3に変更して再度描きなおします
g1 <- ggplot2::ggplot(
dat3,
aes(x = method, y = range)
)
(g2 <- g1 + geom_jitter(
height=0,
width =0.15,
size = 3,
aes(colour = age)
) + theme_test()) + theme(axis.text.x = element_text(size = 12))
これでageとmethodの関係がなんとなく分かりますね。
データを可視化することはとても大切ですね。
分散分析表
fit <- lm(range ~ method*age, data=dat2)
anova(fit)
あるmethodには効果があり、またどちらかのageに効果があることが理解できます
また交互作用もあることから、特定のmethodと特定のageで効果があるようです
結果まででよければ、ここまででOKです
ここからもう少し詳しく勉強してみましょう
各層の平均値
全体の平均と11グループの平均値が必要になります
ageで層別、methodで層別
mean(dat3$range)
psych::describeBy(dat3$range,
group = list(age = dat3$age)
)
psych::describeBy(dat3$range,
group = list(method = dat3$method)
)
不要な値も含まれていますが、縦データの要約には、この方法が一番簡単なようです。他の方法があったらまた紹介します。
age✕methodで層別（6ブロック）
ここでは平均値だけを取り出してみましょう。
dat %>%
dplyr::group_by(age) %>%
dplyr::summarise_if(
is.numeric, mean, na.rm = TRUE
)
注意）dat→横データを使用しています
縦データのdat3を使うのであれば・・・
aggregate(
range ~ method + age,
data=dat3,
FUN=mean
)
（ちょっと見にくいか・・・）
例）下のグラフには全体、60歳代、B群、60歳代Ｂ群の平均値を記載しています
統計モデル
繰り返し測定のない二元配置分散分析とは異なり、ここでは交互作用の効果も検証します
\(y_{ijk}=\mu + \alpha_{i} + \beta_{j} + \alpha_{i}*\beta_{j}+ e_{ijk}\)
\(\quad (i=1,2,3 \quad j=1,2 \quad k=1,2, \cdots, 5)\)
$\quad \alpha_{i}=\bar{y_i..}-\bar{y}=$方法 A, B, AB の効果
$\quad \beta_{j}=\bar{y_{.j.}}-\bar{y}=$年齢の効果
$\quad \alpha_{i}*\beta_{j}=\bar{y_{ij.}}-\bar{y_{i..}}-\bar{y_{.j.}}+\bar{y}=$交互作用の効果
$\quad e_{ijk}=y_{ijk}-\bar{y_{ij.}}=$残りの効果（上述した効果では説明できない効果=各年齢、各方法の組み合わせた群内の残差）
以下のように書けます
\(y_{ijk} = \bar{y} + (\bar{y_i..}-\bar{y}) + (\bar{y_{.j.}}-\bar{y})+(\bar{y_{ij.}}-\bar{y_{i..}}-\bar{y_{.j.}}+\bar{y})+(y_{ijk}-\bar{y_{ij.}})\)
\((y_{ijk}-\bar{y}) = (\bar{y_i..}-\bar{y}) + (\bar{y_{.j.}}-\bar{y})+(\bar{y_{ij.}}-\bar{y_{i..}}-\bar{y_{.j.}}+\bar{y})+(y_{ijk}-\bar{y_{ij.}})\)
平方和の分解
\(\sum\sum\sum(y_{ijk}-\bar{y})^2 = \)
\(\quad \sum\sum\sum(\bar{y_i..}-\bar{y})^2+ \)
\(\quad \sum\sum\sum(\bar{y_{.j.}}-\bar{y})^2+\)
\(\quad \sum\sum\sum(\bar{y_{ij.}}-\bar{y_{i..}}-\bar{y_{.j.}}+\bar{y})^2+\)
\(\quad \sum\sum\sum(y_{ijk}-\bar{y_{ij.}})^2\)
便宜上、上記の偏差平方和を以下のような名称にします
総平方和＝要因A平方和＋要因B平方和＋交互作用平方和＋残差平方和
例題では各グループのサンプルサイズが同じなので以下のように書けます
\(\sum\sum\sum(y_{ijk}-\bar{y})^2 = \)
\(\quad 10*\sum\sum(\bar{y_i..}-\bar{y})^2+ \)
\(\quad 15*\sum\sum(\bar{y_{.j.}}-\bar{y})^2+\)
\(\quad 5*\sum\sum(\bar{y_{ij.}}-\bar{y_{i..}}-\bar{y_{.j.}}+\bar{y})^2+\)
\(\quad \sum\sum\sum(y_{ijk}-\bar{y_{ij.}})^2\)
これで分散分析表ができます
分散分析表
再掲
3つの帰無仮説が全て棄却されました
Rの結果を確認してみましょう
m <- mean(dat3$range)
#要因A平方和
10*(149.5-m)^2 +
10*(142.5-m)^2 +
10*(152.5-m)^2
#要因B平方和
15*(161.67-m)^2 +
15*(134.67-m)^2
#交互作用平方和
5*(162-149.5-161.67+m)^2 +
5*(137-149.5-134.67+m)^2 +
5*(162-142.5-161.67+m)^2 +
5*(123-142.5-134.67+m)^2 +
5*(161-152.5-161.67+m)^2 +
5*(144-152.5-134.67+m)^2
#残差平方和
sum(
(dat3[dat3$age=="age60" & dat3$method=="A", "range"]-162)^2 +
(dat3[dat3$age=="age60" & dat3$method=="B", "range"]-162)^2 +
(dat3[dat3$age=="age60" & dat3$method=="AB", "range"]-161)^2 +
(dat3[dat3$age=="age80" & dat3$method=="A", "range"]-137)^2 +
(dat3[dat3$age=="age80" & dat3$method=="B", "range"]-123)^2 +
(dat3[dat3$age=="age80" & dat3$method=="AB", "range"]-144)^2
)
#もっとスマートな方法がみつかったら更新します
F検定については 二元配置分散分析（繰り返し測定なし）、F分布 をご参照ください
コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください