母分散を知っているか否かで検定が変わる？

　統計の初学者だったころ（今もだけど）、私はある問題に大いに頭を悩ませた。それが、「母分散が既知ならz検定、未知ならt検定」というフレーズだ。実用上、私たちは母分散を知りえないから、選択としては実質t検定の1択で、困りはしないのだが、それでもこのフレーズを聞くたびに頭をかしげたものだ。第一に母分散が既知であることはふつうあり得ない。じゃあ、Z検定がなぜ方法論として確立しているのかわからなかった。第二にそもそも、検定が変わるということは検定統計量が従う帰無分布が変わる、ということだ。標本が従う分布が変わるから検定統計量の分布が変わるのならまだわかる。また、検定統計量の計算が変わるから帰無分布が変わる、というのもわかる。けれども、Z検定もt検定も、標本の分布はともに正規分布に従うと主張するし、ほとんど計算は同じである（実は使われてる指標の名前が似ているせいであり、よく似てるけど別の性質の統計量を計算していることになるのだけど）。それなのに、私が母分散を知っているか否かで、検定統計量が従う分布が変わるというのか。まあ、実際にそうなんだけど、「母分散が既知ならz検定、未知ならt検定」というフレーズは大いに混乱を招く、と思っている。ただ、Z検定とt検定の違いは、実データの母分散がわからない以上、確認することは困難だ。本項では、シミュレーションという武器を使って、これらの検定の違いを理解することを目指す。

再生性と確率変数の1次結合

　さて、一応の前段階として、一つある事実を認めておく必要がある。そしてこの性質は、今後暗に仮定することがあるので覚えていてほしい。例えば、今回データが従う確率分布である正規分布について考える。ある確率変数Xは平均μ1、分散σ1^2の正規分布に従い（これらのパラメータを持つ正規分布から確率変数Xが生成される、と読み替えてもよい）、確率変数Yは平均μ2、分散σ2^2の正規分布に従うとする。これらの確率変数の和X+Yは、同様に正規分布に従い、その平均はμ1+μ2、分散はσ1^2+σ2^2である。図示すると以下のような感じだ。

このような性質を持つことを「確率分布が再生性reproductive propertyを有する」、という。上記の例は、正規分布に限った定義だが、より一般には、もとの確率変数が従う確率分布と、確率変数の和が従う確率分布が、パラメータだけが変わった同一の確率分布に従うことを指す（例えば、正規分布に従う確率変数の和がガンマ分布などに変わったりしないということである）。

　また、再生性のより一般には、正規分布に従う確率変数の1次結合aX+bY（aとbは実数定数）も、やはり正規分布に従い、その平均はaμ1+bμ2、分散は(aσ1)^2+(bσ2)^2である。図示すれば下記のようになるだろう。

これらの性質の証明は補足にて行う。再生性は重要な確率分布ではしばしばみられる性質であり、正規分布をはじめ、ポアソン分布、二項分布、ガンマ分布、負の二項分布などでもみられる。一方、指数分布などでは再生性は見られない。

　なぜ突如、この話をしたかという言うと、今回は正規分布に従う確率変数の差X - Yを考えたい。上記の情報に基づけば、正規分布の確率変数の差は、平均はμ1 - μ2、分散はσ1^2+σ2^2に従う。いったん、このことを覚えておこう。

正規分布に従う確率変数の差のシミュレーション

　では、実際に正規分布に従う確率変数の差をシミュレーションしてみよう。下記のようなスクリプトを実行する。一つは平均1、分散1の正規分布、もう一つは平均3、分散4の正規分布に従う。これらの確率変数の差は、平均1 - 3 = -2、分散1 + 4 = 5の正規分布に従うと期待される。

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

library(plotn)

n <- 100000

X <- rnorm(n = n, mean = 1, sd = 1)#平均1、分散1の正規分布に従う乱数を100000個生成

Y <- rnorm(n = n, mean = 3, sd = 2)#平均3、分散4の正規分布に従う乱数を100000個生成

dif <- X - Y #正規分布から生成された値の差

histn(dif, xlab = "X - Y")#図1の図示

mean(dif)#平均

## [1] -2.008692

var(dif)*(n-1)/n #標本分散

## [1] 4.990047

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

確かに生成された確率変数の差は、上記で議論された性質を満たしている。平均や分散をいろいろ変更して試してみるとよい。

標準正規分布とZ検定

　さて、確率変数Xは平均μ1、分散σ1^2、確率変数Yは平均μ2、分散σ2^2の正規分布に従うとき、確率変数XやYの平均が従う確率分布を考えよう。それぞれのサンプルサイズをn1とn2とすると、中心極限定理より、確率変数Xは平均μ1、分散(σ1^2)/n1、確率変数Yは平均μ2、分散(σ2^2)/n2の正規分布に従うと考えられる。

　上記の議論を総合して、確率変数の平均の差E(X) - E(Y)は、平均はμ1 - μ2、分散は(σ1^2)/n1+(σ2^2)/n2に従う、と考えられる。ではこれを確認しよう。併せて、従うであろう正規分布も描画しておくことにする。下記の設定に基づくと平均の差は、平均1 - 4 = -3、分散(2^2)/4 + (3^2)/9 = 2の正規分布に従う。

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

n <- 100000

n1 <- 4 #標本Xのサンプルサイズ

n2 <- 9 #標本Yのサンプルサイズ

m1 <- 1 #標本Xの平均

m2 <- 4 #標本Yの平均

s1 <- 2 #標本Xの標準偏差

s2 <- 3 #標本Yの標準偏差

dif <- NULL

for(i in 1:n){

  X <- rnorm(n = n1, mean = m1, sd = s1)#標本Xの生成

  Y <- rnorm(n = n2, mean = m2, sd = s2)#標本Yの生成

  dif <- c(dif, mean(X) - mean(Y))#平均値の差

}

histn(dif, xlab = "E(X) - E(Y)", freq = F)#図2の図示

xx <- seq(min(dif), max(dif), length = 200)

overdraw(points(xx, 

                dnorm(xx, mean = m1-m2,sd = sqrt((s1^2)/n1+(s2^2)/n2)), 

                col = "red", type = "l"))#期待される正規分布

mean(dif)

## [1] -2.993786

var(dif)*(n-1)/n

## [1] 2.0009

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

確かに期待通りの正規分布になった。では、この確率変数の差の分布に関して、分散が1になるように調整したいとする。これは、下記のように確率変数の差E(X) - E(Y)を√{(σ1^2)/n1+(σ2^2)/n2}で割ればよい。この検定統計量をzと呼んでいる。

実際にこれで問題ないことは下記のスクリプトで、分散1の正規分布が生成されることから確認できる。

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

z <- dif/sqrt((s1^2)/n1 + (s2^2)/n2)

histn(z, xlab = "z", freq = F)#図3の図示

xx <- seq(min(z), max(z), length = 200)

overdraw(points(xx, 

                dnorm(xx, mean = (m1-m2)/sqrt((s1^2)/n1+(s2^2)/n2), sd = 1),

                col = "red", type = "l"))

var(z)*(n-1)/n

## [1] 1.00045

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

特に、μ1 - μ2 = 0、すなわち、μ1 = μ2であるとき、生成される正規分布は平均0、分散1であり、この正規分布を標準正規分布standard normal distributionと呼ぶ。

　平均値の差が0の場合もシミュレーションしておこう。同時に、後述で使うため累積分布関数cumulative distribution function (CDF)の計算を下記のように定義する（詳しくはリンク先の補足）。確率分布が一致すれば、標準正規分布のCDF（R上ではpnorm関数を使って計算）とヒストグラムから計算されるCDFが一致する。ヒストグラムを確認するよりも、CDFを確認するほうがより正確に判断できる。

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

n <- 100000

n1 <- 4

n2 <- 4

m1 <- 1

m2 <- 1

s1 <- 2

s2 <- 3

z <- NULL

for(i in 1:n){

  X <- rnorm(n = n1, mean = m1, sd = s1)

  Y <- rnorm(n = n2, mean = m2, sd = s2)

  z <- c(z, (mean(X) - mean(Y))/sqrt((s1^2)/n1 + (s2^2)/n2))

}

xx <- seq(min(z), max(z), length = 200)

ylim <- range(dnorm(xx, mean = 0, sd = 1))

histn(z, xlab = "z", freq = F, ylim = ylim)#図4の図示

overdraw(points(xx, 

                dnorm(xx, mean = 0, sd = 1), 

                col = "red", type = "l"))

cdf <- function(sample, x){

  sapply(x, FUN = function(tmp){

    sum(sample < tmp)/length(sample)

  })

}

plotn(pnorm(xx, mean = 0, sd = 1) ~ xx, type = "l", 

      xlab = "X", ylab = "CDF")#図5の図示

overdraw(points(xx, cdf(z, xx), type = "l", col = "red"))

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

確かに、2つのCDFはぴったり一致することがわかる。

　さて、この検定統計量zは、検定統計量tと類似した性質を持つことがわかる。つまり、検定統計量zは、(1) 2つの確率変数の平均の差が大きいほど、(2) 2つの確率変数の分散が小さいほど、(3) サンプルサイズが大きいほど、zは大きな値を示し、2つの確率変数が異なる確率分布から生成されたと期待できるだろう。帰無仮説のもと、つまりμ1 = μ2における検定統計量zは上記の議論から、標準正規分布に従うと言えるので、標準正規分布を使って仮説検定ができそうだ。この方法の検定をZ検定Z-testと呼ぶ。Z検定の発想は下記の通りだ。

　実際、特定の条件のもと、Z検定は確かに機能する。例えば、下記のデータを考える。標本XおよびYは、それぞれ分散1および4の正規分布から生成されたことがわかっている。z.test関数を定義して使ってみよう。z.test関数内で使っているabs関数は、絶対値関数であり、平均の差を正にするように調整している。こうすると場合分けがなくなり、計算が楽になるので使っている。z.test関数は、2つの標本xとyを引数とし、さらにxとyの母標準偏差＝√母分散も引数とする。また、検定結果のP値だけを返すように定義してある。

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

z.test <- function(x, y, sx, sy){

  z <- abs(mean(x)-mean(y))/sqrt((sx^2)/length(x) + (sy^2)/length(y))

  pnorm(z, mean = 0, sd = 1, lower.tail = F) * 2

}

X <- c(1.4, 2.5, 0.7, 3.4, 1.3, 1, 2.8, -0.1, 2.4, 2)

Y <- c(1.7, 0.1, 0.4, -0.9, 2.2, 3.9, 2.3, 3.8, 3.5, 3.1)

z.test(X, Y, 1, 2)

## [1] 0.7025818

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

今回P = 0.70であり、平均値に有意差はなさそうだ。もう一例出そう。標本XおよびYは、それぞれ分散1および4の正規分布から生成されたことがわかっている。

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

X <- c(0.5, -0.6, 3.1, 2, 3.1, 0.9, 1.3, 1.5, 1.6, 3.1)

Y <- c(2.4, 2.6, 4.3, 1.8, 5.1, 2.8, 4.7, 2.1, 4.5, 3.2)

z.test(X, Y, 1, 2)

## [1] 0.01620954

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

今回は、P = 0.016で有意差があった。

　一見すれば、t検定と同様に2標本の検定ができるような気がするのに、なぜ世の中はt検定ばかりでz検定を見かけることが少ないのだろうか。それは、後述するように、実用上の明確な弱点が存在するからだ。z検定のほうがより伝統的な検定方法であった。しかし、より近代にGossetがz検定の弱点を克服し、確立した方法が（Studentの）t検定であり、今日に至るまで使われ続けているのである（普通そこは、Gossetのt検定だろ、という感じがするが、このStudentはGossetのペンネーム。当時勤めていた会社の関係でGossetは実名で論文を出せる立場になかった）。

統計量zは正規分布に従わない！？ー母分散が未知であるという問題

　z検定の明確な問題点は、「実用上、仮に検定統計量zを計算してもzは正規分布に従わない」、これに尽きてしまう。「いや、さっきまで検定統計量zは正規分布に従う」、とそう言っていたじゃないか。ところで思い出してほしいのだが、私たちは実際のところ、標本が従っている正規分布の母平均も母分散も知らない。母平均も母分散も、標本から推定するしかないのであった。これに対して、さっき使ったスクリプトを確認してみよう。

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

n <- 100000

n1 <- 4

n2 <- 4

m1 <- 1

m2 <- 1

s1 <- 2

s2 <- 3

z <- NULL

for(i in 1:n){

  X <- rnorm(n = n1, mean = m1, sd = s1)

  Y <- rnorm(n = n2, mean = m2, sd = s2)

  z <- c(z, (mean(X) - mean(Y))/sqrt((s1^2)/n1 + (s2^2)/n2))

}

xx <- seq(min(z), max(z), length = 200)

ylim <- range(dnorm(xx, mean = 0, sd = 1))

histn(z, xlab = "z", freq = F, ylim = ylim)#図4の図示

overdraw(points(xx, 

                dnorm(xx, mean = 0, sd = 1), 

                col = "red", type = "l"))

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

データの生成過程で私たちは正規分布を仮定したから、当然、その時の平均と分散も知っている。スクリプトを確認すると、検定統計量zを計算するときに、母集団の情報を使ってしまっている。繰り返すが、実務上は母分散はわからない。そこで、母分散をその不偏推定量である不偏分散に変えたらどうなるだろうか？　不偏分散は手元の標本から計算可能であった。母分散を不偏分散に変えて計算した検定統計量をz'としておこう。

　では検定統計量z'はどんな分布になるか、シミュレーションしてみよう。併せてCDFも確認しておく。不偏分散はvar関数を使えば計算できる。

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

n <- 100000

n1 <- 4

n2 <- 4

m1 <- 1

m2 <- 1

s1 <- 2

s2 <- 3

z <- NULL

for(i in 1:n){

  X <- rnorm(n = n1, mean = m1, sd = s1)

  Y <- rnorm(n = n2, mean = m2, sd = s2)

  z <- c(z, (mean(X) - mean(Y))/sqrt(var(X)/n1 + var(Y)/n2))#母分散ではなく不偏分散に置き換えた

}

xx <- seq(min(z), max(z), length = 200)

ylim <- range(dnorm(xx, mean = 0, sd = 1))

histn(z, xlab = "z'", freq = F, ylim = ylim)#図6の図示

overdraw(points(xx, 

                dnorm(xx, mean = 0, sd = 1), 

                col = "red", type = "l"))

plotn(pnorm(xx, mean = 0, sd = 1) ~ xx, type = "l", 

      xlim = c(-5,5), xlab = "X", ylab = "CDF")#図7の図示

overdraw(points(xx, cdf(z, xx), type = "l", col = "red"))

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

すると、検定統計量z'は、標準正規分布に従うことが期待されたのに、なんだか分布の両端の裾が重たい分布になっている。こうなってしまうことが最大の問題点である。つまり、母分散がわからないので不偏分散を代わりに使ってしまうと、この検定統計量zもどきは標準正規分布に従わないのである。

　これによってどんな問題が生じるかというと、z'は標準正規分布から期待されるよりも外れ値が出やすい。これによって、z'が標準正規分布に従うと考えて、標準正規分布を帰無分布として検定を行うと有意になりやすい＝危険率が高まるのである。

　これもシミュレーションで確認することができる。上記のz.test関数のsx、sy引数に不偏分散の平方根を入力すればz'を計算して、標準正規分布を帰無分布とした解析になることに相当する。同時に、母分散を入力して「正しくZ検定した」場合も併せて、P値を収集してみよう。

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

n <- 100000

n1 <- 4

n2 <- 4

m1 <- 2

m2 <- 2

s1 <- 1

s2 <- 2

p1 <- NULL

p2 <- NULL

for(i in 1:n){

  X <- rnorm(n = n1, mean = m1, sd = s1)

  Y <- rnorm(n = n2, mean = m2, sd = s2)

  p1 <- c(p1, z.test(X, Y, s1, s2))

  p2 <- c(p2, z.test(X, Y, sd(X), sd(Y)))

}

histn(p1, xlab = "P value")#図8の図示

histn(p2, xlab = "P value")#図9の図示

sum(p1 < 0.05)#母分散を使ったZ検定の有意なP値の個数

## [1] 4912

sum(p2 < 0.05)#不偏分散を使ったZ検定の有意なP値の個数

## [1] 10860

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

結果、正しくZ検定した場合は危険率が5%程度であったのに対して、検定統計量z'を使った場合は危険率が10.8%程度になってしまった。これでは実用に耐えるとは到底言えない。

　じゃあ、結局、検定統計量z'はどんな分布に従うというのだろうか。そう、それを発見したのがGossetであり、他でもないt分布なのだ。実際にはt検定の項目を見るとわかるように、Gossetが仮定した検定統計量tはちょっと複雑な形であるが、実質、検定統計量z'と同じものである。Gossetが仮定した検定統計量tは、自由度の計算が簡単である。将来的にはWelchのt検定という、Studentのt検定を拡張したより利便性の高い検定を紹介する。Welchのt検定で使う検定統計量tは、検定統計量z'そのものである。その代わり、その自由度の計算はややこしい。

Gossetが小標本問題を解決した

　「検定統計量z'を使ったZ検定は実用に耐えない」、と表現したがある条件を加えれば問題がなくなる。それは、サンプルサイズが十分大きい場合だ。サンプルサイズが十分大きい場合、自由度も大きくなり、t（z'）分布は標準正規分布に近づく。ゆえに、正しくZ検定を行うことができる。けれども、実際のデータで、Z検定が実用に耐えられるほどのサンプルサイズをそろえられるだろうか、それはどれくらいのサンプルサイズだろうか。

　それは以下のように確認できる。標準正規分布において両側検定P値が0.05となる検定統計量を求め、それがt分布においてはP値がどうなるかを計算する。自由度によってt分布は形を変えるが、自由度が大きくなるほど標準正規分布に近づく。自由度を変化させて、t分布におけるP値が0.05に「十分近づく」自由度を明らかにしよう。ここでは、t分布におけるp値が0.051を切る点を十分近づいたとする。

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

z <- qnorm(p = 0.975, mean = 0, sd = 1)

df <- 1:302

p <- NULL

for(i in df){

  p <- c(p, pt(q = z, df = i, lower.tail = F) * 2)

}

ylim <- range(c(0,p))

plotn(p ~ df, ylim = ylim, type = "l")#図10の図示

overdraw(abline(h = 0.05, col = "red"))

plotn(p ~ df, ylim = c(0.048, 0.052), type = "l")#図11の図示

overdraw(abline(h = 0.05, col = "red"),

         abline(h = 0.051, col = "blue"))

sum((p - 0.05) > 0.001)

## [1] 277

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

すると、衝撃的なことに検定統計量z'を使って危険率を5.1%に抑えたZ検定を行おうとすると、自由度は277必要になる（！！）。サンプルサイズに直すと、2標本で合計、279のサンプルサイズが必要で、2標本均等に割れば140は必要だ。周りのデータ解析を見て、どうだろう。それだけのサンプルサイズをそろえた研究なんて、かなりレアだろう。生態学関係なんかでは、まずこんな数字を見かけない。

　Gossetが師事したPearsonなら、「それなら、それだけデータをそろえればいい」と言い放ちそうなものだ（実際そういうことを言ったかは知らない、中村治 著 『基礎から学ぶ統計学』内で語られる逸話をもとに想像）。対して、Gossetは小標本でも頑健な解析をしたいと考えて、Studentのt検定を開発したと伝えられる。以前紹介したように、t検定は小標本でも十分、威力を発揮できる解析である。

母分散が既知ならZ検定、は間違ってないけど……

　統計初学者だった頃、「母分散が既知ならZ検定、未知ならt検定」（あるいは既知ならこう計算できる、未知ならこう計算できる、みたいなのもあった）と習い、私の個人的感想としては「なんで母集団が一緒なのに、母分散が既知かで検定（帰無分布）が変わるんだよ！！！」と思い、理解が至らなかった覚えがある。なんで私の持つ事前情報が、バックに潜む確率分布を変えてしまうのかが理解できなかったのだ。間違ってはないけど、正直言って「母分散が既知ならZ検定/未知ならt検定」という表現は混乱を招くと思う。それよりかは「検定統計量の計算に母分散を使ったらZ検定/不偏分散を使ったらt検定」というように、検定統計量の計算に即した表現のほうが良いのではないかと思う。この表現は暗に、母分散を使う＝母分散の数値を知っている、不偏分散を使う＝母分散の推定している、という状況を含意している。この記事で多少なりとも、理解の助けになったら幸いである。

補足：再生性の証明

　この補足では、上記で仮定した正規分布の再生性を証明する。証明は2段階に分かれる。まず、和の再生性を証明し、次に定数倍の再生性を示す。これらの2つを示せば、任意の1次結合で再生性が示される。

・和の再生性

　確率変数Xは平均μ1、分散σ1^2の正規分布に従い、確率変数Yは平均μ2、分散σ2^2の正規分布に従うとする。この時、X+Yの従う分布を考えたい。ここで一度、X + Y = k（定数）とする。下図をもとに考えると、あるkにおける確率密度の「一つ」は、正規分布の確率密度関数をf(x)とすれば、f(x|μ1, σ1)f(y|μ2, σ2)という積のかたちで与えられる。

しかし、あくまで「一つ」であるので、あるkにおける確率密度全体は、xを全域にわたって積分したものに相当する。また、Y = k - Xであることに注意すれば、積分は下記のように計算できる。

さて、このkに改めて注目すると、得られた結果は平均μ1+μ2、分散σ1^2+σ2^2の正規分布の確率密度関数になっている。k = X + Yだったので、X + Yは平均μ1+μ2、分散σ1^2+σ2^2の正規分布に従うと言える。なお、上記の∫f(x)f(k-x)dxの積分法は畳み込みconvolutionと呼ばれる。これについての視覚的理解は、Youtubeの3Blue1BrownJapanチャンネルに素晴らしい動画があるので、見てみるとよいだろう（原典は3Blue1Brownチャンネル、英語がわかるならこちらもおすすめ）。

・定数倍の再生性

　続いて定数倍の再生性を証明する。まず真っ先に、イメージだけ考えておこう。確率変数Xが平均μ、分散σ^2の正規分布に従っているとき、確率を1にするような正規化定数を無視すると、aXが従う分布は元の正規分布が横方向にa倍に引き伸ばされたような形になっている。

図形的に考えれば、同じ次元単位をもつ標準偏差も同様にa倍になるはずだ。それならば、標準偏差の2乗である、分散は(aσ)^2となる。

　図形的イメージを使ってもう少し説明しよう（もう少し厳密な証明は後述）。以下のような図形変形をたどることを考える。左下が元の確率変数Xが従う確率密度分布である。

左下の確率密度関数の面積は1である。それを縦方向に√(2πσ^2)倍する。さらに、横軸方向にa倍すると、確率変数Xもa倍された位置に移動する。同時に青で示された面積はa倍されるので、aX = Yとして示すと関数の積分と面積の関係は右上のようにあらわせる。求めたい値は、関数の広がりを表すAである。以下のように求めればA = (aσ)^2である。

ゆえに最後、右下のように、右辺を1になるように割れば、左辺の積分の中身は確率密度関数が満たすべき性質を満たしている。この時Y = aXは、平均aμ、 分散(aσ)^2の正規分布に従うと言える。

　最後、数式的に証明しておこう。以下のようにXが従う確率密度分布に関して、sからtまで積分し、x = y/aと置換して積分することを考える。

すると、置換によって確率密度関数のかたちは、横方向にa倍したことに相当することがわかる。得られた値は、平均aμ、 分散(aσ)^2の正規分布の確率密度関数そのものである。