t検定

対応のないt検定

データの読み込みと準備

サンプルデータは、統計教育推進委員会が提供してくださっている運動習慣と脈拍[1]を使います。

R Console
> x <- read.table(pipe("pbpaste"), header=TRUE) # クリップボードからデータを読み込み、xに格納。header=TRUEで先頭行を変数名と設定。ちなみにこちらはMac用のコードでWindowsは、read.table("clipboard", header=TRUE)。 
> str(x)   # xの構造を確認、intは整数を示すが、exerciseの0,1は本当は順序なし因子(0:運動習慣あり、1:運動習慣なし)。
'data.frame':	102 obs. of  2 variables:
 $ HR      : int  59 59 81 70 68 66 65 71 75 61 ...
 $ exercise: int  0 0 0 0 1 1 1 1 0 0 ...
> x$exercise <- as.factor(x$exercise) # exerciseのデータ型を順序なし因子に変換し、上書き。 
> str(x)   # 再びxの構造を確認。exerciseがfactorに変換されている。
'data.frame':	102 obs. of  2 variables: 
 $ HR      : int  59 59 81 70 68 66 65 71 75 61 ...
 $ exercise: Factor w/ 2 levels "0","1": 1 1 1 1 2 2 2 2 1 1 ...
> summary(x)   # xの要約を算出
       HR         exercise
 Min.   : 40.00   0:55    
 1st Qu.: 61.25   1:47    
 Median : 70.00           
 Mean   : 70.21           
 3rd Qu.: 76.75           
 Max.   :122.00
[1]  統計教育推進委員会 – 統計活用提案のための教材サイト

 

基本統計量の算出

運動習慣あり(0)、運動習慣なし(1)別に、平均値、標準偏差、標準誤差を求める。

R Console
> l0 <- nrow(subset(x, exercise=="0")) # 運動習慣ありの対象者数をl0に代入。
> l0
[1] 55
> l1 <- nrow(subset(x, exercise=="1")) # 運動習慣なしの対象者をl1に代入。
> l1
[1] 47
> m0 <- mean(subset(x, exercise=="0")$HR) # 運動習慣ありの心拍数の平均値。
> m0
[1] 69.63636
> m1 <- mean(subset(x, exercise=="1")$HR) # 運動習慣なしの心拍数の平均値。
> m1
[1] 70.87234
> d0 <- sd(subset(x, exercise=="0")$HR) # 運動習慣ありの心拍数の標準偏差。
> d0
[1] 9.972859
> d1 <- sd(subset(x, exercise=="1")$HR) # 運動習慣なしの心拍数の標準偏差。
> d1
[1] 13.94651
> e0 <- d0/sqrt(l0) # 運動習慣ありの心拍数の標準誤差。
> e0
[1] 1.34474
> e1 <- d1/sqrt(l1) # 運動習慣なしの心拍数の標準誤差。
> e1
[1] 2.034307

 

以下の方法でも同様に算出できます。こちらの方が簡便ですね。

R Console
> l <- table(x$exercise) # 2群の対象者数 
> l
 0  1 
55 47 
> m <- tapply(x$HR, x$exercise, mean) # 2群の平均値
> m
       0        1 
69.63636 70.87234 
> d <- tapply(x$HR, x$exercise, sd) # 2群の標準偏差
> d
        0         1 
 9.972859 13.946508 
> e <- d/sqrt(l) # 2群の標準誤差
> e
       0        1 
1.344740 2.034307 

 

等分散性の検定

等分散性の検定としてF検定の一種であるlevene検定を行います。levene検定は等分散であることを帰無仮説としますので、p値が0.05以上であれば帰無仮説は棄却されず、等分散であると言えます。levene検定にはcarパッケージを使います。

R Console
> library(car) # carを読み込む。
> leveneTest(x$HR, x$exercise) #levene検定の実施。
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   1  1.0037 0.3188
      100               

p値が0.05より大きいので、運動習慣あり群と運動習慣なし群は等分散といえます。

等分散の場合はStudentのt検定、等分散でない場合はWelchのt検定

を行います。

MEMO
ちなみにStudentのt検定もWelchのt検定も正規分布を仮定した検定です。正規分布を仮定できないノンパラメトリックなデータにはMann-WhitneyのU検定を行うことになります。しかし、Mann-WhitneyのU検定も実は不等分散では上手く検定できないらしく、そういう時には、Brunner-Munzel検定という手法が有効だとのことです[1],[2]。
MEMO
ちなみにが多くて申し訳ありませんが、このF検定からのStudentもしくはWelch検定の手順については、二段階の多重検定であり、問題視されることも増えているようです[3]。
[1] 奥村晴彦 – Brunner-Munzel検定

[2] ほくそ笑む – マイナーだけど最強の統計的検定 Brunner-Munzel検定

[3] 生物科学研究所 Laboratory of Biology 井口研究室ホームページ – 等分散検定から t検定,ウェルチ検定,U検定への問題点

 

対応のないt検定の実行

t.test(変数~群分け水準)の式になります。 var.equal=を記述しなければデフォルトではFALSEになっており、Welchのt検定が行われます。

t.test(A群の変数, B群の変数, paired=FALSE)でも同様の結果になります。

R Console
> t.test(x$HR~x$exercise, var.equal=TRUE) # var.equal=TRUEで等分散を仮定したStudentのt検定を行うよう指示。

	Two Sample t-test

data:  x$HR by x$exercise
t = -0.52, df = 100, p-value = 0.6042
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5.951673  3.479720
sample estimates:
mean in group 0 mean in group 1 
       69.63636        70.87234 

今回は残念ながら運動経験あり群となし群の間で心拍数に有意差は認められませんでした。

 

効果量の計算

パッケージcompute.esを用いて計算します。

mes(平均1, 平均2, 標準偏差1, 標準偏差2, サンプルサイズ1, サンプルサイズ2)

で計算できます。

MEMO

効果量には数十種類もあるそうで、t検定では群間差についての効果量であるd族が用いられる。有名な効果量としてCohenのdやHedgesのgがあり、Cohenのdは平均値の差を標本標準偏差で割った値、Hedgesのgは平均値の差を不偏標準偏差で割った値である。従って、Cohenのdは記述統計的な意味合い、Hedgesのgは推測統計的な意味合いで使われる。推測統計であるt検定ではHedgesのgが推奨されている(Cohen’s dHedges’ g)[1],[2]。

R Console
> library(compute.es) #compute.esパッケージを読み込む。
> mes(m0, m1, d0, d1, l0, l1) # 上で定めた変数を使う。
Mean Differences ES: 
 
 d [ 95 %CI] = -0.1 [ -0.5 , 0.29 ] 
  var(d) = 0.04 
  p-value(d) = 0.6 
  U3(d) = 45.89 % 
  CLES(d) = 47.09 % 
  Cliff's Delta = -0.06 
 
 g [ 95 %CI] = -0.1 [ -0.49 , 0.29 ] 
  var(g) = 0.04 
  p-value(g) = 0.6 
  U3(g) = 45.92 % 
  CLES(g) = 47.11 % 
 
 Correlation ES: 
 
 r [ 95 %CI] = -0.05 [ -0.25 , 0.15 ] 
  var(r) = 0.01 
  p-value(r) = 0.61 
 
 z [ 95 %CI] = -0.05 [ -0.25 , 0.15 ] 
  var(z) = 0.01 
  p-value(z) = 0.61 
 
 Odds Ratio ES: 
 
 OR [ 95 %CI] = 0.83 [ 0.41 , 1.7 ] 
  p-value(OR) = 0.6 
 
 Log OR [ 95 %CI] = -0.19 [ -0.9 , 0.53 ] 
  var(lOR) = 0.13 
  p-value(Log OR) = 0.6 
 
 Other: 
 
 NNT = -36.17 
 Total N = 102
[1] 統計Web – Cohen’s d

[2] 統計Web – Hedges’ g

 

対応のあるt検定

データの読み込みと準備

サンプルデータは、統計教育推進委員会が提供してくださっているいろいろな樹木の葉の長さと幅のデータ[1]を使います。ダウンロードしたエクセルファイルにはいろいろな樹木のデータがシート別にまとめられています。今回はハナミヅキの葉の長さと幅に差があるかどうか比較することにします。葉の長さをleaf_l、葉の幅をleaf_wと名前を付けてデータを読み込みます。

R Console
> x <- read.table(pipe("pbpaste"), header=TRUE) # クリップボードからデータを読み込み、xに格納。header=TRUEで先頭行を変数名と設定。ちなみにこちらはMac用のコードでWindowsは、read.table("clipboard", header=TRUE)。 
> str(x) # xの構造を確認。numは実数のこと。
'data.frame':	50 obs. of  2 variables:
 $ leaf_l: num  8 7 7.5 7.5 7 7.5 6.5 7.5 6.5 5.5 ...
 $ leaf_w: num  6.5 5 6.5 7 4 5.5 3.5 5 4.5 3 ...
> summary(x) # xの要約を算出。
     leaf_l          leaf_w     
 Min.   :4.300   Min.   :2.200  
 1st Qu.:5.625   1st Qu.:3.425  
 Median :6.500   Median :4.000  
 Mean   :6.378   Mean   :4.164  
 3rd Qu.:7.000   3rd Qu.:5.000  
 Max.   :8.000   Max.   :7.000  
[1]  統計教育推進委員会 – 統計活用提案のための教材サイト

 

 

基本統計量の算出

leaf_lとleaf_wのサンプルサイズ、平均値、標準偏差、標準語さを計算します。

R Console
> attach(x) # 「x$」の部分を記述しなくて良いようにするおまじない。
> l <- c(length(leaf_l), length(leaf_w)) # length関数でleaf_lとleaf_wのサンプルサイズ(ベクトルの長さ)を求める。c()で1つのベクトルにくくる。このデータは1つの葉に対して長さと幅を測定しているので、当然ながらサンプルサイズは一致する。 
> l
[1] 50 50
> m <- c(mean(leaf_l), mean(leaf_w)) # 2変数の平均値 
> m
[1] 6.378 4.164
> d <- c(sd(leaf_l), sd(leaf_w)) # 2変数の標準偏差 
> d
[1] 0.8572119 1.1385203
> e <- d/sqrt(l) # 2群の標準誤差  
> e
[1] 0.1212281 0.1610111

 

対応のあるt検定の実行

R Console
> t.test(leaf_l, leaf_w, paired = TRUE) # paired=TRUEで対応のあるt検定になる。

	Paired t-test

data:  leaf_l and leaf_w
t = 16.091, df = 49, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.937505 2.490495
sample estimates:
mean of the differences 
                  2.214 

p値がとても小さいので有意差が認められました。平均値からハナミヅキ の葉は縦長だと分かります。

 

効果量の計算

t検定の効果量gの計算方法は、

g = t検定のt値 * sqrt( 1/(サンプルサイズ))

です。

R Console
> t <- t.test(leaf_l, leaf_w, paired = TRUE) # 上述の対応あるt検定をt定める。 
> t$statistic # tに含まれるt値を確認。
       t 
16.09142 
> g <- t$statistic * sqrt( 1/ length(leaf_l)) # gを計算。サンプルサイズはleaf_lとleaf_wは同数なので、ここではleaf_lを指定。 
> g
      t 
2.27567 

効果量は非常に大きいことがわかります。

[1] 図の学のWebPage 総合案内所 目次 梶山喜一郎 – 対応のあるt検定の効果量計算

[2] 水本篤 – 効果量と検定力分析入門 ―統計的検定を正しく使うために―(PDF)

Print Friendly, PDF & Email

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください