文書の過去の版を表示しています。
1要因の分散分析 対応あり
対応ありの1要因分散分析。繰り返しあり、被験者内要因ともいわれる。同じ人について、複数の計測値があるような場合。
データの準備
ID | テストA | テストB | テストC |
---|---|---|---|
1 | 70 | 80 | 68 |
2 | 75 | 85 | 66 |
3 | 68 | 78 | 68 |
4 | 71 | 79 | 69 |
5 | 72 | 81 | 70 |
この場合も、このままデータを入力すると分析できないので(SPSSではできる、、がかえって理解しにくくややこしくしている要因のようにも思える)、ID(個人番号)、テストの種類、スコアの変数に再編する。
ID | スコア | テスト |
---|---|---|
1 | 70 | A |
2 | 75 | A |
… | ||
1 | 80 | B |
2 | 85 | B |
ID+テストが主キーになるようなテーブルを作るイメージ。IDは、数値ではなくカテゴリ変数にしておく。
> ID <- factor(rep(1:5, 3)); > SCORE <- c(70,75,68,71,72,80,85,78,79,81,68,66,68,69,70); > TEST <- c(rep("A",5), rep("B",5), rep("C",5)); > df <- data.frame(ID=ID, SCORE=SCORE, TEST=TEST); > df ID SCORE TEST 1 1 70 A 2 2 75 A 3 3 68 A 4 4 71 A 5 5 72 A 6 1 80 B 7 2 85 B 8 3 78 B 9 4 79 B 10 5 81 B 11 1 68 C 12 2 66 C 13 3 68 C 14 4 69 C 15 5 70 C
では、このデータで対応ありの分散分析。いろいろやり方がある。
> res.anova <- lm(SCORE ~ TEST + ID, df); > anova(res.anova); Analysis of Variance Table Response: SCORE Df Sum Sq Mean Sq F value Pr(>F) TEST 2 418.53 209.267 46.3321 3.989e-05 *** ID 4 28.67 7.167 1.5867 0.2678 Residuals 8 36.13 4.517
まずは、lm()関数の結果をanova()関数に渡す方法。
ここでのポイントは、モデル式に「ID」を含めること。個人による要因を統制した分散分析ということがわかる。結果を用いる際には、IDの推定結果を無視すればよい。
別の方法として、aov()関数。
> res.anova <- aov(SCORE ~ TEST + Error(ID/TEST), df); > summary(res.anova); Error: ID Df Sum Sq Mean Sq F value Pr(>F) Residuals 4 28.67 7.167 Error: ID:TEST Df Sum Sq Mean Sq F value Pr(>F) TEST 2 418.5 209.27 46.33 3.99e-05 *** Residuals 8 36.1 4.52
モデル式にて、個人の要因をError()で表現する。結果は同じだが、ID行は出力されない。
分散分析だけならよいが、Mauchlyの球面性検定を別に行う必要がある。この手順が結構面倒。
まず、データ形式を縦長形式から横長形式(上記の上のテーブル)に変換する。
> df_w <- reshape(df, timevar="TEST", idvar=c("ID"), direction="wide"); > df_w ID SCORE.A SCORE.B SCORE.C 1 1 70 80 68 2 2 75 85 66 3 3 68 78 68 4 4 71 79 69 5 5 72 81 70
reshape()関数が結構便利。timevarには横に展開する変数、idvarには個人を特定する変数を指定する。
次に、各列をベクトルとして取り出して結合して、行列を作成する。
> df_m <- cbind(df_w$SCORE.A, df_w$SCORE.B, df_w$SCORE.C); > colnames(df_m) <- c("A","B","C"); > df_m A B C [1,] 70 80 68 [2,] 75 85 66 [3,] 68 78 68 [4,] 71 79 69 [5,] 72 81 70
次に、行列をlm()関数で処理する。
> res.mlm <- lm(df_m ~ 1); > res.mlm Call: lm(formula = df_m ~ 1) Coefficients: A B C (Intercept) 71.2 80.6 68.2
きちんと調べられていないが、行列を目的変数にlm()関数を適用すると、各列をベクトルとして線形モデル(回帰)の推定が行われるようだ。そして、返されるオブジェクトは普通のlm()のオブジェクトが複数結合したような形になるみたい。推定の内容は、ここでは切片のみのモデルとなっているので、各列の平均値を推定していることになるはず。実際に結果をみても、そうなっている。
Mauchlyの球面性検定を単体で行う関数もあるが、carパッケージのAnova()関数の方がεも同時に計算してくれたり、分散分析も一気にやってくれるので便利。以下でTEST.factorに水準数を設定したベクトルを作ったりしているが、これの意味はまだ調べていない。とりあえず、適当に3要因のベクトルを作成しておく。
> TEST.factor <- factor(c("A","B","C")); > res.anova <- Anova(res.mlm, idata=data.frame(TEST.factor), idesign=~TEST.factor, type=3); > summary(res.anova, multivariate=F); Univariate Type III Repeated-Measures ANOVA Assuming Sphericity SS num Df Error SS den Df F Pr(>F) (Intercept) 80667 1 28.667 4 11255.814 4.733e-08 *** TEST.factor 419 2 36.133 8 46.332 3.989e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Mauchly Tests for Sphericity Test statistic p-value TEST.factor 0.14767 0.056746 Greenhouse-Geisser and Huynh-Feldt Corrections for Departure from Sphericity GG eps Pr(>F[GG]) TEST.factor 0.53986 0.001744 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 HF eps Pr(>F[HF]) TEST.factor 0.5818973 0.001229061
Anova()関数のidataには、先程のTEST.factorをデータフレームにしたものを指定し、idesignにはTEST.factorを指定する。idesignは「=」ではなく「=~」である点に注意。平方和のタイプを3にする場合は、contrastsをcontr.sumに変更しておく必要があるかもしれない。ここでは1要因なので、結果には影響しないので特に設定していない。
出力の下半分で、Mauchlyの球面性検定が行われている。MauchlyのWは「0.14767」、p値は「0.056746」となり、SPSSと一致。イプシロンのGreenhouse-Geisserは「GG eps」の値で「0.53986」、Huynh-Feldtは「HF eps」の値で「0.5818973」。
Mauchlyの球面性検定が有意なら、球面性が仮定されないので、自由度の補正が必要なようだ。ここでは、分散分析表のSSの「419」がSPSSでいう「タイプIII平方和」で要因の平方和、その右側のError SSが誤差の平方和で、自由度がそれぞれ2,8となっている。球面性が仮定されなければ、自由度に「GG eps」および「HF eps」をかけて自由度を調整し検定をすればよい。オブジェクトから直接値を抜き出す方法がまだわからないので、いまのところ手計算。とりあえずSPSSと同じ結果が得られたので今のところはOKとする。あまり勉強もしないままソースまで追いかけて、半日かかって単純なことに気がつく。値を抜き出す方法が分かった。
summary関数は結果を表示するものだと思っていたが、これをオブジェクトに入れておくと後で個別の値にアクセスできる。また、GG eps、HF epsで自由度を調整した有意確率についても、実はGG、HFの値の後にすでに計算済みの値が表示されていた。
> summary.anova <- summary(res.anova); #要因の平方和 > summary.anova$univariate.tests[2,1] #要因の自由度 > summary.anova$univariate.tests[2,2] #誤差の平方和 > summary.anova$univariate.tests[2,3] #誤差の自由度 > summary.anova$univariate.tests[2,4] #Greenhouse-Geisserのイプシロン > summary.anova$pval.adjustments[1] #Greenhouse-Geisserのイプシロンによって自由度を調整した後の有意確率(計算しなくてもここを見れば良かった) > summary.anova$pval.adjustments[2] #Huynh-Feldtのイプシロン > summary.anova$pval.adjustments[3] #Huynh-Feldtのイプシロンによって自由度を調整した後の有意確率(計算しなくてもここを見れば良かった) > summary.anova$pval.adjustments[4]
遠回りしたが、Rについて幾つかわかったことがあるし、球面性の仮定が成り立たないときの有意確率の計算方法もわかったので、無駄ではなかったということにしよう。
最後に、多重比較。ここではボンフェローニの方法で検討されているので、以下のようになる。
> pairwise.t.test(df$SCORE, df$TEST, paired=T, p.adjust.method="bonferroni"); Pairwise comparisons using paired t tests data: df$SCORE and df$TEST A B B 5.8e-05 - C 0.3746 0.0055 P value adjustment method: bonferroni
最初の引数は平均を求めるデータ、次がグループ、pairedは対応のある場合、p.adjust.methodに有意確率の調整方法としてボンフェローニの方法。