Edited at

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

More than 1 year has passed since last update.


はじめに


実験デザインについて(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.