目次

・Studentのt検定が使える条件

・独立な標本とは？

・対応のあるt検定

・独立性と検出力

・独立性と危険率

・Rで行う対応のあるt検定

Studentのt検定が使える条件

　2標本のStudentのt検定は、母集団に次のような過程が成り立つときに有効に使うことができる。

1. 等分散性……2標本の分散が等しい

2. 独立性……2標本は独立にサンプリングされている

3. 正規性……2標本は正規分布に従う

では上記の3条件が崩れたときには、どれくらい信頼をおけない解析となってしまうのだろうか？　本稿では、シミュレーションを通して、t検定の頑健性と、仮定が満たされないときの対策を紹介する。

独立な標本とは？

　前回は等分散性について検討したが、今回は独立性について検討しよう。これは母集団の分布もそうだが、サンプリングの性質にも依存した仮定となっている。そもそも独立な標本とはどういうことなのだろうか。一言で言えば、2標本は互いに数値の影響がない状態である。例えば、標本1で高い値が出るときは、標本2でも高い値が出やすいような関係があるときは独立でない可能性が高い。もし、独立であれば標本1で高い値が出ようが、標本2には影響がなく高い値も低い値も出る。独立でない状態を相関がある、と表現することもある。実践的なサンプリングの例ならば、ある植物を60個体用意し、30個体はそのまま、30個体は肥料を与えて育てたとき、生育量に差があるかを検定するとしよう。肥料を与えた30個体と与えなかった30個体は、全く別の個体なので、互いに影響を与えるような関係にはない。ゆえに独立と考えてもよいだろう。

　一方で、ある植物30個体を用意し、ある時に肥料を与えて、肥料を与える前と後で生育量を比較する状況を考える。この時は、肥料に生育を促進する効果があるのだとしたら、肥料を与える前から大きく育っていた個体はより大きく育つし、小さかった個体も大きく育つだろう。しかし、もともと小さかった個体が大きかった個体を追い抜いて生育するとは考えにくく、肥料を与える前と後で大きさの順序関係はほとんど変わらない、つまり肥料を与える前の大きさが、肥料を与えた後の大きさに影響を与えてしまっている。このような時は、独立な標本とは言えない。

　シミュレーションにより独立ではない2標本を独立を仮定したStudentのt検定にかけてしまった場合、どんな影響があるかを確認してみよう。結論から言うと、独立ではない標本の場合は、対応関係を考えた解析を行う方が間違いを犯す確率がぐっと低下する。以下のようなコードで、独立な標本と独立ではない標本を作る。標本m11は母平均50、標準偏差20の正規分布から生成する。独立な標本m21は母平均60、標準偏差20の正規分布から生成する。一方、独立ではない標本m22は標本m11に平均10、標準偏差2の正規分布から生じた値を足して作っている。

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

m11 <- rnorm(10000, mean = 50, sd = 20)

m21 <- rnorm(10000, mean = 60, sd = 20)

m22 <- m11 + rnorm(length(m11), mean = 10, sd = 2)

mean(m11)

## [1] 50.12232

sd(m11)

## [1] 19.75059

mean(m21)

## [1] 59.84506

sd(m21)

## [1] 20.11972

mean(m22)

## [1] 60.13049

sd(m22)

## [1] 19.82207

m <- data.frame(m = c(m21, m22), g = c(rep("i", 10000), rep("d",10000)))

histn(m$m ~ m$g, xlab = "value", ylab = "頻度") #図1の描画

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

平均や標準偏差、分布を比較しても標本m21とm22はほとんど違いがなく、一見区別ができない。だが、m11との関係を見れば、全く性質の異なる標本であることがわかる。

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

plotn(m21 ~ m11, col.dot = "#00000020") #図2の描画

plotn(m22 ~ m11, col.dot = "#00000020") #図3の描画

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

対応のあるt検定

　では、独立ではない標本に対して検定を行う場合を考えてみよう。平均値の差は10に設定している。このとき、2標本のt検定を行うことは良い戦略ではない。独立ではなく対になるデータがある場合、いわゆる対応関係のあるStudentのt検定を行うことになる。これは、対となるデータで差をとり、差の標本を作る。この差の標本に関して、1標本のt検定（平均値のt検定）に帰着させる。このときの差の標本の統計量tは以下のように計算する。

XDは差の標本の平均、μは帰無仮説に相当する平均（差がないなら0）、UDは差の標本の不偏分散、nDは差の標本のサンプルサイズである。2標本のt検定では、各標本で平均値をとり、その差を使っていた点に違いがある。

独立性と検出力

　独立ではない標本を10000回生成し、2標本のt検定と対応関係のあるt検定を行うコードは以下のようになる。

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

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(30, mean = 50, sd = 20) #標本1作成

  m2 <- m1 + rnorm(length(m1), mean = 10, sd = 2) #標本2作成、標本1と独立ではない点に注意

  ue <- ((30 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(30 + 30 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/30 + 1/30))) #2標本のt検定の統計量tの算出

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value) #2標本のt検定を行い、p値だけ取り出す

  m <- m1 - m2

  t3 <- c(t3, (mean(m) - 0)/(sd(m)/sqrt(30))) #対応関係のあるt検定の統計量tの算出

  p3 <- c(p3, t.test(m1 - m2)$p.value) #2標本の差について、1標本のt検定を行い、p値だけ取り出す

}

histn(t2, xlab = "統計量t", ylab = "頻度") #図4の描画

overdraw(abline(v = qt(df = 58, p = 0.975), col = "red"),

         abline(v = qt(df = 58, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図5の描画

histn(t3, xlab = "統計量t", ylab = "頻度") #図6の描画

overdraw(abline(v = qt(df = 29, p = 0.975), col = "red"),

         abline(v = qt(df = 29, p = 0.025), col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図7の描画

sum(p2 < 0.05)

## [1] 4249

sum(p3 < 0.05)

## [1] 10000

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

2標本のt検定では、検出力は43％程度となる。対応関係に気付かなかった場合、独立な標本と独立ではない標本はほとんど同じような母集団の分布になるので、検出力は独立な標本のときと大差がない。ところが、対応関係のあるt検定を行えば、検出力は100%となる。有意水準は両検定で5%なので、これは大幅な解析の改善である。逆に言えば、適切な検定を行えばほぼ確実に差を検出できるのに、不適切な検定を行ったために第2種の過誤を犯す確率が不必要に高まっている、と考えることができる。

　なぜこのような劇的な検出力の改善が起こったかというと、対となるデータで差をとったことにより、より標準偏差の小さい母集団の分布に帰着できたためである。独立ではない標本は、独立な標本と同様に平均60、標準偏差20であることに変わりはない。しかし、独立ではない標本2を生成した過程では、平均50、標準偏差20の標本1に平均10、標準偏差2の正規分布のデータを足していた。逆に標本1と標本2の差をとると、その差は平均10、標準偏差2の正規分布になってしまうということだ。統計量tの性質から、標準偏差が小さくなれば、検出力の大幅な改善が期待される。このように独立性を念頭に置いた解析は非常に重要なのだ。

　もともと、独立な2標本を対応のあるt検定にかけることはまれだろう。2標本のサンプルサイズがそろっていないことも多いので、気づきやすい。問題は、独立ではないのに、普通の2標本のt検定にかけてしまう場合である。自分の実験が、同じ個体から複数回、データを得ている、というような状態になっていないか常に確認しておこう。

　ちなみに、前回までのt検定では顕著ではなかったが、対立仮説が正しいときのt分布がゆがんだ形になっている。これは間違いではなく、このようなt分布を非心t分布noncentric t-distributionと呼ぶ（前回までの条件でも母集団の分布に対して平均値の差が大きければ歪みが顕著になる）。なぜ、非心t分布となるかは、カイ二乗分布の項を参照。

独立性と危険率

　次に危険率について検討しよう。今回は2標本に差がない状況を作るので、独立ではない標本を作るときは、平均0の正規分布を足せばよい。

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

t2 <- NULL

t3 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  m1 <- rnorm(30, mean = 50, sd = 20) #標本1作成

  m2 <- m1 + rnorm(length(m1), mean = 0, sd = 2) #標本2作成、標本1と独立ではないが、差はない

  ue <- ((30 - 1)*(sd(m1)^2) + (30 - 1)*(sd(m2)^2))/(30 + 30 - 2)

  t2 <- c(t2, (mean(m1) - mean(m2))/sqrt(ue*(1/30 + 1/30))) #2標本のt検定の統計量tの算出

  p2 <- c(p2, t.test(m1,m2,var.equal = T)$p.value) #2標本のt検定を行い、p値だけ取り出す

  m <- m1 - m2

  t3 <- c(t3, (mean(m) - 0)/(sd(m)/sqrt(30))) #対応関係のあるt検定の統計量tの算出

  p3 <- c(p3, t.test(m1 - m2)$p.value) #2標本の差について、1標本のt検定を行い、p値だけ取り出す

}

histn(t2, xlab = "統計量t", ylab = "頻度") #図8の描画

overdraw(abline(v = qt(df = 58, p = 0.975), col = "red"),

         abline(v = qt(df = 58, p = 0.025), col = "red"))

histn(p2, xlab = "P value", ylab = "頻度") #図9の描画

histn(t3, xlab = "統計量t", ylab = "頻度") #図10の描画

overdraw(abline(v = qt(df = 29, p = 0.975), col = "red"),

         abline(v = qt(df = 29, p = 0.025), col = "red"))

histn(p3, xlab = "P value", ylab = "頻度") #図11の描画

sum(p2 < 0.05)

## [1] 0

sum(p3 < 0.05)

## [1] 471

差のある仮説が正しいときとは、また違った様相のtとpの分布になることが見て取れる。今回は差のない仮説が正しいので、pが有意水準を下回った場合、まちがった判断を下すこととなる。ゆえに、ここでのp値は検出力ではなく、危険率の指標となる。普通のStudentのt検定の危険率は0%であり、一方、対応関係のあるt検定の危険率は約5%だ。一見すれば、間違いを犯す確率が低いので、普通のt検定のほうが優れているように見えなくもないが、対応のあるt検定の危険率5%は、正しく検定を行ったために生じていると考えることもできる。むしろ普通のt検定の場合は、危険率が低すぎるゆえに検出力が低くなっているわけなので、対応のある標本において普通の2標本のt検定を行うことは、決して得策ではないだろう。

　特に、危険率を確認する場合、p値の分布に関する累積分布関数(CDF)を確認するのが良い。標本に差がないとき、期待されるp値のCDFはy=xの直線となる。それに対するt検定とU検定のCDFも確認しよう。 

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

cdf <- function(p, v){

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

}

v <- seq(0, 1, length = 200)

p2c <- cdf(p2,v)

p3c <- cdf(p3,v)

cs <- data.frame(v,v,p2c,p3c)

plotn(cs[,1],cs[,2:4], mode = "m", xlab = "P-value", ylab = "CDF", type = "l",

      col.line = c("black","red","blue"), 

      pt.col.leg = c("black","red","blue"), 

      legend = T, line = T, pt.cex.leg = 0, 

      leg.lab = c("expected", "unpaired","paired"), leg.sp = 6) #追加図1の描画

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

すると、対応のあるt検定ではy = xの直線状に乗るが、対応がないt検定では明らかに下側になっている。これは、極めて対応のないt検定での検出力が劣ることを示している。

Rで行う対応のあるt検定

　上記のような背景のもと、対応のあるt検定を行ってみよう。以下のようなサンプルサイズ30のデータがあるとする。また、これらのデータに対応関係があることがわかっているとする。有意水準は5%としよう。まずは、データを描画してみて様子を確認する。

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

m1 <- c(51, 41, 51, 80, 52, 50, 58, 8, 29, 62,

        45, 53, 62, 22, 45, 27, 68, 43, 40, 37, 

        78, 44, 63, 17, 56, 42, 33, 33, 72, 30)

m2 <- c(54, 40, 51, 83, 51, 52, 57, 10, 31, 65, 

        46, 54, 60, 22, 48, 28, 73, 44, 41, 37,

        79, 42, 63, 18, 57, 44, 33, 38, 72, 32)

d <- data.frame(d = c(m1, m2), g = c(rep("m1", length(m1)), rep("m2", length(m2))))

boxplotn(d$d ~ d$g, ylab = "data", jitter.method = "swarm") #図12の描画

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

一見するとこれらの2標本の平均値の差は極めて小さそうだ。これらのデータには対応関係があるので、箱ひげ図だけではなく2次元のプロットも確認してみよう。

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

plotn(m2 ~ m1) #図13の描画

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

図13を確認すると、明らかに2標本は独立の関係にはない。今回は標本が具体的に何かを考えずに架空データを用いているが、実際の実験ではこのようなデータは同じ個体から複数回サンプリングを行うと生じやすい。このような独立ではない標本では2標本で差をとり、1標本のt検定に帰着する。

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

m <- m1 - m2

m

##  [1] -3  1  0 -3  1 -2  1 -2 -2 -3 -1 -1  2  0 -3 -1 -5 -1 -1  0 -1  2  0 -1 -1

## [26] -2  0 -5  0 -2

n <- length(m)

fit <- t.test(m, mu = 0)

fit

## 

##  One Sample t-test

## 

## data:  m

## t = -3.445, df = 29, p-value = 0.001761

## alternative hypothesis: true mean is not equal to 0

## 95 percent confidence interval:

##  -1.7530467 -0.4469533

## sample estimates:

## mean of x 

##      -1.1

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

fitの中身を確認してみると、いくつかのデータが並んでいる。tは検定統計量tのことで以下のコマンドで計算できるものと同じだ。

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

t0 <- (mean(m) - 0)/(sd(m)/sqrt(n))

t0

## [1] -3.44501

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

次にdfは自由度（degree of freedom）のことで、今回は以下のように計算する。次にp-valueで、あらかじめ定めた有意水準を下回っていれば、帰無仮説を棄却する。p-valueはt分布を使って次のように計算できる。

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

df <- n - 1 #自由度

df

## [1] 29

pt(df = df, q = t0, lower.tail = T) * 2

## [1] 0.001760914

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

今回は0.0018なので、0.05を下回っており、帰無仮説を棄却する。次の行の”alternative hypothesis: true mean is not equal to 0”は「対立仮説(alternative hypothesis)は平均が0ではないこと」を示している。95 percent confidence intervalは平均に関する95%信頼区間である。これらの2標本をもとにしたt分布は、平均0のt分布からt0 = -3.45ずれた形となっている。このときの、95％のtが含まれる範囲は、下記のように計算できる。

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

t0 - qt(df = df, p = 0.975)

## [1] -5.49024

t0 + qt(df = df, p = 0.975)

## [1] -1.39978

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

これはt分布における範囲なのでもとの平均の分布に戻すときは下記のように計算できる。

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

(t0 - qt(df = df, p = 0.975)) * (sd(m)/sqrt(n))

## [1] -1.753047

t0 * (sd(m)/sqrt(n))

## [1] -1.1

(t0 + qt(df = df, p = 0.975)) * (sd(m)/sqrt(n))

## [1] -0.4469533

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

確かにt検定を行った時の95%信頼区間と一致している。また、t0を平均の分布に戻したときの値はちょうど、差の標本の平均値となっている。これはmean of xの値でもある。

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

mean(m)

## [1] -1.1

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