二元配置分散分析（繰り返し測定あり）

例題
地域在住の60代と80代の女性を対象として、自主的に毎日体操しているA群、運動指導を受けているB群、運動指導を受けてかつ毎日自主的に体操しているAB群の肩関節屈曲可動域を計測します。肩関節屈曲可動域に年齢や運動により可動域の差はあるか、また年齢と運動の交互作用の効果について有意水準5％で検証します。

僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。平方和の分解などは、以下の二元配置分散分析（繰り返し測定なし）のページをご参照ください。

データの準備と要約
縦のデータセットに変換（縦持ち）
分散分析表
各層の平均値
統計モデル
平方和の分解
分散分析表

データの準備と要約

使用するパッケージ（パッケージのインストール）

library(psych)
library(tidyr)
library(ggplot2)
library(dplyr)

データセット anova03 をdatに格納します（ファイルの読み込み）

anova03

要因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 
)
#もっとスマートな方法がみつかったら更新します