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

統計学備忘録 since2016

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

2標本平均の差の検定

2標本問題のまとめ
以前投稿した記事を全面改修して再投稿です
投稿日2017.9.6


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

ここでは2標本の平均の差の検定のみを扱います.

まずは記号の定義・・・これ大事
母平均μx, 標本x, サイズm, 標本平均xa
母平均μy, 標本y, サイズn, 標本平均ya


母分散が既知の場合

標本x, y はそれぞれ正規分布 N (μx, σx^2) , N (μy, σy^2) に従います.標本 x の平均 xa~N(μx,\frac{σx^2}{m}) 、標本 y の平均 ya~N(μy,\frac{σy^2}{n})となり、μx-μy も正規分布に従います.ここで、まず標本 x の平均 xaがなぜ正規分布 N(μx,\frac{σx^2}{m}) に従うのか考えてみましょう.

前述の標本xを例にとります.
E(xa)=μx       期待値:標本平均は常に母平均の不偏推定量となります
V(xa)が問題です!   xaの分散:「標本の平均値」の分散
ここで多くの標本を x1, x2, ..., xn と仮定して考えてみます.(ここがポイント
V( x1 + x2 + ... + xm ) = V(x1) + V(x2) + ... + V(xm)
V(x1) = V(x2) = ... = V(xm) = σx^2 なので
V( x1 + x2 + ... + xm ) = m×σx^2
xaの分散:「標本の平均値」の分散は
V(xa) =  V(\frac{x1 + x2 + ... + xm}{m}) = \frac{V(x1 + x2 + ... + xm)}{m^2} となります(ここがポイント
したがって
V[xa] = \frac{σx^2}{m}
標本xの母集団の分散の推定値は\frac{σx^2}{m} となり、標本 x の平均 xaは正規分布N(μx,\frac{σx^2}{m}) に従うことになります.
また標準誤差はその平方根\frac{σx}{\sqrt{m}}となります.

分散の加法性

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{ (A-μa)*(B-μb) }=0 相関なし
したがって
V(A-B) = V(A) + V(B)

これで母分散が既知の場合の2標本検定の準備が整いました.xa-yaは正規分布N(E[xa-ya], V[xa-ya])に従います.
E[xa-ya] = μx - μy
V[xa-ya] = \frac{σx^2}{m}+\frac{σy^2}{n}\ \ \ 分散の加法性より
つまり、xa-yaは正規分布N(μx - μy, \frac{σx^2}{m}+\frac{σy^2}{n}) に従うことになり、
標準化した結果、N(0,1)に従うことになり検定が可能となります.

\ \ \ \ Z\ =\ \frac{(xa-ya)-(μx - μy)}{\sqrt{(\frac{σx^2}{m}\ +\ \frac{σy^2}{n})}}

Rを使って標準正規分布確率密度関数、累積分布関数のグラフを描いてみます

par(mfrow = c(1, 2))
curve(dnorm(x, 0, 1), xlim = c(-5, 5), xlab = "確率密度関数", ylab = "")
curve(pnorm(x, 0, 1), xlim = c(-5, 5), xlab = "累積分布関数", ylab = "")
dev.off()

f:id:yoshida931:20170905105517j:plain:w600
母分散\sigma^2が未知で、標本サイズが大きい場合(n≧30と書いてあることが多い)には、不偏分散を母分散として代用します(一致性
例題)帰無仮説xとyの母集団には差がない (μx=μy)ことを有意水準5%で検定せよ
標本x:xa=22.3、m=50、Sx=3.8
標本y:ya=21.4、n=40、Sy=3.1

z <- (22.3-21.4)/sqrt((3.8^2/50)+(3.1^2/40))
= 1.237355
qnorm(0.975)
= 1.959964
qnorm(0.975) > z
図で理解しましょう
curve(dnorm(x, 0, 1), xlim = c(-5, 5), xlab = "", ylab = "",xaxt="n")
name<-c("z","1.96")
axis(side=1,at=c(1.24, 1.96),labels=name,cex.axis=1.3)

f:id:yoshida931:20170907094506j:plain:w300
したがって帰無仮説は棄却できません(採択).この検定結果では両標本間には差がないと判断します.

母分散が未知で等しい場合
まずは記号の定義・・・これ大事
母平均μx, 標本x, サイズm, 標本平均xa, 不偏分散 Sx^2
母平均μy, 標本y, サイズn, 標本平均ya, 不偏分散 Sy^2

準備として合併した分散S^2を勉強しましょう.
母分散が等しい場合には、推定量はプールした分散 (合併した分散)になります.
まずそれぞれの標本の不偏分散を考えてみます
Sx^2=\frac{\sum(x-xa)^2}{m-1}\ ,\ Sy^2=\frac{\sum(y-ya)^2}{n-1}
合併した場合の分母はm+n-1ではなく(m-1)+(n-1)=m+n-2になることに注意
S^2=\frac{\sum(x-xa)^2+\sum(y-ya)^2}{m+n-2}=\frac{(m-1)Sx^2+(n-1)Sy^2}{m+n-2}
ここで
Z\ =\ \frac{(xa-ya)-(μx - μy)}{\sqrt{(\frac{σ^2}{m}\ +\ \frac{σ^2}{n})}}σ^2が未知であるので、S^2で代用して

t\ =\frac{(xa-ya)-(μx - μy)}{\sqrt{(\frac{S^2}{m}\ +\ \frac{S^2}{n})}}=\frac{(xa-ya)-(μx - μy)}{S(\frac{1}{m}+\frac{1}{n})}
は自由度m+n-2のt分布に従うことがわかります.

標本xと標本yの差の検定を行う場合のt統計量を求めてみます.
帰無仮説はμx=μyなので\ t=\frac{(xa-ya)}{\sqrt{(\frac{S^2}{m}\ +\ \frac{S^2}{n})}}\ に代入します.

x <- c(39, 47, 48, 59, 65, 68, 74, 79, 80, 82, 87, 88, 89, 94, 99)
y <- c(25, 27, 32, 37, 38, 44, 45, 52, 54, 60, 65, 67, 68, 69, 71)
pool <- ((length(x) - 1) * var(x) + (length(y) - 1) * var(y)) / (length(x) + length(y) - 2)
t <- (mean(x) - mean(y)) / sqrt((pool / length(x)) + (pool / length(y)))
= 3.6449

例)Rを使って標本xと標本yの差の検定を有意水準5%で検定してみます.

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

data: x and y
t = 3.6449, df = 28, p-value = 0.001079
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
10.04508 35.82159
sample estimates:
mean of x mean of y
73.20000 50.26667

母分散が未知で等しくない場合(等分散が否定された場合)
ウェルチの近似法

自由度 = v=\frac{(\frac{Sx^2}{m}+\frac{Sy^2}{n})^2}{\frac{(\frac{Sx^2}{m})^2}{m-1}+\frac{(\frac{Sy^2}{n})^2}{n-1}}

母分散が未知で等しくない場合には、t統計量は、近似的に上記の自由度vにもっとも近い整数v'のt分布t(v')に従うことが知られています.
例)次のx、yの標本を使用して2群の差の検定を実施します.まずLeveneの検定にて等分散を検定(帰無仮説:xの分散=yの分散)します.

install.packages("lawstat")
library(lawstat)

x <- c(100,100,101,101,106,107,108,111,116,120)
y <- c(89,100,145,85,106,107,90,149,116,78)
xy <- c(x,y)
group <- factor(rep(c("x", "y"), c(10,10)))
levene.test(xy, group)

modified robust Brown-Forsythe Levene-type test based on the absolute deviations from the median
data: xy
Test Statistic = 6.3953, p-value = 0.02101

帰無仮説は棄却され、母分散が等しくないことが判明したためウェルチの検定を行います.
自由度vに最も近い整数v'を自由度として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