★手計算とRで学ぶ統計学: t検定

まずは，どれだけのデータを集めればよいのかを調査の前に把握しておきましょう。

以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする
install.packages("pwr")
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む
library(pwr)

# 1標本t検定
# サンプルサイズ設計
pwr.t.test(n = NULL, 
           d = 0.5, 
           sig.level = 0.05, 
           power = 0.8,
           type = "one.sample", 
           alternative = "greater")

最低でも27本のスナックがサンプルとして必要であることがわかりました。今回は30本のスナックを使って仮説を検証することにしましょう。

     One-sample t test power calculation
               n = 26.13753
               d = 0.5
       sig.level = 0.05
           power = 0.8
     alternative = greater

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まずは作業ディレクトリを設定します。次にデータを読み込んで変数に代入し，列名でデータを指定できるようにします。

# 作業ディレクトリを設定 
setwd("csvファイルが入っているフォルダ名")
# 使用するデータの設定
data_snack <- read.csv("t_test_onesample.csv")
# 列名でデータを指定できるようにする
attach(data_snack)

データを読み込んだら，必ずデータの中身を確認しましょう。解析を進めた後にデータの中身が違うと，それまでの作業が徒労に終わります。

# データ内容の確認
head(data_snack)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります (※1)。

> head(data_snack)
  ID length 
1  1  74.64 
2  2  68.13 
3  3  72.66 
4  4  71.82 
5  5  68.89 
6  6  71.61

データを読み込んだあとは， ヒストグラムを表示して分布の形状を確認します。

# ヒストグラムの表示
hist(length)

ヒストグラムを見る限りでは，このデータに2つ以上のピークがないことがわかります。
統計学の教科書だと，この後にデータの正規性や等分散性を検討するための検定を実施してから手法を決めるように記述があるものがあります。しかし，多重検定の問題 (同じデータに対して統計的仮説検定を複数回行うと第一種の誤りが発生しやすくなる) が指摘されています。
そのため，最近はこれらの検定を行わずに目的の検定に進むことがデフォルトになりつつあります。

(※1) 列数が多いデータの場合はhead()で表示できない箇所に誤りがある可能性もあります。
以前質問紙調査で1つ設問が重複していたため，列名を変更するときにエラーが出てしまいました。慌ててデータを見直し，顔面蒼白になって指導教員に相談したことがあります。

実際にRで1標本\(t\)検定を行い，抽出したポテトスナックの母平均\(\mu_{0}\) が\(70.0\;mm\)であるかどうかを検討しましょう。

以下のコードをR上で実行すると

# 1標本t検定
t.test(length, #データの指定
       mu = 70.00, #比較対象となる値
       alternative = "greater", #「ポテトスナックの長さは70.00 mmより長い」というのが対立仮説
       conf.level = 0.95)

たちまちに結果が得られます。
なお，95 percent confidence intervalは95%信頼区間 (以後，95% CIと表記します) といい，「この区間に母集団の平均値が存在する」ということを示します。

    One Sample t-test data:
 length t = 4.6299, df = 29, p-value = 3.544e-05
 alternative hypothesis: true mean is greater than 70
 95 percent confidence interval:  71.01492      Inf
 sample estimates:
 mean of x
  71.60333 

\(t(29) = 4.63,\;p < .001\) という結果から，5%有意水準で帰無仮説が棄却され，抽出したポテトスナックの長さの母平均は\(70.0mm\)とは異なると示されました。

近年は，統計検定量と\(p\)値に加えて，「どの程度の差があったのか」を示す値である効果量を論文中で報告することを義務付けるジャーナルが増えてきています。もちろん，投稿規定で義務付けられていなくても「この研究において得られた効果はどの程度だったのか」を読者にわかりやすく示す上では効果量を示すことが望ましいです。
効果量を計算する以下のコードをR上で実行すると，\(Cohen's\;d = 0.845293...\) と非常に大きな値であったことがわかりました。

# 効果量 (Cohen's d) を求める 
# (初めてパッケージを使用するときのみ行う) パッケージをインストールする
install.packages("effsize")
# (RStudioを新しく起動したら毎回行う) パッケージを読み込む
library(effsize) 

# 効果量の計算 
# muは比較対象とする値
# hedges.correctionをTRUEとすると，計算結果はHedge's gになる
cohen.d(data_snack$length,f=NA, mu=70, hedges.correction = FALSE)

計算結果は以下のように出力されます。

> cohen.d(data_snack$length,f=NA, mu=70, hedges.correction = FALSE) 
Cohen's d (single sample) 

d estimate: 0.845293 (large) 
Reference mu: 70 
95 percent confidence interval: 
     lower      upper 
0.06584309 1.62474285 

サンプルサイズ設計の段階では \(Cohen's\;d = 0.8\) (大きな効果を検出できることを想定) と設定しましたが，実際にはそれを上回る大きな効果であったことがわかります。

どのような仮説のもとでt検定を行なったか，その結果はどうだったか，を論文やレポートに記述する際は，今回の解析例の場合は以下のように記載します。

(※) この書き方は心理学をはじめとする社会科学分野で一般的なものです。実際にレポートや論文を書く際には，担当教員や指導教員の指示，それぞれの分野の慣例に従ってください。  

まずは，実験や調査の前にどれだけのデータを集めればよいのかを把握しておきましょう。

以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする 
install.packages("pwr") 
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(pwr) 

# 対応のないt検定
# サンプルサイズ設計
pwr.t.test(n = NULL, # サンプルサイズを調べたいので設定しない
           d = 0.5, # 中程度の効果量
           sig.level = 0.05, # 有意水準
           power = 0.8, # 検定力
           type = "two.sample", # 対応のないt検定
           alternative = "two.sided" # 「2群の平均値に差がある」という対立仮説を設定。差の正負は問わない
　　　　　　)

最低でも各群64人の参加者が必要であることがわかりました。

Two-sample t test power calculation
               n = 63.76561
               d = 0.5
       sig.level = 0.05
           power = 0.8
     alternative = two.sided 
NOTE: n is number in *each* group # 各群の参加者数であることに留意

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まずは作業ディレクトリを設定します。次にデータを読み込んで変数に代入し，列名でデータを指定できるようにします。

# 作業ディレクトリを設定 
setwd("csvファイルが入っているフォルダ名")
# 使用するデータの設定
data_ttest_indep <- read.csv("t_test_twosample.csv")
# 列名でデータを指定できるようにする
attach(data_ttest_indep)

データを読み込んだら，必ずデータの中身を確認しましょう。

# データ内容の確認
head(data_ttest_indep)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります。

> head(data_ttest_indep)
  ID school_A school_B 
1  1       80       76 
2  2       88       61 
3  3       93       68 
4  4       76       77 
5  5       79       63 
6  6       98       61

データを読み込んだあとは， ヒストグラムを表示して分布の形状を確認します。

# ヒストグラムの表示
hist(school_A)
hist(school_B)

ヒストグラムを見る限りでは，いずれのデータにも2つ以上のピークがないことがわかります。

事前に計算したサンプルサイズをもとにそれぞれの高校でテストを行ったところ，A高校の平均点は85.97点 \(\left(SD = 8.63\right)\), B高校の平均点は74.42点 \(\left(SD = 9.53\right)\) でした。
実際にRで対応のない\(t\)検定を行い，2校の平均点に差があるかどうかを検討しましょう。

以下のコードをR上で実行すると

# 対応のない2標本t検定の実施
t.test(school_A, school_B, #データの指定
       alternative = "two.sided", #対立仮説の指定 (今回は「平均値に差がある」だからtwo.sided)
       paired = FALSE, #対応の有無。対応なしなのでFALSE
       var.equal = TRUE, #等分散性の仮定。等分散を仮定しないWelchのt検定を行うことが近年は一般的。
       conf.level = 0.95)

たちまちに結果が得られます。
なお，95 percent confidence intervalは95%信頼区間 (以後，95% CIと書きます) といい，「この区間に真の母平均が存在する」ということを示します。

Two Sample t-test 
data:  school_A and school_B 
t = 7.1836, df = 126, p-value = 5.238e-11 
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:   8.365889 14.727861 
sample estimates: 
mean of x     mean of y  
85.96875      74.42188  

\(t(126) = 7.18,\;p < .001\) という結果から，5%有意水準で帰無仮説が棄却され，A高校とB高校で数学のテストの平均点に有意な差があった (正確には，差がないとはいえない) ことが示されました。

効果量を計算する以下のコードをR上で実行すると，\(Cohen's\;d = 1.26989...\) と非常に大きな値であったことがわかりました。

# 効果量 (Cohen's d) を求める
# (初めてパッケージを使用するときのみ行う) パッケージをインストールする 
install.packages("effsize") 
# (RStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(effsize) 

# 効果量の計算
# hedges.correctionをTRUEにするとHedges's gという値になる
cohen.d(school_A, school_B, hedges.correction = FALSE)

計算結果は以下のように出力されます。

> cohen.d(school_A, school_B, hedges.correction = FALSE) 
Cohen's d d estimate: 1.269893 (large) 
95 percent confidence interval: #効果量の95%信頼区間    
lower     upper 
0.8864147 1.6533706 

サンプルサイズ設計の段階では \(Cohen's\;d = 0.5\) (中程度の効果を検出できることを想定) と設定しましたが，実際にはそれを上回る大きな効果であったことがわかります。

どのような仮説のもとでt検定を行なったか，その結果はどうだったか，を論文やレポートに記述する際は，今回の解析例の場合は以下のように記載します。

(※) この書き方は心理学をはじめとする社会科学分野で一般的なものです。実際にレポートや論文を書く際には，担当教員や指導教員の指示，それぞれの分野の慣例に従ってください。  

まず，対応のある\(t\)検定のサンプルサイズ設計を行いましょう。以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする 
install.packages("pwr") 
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(pwr) 

# 対応のあるt検定
# サンプルサイズ設計
pwr.t.test(n = NULL, #サンプルサイズを調べたいので設定しない
           d = 0.8, #大きな効果量を得たい
           sig.level = 0.05, #有意水準
           power = 0.8, #検定力
           type = "paired", #対応のあるt検定
           alternative = "greater" #「どちらかが大きい」という対立仮説を採用。差があるだけじゃだめ
　　　　　　)

最低でも12人の参加者が必要であることがわかりました。

Two-sample t test power calculation
               n = 11.14424
               d = 0.5
       sig.level = 0.05
           power = 0.8
     alternative = greater 
NOTE: n is number of *pairs* #検定に使用するペアの数 (=参加者数) であることに留意

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まずは作業ディレクトリを設定します。次にデータを読み込んで変数に代入し，列名でデータを指定できるようにします。

# 作業ディレクトリを設定 
setwd("csvファイルが入っているフォルダ名")
# 使用するデータの設定
data_ttest_paired <- read.csv("t_test_paired.csv")
# 列名でデータを指定できるようにする
attach(data_ttest_paired)

データを読み込んだら，必ずデータの中身を確認しましょう。解析を進めた後にデータの中身が違うと，それまでの作業が徒労に終わります。

# データ内容の確認
head(data_ttest_paired)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります。

> head(data_ttest_paired)
  before after
1    389   431
2    424   430
3    415   463
4    395   294
5    445   452
6    349   349

データを読み込んだあとは， ヒストグラムを表示して分布の形状を確認します。

# ヒストグラムの表示
hist(before)
hist(after)

ヒストグラムを見る限りでは，いずれのデータにも2つ以上のピークがないことがわかります。

対応のある\(t\)検定を実施するコードを以下に示します。
以下のコードをR上で実行すると

# 対応のあるt検定の実施
t.test(before, after, #データの指定
       alternative = "greater", #対立仮説の指定 (今回はBefore < After (差の方向性が重要) なのでgreater)
       paired = TRUE, #対応の有無。今回は対応ありなのでTRUE
       var.equal = FALSE, #等分散性の仮定。
       conf.level = 0.95)

t検定の結果がたちまちに得られます。

Paired t-test
data:  before and after # 使用したデータ
t = -1.4358, df = 11, p-value = 0.9106 # 左から順に統計検定量 (t値)，自由度，p値
alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:   
 -91.90869       Inf # 母平均の差の95%信頼区間

sample estimates: 
mean of the differences 
              -40.83333 

\( t (11) = -1.44, p = .911\) という結果から，5%有意水準のもとで帰無仮説に関する判断が保留され，ドリンクを投与する前後での平均睡眠時間には有意な差があるとはいえないことがわかりました。
 

効果量を計算する以下のコードを走らせると，\(Cohen's\,d = -0.66873...\) と (値としては) 中程度の効果がありそうです。しかし，効果量の95%信頼区間の中に0が含まれる ( =before, afterの差の母平均が0の可能性がある) ために，「統計的仮説検定では有意ではなかったが，効果量をみると差はありそう」という結論を出すことも難しそうです。

# 効果量 (Cohen's d) を求める
# (初めてパッケージを使用するときのみ行う) パッケージをインストールする 
install.packages("effsize") 
# (RStudioを新しく起動したら毎回行う) 
パッケージを読み込む library(effsize) 

# 効果量の計算
# hedges.correctionをTRUEにするとHedges's gという値になる
cohen.d(before, after, hedges.correction = FALSE)

計算結果は以下のように出力されます。

> cohen.d(before, after, hedges.correction = FALSE) 
Cohen's d d estimate: -0.6687328 (medium)
95 percent confidence interval: # 効果量の95%信頼区間    
lower     upper 
-1.5387303  0.2012648

どのような仮説のもとで手元のデータに対して統計的仮説検定を行なったか，その結果はどうだったか，を論文やレポートに記述する際は以下のように記載します。