実験デザイン別の統計的検定Rスクリプト

はじめに

実験デザインについて(TODO)

• 要因数
• 水準数
• 参加者間デザイン（対応なし）
• 参加者内デザイン（対応あり）
• 混合デザイン（対応あり x 対応なし）

スクリプトの使い方

Rscript --vanilla sample.R sample.csv
--vanilla: 前に実行したオブジェクトを破棄
--slave: Rの標準出力をしない

ソースコード

1要因2水準参加者間デザイン

例題：t_test.csv

sample_id	factor:food	q1	q2
1	banana	3	4
2	banana	3	4
3	banana	5	4
4	apple	7	7
5	apple	7	7
6	apple	6	6

プログラム：t_test.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# 独立変数に名前をつける
x = factor(data[, 2])


# 分散の等質性の検定とt検定
for ( i in 3:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- summary ---\n")
  print(by(y, x, summary))
  cat("\n--- standard Deviation ---\n")
  print(by(y, x, sd))

  # 分散の等質性の検定
  cat("\n--- 分散の等質性（ p-value < .05 の場合はt検定はできない） ---\n")
  print( var.test(y ~ x) )

  # t検定
  cat("\n--- t検定 ---\n")
  print( t.test(y ~ x, var.equal=TRUE) )
}

実行結果

$ Rscript --vanilla t_test.R t_test.csv

### Data Format ###
'data.frame':   6 obs. of  4 variables:
 $ sample_id  : int  1 2 3 4 5 6
 $ factor.food: Factor w/ 2 levels " apple"," banana": 2 2 2 1 1 1
 $ q1         : int  3 3 5 7 7 6
 $ q2         : int  4 4 4 7 7 6
NULL

（省略）

########################
q1
########################

（省略）

--- 分散の等質性（ p-value < .05 の場合はt検定はできない） ---

    F test to compare two variances

data:  y by x
F = 0.25, num df = 2, denom df = 2, p-value = 0.4
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.006410256 9.750000000
sample estimates:
ratio of variances
              0.25


--- t検定 ---

    Two Sample t-test

data:  y by x
t = 4.0249, df = 4, p-value = 0.0158
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.93056 5.06944
sample estimates:
 mean in group  apple mean in group  banana
             6.666667              3.666667

（省略）

実行結果の確認

• 分散の等質性の検定で有意差がなければ (p > .05)、t検定をできる
• t値 = 4.025
• 自由度（df）= 4
  • このとき、df が各条件の「参加者数-１」の和と等しいことを確認する
• t検定で p < .05 なら危険水準.05で有意差あり、 p < .01 なら危険水準.01で有意差あり

結果の報告（注：t や p などの統計量はイタリックにする）

• 対応のない t 検定を行なった結果, 条件Aが条件Bよりも有意に得点が高かった（ t (4) = 4.025, p = 0.016）.

1要因2水準参加者内デザイン

paired_t_test.csv

sample_id	factor:food	q1	q2
1	banana	3	4
2	banana	3	4
3	banana	5	4
1	apple	7	7
2	apple	7	7
3	apple	6	6

paired_t_test.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# factorのカラム（独立変数）に名前をつける
x = factor(data[, 2])


# t検定
for ( i in 3:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- summary ---\n")
  print(by(y, x, summary))
  cat("\n--- standard Deviation ---\n")
  print(by(y, x, sd))

  # 2群で同じ参加者なので分散の等質性の検定は不要

  # t検定
  cat("\n--- t検定 ---\n")
  print( t.test(y ~ x, paired=TRUE) )
}

実行結果を確認
- dfは「参加者数-1」

$ Rscript --vanilla paired_t_test.R paired_t_test.csv

（省略）

--- t検定 ---

    Paired t-test

data:  y by x
t = 8, df = 2, p-value = 0.01527
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.232449 4.100884
sample estimates:
mean of the differences
               2.666667

1要因3水準参加者間デザイン

one_way_between_design_anova.csv

sample_id	factor:food	q1	q2
1	banana	3	4
2	banana	3	4
3	banana	5	4
4	banana	7	7
5	peach	7	7
6	peach	6	6
7	peach	6	4
8	peach	5	6
9	apple	2	5
10	apple	3	5
11	apple	4	5
12	apple	5	2

one_way_between_design_anova.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# 参加者のカラムに名前をつける
# subject = factor(data[, 1])

# factorのカラム（独立変数）に名前をつける
x = factor(data[, 2])


# 分散の等質性の検定とt検定
for ( i in 3:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- Summary: factor ---\n")
  print(by(y, x, summary))

  # ANOVA
  cat("\n-- ANOVA ---\n")
  print(summary(aov(y ~ x)))

  # 多重検定
  cat("\n-- 多重検定 ---\n")
  print(pairwise.t.test(y, x, p.adjust.method="bonferroni"))
}

実行結果

$ Rscript --vanilla one_way_between_design_anova.R one_way_between_design_anova.csv

（省略）

-- ANOVA ---
            Df Sum Sq Mean Sq F value  Pr(>F)
x            2  28.50  14.250   10.91 0.00392 **
Residuals    9  11.75   1.306
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-- 多重検定 ---

    Pairwise comparisons using t tests with pooled SD

data:  y and x

         apple  banana
 banana 0.0037 -
 peach  0.2890 0.0637

結果の確認（Q1とQ2は別にみる）

• 全体の自由度＝データ数ー１＝１１
• 要因の自由度＝条件数ー１＝２
• 残差（誤差）の自由度＝全体の自由度ー要因の自由度＝２２

結果の報告
- F値は「要因の自由度」と「残差の自由度」と一緒に報告すること（F(2, 11)=10.91）

1要因3水準参加者内デザイン

one_way_within_design_anova.csv

sample_id	factor:food	q1	q2
1	banana	3	4
2	banana	3	4
3	banana	5	4
4	banana	7	7
1	peach	7	7
2	peach	6	6
3	peach	6	4
4	peach	5	6
1	apple	2	5
2	apple	3	5
3	apple	4	5
4	apple	5	2

one_way_between_design_anova.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# 参加者のカラムに名前をつける
subject = factor(data[, 1])

# factorのカラム（独立変数）に名前をつける
x = factor(data[, 2])


# 分散の等質性の検定とt検定
for ( i in 3:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- Summary: factor ---\n")
  print(by(y, x, summary))

  # ANOVA
  cat("\n-- ANOVA ---\n")
  print(summary(aov(y ~ x + Error(subject/x))))

  # 多重検定
  cat("\n-- 多重検定 ---\n")
  print(pairwise.t.test(y, x, p.adjust.method="bonferroni", paired=TRUE))
}

実行結果

$ Rscript --vanilla one_way_within_design_anova.R one_way_within_design_anova.csv

（省略）

-- ANOVA ---

Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  3   6.25   2.083

Error: subject:x
          Df Sum Sq Mean Sq F value  Pr(>F)
x          2   28.5  14.250   15.54 0.00423 **
Residuals  6    5.5   0.917
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-- 多重検定 ---

    Pairwise comparisons using paired t tests

data:  y and x

         apple  banana
 banana 0.066  -
 peach  0.308  0.055

結果の確認

• 全体の自由度＝データ数-１＝１１
• 要因の自由度＝条件数-１＝２
• 残差の自由度＝（参加者数－１）ｘ（条件数ー１）= ６

結果の報告
- F(2, 6)=15.54

2要因参加者間デザイン

two_way_between_design_anova.csv

sample_id	factor:food	factor:drink	q1	q2
1	banana	coke	3	5
2	banana	coke	3	7
3	banana	coke	5	4
4	banana	milk	7	7
5	banana	milk	7	4
6	banana	milk	6	4
7	apple	milk	6	2
8	apple	milk	7	4
9	apple	milk	6	2
10	apple	coke	4	3
11	apple	coke	6	2
12	apple	coke	6	1

two_way_between_design_anova.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# factorのカラム（独立変数）に名前をつける
x1 = factor(data[, 2])
x2 = factor(data[, 3])


# 分散の等質性の検定とt検定
for ( i in 4:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- Summary: factor 1 ---\n")
  print(by(y, x1, summary))
  cat("\n-- Summary: factor 2 ---\n")
  print(by(y, x2, summary))

  # ANOVA
  cat("\n-- ANOVA ---\n")
  print(summary(aov(y ~ x1 * x2)))

  # 多重検定
  cat("\n-- 多重検定 ---\n")
  print(pairwise.t.test(y, x1:x2, p.adjust.method="bonferroni"))
}

2要因参加者内デザイン

two_way_within_design_anova.csv

sample_id	factor:food	factor:drink	q1	q2
1	banana	coke	3	5
2	banana	coke	3	7
3	banana	coke	5	4
1	banana	milk	7	7
2	banana	milk	7	4
3	banana	milk	6	4
1	apple	milk	6	2
2	apple	milk	7	4
3	apple	milk	6	2
1	apple	coke	4	3
2	apple	coke	6	2
3	apple	coke	6	1

two_way_within_design_anova.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# 参加者のカラムに名前をつける
subject = factor(data[, 1])

# factorのカラム（独立変数）に名前をつける
x1 = factor(data[, 2]) # within
x2 = factor(data[, 3]) # within


# 分散の等質性の検定とt検定
for ( i in 4:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- Summary: factor 1 ---\n")
  print(by(y, x1, summary))
  cat("\n-- Summary: factor 2 ---\n")
  print(by(y, x2, summary))

  # ANOVA
  cat("\n-- ANOVA ---\n")
  print(summary(aov(y ~ x1 * x2 + Error(subject/(x1*x2)))))

  # 多重検定
  cat("\n-- 多重検定 ---\n")
  print(pairwise.t.test(y, x1:x2, paired=TRUE, p.adjust.method="bonferroni"))
}

2要因混合デザイン

two_way_mixed_design_anova.csv

sample_id	factor:food	factor:drink	q1	q2
1	banana	milk	3	5
2	banana	milk	3	7
3	banana	milk	5	4
4	banana	coke	7	7
5	banana	coke	7	4
6	banana	coke	6	4
1	apple	milk	6	2
2	apple	milk	7	4
3	apple	milk	6	2
4	apple	coke	4	3
5	apple	coke	6	2
6	apple	coke	6	1

two_way_mixed_design_anova.R

# 引数取得
args <- commandArgs(trailingOnly = T)

# csvファイルの読み込み
data = read.csv(args[1])

# データの確認
cat("\n### Data Format ###\n")
print( str(data) )

# データの要約
cat("\n### Data Summary ###\n")
print( summary(data))

# カラム名を取り出す
columns = colnames(data)

# 参加者のカラムに名前をつける
subject = factor(data[, 1])

# factorのカラム（独立変数）に名前をつける
x1 = factor(data[, 2]) # within
x2 = factor(data[, 3]) # between


# 分散の等質性の検定とt検定
for ( i in 4:ncol(data) )
{
  cat("\n########################\n")
  cat(columns[i])
  cat("\n########################\n")

  # データカラム（従属変数）に名前をつける
  y = data[, i]

  # 要約
  cat("\n--- Summary: factor 1 ---\n")
  print(by(y, x1, summary))
  cat("\n-- Summary: factor 2 ---\n")
  print(by(y, x2, summary))

  # ANOVA (within: x1, between: x2)
  cat("\n-- ANOVA ---\n")
  print(summary(aov(y ~ x1*x2 + Error(subject/x1) + x2 )))

  # 多重検定
  cat("\n-- 多重検定 ---\n")
  print(pairwise.t.test(y, x1:x2, paired=TRUE, p.adjust.method="bonferroni"))
}

Ref.

• 実吉 綾子, フリーソフト「R」ではじめる 心理学統計入門 (知識ゼロでもわかる統計学), 技術評論社 (2013/1/17)
• コピペで学ぶ Rでテクニカルデータプレゼンテーション
• Quick-R: ANOVA