Help us understand the problem. What is going on with this article?

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした