なぜポアソン/二項回帰ではz値で、ガンマ回帰ではt値なの？

目次

・ポアソン/二項回帰ではz値、ガンマ回帰ではt値？

・「標準化した平均値」の従う分布は標準正規分布？t分布？

・ポアソン/二項分布は平均と分散が独立ではない

　　・ポアソン分布

　　・二項分布

・ガンマ分布は平均と分散が独立に決められる

・まとめ

・補足：ポアソン分布で分散を独立に推定した場合

・補足：二項分布で分散を独立に推定した場合

ポアソン/二項回帰ではz値、ガンマ回帰ではt値？

　これまでに、様々な確率分布の一般化線形モデルを扱ってきた。その中で、係数が有意か判断するのに用いられたWald統計量と呼ばれる値があった。具体的に見てみよう。例えば、ポアソン回帰では以下。

------------------------------------------------------

library(plotn)

y <- c(3, 3, 3, 2, 3, 3, 3, 4, 1, 1, 1, 5, 3, 1, 3, 5, 4, 5,

       4, 5, 2, 1, 1, 2, 4, 1, 4, 2, 0, 2, 3, 2, 1, 2, 3, 1,

       2, 3, 1, 6, 1, 2, 4, 3, 1, 4, 1, 2, 0, 3, 6, 1, 2, 4,

       6, 2, 1, 4, 2, 3, 1, 0, 0, 3, 1, 8, 3, 3, 1, 3, 1, 2, 

       2, 3, 3, 1, 1, 1, 0, 1, 4, 1, 4, 5, 0, 4, 2, 0, 8, 4,

       3, 2, 1, 5, 2, 2, 3, 3, 3, 2)

fit1 <- glm(y ~ 1, family = poisson(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##    Min      1Q  Median      3Q     Max  

## -2.249  -1.097  -0.346   0.287   2.735  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  0.92822    0.06287   14.76   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 115.76  on 99  degrees of freedom

## Residual deviance: 115.76  on 99  degrees of freedom

## AIC: 374.02

## 

## Number of Fisher Scoring iterations: 5

------------------------------------------------------

例えば、二項回帰では以下。

------------------------------------------------------

n <- c(18, 10, 16, 15, 11, 17, 20, 14, 18, 16, 15, 11, 

       19, 17, 18, 19, 17, 19, 16, 15, 16, 14, 13, 12, 

       10, 14, 20, 14, 10, 20, 19, 15, 11, 17, 12, 11, 

       17, 11, 16, 13, 14, 13, 11, 17, 17, 14, 13, 13, 

       19, 12, 11, 17, 18, 13, 12, 18, 19, 15, 15, 11, 

       14, 20, 12, 19, 13, 12, 17, 14, 17, 20, 10, 20, 

       18, 15, 13, 18, 14, 16, 14, 17, 11, 11, 13, 11, 

       20, 20, 17, 14, 14, 20, 14, 17, 20, 14, 13, 17, 

       15, 15, 20, 13)

y <- c(2, 2, 2, 2, 4, 8, 2, 4, 5, 1, 5, 1, 6, 4, 8, 8, 

       5, 4, 3, 6, 5, 5, 3, 4, 2, 6, 4, 2, 5, 5, 6, 2, 

       1, 3, 2, 3, 5, 2, 3, 6, 7, 3, 3, 2, 4, 3, 2, 4, 

       7, 4, 3, 4, 7, 4, 5, 5, 4, 2, 5, 6, 3, 7, 4, 6, 

       5, 5, 3, 5, 5, 6, 2, 4, 2, 1, 2, 5, 1, 7, 3, 5, 

       2, 5, 4, 3, 0, 5, 6, 4, 6, 8, 3, 6, 5, 3, 4, 3, 

       5, 3, 4, 3)

fit1 <- glm(cbind(y, n-y) ~ 1, family = binomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = cbind(y, n - y) ~ 1, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.5103  -0.6798   0.0569   0.5841   1.9629  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept) -1.01947    0.05811  -17.54   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for binomial family taken to be 1)

## 

##     Null deviance: 101.97  on 99  degrees of freedom

## Residual deviance: 101.97  on 99  degrees of freedom

## AIC: 387.3

## 

## Number of Fisher Scoring iterations: 4

------------------------------------------------------

そして、ガンマ回帰では以下。

------------------------------------------------------

y <- c(4.78, 1.41, 0.24, 2.56, 2.74, 0.05, 2.42, 0.74, 2.1, 

       0.08, 1.39, 1.53, 0.39, 7.54, 0.31, 3.01, 2.85, 1.22, 

       3.19, 2.96, 4.03, 1.29, 2.59, 3.72, 3.04, 0.001, 1.15, 0.4,

       1.13, 1.11, 1.32, 0.83, 1.47, 0.13, 2.5, 1.52, 0.54,

       4.84, 4.43, 3.45, 2.74, 3.02, 0.9, 0.69, 0.36, 1, 6.39, 

       1.08, 0.38, 1.13, 0.32, 4.23, 0.86, 1.1, 2.18, 6.48, 0.22, 

       4.27, 7.49, 7.33, 0.86, 2.32, 3.28, 2.2, 3.44, 8.45, 1.68, 

       0.61, 2.07, 0.59, 2.76, 2.23, 0.8, 1.65, 0.01, 1.47, 0.68, 

       3.65, 0.58, 0.48, 6.05, 0.8, 0.51, 2.37, 0.51, 0.81, 1.42, 

       1.29, 0.42, 1.57, 0.64, 2.02, 4.1, 0.1, 1.44, 1.35, 5.14, 3.29, 0.26, 1.78)

fit1 <- glm(y ~ 1, family = Gamma(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = Gamma(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.6455  -0.9384  -0.3410   0.3750   1.8156  

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  0.73645    0.09019   8.166 1.06e-12 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for Gamma family taken to be 0.8133505)

## 

##     Null deviance: 107.33  on 99  degrees of freedom

## Residual deviance: 107.33  on 99  degrees of freedom

## AIC: 352.16

## 

## Number of Fisher Scoring iterations: 5

------------------------------------------------------

以上のように、ポアソン/二項回帰では、Wald統計量がzと表示され、ガンマ回帰ではtと表示されている。これは、ポアソン/二項回帰ではWald統計量が標準正規分布に従い（つまりz検定）、ガンマ回帰ではt分布に従う（つまりt検定）と考えてp値を計算しているということである。発想は、線形回帰の時にt分布を使ったことと同じである。

　しかし、気になるのは、なぜ、ポアソン/二項回帰では、標準正規分布を使うのだろうか。普通に考えたらt分布に従うと考えるべきではないのか。あるいは、GLMだからというならすべてz値に統一するべきではないのだろうか。今回は、これらのテーマを深堀りしよう。

「標準化した平均値」の従う分布は標準正規分布？t分布？

　復習になるが、母集団の確率分布が”だいたいどんなものでも”、母集団の平均がμ、分散がσ^2であれば、そこから得られたサンプルサイズnの確率変数Xの平均値‾Xは、平均μ、分散(σ^2)/nの正規分布に従う。これが中心極限定理である。シミュレーションによって、多数の平均値を生成することで、これが正しいとわかる。平均値‾Xは、平均μ、分散(σ^2)/nの正規分布に従うということは、標準化した平均値z = (‾X - μ)/(σ/√n)は、標準正規分布に従うということである。

　実際はデータは1つ、つまりzも1つだ。そこから、この平均値がどんな分布に従うかを、予測する必要がある。これも、もとの母集団の分散さえ分かっていれば、そのzの信頼区間を求めることができ、自動的に平均値の信頼区間も求めることができる。ところが、実際問題として母分散はわからない。そこで、z値の母分散σ^2の部分を母分散の不偏推定量である不偏分散Uで置き換えることを考える。不偏分散への置き換え後の値をt値とよび、これは標準正規分布ではなく自由度n-1のt分布に従うのだった。母分散は唯一の値だが、不偏分散は標本ごとにばらつくので、それが標準正規分布よりもt分布のばらつきが大きい状態を生み出している。

　一応、シミュレーションを使って、z値が従う分布を確認してみよう。以下のように、母分散を使って平均値を標準化すると、確かに標準正規分布に従うことが確認できる。ヒストグラムと累積密度関数を使って確認する。

------------------------------------------------------

cdf <- function(p, v){

  sapply(v, function(x) sum(p < x))/length(p)

}

m <- 1 #平均

s <- 2 #標準偏差

n <- 5 #サンプルサイズ

d <- NULL

for (i in 1:10000){

  d <- c(d, (mean(rnorm(n, mean = m, sd = s)) - m)/(s/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "z-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax)) #図1の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "z-value", ylab = "cdf")

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red")) #図2の図示

------------------------------------------------------

　z値の母分散を不偏分散に置き換えたt値は、下記のように、標準正規分布ではなく、自由度5 - 1 = 4のt分布に従うことがわかる。

------------------------------------------------------

d <- NULL

for (i in 1:10000){

  x <- rnorm(n, mean = m, sd = s)

  d <- c(d, (mean(x) - m)/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax)) #図3の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図4の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

　線形回帰では、(係数推定値)/(推定値の分散推定値)はt分布に従うことがわかっているので、Wald統計量にあたる値はtとして表記されている。では、ポアソン分布、二項分布、ガンマ分布を使った回帰では、(係数推定値)/(推定値の分散推定値)はどんな分布に従うのだろうか。そんなの、中心極限定理の発想から考えれば、t分布になると考えるだろう。しかし、分布の推定を行う上で、ある性質を考慮すると、必ずしもそうなるとは限らないのである。

ポアソン/二項分布は平均と分散が独立ではない

　ポアソン分布は平均＝分散が成り立つ確率分布だった。二項分布は、試行回数をnとしたとき、平均がμであれば、成功確率p = μ/nが成り立ち、分散はμ(1 - μ/n)が成り立つ。これらの確率分布に共通することは、平均が決まれば分散も決まるという性質である。この性質がGLMにおけるWald統計量の検定に標準正規分布を使う理由となっている。

・ポアソン分布の場合

　まず、ポアソン分布の場合を確認しよう。統計モデリングにおいて、分布がポアソン分布に従うと仮定した場合、平均を求めたら、自動的にポアソン分布の分散を求めたことになり、独立に分散を計算する必要がない。この条件の下で標本平均を分散で標準化した統計量は以下のようになる。どんな分布に従うかわからないし、分散未知の場合を想定しているはずなのでいったんt?（= Wald統計量）とでも置いておこう。

この条件の元、t?をたくさん求めて、どんな分布に従うか確認しよう。

------------------------------------------------------

n <- 5

l <- 2

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/sqrt(mean(x/n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図5の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図6の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

一概にどちらの分布に近いとか、そういうわけでない。平均値が小さい場合はこのような状況になる。では、平均値を大きくしてみよう。

------------------------------------------------------

l <- 5

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/sqrt(mean(x/n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図7の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図8の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

すると、徐々に赤い線（標準正規分布）のほうにt?の分布が近づいていることがわかる。さらに平均値を大きくする。

------------------------------------------------------

l <- 10

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/sqrt(mean(x/n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図9の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図10の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

------------------------------------------------------

l <- 30

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/sqrt(mean(x/n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 50)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 50, ylim = c(0, ymax)) #図11の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図12の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

すると、t?値は、t分布よりも標準正規分布により近づいた。すなわち、ポアソン分布に従うと仮定した場合、標準化した平均値 = Wald統計量を検定するのであれば標準正規分布を使う方がより適切であることがわかる。もし、t分布に従うと考えた場合、t分布は裾が重いので、検出力が低下することになる。

　これまで、サンプルサイズ5で検証してきたが、サンプルサイズが大きい場合がどうなるか確認する。

------------------------------------------------------

l <- 2

n <- 30

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/sqrt(mean(x/n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図13の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図14の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

すると、ほとんど正規分布になることがわかる。そもそも、サンプルサイズが大きければ、t分布自体も標準正規分布に近づくので、どちらに近いということもなくなる。

　以上をまとめると、Wald統計量が標準正規分布に近似できるのは以下の2パターンである。

・母集団のポアソン分布の平均値が十分大きい場合

・サンプルサイズが十分大きい場合

前者は、つまり母集団が正規分布に近づくということである。後者は、ポアソン分布の再生性に関連している。平均λ1およびλ2のポアソン分布からそれぞれ確率変数X1とX2を得たとき、その和X1 + X2は平均λ1 + λ2のポアソン分布に従う。もしたくさんの変数を足せば、その分だけポアソン分布の平均は大きくなるから、変数の和の従うポアソン分布は正規分布に近づく。これがポアソン分布における中心極限定理であり、変数の平均も同様に正規分布に従うようになる、ということである。

　では、もし、独立に分散を求めたらどうなるだろうか。その場合は、今度はt分布に収束していくことがわかる（追記）。

・二項分布の場合

　次に、二項分布の場合を確認しよう。統計モデリングにおいて、分布が二項分布に従うと仮定した場合、平均を求めたら、自動的に二項分布の分散を求めたことになり、独立に分散を計算する必要がない。この条件の下で標本平均を分散で標準化した統計量は以下のようになる。サンプルサイズnと、二項分布における試行回数nbは違う値であることに注意しよう。どんな分布に従うかわからないし、分散未知の場合を想定しているはずなのでいったんt?（= Wald統計量）とでも置いておこう。

この条件の元、t?をたくさん求めて、どんな分布に従うか確認しよう。今回、成功確率は1/5に固定してある。

------------------------------------------------------

n <- 5

s <- 10

p <- 1/5

d <- NULL

for (i in 1:10000){

  x <- rbinom(n, size = s, prob = p)

  p_hat <- mean(x)/s

  d <- c(d, (mean(x) - s*p)/(sqrt(s*p_hat*(1-p_hat)/n)))

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図15の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図16の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

ポアソン分布の時と同様に、一概にどちらの分布に近いとか、そういうわけでない。平均値＝試行回数×成功確率が小さい場合はこのような状況になる。では、試行回数を大きくすることで平均値を大きくしてみよう。

------------------------------------------------------

s <- 20

p <- 1/5

d <- NULL

for (i in 1:10000){

  x <- rbinom(n, size = s, prob = p)

  p_hat <- mean(x)/s

  d <- c(d, (mean(x) - s*p)/(sqrt(s*p_hat*(1-p_hat)/n)))

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図17の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図18の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

------------------------------------------------------

s <- 50

p <- 1/5

d <- NULL

for (i in 1:10000){

  x <- rbinom(n, size = s, prob = p)

  p_hat <- mean(x)/s

  d <- c(d, (mean(x) - s*p)/(sqrt(s*p_hat*(1-p_hat)/n)))

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図19の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図20の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

すると、t?値は、t分布よりも標準正規分布により近づいた。すなわち、二項分布に従うと仮定した場合、標準化した平均値 = Wald統計量を検定するのであれば標準正規分布を使う方がより適切であることがわかる。もし、t分布に従うと考えた場合、t分布は裾が重いので、検出力が低下することになる。

　これまで、サンプルサイズ5で検証してきたが、サンプルサイズが大きい場合がどうなるか確認する。試行回数も5まで減らす。

------------------------------------------------------

n <- 30

s <- 5

p <- 1/5

d <- NULL

for (i in 1:10000){

  x <- rbinom(n, size = s, prob = p)

  p_hat <- mean(x)/s

  d <- c(d, (mean(x) - s*p)/(sqrt(s*p_hat*(1-p_hat)/n)))

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value?", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図21の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value?", ylab = "cdf") #図22の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

やはり、ほとんど正規分布になることがわかる。そもそも、サンプルサイズが大きければ、t分布自体も標準正規分布に近づくので、どちらに近いということもなくなる。

　以上をまとめると、Wald統計量が標準正規分布に近似できるのは以下の2パターンである。

・母集団の二項分布の平均値と分散が十分大きい場合

・サンプルサイズが十分大きい場合

前者は、つまり母集団が正規分布に近づくということである。後者は、二項分布の再生性に関連している。成功確率pで、試行回数n1およびn2二項分布からそれぞれ確率変数X1とX2を得たとき、その和X1 + X2は成功確率p、試行回数n1 + n2の二項分布に従う。もしたくさんの変数を足せば、その分だけ二項分布の平均と分散は大きくなるから、変数の和の従う二項分布は正規分布に近づく。これが二項分布における中心極限定理であり、変数の平均も同様に正規分布に従うようになる、ということである。

　では、もし、独立に分散を求めたらどうなるだろうか。その場合は、今度はt分布に収束していくことがわかる（追記）。

ガンマ分布は平均と分散が独立に決められる

　ポアソン分布や二項分布と異なり、ガンマ分布は平均と分散が独立に決めることができる。ガンマ分布の形状パラメータk、尺度パラメータθがあった時、平均はkθ、分散はkθ^2である。同じ文字が使われているじゃないかと感じるかもしれないが、未知の文字は2つあるので、平均と分散は独立に決めることが可能である。このとき、平均がわかったとしても、分散はわからない。ゆえに両方とも推定が必要になる。分散を推定する必要があるので、標準化した平均値はt分布になる。

　ガンマ回帰の元、tをたくさん求めて、どんな分布に従うか確認しよう。先に述べておくと、ガンマ分布の中心極限定理の性質から、形状パラメータが大きくなるほどWald統計量はt分布に近づくが、尺度パラメータが大きくなってもt分布には近づかない。まずは、尺度パラメータを固定した場合。

------------------------------------------------------

n <- 5

sh <- 2

sc <- 2

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax))

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue")) #図23の図示

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図24の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

ポアソン分布や二項分布の時と同様に、一概にどちらの分布に近いとか、そういうわけでない。形状パラメータが小さい場合はこのような状況になる。では、形状パラメータを大きくしてみよう。

------------------------------------------------------

sh <- 4

sc <- 2

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax))

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue")) #図25の図示

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図26の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

----------------------------------------------------

sh <- 32

sc <- 2

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax))

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue")) #図27の図示

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図28の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

すると、t値は、標準正規分布よりもt分布により近づいた。すなわち、ガンマ分布に従うと仮定した場合、標準化した平均値 = Wald統計量を検定するのであればt分布を使う方がより適切であることがわかる。もし、標準正規分布に従うと考えた場合、標準正規分布は裾が軽いので、危険率が増加することになる。

　もし、尺度母数を大きくするとどうなるか、確認してみよう。

------------------------------------------------------

sh <- 2

sc <- 4

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax))

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue")) #図29の図示

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図30の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

-----------------------------------------------------

sh <- 2

sc <- 32

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax))

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue")) #図31の図示

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図32の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

形状パラメータが尺度母数に比較して相対的に小さい場合、サンプルサイズが大きくなればt分布（実質上、標準正規分布に近づく。）

------------------------------------------------------

n <- 500

sh <- 1/2

sc <- 2

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax)) #図33の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図34の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

------------------------------------------------------

n <- 500

sh <- 2

sc <- 32

d <- NULL

for (i in 1:10000){

  x <- rgamma(n, shape = sh, scale = sc)

  d <- c(d, (mean(x) - (sh*sc))/(sd(x)/sqrt(n))) 

}

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 100)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 100, ylim = c(0, ymax)) #図35の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図36の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

ガンマ分布の正規分布への収束が遅いため、サンプルサイズ500でも歪みがまだ見えるが、ほとんどt分布になることがわかる。そもそも、サンプルサイズが大きければ、t分布自体も標準正規分布に近づくので、どちらに近いということもなくなる。

　以上をまとめると、Wald統計量がt分布に近似できるのは以下の2パターンである。

・母集団のガンマ分布の形状パラメータが十分大きい場合

・サンプルサイズが十分大きい場合

前者は、つまり母集団が正規分布に近づくということである。後者は、ガンマ分布の再生性に関連している。尺度パラメータθで、形状パラメータk1およびk2のガンマ分布からそれぞれ確率変数X1とX2を得たとき、その和X1 + X2は尺度パラメータθ、形状パラメータk1 + k2のガンマ分布に従う。もしたくさんの変数を足せば、その分だけガンマ分布の形状パラメータは大きくなるから、変数の和の従うガンマ分布は正規分布に近づく。これがガンマ分布における中心極限定理であり、変数の平均も同様に正規分布に従うようになる、ということである。

まとめ

　GLMにおいて、Wald統計量がポアソン/二項回帰ではz値、ガンマ分布ではt値で表示される理由について、標準化された平均値の収束先に注目して紹介した。これらの違いは、平均と分散が互いに独立に推定されるかどうかに依存していることがわかっただろう。

　他の統計解析においても、zと表記されていれば、その統計量が標準正規分布に近づくことを示している。例えば、二項検定において漸近検定を行う場合は統計量はしばしばzと表示される。この理由は、二項分布に従う確率変数の挙動を考えれば納得だろう。二項検定同じような考え方で、ポアソン分布を使った検定、ポアソン検定というものがある。ポアソン検定で使われる漸近検定でも同様に統計量はzである。

補足：ポアソン分布で分散を独立に推定した場合

　試しにポアソン分布で、分散を独立に推定した場合にどうなるかを紹介しよう。以下のようにt分布に近づくことが見て取れる。

------------------------------------------------------

n <- 5

l <- 2

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/(sd(x)/sqrt(n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図37の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図38の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

------------------------------------------------------

l <- 30

d <- NULL

for (i in 1:10000){

  x <- rpois(n, lambda = l)

  d <- c(d, (mean(x) - l)/(sd(x)/sqrt(n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図39の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図40の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

補足：二項分布で分散を独立に推定した場合

　試しに二項分布で、分散を独立に推定した場合にどうなるかを紹介しよう。以下のようにt分布に近づくことが見て取れる。

------------------------------------------------------

n <- 5

s <- 10

p <- 1/5

d <- NULL

for (i in 1:10000){

  x <- rbinom(n, size = s, prob = p)

  d <- c(d, (mean(x) - s*p)/(sd(x)/sqrt(n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図41の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図42の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

-----------------------------------------------------

n <- 5

s <- 50

p <- 1/5

d <- NULL

for (i in 1:10000){

  x <- rbinom(n, size = s, prob = p)

  d <- c(d, (mean(x) - s*p)/(sd(x)/sqrt(n))) 

}

d <- d[!is.infinite(d)]

d <- d[!is.nan(d)]

xx <- seq(-10,10, length = 200)

ymax <- max(hist(d, plot = F, breaks = 20)$density, dnorm(xx, mean = 0, sd = 1))

histn(d, xlab = "t-value", ylab = "頻度", freq = F, breaks = 20, ylim = c(0, ymax)) #図43の図示

overdraw(points(xx, dnorm(xx, mean = 0, sd = 1), type = "l", col = "red"),

         points(xx, dt(xx, df = n-1), type = "l", col = "blue"))

v <- seq(min(d), max(d), length = 200)

cd <- cdf(d,v)

plotn(cd ~ v, type = "l", xlab = "t-value", ylab = "cdf") #図44の図示

overdraw(points(v, pnorm(v, mean = 0, sd = 1), type = "l", col = "red"),

         points(v, pt(v, df = n-1), type = "l", col = "blue"))

------------------------------------------------------

参考

https://stats.stackexchange.com/questions/56066/wald-test-in-regression-ols-and-glms-t-vs-z-distribution