僕が作成した架空の例題です。リハビリテーションの実施期間で運動機能に差はあるのか？ただし、効果判定は別問題ですので、ここでは議論しません。
データプレパレーション
例題 この研究には5施設（施設ID=A, B, C, D, E）の回復期病棟が参加した。乱塊法（Random Block method）に従い、理学療法を１週間受けている患者１名、２週間受けけている患者１名、３週間受けている患者１名を５施設からそれぞれランダムに抽出し（n=15）、TUG（time up & go test）を測定して各期間の差を検証した（患者情報は除外）。
使用するパッケージ
library(tidyr)
library(ggplot2)
以下のファイルを読み込んでください（ファイルの読み込み方）
(dat <- read.csv(
"anova02.csv",
header=T,
fileEncoding = "Shift-JIS"
))
rowMeans(dat[,-1]) #行平均(各患者の平均)
colMeans(dat[,-1]) #列平均(各期間での平均)
PT実施期間を要因A、施設を要因Bとします
- 要因Aに着目するので、要因Bの施設水準はブロック因子となります
- １～３週間のことを要因Aの水準、施設A～Dのことを要因Bの水準と言うこともあります
- 繰り返し測定なし＝1つの水準に含まれるデータが1つ
縦のデータセットに変換（縦持ち）
dat2 <- tidyr::gather(
dat,
key=point,
value=data,
-ID)
head(dat2)
全体の平均
mean(dat2$data)
名義変数はfactorに変更
str(dat2)
pointはcharacter変数ですが、解析する場合にはレベルが付与されるfactor変数に変換しておきましょう（character、factorはともに名義変数です）
dat2$ID <- as.factor(dat2$ID)
dat2$point <- as.factor(dat2$point)
str(dat2)
- IDとpointはFactor変数になり、レベルも付いてます
- 基準レベルは、ID＝1、point=PT1wとなります（先頭の文字）
グラフ
#グラフの描き方
g1 <- ggplot2::ggplot(
dat2,
aes(x = point, y = data)
)
g2 <- g1 + geom_point(
aes(x = point, y = data, color=ID),
size = 3) +
theme_test()
y <- mean(dat2$data)
(g3 <- g2 + geom_hline(
aes(
yintercept = y
),
color = "red"
))
対応のある一元配置分散分析、1要因ランダムブロックデザイン
通常行われている二元配置分散分析は、要因Aか要因Bのどちらかの効果を検証します。ここでは、理学療法の実施期間（要因A）が解析の対象となります。施設要因は、効果を検証することが目的ではなく、誤差を減らすために取り込みます。そのような場合には、施設A、B、C、D、Eはブロック因子となります。
もし施設A、B、C、D、Eを無視して、以下のように要因Bを誤差と見なせば、繰り返し測定と同じ扱いになり、一元配置分散分析の適応となります
しかし今回の例題では、PT実施期間（1週間、２週間、３週間）である要因Aに加えて要因Bも結果に影響します
このような場合には、要因が2つあるので二元配置分散分析となります
したがって 全体の誤差（ばらつき）＝ 要因Aの効果 ＋ 要因Bの効果 ＋ 誤差 となります
要因Aの効果、要因Bの効果は主効果と呼ばれます
回帰分析（カテゴリカルデータ）
要因A間の差：PT1w(平均13.66)を基準として他の期間との差を求めてみます
- PT2wの平均ーPT1wの平均\(=10.44-13.66=-3.22\)
- PT3wの平均ーPT1wの平均\(=9.22-13.66=-4.44\)
要因B間の差：施設A(平均14.1)を基準として他の施設との差を求めてみます
- 施設Bの平均ー施設Aの平均\(=10.43-14.1=-3.67\)
- 施設Cの平均ー施設Aの平均\(=9.76-14.1=-4.34\)
- 施設Dの平均ー施設Aの平均\(=11.4-14.1=-2.7\)
- 施設Eの平均ー施設Aの平均\(=9.83-14.1=-4.27\)
２変数の重回帰分析を実行（ただし両変数とも名義変数）
fit <- lm(data ~ point + ID, data=dat2)
summary(fit)
説明変数がカテゴリカルデータの場合には解釈に気を付けましょう
- 回帰分析の係数の推定値は、上で求めた要因A間の差、要因B間の差になります
- ここのF値は要因Aと要因Bの平方和の総和（自由度=2+4）から算出されています
- p<0.05となり両要因の効果があることが示されました
分散分析表
回帰分析の結果をanova関数に渡すことで分散分析表が算出されます
anova(fit)
- このF値は要因Aと要因Bそれぞれの平均平方和から算出されます
- pointのp値<0.05なので、PTの実施期間の効果があることが言えます
- 要因Bはブロック因子となるので、分散分析表には以下のように要因Aの結果のみを記載します（当たり前ですが、結果は二元配置分散分析と全く同じものです）
統計モデルと仮説
\( y_{ij}-\bar{y}=(\bar{y_{i\cdot}}-\bar{y})+(\bar{y_{\cdot j}}-\bar{y})+\)残差
残差=\(y_{ij}+\bar{y}-\bar{y_{i\cdot}}-\bar{y_{\cdot j}}\)
\( y_{ij}-\bar{y}=(\bar{y_{i\cdot}}-\bar{y})+(\bar{y_{\cdot j}}-\bar{y})+(y_{ij}+\bar{y}-\bar{y_{i\cdot}}-\bar{y_{\cdot j}})\)
二元配置分散分析の統計モデル
\(y_{ij}=\mu + \alpha_{i} + \beta_{j} + e_{ij}\) \(\quad (i=1,2,3, \quad j=1,2,3,4,5)\)
\(\alpha_{i}：\)要因AのPT期間 \(i\) の効果
\(\beta_{j}：\)要因Bのブロック因子 \(j\) の効果
仮説
要因Aの効果を検証することが目的となるので、仮説は以下のようになります
帰無仮説A：\(\alpha_{1}=\alpha_{2}=\alpha_{3}=0\)
対立仮説：要因Aに効果がある
理学療法期間とTUGの関連性を検証するためには多重比較が必要になります。ただし因果関係を証明したいのであれば、単に比較するだけでは不十分です。
残差
理学療法の実施期間にともなうTUGの変化が施設間で異なる（要因A×要因Bの交互作用がある）場合でも、その効果は残差 \(e_{ij}\)に含まれることになります
繰り返し測定のある二元配置分散分析のように交互作用と誤差項の分離ができないので、ここでは要因A×要因Bの交互作用を誤差とみなしています。つまり繰り返しのない二元配置分散分析では、たとえ要因A×要因Bの交互作用があったとしても誤差として捉えて、研究目的である要因Aの主効果を評価しているのです。
平方和の分解
一元配置分散分析と同様の方法で平方和を分解します（分かりやすくするために例題の数値を当てはめて記載します）
全体平方和\(\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{5} (y_{ij}-\bar{y})^2\)
要因Aの平方和\(\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{5} (\bar{y_{i\cdot}}-\bar{y})^2=5\sum_{i = 1}^{3} (\bar{y_{i\cdot}}-\bar{y})^2\)
自由度：３－１＝２
要因Bの平方和\(\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{5} (\bar{y_{\cdot j}}-\bar{y})^2=3\sum_{j = 1}^{5} (\bar{y_{\cdot j}}-\bar{y})^2\)
自由度：５－１＝４
残差平方和\(\displaystyle= \sum_{i = 1}^{3} \sum_{j = 1}^{5} (y_{ij}+\bar{y}-\bar{y_{i\cdot}}-\bar{y_{\cdot j}})^2\)
自由度：（３－１）×（５－１）＝８
全体平方和 = 要因Aの平方和 + 要因Bの平方和 + 残差平方和
証明
誤差の合計は0なので
\(\sum e_{ij}=0\)
\(\bar{y_{i\cdot}}-\bar{y}=(y_{ij}-\bar{y})-(y_{ij}-\bar{y_{i\cdot}})\)
\(\bar{y_{\cdot j}}-\bar{y}=(y_{ij}-\bar{y})-(y_{ij}-\bar{y_{\cdot j}})\)より
\(\sum \sum \alpha_{i} \beta_{j}=0\)
以上のことから
\(\sum(y_{ij}-\mu)^2= \sum(\alpha_{i})^2 + \sum(\beta_{j})^2 + \sum(e_{ij})^2\)
検定
平均平方和
勉強のために要因Bの平均平方和も求めてみましょう
再掲
分散分析表の値を確認してみましょう
#全体の平均
ybar <- mean(dat2$data)
#A群平方和
(Agun <- sum(
(mean(dat2[dat2$point=="PT1w", "data"]) - ybar)^2*5 +
(mean(dat2[dat2$point=="PT2w", "data"]) - ybar)^2*5 +
(mean(dat2[dat2$point=="PT3w", "data"]) - ybar)^2*5
))
#B群平方和（同じ計算の繰り返しはfor関数が便利です）
vec1 <- c()
vec2 <- c()
for(a in c("A","B","C","D","E")){
p <- (mean(dat2[dat2$ID==a, "data"]) - ybar)^2*3
vec1 <- c(p)
vec2 <- c(vec2, vec1)
}
(Bgun <- sum(vec2))
#残差平均平方和(全体平方和ーA群平方和ーB群平方和 )
(Zan <- sum((dat2$data - ybar)^2) -
Agun -
Bgun
)
Rの答えと同じになりました
F値
あとは簡単です
平均平方和＝平方和÷自由度
F値＝群平均平方和÷残差平均平方和
要因AのF値\(=\dfrac{Agun/2}{Zan/8}\)
要因BのF値\(=\dfrac{Bgun/4}{Zan/8}\)
#要因AのF値
(Fa <- (Agun/2)/(Zan/8))
#要因BのF値
(Fb <- (Bgun/4)/(Zan/8))
F検定
要因A（PT週）の効果
自由度(2, 8)のF分布で検定します
pf(Fa, 2, 8, lower.tail=F)
コメント