[R][メモ] 卒論で使った分析

[R][メモ] 卒論で使った分析

卒論のためにRを使って行った分析をメモ。

やったことは、

  1. 3要因分散分析 (被験者間計画)
  2. 検出された交互作用を検定するために単純主効果検定
  3. 単純主効果の強さを測定するために効果量を計算 (イータ2乗)

統計のこともRのこともあまりよくわかっていないけど、少なくとも手元のデータではSPSSと同じ数値が出た。

分散分析

コントラスト配列の作成方法を指定するのがポイント?

> # Type III Anova
> library(car) # load Anova()
> options(contrasts = c("contr.sum", "contr.sum"))
> model <- lm(response ~ factor.A * factor.B * factor.C, data = results)
> tta   <- Anova(model, type = "III")
> tta
Anova Table (Type III tests)

Response: response
                           Sum Sq Df   F value    Pr(>F)
(Intercept)                 53179  1 1469.5196 < 2.2e-16 ***
factor.A                       39  1    1.0730   0.30354
factor.B                     2349  1   64.9153 8.563e-12 ***
factor.C                        7  1    0.2037   0.65302
factor.A:factor.B             152  1    4.1995   0.04389 *
factor.A:factor.C              43  1    1.1757   0.28166
factor.B:factor.C               1  1    0.0237   0.87809
factor.A:factor.B:factor.C      2  1    0.0452   0.83212
Residuals                    2750 76
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

単純主効果検定

残差項を指定するのがポイント? (参考サイト)

> # Simple Main Effects
> sme <- by(results, results$factor.A,
+           function(x) Anova(aov(response ~ factor.B, data = x),
+                             type = "II", error = model))
> sme
results$factor.A: Level.A
Anova Table (Type II tests)

Response: response
          Sum Sq Df F value    Pr(>F)
factor.B  1797.2  1  49.661 7.158e-10 ***
Residuals 2750.3 76
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------------------------------------------
results$factor.A: Level.B
Anova Table (Type II tests)

Response: response
           Sum Sq Df F value    Pr(>F)
factor.B   672.82  1  18.592 4.801e-05 ***
Residuals 2750.30 76
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

注: factor.A, factor.B, factor.Cの水準数は全て2なので、多重比較は無し。

イータ2乗値

分散分析のイータ2乗値は、 η2 = SSeffect / SStotal で計算できるらしい。ちなみに、偏イータ2乗値は、2 = SSeffect / (SSeffect + SSerror)らしい。

> # Effect Sizes (Eta-Squared)
> es <- lapply(sme, function(x) x$`Sum Sq` / sum(x$`Sum Sq`))
> es
$Level.A
[1] 0.3951999 0.6048001

$Level.B
[1] 0.1965528 0.8034472

まとめ

  • 単純主効果検定は面倒
  • 分散分析の結果からイータ2乗値を出すのは簡単っぽい
  • RのコードとGoogle Code Prettifyは相性がいい
  • Tsukuba.R#4に行けなくなって残念><
スポンサーサイト

関連記事

トラックバック URL

http://liosk.blog103.fc2.com/tb.php/177-58653ba3

トラックバック

コメント

コメントの投稿

お名前
コメント
編集キー