R
心理学

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

はじめに

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

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

スクリプトの使い方

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

ソースコード

https://github.com/DaiHasegawa/R_Scripts

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 が各条件の「参加者数-1」の和と等しいことを確認する
  • t検定で p < .05 なら危険水準.05で有意差あり、 p < .01 なら危険水準.01で有意差あり

結果の報告(注:tp などの統計量はイタリックにする)

  • 対応のない 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は別にみる)

  • 全体の自由度=データ数ー1=11
  • 要因の自由度=条件数ー1=2
  • 残差(誤差)の自由度=全体の自由度ー要因の自由度=22

結果の報告
- 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

結果の確認

  • 全体の自由度=データ数-1=11
  • 要因の自由度=条件数-1=2
  • 残差の自由度=(参加者数-1)x(条件数ー1)= 6

結果の報告
- 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.