実験デザイン別の統計的検定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`

``` 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) )
}
```

実行結果

``` sh
$ 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`

``` 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」

``` sh
$ 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`

``` 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"))
}
```

実行結果

``` sh
$ 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`

``` 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`

``` 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`

``` 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`

``` 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でテクニカルデータプレゼンテーション](http://monge.tec.fukuoka-u.ac.jp/R_analysis/0r_analysis.html)
- [Quick-R: ANOVA](http://www.statmethods.net/stats/anova.html)