R 対応のないt検定

数学いらずの医科統計学 第2版

対応のないt検定について、気づいたことがあったので備忘も兼ねて書いておく。扱うデータは本にも記載されていた「高齢ラットと若齢ラットの膀胱筋片における高濃度ノルアドレナリン刺激に伴う最大弛緩」のデータ(単位は%Emax)。

old <- c(20.8,2.8,50.0,33.3,29.4,38.9,29.4,52.6,14.3)
young <- c(45.5,55.0,60.7,61.5,61.1,65.5,42.9,37.5)
boxplot(list(old, young))

boxplotしてみる。

お馴染みのRの関数、t.testで検定をする。

test <- t.test(young,old)

…で、結果がこちら。

	Welch Two Sample t-test
data:  young and old 
t = 3.6242, df = 13.778, p-value = 0.002828
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
  9.590586 37.501081 
sample estimates:
mean of x mean of y 
 53.71250  30.16667 

表示されたp値0.002828と本に記載されているp値(0.0030)の数値が合わないことが気になっていた。群馬大青木先生によると、t.testは1次データ（つまりは母集団ということか？）の場合のt検定向けなので、t値計算時は標本分散を用いていると思われる。青木先生の自作関数を用いると本と不偏分散を用いているため結果が一致する。

sa <- sd(young)
sb <- sd(old)
fa <- length(young)-1
fb <- length(old)-1
na <- length(young)
nb <- length(old)

s <- sqrt{(fa*(sa^2)+fb*(sb^2))/(fa+fb)} #文字化け対策で()を{}にしてます
ma <- mean(young)
mb <- mean(old)
t <- (ma-mb)/(s*sqrt(1/length(young)+1/length(old)))

source("http://aoki2.si.gunma-u.ac.jp/R/src/my_t_test.R", encoding="euc-jp")　#関数の使い方はこちら

test <- my.t.test(length(young),mean(young), var(young), 
        length(old), mean(old), var(old), var.equal = TRUE)

結果は、確かに本と一致している。

	等分散を仮定した，二群の平均値の差の検定
data:  
n1 = 8, mean1 = 53.7125, variance1 = 107.406964285714
n2 = 9, mean2 = 30.1666666666667, variance2 = 259.0375 
t = 3.5315, df = 15, p-value = 0.003022