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

統計学備忘録 since2016

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

Z検定

Z検定
投稿日2017.9.19

記号の定義

母平均μx, 母分散σx, 標本x, サイズm, 標本平均xa
母平均μy, 母分散σy, 標本y, サイズn, 標本平均ya

平均値のZ検定

2標本平均の差の検定 - 統計学備忘録 since2016より
平均値の標準誤差=\sqrt{V(xa)=\frac{σx^2}{m}}\ =\ \frac{サンプルの標準偏差}{\sqrt{サンプルサイズ}}
Z=\frac{(xa-μx)}{\sqrt{\frac{σx^2}{m}}}

問題 xa=38,σx=15, n=100, α=0.05 で次の帰無仮説を検定せよ!
帰無仮説H0:μx = 40 対立仮説H1:μx ≠ 40
(出典:村上 正康,安田 正実; 統計学演習, 培風館,1989,p119)
上記のような問題は次のように書き換えて記憶しておくとよいと思います.(現実的ではないのですが)十分大きなサンプル(n=100)の平均が38、分散が15だった場合に、その母集団の平均は40であることを有意水準0.05で検定せよ.
以下Rで計算してみます

z <- (38 - 40) / sqrt(15^2/100)
- 1.959964 = -1.33
qnorm(0.025,lower.tail = T)
= - 1.959964
z>- 1.959964のため帰無仮説は採択される
pnorm(z,lower.tail = T)
= 0.09121122

f:id:yoshida931:20170919125119j:plain:w500
水色の部分の面積が0.09121122ということになります.したがって、有意水準0.05で母集団の平均は40であると言えます.またサンプル数が少ない場合にはt分布を使用した検定になります.ちなみにこの図の描き方は・・・polygonという関数を使うらしいです.http://cse.naro.affrc.go.jp/takezawa/r-tips/r/51.html

plot(dnorm, -4, 4, xaxt="n")
xvals <- seq(-4, -1.33, length=10)        # -4以上-1.33以下 領域をx軸方向に10個の多角形(台形)に等分割
dvals <- dnorm(xvals)             # 対応するグラフの高さ
polygon(c(xvals,rev(xvals)),
        c(rep(0,10),rev(dvals)),col=5)            # 塗りつぶす
name<-c("-1.96","-1.33","0")
axis(side=1,at=c(-1.96,-1.33,0),labels=name)         #指定した場所にnameを挿入

平均値の差のZ検定

平均値の差の標準誤差
V(xa-xy)=V(xa)+V(xy)   分散の加法性
V(xa)=\frac{σx^2}{m}, V(ya)=\frac{σy^2}{n}

平均値の標準誤差 = {\sqrt{\frac{σx^2}{m}+\frac{σy^2}{m}}}

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

問題 xa=22.31, xy=21.54, σx=3.8, σx=3.2, m=50, n=40, α=0.05 で次の帰無仮説を検定せよ.また2群の差の95%信頼区間を求めてよ.
帰無仮説H0:μx = μy 対立仮説H1:μx ≠ μy
(出典:村上 正康,安田 正実; 統計学演習, 培風館,1989,p119)
上記のような問題は次のように書き換えて記憶しておくとよいと思います.
m=50、平均22.31、分散3.8 と n=40、平均21.54、分散3.2 の2群には差がないことを有意水準0.05で検定せよ.

z <- (22.31-21.54)/sqrt(3.8^2/50 + 3.2^2/40)
z = 1.043211
qnorm(0.025,lower.tail = F)
= 1.959964
  #  z<1.959964のため、帰無仮説は採択され、2群には差がないという検定結果になります
2群の差の95%信頼区間は
(22.31-21.54) - 1.959964 * sqrt(3.8^2/50 + 3.2^2/40)
=-0.6766606
(22.31-21.54) + 1.959964 * sqrt(3.8^2/50 + 3.2^2/40) 
= 2.216661
  # -0.6766606 < 2群の差 < 2.216661 より2群には差がないことが推測できます.

比率のZ検定

確率変数Xが二項分布B(n,p)に従う場合、期待値E[X]=np, 分散V[X]=np(1-p)となります.nがある程度大きい場合に、比率pの推定値としてp'=x/n を用います. n=5の場合yesが3回でたらP(YES)=3/5.
p'の期待値はE(p')=p
p'の分散はV(p')=V(x/n)=\frac{V(x)}{n^2}=\frac{np(1-p)}{n^2}=\frac{p(1-p)}{n}

比率の標準誤差={\sqrt{\frac{p(1-p)}{n}}}

z=\frac{p'-p}{\sqrt{\frac{p(1-p)}{n}}}

問題  ある月に1020人を調査したところ内閣支持率は57%でした.内閣支持率の95%信頼区間を推定せよ.

z <- (0.57 - p)/sqrt(0.57*0.43/1020) # ~N(0,1)
pmi <- 0.57-z*sqrt(0.57*0.43/1020) # 下限
pma <- 0.57+z*sqrt(0.57*0.43/1020) # 上限
z <- qnorm(0.025,lower.tail = F) # p値0.025のときのZ値
0.5396178 < p < 0.6003822

比率の差のZ検定

比率の標準誤差={\sqrt{\frac{p(1-p)}{n}}} より比率の分散は\frac{p(1-p)}{n}
例としてYes, Noの割合を考えます.サイズmのA群のyesの比率をp1, サイズnのB群のyesの比率をp2とします.p1-p2の期待値はE(p1-p2), p1-p2の分散はV(p1-p2)=V(p1)+V(p2)・・・分散の加法性より.

比率の差の標準誤差=\sqrt{{\frac{p1*(1-p1)}{m}+\frac{p2*(1-p2)}{n}}}

A群、B群の母集団の比率をP1、P2とした場合

z = \frac{(p1-p2)-(P1-P2)}{\sqrt{\frac{p1*(1-p1)}{m}+\frac{p2*(1-p2)}{n}}}

どんなときに使用するか?まずはクロス表から考えてみます.
問題 男性530人、女性580人にある質問した場合に、yesと答えた割合は男性0.24、女性0.30だった.慌て者は迷わず女性が多いと考えると思います.この回答結果に有意差(α=0.05)はあるのでしょうか?帰無仮説H0:P1=P2

z <- (0.25-0.30)/sqrt(0.25*0.75/530+0.30*0.70/580)
z = -1.868793
qnorm(0.025, lower.tail = T)
= -1.959964
z > -1.959964
-0.1024393 < P1-P2 < 0.002439296
したがって帰無仮説は採択され、男女の割合に差は認められない.

仮説検定のP値や信頼区間は「自分があわて者である」かどうかを教えてくれる.それをどう活かすかというところにこそ、あなたの経験と勘を活かせばいいのである.

西内 啓; 統計学が最強の学問である[実践編]---データ分析のための思想と方法, ダイヤモンド社, 2014

独立性の検定

独立性の検定
投稿日2016.11.4
更新日2017.9.13

f:id:yoshida931:20170913173322j:plain:w500
  
この分割表において独立とは  Ai ∩ Bj の各確率に対して
全ての i j(どの i j でも)に対して  P(Ai ∩ Bj )=P( Ai ) * P( Bj )  となる

AとBが完全に独立な場合の例
  A1 : A2 : A3  がどのBでも  1 : 3 : 6 
  B1 : B2 : B3  がどのAでも  1 : 3 : 5
f:id:yoshida931:20170913173555j:plain:w250

帰無仮説の表現いろいろ
事象Ai とBjは独立である 
A1 , A2 , ... , Ak の起こり方はどのBjに対しても共通である
全ての i jに対して  P(Ai ∩ Bj )=P( Ai ) * P( Bj )

期待度数(理論度数)
Ai の出現確率の推定値は  P( Ai ) = 行の周辺度数 / 総度数 = ni・ / n
Bj の出現確率の推定値は  P( Bj ) = 列の周辺度数 / 総度数 = nj・ / n
したがって セルOij の期待度数Eiは、総度数×セルの確率
Ei  =  n * P( Ai ∩ Bj )  = n * P( Ai ) * P( Bj )  =  n * ( ni・ / n ) * ( nj・ / n )

カイ2乗値
=  ∑∑ ( ( ( 観測度数 - 期待度数 )^2 ) / 期待度数 )

2つの異なる属性A,Pについて度数を集計する.
rn <- c ( "P1" , "P2" , "P3" )
A1 <- c ( 12, 10, 13 )
A2 <- c ( 20, 12, 16 )
A3 <- c ( 20, 15, 19 )
y <-  data.frame (A1,A2,A3)
rownames ( y ) <- paste0 ( rn )
f:id:yoshida931:20170913194511j:plain:w200
周辺度数を算出します

ym <- matrix(c(A1,A2,A3),3,3)
addmargins(ym)

f:id:yoshida931:20170913195030j:plain
次は期待値を求めます

E1A <- (52*35)/134
E1B <- (52*48)/137
E1C <- (52*54)/137
E2A <- (37*35)/137
E2B <- (37*48)/137
E2C <- (37*54)/137
E3A <- (48*35)/137
E3B <- (48*48)/137
E3C <- (48*54)/137

カイ2乗値を求めてみます  自由度 = ( 3-1 ) * ( 3-1 ) = 4

( 12 - E1A )^2 / E1A +
  ( 20 - E1B )^2 / E1B +
  ( 20 - E1C )^2 / E1C +
  ( 10 - E2A )^2 / E2A +
  ( 12 - E2B )^2 / E2B +
  ( 15 - E2C )^2 / E2C +
  ( 13 - E3A )^2 / E3A +
  ( 16 - E3B )^2 / E3B +
  ( 19 - E3C )^2 / E3C
=0.5099473

P値

pchisq(0.5099473, 4, lower.tail = F)      # 累積分布 pchisq ( カイ二乗値 , 自由度 ) 
=0.9725254

帰無仮説は棄却できないため、P1~3におけるAの比率に明確な違いがあるとは言えません.

Rの関数で確認します

chisq.test ( y , correct=F )

f:id:yoshida931:20170913200859j:plain

Yates(イェ-ツ)補正
chisq.test(x)
2行×2列のクロス集計表のデータに対して行われる補正分割表から得られる.分割表から得られるカイ二乗値は跳び跳びの値しかとらないがχ2 分布は連続分布である。 このため,2×2 分割表の場合には連続性の補正をしたほうがよい。検出力は低下するが、より正確な検定が可能になる。

フィッシャーの正確確率検定
fisher.test(x)
分割表において,期待値が 5 以下の桝目が全体の桝目の 20% 以上あるか,期待値が 1 以下の桝目が 1 つでもある場合には,χ2 分布を利用する「独立性の検定」は不適当である。そのような場合に実施する検定。

http://aoki2.si.gunma-u.ac.jp/lecture/Cross/Fisher.html


例)HairEyeColorを使って練習してみましょう!

HairEyeColor
str(HairEyeColor)

table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
table [行数, 列数, 表数]

y<-addmargins(HairEyeColor) #周辺度数付き
#HairEyeColorの第一表(Sex = Male)にアクセスする
y[,,1]

f:id:yoshida931:20170913184150j:plain
それでは第一表(Sex = Male)のカイ二乗値を求めてみます

addmargins(HairEyeColor) #これで周辺度数を確認します
mh<-HairEyeColor[,,1] #第1表のみ抜き出します
xx<-matrix(c(56,143,34,46),4,4) #行合計マトリックス
yy<-t(matrix(c(98,101,47,33),4,4)) #列合計マトリックス
xy<- xx*yy/279 #各データの期待度数マトリックス
sum((mh-xy)^2/xy) #t値の算出
= 41.28029

行合計マトリックス
f:id:yoshida931:20170913193309j:plain
列合計マトリックス
f:id:yoshida931:20170913193322j:plain
各データの期待度数マトリックス
f:id:yoshida931:20170913193423j:plain
Rで確認してみます

chisq.test ( HairEyeColor[,,1] )

f:id:yoshida931:20170913193701j:plain

テキストファイルに検定結果を出力

テキストファイルに検定結果を出力
投稿日2017.9.8

ちょっとした小技

テキスト形式のファイルに検定結果を書き込む練習をします
sprintf:書式指定変換した出力を文字列に格納
cat:文字列を表示する基本的な関数, file=でファイルに書き込み

備忘録
\n 改行
%.3f 小数点3桁まで出力
append = F:ファイルに上書き
append = T:ファイルに追加

Rのサンプルirisを使います

x <- iris$Sepal.Length[51:75]
y <- iris$Petal.Length[101:125]
#noteというファイルに、xとyの平均を書き込みます
cat("x平均", mean(x), "y平均", mean(y), file = "note.txt") #noteというファイル名

f:id:yoshida931:20170908160200j:plain

sprintfの使い方
xの平均を小数点以下1桁まで、項目名「Xの平均」と付いた文字列にして取り出します

sprintf("xの平均 %.1f", mean(x))
"xの平均 6.0"

t検定を実施して、strで中身を確認してみます

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

Two Sample t-test

data: x and y
t = 2.1954, df = 48, p-value = 0.033
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.03131358 0.71268642
sample estimates:
mean of x mean of y
6.012 5.640


List of 9
$ statistic : Named num 2.2
..- attr(*, "names")= chr "t"
$ parameter : Named num 48
..- attr(*, "names")= chr "df"
$ p.value : num 0.033
$ conf.int : atomic [1:2] 0.0313 0.7127
..- attr(*, "conf.level")= num 0.95

t値、p値、95%信頼区間のみ取り出してテキストファイル「note」に書き出してみます

tv <- sprintf(" t値 %.4f", t[[1]]) 
pv <- sprintf("p値 %.4f", t[[3]]) 
cin <- "95%信頼区間"
ci <- sprintf("%.4f", t[[4]])
cat(tv,"\n", pv,"\n", cin, ci, file = "note.txt", append = F)

f:id:yoshida931:20170908162051j:plain
なぜか1行目が揃わないので、スペース入れてます?

参考と引用 Rプログラム (TAKENAKA's Web Page)

対応のないt検定の検定力分析(事後の分析)

検定力(1-β):帰無仮説が「偽」であるとき(母集団に差があるとき)にサンプルから有意差を得る確率
効果量(effect size):標準化された平均値差(検定により色々あります...d family, r family)
標本数:研究者が決定しなければならない
α(有意水準):慣例的に5%、1%で固定されている

検定力分析を実行する場合に、効果量の問題を解決しなければなりません.効果量を表す指標は数多く存在するので、まずどれを使うかが問題となります.今回はt検定の検定力分析に利用されている、次の二つの効果量について勉強していきます.

効果量=効果の大きさ
実験的操作の効果、変数間の関係性の強さを表す指標(大 中 小)

   Cohoen's d : グループごとの平均値の差を標準化したd family
   r : 変数間の関係の強さを示すr family

効果量の目安
   小  中  大
d 0.2 0.5 0.8
r 0.1 0.3 0.5

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

標本効果量(cohensのd効果量)の求め方
d = xの平均値とyの平均値の差の絶対値÷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}

例
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

d <- (abs(mean(x) - mean(y)))/pool
=1.33 (dは1を超える場合もある)

2標本平均の差の検定 - 統計学備忘録 since2016

Rを使ってt値から算出してみます

t.test(x,t,var.equal = T)
d = t*sqrt((length(x)+length(y))/length(x)*length(y))
  = 3.6449*sqrt(30/(15*15))
  = 1.33

Rを使って検定力を算出します

library(pwr)
pwr.t2n.test(n1=15,n2=15, d =1.33 , sig.level =0.05 , power = NULL,alternative = "greater")

t test power calculation
n1 = 15
n2 = 15
d = 1.33
sig.level = 0.05
power = 0.971769
alternative = greater

power = 0.97より、帰無仮説が偽の場合に有意差が検出される確率は97%です.高い検定力であると言えます.

標本効果量(r)の求め方
r= \sqrt{\frac{t^2}{t^2+自由度}}

t.test(x,y,var.equal = T)
t = 3.6449
r <- sqrt(t^2/(t^2+18))
=0.6516516

効果量は「大」でした

参考)水本篤; 竹内理. 研究論文における効果量の報告のために―基本的概念と注意点―. 2008.

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

指数分布

指数分布
投稿日2017.9.4


今回はRを使って苦手な指数関数に挑戦します

定義
f(x)を確率変数X確率密度関数とします
F(X)は確率変数Xの累積分布関数とします
F(X)=\int_{-\infty}^{\infty}f(x)dx\ =\ 1

離散型の変数で確率密度関数の復習をします
例)正確なサイコロ
確率変数X\ \ (1, 2, 3, 4, 5, 6)
f(x) = \frac{1}{6}\ \ \ \ (離散型一様分布)

y<-c(rep(1/6,6))
plot(y, type = "h")

f:id:yoshida931:20170904101052j:plain:w250

指数分布の密度関数
f(x) = λe^{-λx}\ \ \ (x ≧ 0)
f(x) = 0\ \ \ (x < 0)
f(x)は確率密度関数なのでx ≧ 0の範囲で積分したら1になります.

Rを使ってグラフを描いてみます(dexp=指数分布の密度関数)

curve(dexp(x, 1.0), from = 0, to = 10)

f:id:yoshida931:20170904095802j:plain:w250
指数分布の累積分布関数
 F(X)=\intλe^{-λx}dx
公式  ∫e^{cx}dx = \frac{1}{c}e^{-λx} より
 F(X)= 1 - e^{-λx}
Rを使ってグラフを描いてみます(dexp=指数分布の密度関数)
f:id:yoshida931:20170904132809j:plain:w250

指数分布の期待値
確率変数Xの期待値をE(X)、分散をV(X)とします
E(X) は確率の重み付き平均です.客観的な予想値になります.
サイコロ(離散型)の例で考えると・・・
期待値= 1*f(x) + 2*f(x) + ・・・ + 6*f(x) = 3.5
連続変数では、期待値を求めるために積分します.
E(X) = \int xf(x) dx
公式より \ \ \int xe^{cx}dx = \frac{e^{cx}}{c^2}(cx-1)
E(X) =[\frac{e^{-λx}}{λ^2}(-λx-1)]=\ \frac{1}{λ}

指数分布の分散
V(X) = E(X^2) - (E(X))^2
E(X^2)=\int x^2λe^{-λx}dx
公式より \ \ \int x^2e^{cx}dx\ =\ e^{cx}(\frac{x^2}{c}-\frac{2x}{c^2}+\frac{2}{c^3})
E(X^2)=[λe^{-λx}(\frac{x^2}{-λ}-\frac{2x}{-λ^2}+\frac{2}{-λ^3})]\ =\ \frac{2}{λ^2}
V(X)=\frac{2}{λ^2}-\frac{1}{λ^2}\ =\ \frac{1}{λ^2}

指数分布の確率密度関数と累積分布関数のグラフをまとめてRで描いてみます
λ:0.5, 1.0, 1.5
軸範囲を統一して同じグラフに書き加えるのでpar(new=T)を使用します
左に確率密度関数、右に累積分布関数を表示するように指示します

par(mfrow=c(1,2))
curve(dexp(x, 0.5), from = 0, to = 10, ylim = c(0, 1.5), col = 1,ylab ="確率密度")
par(new=T)
curve(dexp(x, 1.0), from = 0, to = 10, col = 2, ylab ="")
par(new=T)
curve(dexp(x, 1.5), from = 0, to = 10, col = 3,ylab ="")
curve(pexp(x, 0.5), from = 0, to = 10, ylim = c(0, 1.0), col = 1,ylab ="累積分布")
par(new=T)
curve(pexp(x, 1.0), from = 0, to = 10, col = 2, ylab ="")
par(new=T)
curve(pexp(x, 1.5), from = 0, to = 10, col = 3, ylab ="")
labels<-c("黒:λ=0.5","赤:λ=1.0","緑:λ=1.5")
legend("bottomright", legend = labels)
dev.off()

f:id:yoshida931:20170904201718j:plain:w600

1標本t検定のまとめ

t検定のまとめ
投稿日2017.1.27
最終更新日2017.9.1

t検定のまとめを復習していきます
Rでは次のような関数を使用してt検定を実行します.

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は不偏分散

例)次のデータを図で概観し,帰無仮説:母平均=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) #箱ひげ図
plot(x)    #散布図
#縦に散布
y <- c(rep("1", length(x)))
plot(x~y, ylim = c(145, 180), xlab = "", xaxt= "n")  #xaxt_x軸目盛り消去
#各図を横一列に並べて表示
par(mfrow = c(1, 3))
boxplot(x);plot(x);plot(x~y, ylim = c(145, 180), xlab = "", xaxt = "n") 
#作業終了後に
dev.off()

f:id:yoshida931:20170901124139j:plain:w550
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%の危険性(この結論が間違える危険性)があるのですが、標本xの母集団の平均が160であることが言えます.

対応のある2標本の場合
1標本問題と同じ原理です.治療前後の比較を考えてみましょう.
治療前のデータpre, 治療後のデータposとして治療の前後差を有意水準5%で検定せよ.

pre <- c(25, 27, 32, 37, 38, 39, 44, 45, 47, 48, 52, 54, 59, 60, 65, 65, 67, 68, 68)
pos <- c(41, 28, 32, 45, 40, 48, 44, 50, 53, 48, 60, 59, 60, 68, 66, 70, 66, 60, 65)

対応あるt検定は「前後差 = 0」つまり「 posの母平均 - preの母平均 = 0 」を帰無仮説とした1標本問題として分析します.2群の差の母平均が0になる確率が高ければ、2群の母平均には「差がない」という結論を得ます.

def <- pre - pos
t.test(def)

One Sample t-test

data: def
t = -2.7097, df = 18, p-value = 0.01435
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-5.8866381 -0.7449409
sample estimates:
mean of x -3.315789

95%信頼区間に「0」が含まれていません.また差の平均は-3.31となっています.有意水準0.05のt検定の結果は「有意差あり」ということになります.しかし、この平均-3.31(95%CI:-5.88~-0.74)という数値は検査内容により実質的な意味合いが異なります.最終的な判断は研究目的に応じてそれぞれの研究者が判断することになります.有意水準0.01のt検定を実行してみます.

def <- pre - pos
t.test(def, conf.level = 0.99)

One Sample t-test
data: def
t = -2.7097, df = 18, p-value = 0.01435
alternative hypothesis: true mean is not equal to 0
99 percent confidence interval:
-6.8380690 0.2064901
sample estimates:
mean of x -3.315789

有意差のみで結果を述べると、今度は「有意差あり」ということになります.αの設定は重要ですね.