理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

t検定のまとめ

投稿日2017.1.27  最終更新日2017.6.28

Rでは次の式を使用します
t.test ( x, alternative=c("two.sided","less","greater") )
mu=0, paired=FALSE, var.equal=FALSE, conf.level=0.95)

引数
alternative:両側検定("two.sided")か片側検定("less","greater")か
mu=0:母平均が0であるかどうかを検定する
paired=TRUE:対応のあるt検定
var.equa=TRUE:等分散を仮定する
conf.level:信頼区間

統計量
t=(標本平均-母平均)/√(不偏分散/n) 
=(標本平均-母平均)/標準誤差  
自由度n-1のt分布t(n-1)に従う 

 

1標本問題(1つの母集団の母数に関する検定)

確率変数Xからのn個の標本x, μは母平均、S^2は不偏分散

     f:id:yoshida931:20170619141944p:plain
例)次のデータを図で概観し,母平均=160.0を有意水準5%で検定せよ
x = c(153.0,148.0,161.0,166.0,160.0,173.0,173.0,165.0,153.0,168.0,171.0)

箱ひげ図
boxplot(x)
ドット図
stripchart  (x, pch=16, at=0, method="stack")
ヒストグラム
hist  (x, right=FALSE, col="gray")

t検定を実行

t.test (x)

One Sample t-test
data: x
t = 62.79, df = 10, p-value = 2.554e-14
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
157.0405 168.5959
sample estimates:
mean of x
162.8182

 

標本xの平均は162.82
95%信頼区間[157.04, 168.60]で「160.0」が含まれている
したがって帰無仮説は5%有意水準で採択される(p値<0.01)

もし、この検定の有意水準が1%であった場合は

t.test(x,conf.level=0.01)

One Sample t-test
data: x
t = 62.79, df = 10, p-value = 2.554e-14
alternative hypothesis: true mean is not equal to 0
1 percent confidence interval:
162.7849 162.8515
sample estimates:
mean of x
162.8182

 

99%信頼区間[162.78, 162.85]で「160.0」が含まれていません
したがって帰無仮説は1%有意水準では棄却されることになります
研究前のαの設定は重要です

 

対応のある2標本の場合
1標本問題と同じ原理です
治療前後の比較を考えてみましょう
治療前のデータPRE,治療後のデータPOSTとします
治療の前後差を有意水準1%で検定せよ
PRE<-c(25,27,32,37,38,39,44,45,47,48,52,54,59,60,65,65,67,68,68)
POST<-c(41,28,32,60,50,48,44,55,95,48,61,67,60,83,66,77,70,81,87)
対応あるt検定の帰無仮説は「前後差=0」と設定します
def<-POST-PRE
つまり1標本となります

t.test(def)

One Sample t-test
data: def
t = 4.1537, df = 18, p-value = 0.0005965
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
5.540338 16.880715
sample estimates:
mean of x
11.21053

 

差の平均は11.21となり有意水準1%で差を認める
(当然この11.21という数値は検査内容により実質的な意味合いが異なります)

 

2標本問題(2つの母集団の母数に関する検定)
ここでは2標本の平均の差の検定のみを扱います
確率変数Xからのm個の標本x, 確率変数Yからのn個の標本y
標本x,y   標本数m,n   標本平均μx,μy   標本分散sx, sy  

 

母分散が既知の場合
標本 x (m個) , 標本 y (n個)
 ・分散はσ^2
 ・正規分布 N (μ1,σ^2) , N (μ2,σ^2) に従う

標本 x の平均 μx~N(μ1,σ^2/m) 、標本 y の平均 μy~N(μ2,σ^2/n)
μxμy も正規分布に従う
E[μxμy]=μ1-μ2、V[μxμy] = V(μx)+V(μy) = (σ^2)/m + (σ^2)/n

ここをすぐに忘れてしまうので、ちょいと書き加えておきます.
AとBを互いに独立な確率変数として、それぞれの期待値をμa、μbとします.
V(A-B) = E { (A-B) - (μa-μb) }^2 = E { (A-μa) - (A-μb) } ^2
            = V(A) + V(B) - 2E{ (A-μa)*(B-μb) }
共分散は確率変数X,Yが互いに独立な場合には0になるので、
E{ (Aa)*(Bb) }=0
したがって
V(A-B) = V(A) + V(B)

 母分散が既知の場合には、等しくなくても xme-yme を利用して検定可能

            Z = { (μxμy) - (μ1-μ2) } ÷ √(σ^2/m+σ^2/n)

 

母分散が未知で等しい場合
標本 x (m個) , 標本 y (n個)
 ・母分散が等しい・・・推定量プールした分散 (合併した分散)を使用
 ・正規分布 N (μ1,σx^2) , N (μ2,σy^2) に従う

f:id:yoshida931:20170619134940p:plain   

f:id:yoshida931:20170628194132p:plainf:id:yoshida931:20170619140210p:plain


プールした分散 s^2

プールした分散
∑( x - μx )^2   +  ∑( y -μy )^2 ÷ ( m+n-2 )
標本分散sx, sy は
sx=( 1/m-1 )*∑( x - μx )^2   ,   sy=( 1/n - 1 )*∑( y - μy )^2

したがってプールした分散s^2
s^2 = { (m-1 )*sx + ( n-1 )*sy }  ÷  ( m+n-2 )  


f:id:yoshida931:20170619134940p:plain
t=  (μx - μy) - ( μ1 - μ2 )  / √( s^2/m + s^2/n )

t= (μx - μy) - ( μ1 - μ2 ) / {s*( √1/m + 1/n ) }
自由度  m+n-2 のt分布 t (m+n-2) にしたがう

 

例)x若年群の片足立ち(sec)、y高齢群の片足立ち(sec)
x<-c(39,47,48,59,65,68,74,79,80,82,87,88,89,94,99,
100,101,106,107,108,111,120,122,129,130,136,138,139,
147,157,158,179,189,197,213)
y<-c(25,27,32,37,38,44,45,52,54,60,65,67,68,69,71,73,
82,88,89,92,99,100,101,116,126,131,151,156,174)

t.test(x, y, var.equal=T)

Two Sample t-test
data: x and y
t = 2.9234, df = 62, p-value = 0.004829
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
9.671654 51.500760
sample estimates:
mean of x mean of y
111.00000 80.41379

有意水準1%で差があるという結果になるのだが、
95%信頼区間は[9.67,51.50]と範囲が大きい
約9秒の差、50秒の差をどのようにとらえるかは研究者が十分に考察しなければならない

 

検定力分析(事後の分析)
検定力1-βとは、帰無仮説が「偽」であるときに、有意差が検出される確率.
効果量(effect size)とは、母集団における効果の大きさ

標本効果量(cohensのd効果量)の求め方
cohensのd効果量=
  xの平均値とyの平均値の差の絶対値  ÷  2群の共通の標準偏差

cohens_d <- function(x, y) {    
lx <- length(x)- 1         #n‐1を算出
ly <- length(y)- 1         #n‐1を算出
md <- abs(mean(x) - mean(y))   # x群とy群の平均値の差の絶対値
var2 <- lx * var(x) + ly * var(y)   # ( n‐1) × 不偏分散
pov <- var2  / ( lx + ly )      
                #プールした分散 s^2 = { (m-1)*sx + (n-1)*sy } ÷ (m+n-2) 
cd <- md/ sqrt(pov)   
     # sqrt(pov) = 2群に共通の標準偏差
}
res <- cohens_d(x, y)
res
[1] 0.7340753

 http://tjo.hatenablog.com/entry/2014/02/24/192655

 
次にt値から標本効果量を算出
A=length(x)+length(y) 
B=length(x)*length(y) 

標本効果量(cohensのd効果量) = t値 × A / B            
2.9234*sqrt(64/(35*29))
[1] 0.7400945

 

検定力
次に5%水準で検定力を求めてみます
library(pwr)
pwr.t2n.test(n1=35,n2=29, d =0.7752781 , sig.level =0.05 , power = NULL,alternative = "greater")

t test power calculation

n1 = 35
n2 = 29
d = 0.7752781
sig.level = 0.05
power = 0.9205285
alternative = greater

検定力は0.92となり
検定力の高い検定であったといえます


母分散が未知で等しくない場合(ウェルチの検定)
まず等分散性を検定する
Leveneの検定
帰無仮説:xの分散=yの分散
library(car)

score<-c(100,100,101,101,106,107,108,111,116,
120,89,100,145,85,106,107,90,149,116,78)
group=factor(rep(c("x", "y"), c(10,10)))
leveneTest(score ~ group)

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 6.3953 0.02101 *
18

帰無仮説は棄却され、母分散が等しくないことが判明したため
ウェルチの検定を行う

t=(標本xの平均-標本yの平均)/√( s1^2/m + s2^2/n )
このtは帰無仮説の下で近似的に次の自由度fをもつt分布に従う
g1=s1^2/m ,g2=s2^2/n
f=(g1+g2)^2/{(g1^2/(m-1))-(g2^2/(n-1))}
このtを用いる方法はウェルチの検定(Welch's test)とよばれる

ベクトルscoreよりデータを抽出する
x<-score[1:10]
y<-score[11:20]

ウェルチの補正による2標本のt検定
t.test(x,y,var=F)

Welch Two Sample t-test
data: x and y
t = 0.062759, df = 10.485, p-value = 0.9511
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.14089 18.14089
sample estimates:
mean of x mean of y
107.0 106.5