『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第6章 2群の母平均の差のt検定

Author

SUGINO Isamu, Build with R4.3.2

Published

October 20, 2024

0 全体の章構成

1 2群の母平均差の検定の基本

(二群の母平均の差の検定)

2群の母分散が既知である場合は,

\[z=\frac{(\bar{y}_{1}-\bar{y}_{2})-(\color{blue}{\mu_{1}}-\color{blue}{\mu_{2}})}{\sqrt{\frac{{\sigma_{1}}^2}{n_{1}}+\frac{{\sigma_{2}}^2}{n_{2}}}}\;\sim\;N(0, 1^2)\]

更に等分散であれば,

\[z=\frac{(\bar{y}_{1}-\bar{y}_{2})-(\color{blue}{\mu_{1}}-\color{blue}{\mu_{2}})}{\sigma\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}=\frac{(\bar{y}_{1}-\bar{y}_{2})-(\color{blue}{\mu_{1}}-\color{blue}{\mu_{2}})}{\sigma}\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}\;\sim\;N(0, 1^2)\]

1-1 2群の母分散が等しいかどうかのF検定

 演習用データの男女別幸福度の平均と不偏分散,有効ケース数の表が本文にある。
これを求める方法を示しておく。まずはデータを準備する部分を再掲しておく。
本文では簡便化の為に幸福度の変数名を元のq1700のままで書いているが,実際には11という値を欠損値NAにしなければならないので,その欠損値処理をした変数を新たにq1700bとしている。

データ・フレイムと変数の準備

データファイルはworking directoryに置く

data01 <- read.csv("practice.csv")
with(data01, table(sex,   useNA = "always"))
sex
   1    2 <NA> 
 165  219    0 
with(data01, table(q1700, useNA = "always"))
q1700
   0    1    2    3    4    5    6    7    8    9   10   11 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    1    5 

q1700 = 11は「わからない」で欠損値NAとすべきなので,新変数q1700bを作成。

data01$q1700b <- c(0:10)[data01$q1700 + 1]    # 変数のリコード法にやや工夫
with(data01, table(q1700b, useNA = "always")) # 変換の結果を確認
q1700b
   0    1    2    3    4    5    6    7    8    9   10 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    6 
基本的な要約統計
with(data01, by(q1700b, sex, mean, na.rm = T))
with(data01, by(q1700b, sex, var,  na.rm = T))
with(data01, by(q1700b, sex, function(x) sum(!is.na(x))))
sex: 1
[1] 6.36646
------------------------------------------------------------ 
sex: 2
[1] 6.917051
sex: 1
[1] 4.083618
------------------------------------------------------------ 
sex: 2
[1] 3.08568
sex: 1
[1] 161
------------------------------------------------------------ 
sex: 2
[1] 217
with(data01, tapply(q1700b, sex, mean, na.rm = T))
with(data01, tapply(q1700b, sex, var,  na.rm = T))
with(data01, tapply(q1700b, sex, function(x) sum(!is.na(x))))
       1        2 
6.366460 6.917051 
       1        2 
4.083618 3.085680 
  1   2 
161 217 

等分散性の検定は次のように行う。本文のF分布グラフを描くスクリプトをその後に挙げる。

等分散性の検定
with(data01, var.test(q1700b ~ sex))

    F test to compare two variances

data:  q1700b by sex
F = 1.3234, num df = 160, denom df = 216, p-value = 0.05555
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.9933775 1.7743088
sample estimates:
ratio of variances 
          1.323409 
F分布のグラフ

dofは自由度ヴェクトル(degree of freedom)

dof <- with(data01, by(q1700b, sex, function(x) sum(!is.na(x)))) - 1
p0  <- .05 # 有意水準
plot(x <- seq(0, 2, by = .01), df(x, dof[1], dof[2]), type = "l", bty = "n",
    xlab = "", ylab = "", family = "serif",
    main = sprintf("自由度%d, %dのF分布", dof[1], dof[2]), 
    sub  = sprintf("両側%d%%棄却域", p0 * 100))
segments(x1 <- c(seq(0,     qf(p0/2, dof[1], dof[2]), by = .01),
                 seq(qf(1 - p0/2, dof[1], dof[2]), 2, by = .01)), 0,
         x1, df(x1, dof[1], dof[2]), col = "red")
axis(side = 1, at = round(c(qf(p0/2, dof[1], dof[2]),
                            qf(1 - p0/2, dof[1], dof[2])), 2),
    las = 2, family = "serif", font = 2, col.axis = "red")

 共通の母分散の推定値を計算するには以下の様にすると良い。不偏分散と自由度をいずれも長さ2のヴェクトルとして格納すると,計算式はすっきりとしたものになる。

var0 <- with(data01, by(q1700b, sex, var, na.rm = T))
sum(var0 * dof)/sum(dof)
[1] 3.510335

1-2 母分散が等しい2群の母平均差のt検定

 テキストの\(pp. 086, 089-091\)に,下記の様に\(\hat{y}_{1},\;\hat{y}_{2}\)とあるが,これらは全て\(\bar{y}_{1},\;\bar{y}_{2}\)の間違いです。お詫びして訂正します。

誤植の訂正

\(Student\)(William S. Gosset)のt検定統計量 \[t=\frac{(\bar{y}_{1}-\bar{y}_{2})-(\color{blue}{\mu_{1}}-\color{blue}{\mu_{2}})}{\hat{\sigma}\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}\;\sim\;t_{(n_{1}+n_{2}-2)}\] 但し,\(s_{1}^2\)は群1の標本分散,\(s_{2}^2\)は群2の標本分散とする(不偏分散はそれぞれ\(\hat{\sigma}_{1}^2\), \(\hat{\sigma}_{2}^2\))。

\[\hat{\sigma}^2=\frac{n_{1}s_{1}^2+n_{2}s_{2}^2}{n_{1}-1+n_{2}-1}=\frac{(n_{1}-1)\hat{\sigma}_{1}^2+(n_{2}-1)\hat{\sigma}_{2}^2}{n_{1}-1+n_{2}-1}=\frac{SS_{1}+SS_{2}}{n_{1}+n_{2}-2}\] Null Hypothesis, \(H_{0}: \mu_{1}=\mu_{2}\)が正しいとすると, \[t=\frac{\bar{y}_{1}-\bar{y}_{2}}{\hat{\sigma}\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}=\frac{\bar{y}_{1}-\bar{y}_{2}}{\hat{\sigma}}\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}\;\sim\;t_{(n_{1}+n_{2}-2)}\]

Studentのt検定

 

test1 <- with(data01, t.test(q1700b ~ sex, var.equal = T))
test1

    Two Sample t-test

data:  q1700b by sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 -0.9337912 -0.1673909
sample estimates:
mean in group 1 mean in group 2 
       6.366460        6.917051 
箱ひげ図

本文のような箱ひげ図を描くスクリプト。少し変えた部分もある。

with(data01, 
  boxplot(q1700b ~ sex, las = 1, cex.main = 1, varwidth = T,
          main = "男女別幸福度", 
          sub  = "実線は中央値,赤い線は算術平均", 
          names = c("男性","女性"), 
          par(family = "serif"), 
          horizontal = T)
)

m0 <- with(data01, by(q1700b, sex, mean, na.rm = T))
segments(m0[1], 0.65, m0[1], 1.35, col = "red")
segments(m0[2], 1.60, m0[2], 2.40, col = "red")

この例題をベイズ統計学的手法で計算した例を,このサポートウェブの〔付録〕2-1に示しています。 Bayes

1-3 母分散が異なる2群の母平均差のウェルチ検定(Welch test)

\(Welch\)の検定統計量は, \[t=\frac{\bar{y}_{1}-\bar{y}_{2}}{\sqrt{\frac{\hat{\sigma}_{1}^2}{n_{1}}+\frac{\hat{\sigma}_{2}^2}{n_{2}}}}\;\sim\;t_{df}\] 自由度の式は複雑なので省略する。自由度が小数になる事が通常である。

 本文では省略している部分も一つ一つ確認しながら説明してゆこう。
まずは世帯年収fincomeの度数分布表,男女別の平均や分散,有効ケース数の出力,そして等分散性のF検定を行おう。

男女別の度数分布、要約統計、等分散性検定
with(data01, table(sex, fincome, useNA = "always"))
      fincome
sex     0 25 75 125 200 300 400 500 600 700 800 925 1125 1375 1750 2500 <NA>
  1     2  0  2   4   6   4   9  20  19   9  14  18   12    5    2    1   38
  2     3  1  5   4  11  15  19  16  13  18   9  15   19    8    9    5   49
  <NA>  0  0  0   0   0   0   0   0   0   0   0   0    0    0    0    0    0
with(data01, by(fincome, sex, mean, na.rm = T))
sex: 1
[1] 697.4409
------------------------------------------------------------ 
sex: 2
[1] 737.0588
with(data01, by(fincome, sex, var,  na.rm = T))
sex: 1
[1] 142577.7
------------------------------------------------------------ 
sex: 2
[1] 276015
dof <- with(data01, by(fincome, sex, function(x) sum(!is.na(x))))
dof
sex: 1
[1] 127
------------------------------------------------------------ 
sex: 2
[1] 170
with(data01, 
     var.test(fincome ~ sex)
)

    F test to compare two variances

data:  fincome by sex
F = 0.51656, num df = 126, denom df = 169, p-value = 0.000114
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.3737504 0.7196063
sample estimates:
ratio of variances 
         0.5165579 
p0 <- .05 # 有意水準
plot(x <- seq(0, 2, by = .01), df(x, dof[1] - 1, dof[2] - 1), 
     type = "l", bty = "n",
     xlab = "", ylab = "", family = "serif",
     main = sprintf("自由度%d, %dのF分布", dof[1] - 1, dof[2] - 1), 
     sub  = sprintf("両側%d%%棄却域", p0*100))

segments(x1 <- c(seq(0, qf(p0/2, dof[1] - 1, dof[2] - 1), by = .01),
                 seq(qf(1 - p0/2, dof[1] - 1, dof[2] - 1), 2, by = .01)), 0,
                 x1, df(x1, dof[1] - 1, dof[2] - 1), col = "red")
axis(side = 1, las = 2, family = "serif", font = 2, col.axis = "red",
     at = round(c(qf(p0/2,     dof[1] - 1, dof[2] - 1), 
                  qf(1 - p0/2, dof[1] - 1, dof[2] - 1)), 2))

 男女の世帯年収の分布を,上の1-2でも紹介した箱ひげ図を用いて視覚的に確認しておこう。

with(data01, 
     boxplot(fincome ~ sex, horizontal = T, las = 1,
             names = c("男性", "女性"), main = "男女別世帯年収")
)

ウェルチの検定

 等分散性が棄却されたので,ウェルチの検定を行う。
等分散ではないというオプション “var.equal = FALSE” を指定してt.test( )関数を用いるのだが,このオプションを単に省略すれば,Rのt.test( )関数はデフォルトでウェルチの検定を行うようになっている。

with(data01, 
     t.test(fincome ~ sex)
)

    Welch Two Sample t-test

data:  fincome by sex
t = -0.756, df = 294.6, p-value = 0.4503
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 -142.7534   63.5176
sample estimates:
mean in group 1 mean in group 2 
       697.4409        737.0588 

ウェルチ検定のt統計量の絶対値は.756と小さく,有意確率は45%もある。
つまり,ゼロ仮説が正しく,母平均は男女で等しいとしても,標本でたまたまこの程度の差が生じることはよくあることであり,(ゼロ仮説を含めて)何も疑うべきものはないということである。
母平均の差の95%信頼区間はマイナスからプラスに広がっており,0を間に含んでいる。
ゼロ仮説の検定が5%で有意にならない(否定できない)ことと一致している。

〔付録〕3-2では,RとStanでWelch検定と等価な分析を行う方法を紹介しています。 R & Stan

2 区間推定と効果サイズ

2-1 母平均の差の区間推定

 本文同様,男性平均-女性平均の母平均差の95%信頼区間を式に従って求め,それと検定結果の出力が一致する事を確認しよう。

m0   <- with(data01, by(q1700b, sex, mean, na.rm = T))
var0 <- with(data01, by(q1700b, sex, var,  na.rm = T))
n0   <- with(data01, by(q1700b, sex, function(x) sum(!is.na(x))))

# 男女の幸福度の標本平均差
m_diff <- m0[1] - m0[2]

# 共通の母標準偏差の推定値
sigma <- sqrt(sum(var0*(n0 - 1))/(sum(n0) - 2))
母平均差の95%信頼区間下限値と上限値
m_diff - qt(c(.975, .025), df = sum(n0) - 2) * sigma * sqrt(sum(1/n0))
[1] -0.9337912 -0.1673909
Studentのt-test
test1 <- with(data01, t.test(q1700b ~ sex, var.equal = T))
test1

    Two Sample t-test

data:  q1700b by sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 -0.9337912 -0.1673909
sample estimates:
mean in group 1 mean in group 2 
       6.366460        6.917051 

 1-3のウェルチの検定はどうなるだろうか。
やや難しく見えるが,以下の計算で求められる。
最後にウェルチの検定の結果を示して,信頼区間が一致している事を確認する。

m0  <- with(data01, by(fincome, sex, mean, na.rm = T))
var0 <- with(data01, by(fincome, sex, var,  na.rm = T))
n0  <- with(data01, by(fincome, sex, function(x) sum(!is.na(x))))

m_diff <- m0[1] - m0[2]   # 男女の世帯収入の標本平均差
se0    <- sqrt(sum(var0/n0)) # 標準誤差推定値
Welchの自由度と区間推定

ウェルチ検定の自由度

df1 <- prod(n0 - 1) * sum(var0/n0) ^ 2 /
       ( (n0[2] - 1) * (var0[1] / n0[1]) ^ 2 + 
         (n0[1] - 1) * (var0[2] / n0[2]) ^ 2   )

母平均差の95%信頼区間下限値と上限値

m_diff - qt(c(.975, .025), df = df1) * se0
[1] -142.7534   63.5176
Welchのt-test
test2 <- with(data01, t.test(fincome ~ sex))
test2

    Welch Two Sample t-test

data:  fincome by sex
t = -0.756, df = 294.6, p-value = 0.4503
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 -142.7534   63.5176
sample estimates:
mean in group 1 mean in group 2 
       697.4409        737.0588 

2-2 効果サイズ(effect size)と検定力

\[t=\frac{\bar{y}_{1}-\bar{y}_{2}}{\hat{\sigma}\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}=\frac{\bar{y}_{1}-\bar{y}_{2}}{\hat{\sigma}}\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}=\hat{d}\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}\;\sim\;t_{(n_{1}+n_{2}-2)}\]

 1-1~1-2の幸福度の男女差の例で効果サイズを考える。
テキスト本文や上記スクリプトとは少し違うスクリプトで計算してみる。
様々な方法で同じ様な事が出来ると云う事を読み取って欲しい。

m0 <- with(data01, by(q1700b, sex, mean, na.rm = T)); m0
sex: 1
[1] 6.36646
------------------------------------------------------------ 
sex: 2
[1] 6.917051
v0 <- with(data01, by(q1700b, sex, var,  na.rm = T)); v0
sex: 1
[1] 4.083618
------------------------------------------------------------ 
sex: 2
[1] 3.08568
n0 <- with(data01, by(q1700b, sex, function(x) sum(!is.na(x)))); n0
sex: 1
[1] 161
------------------------------------------------------------ 
sex: 2
[1] 217
m_diff <- diff(rev(m0)); m_diff # 男女の幸福度の標本平均差
         1 
-0.5505911 
sigma  <- sqrt(sum(v0*(n0 - 1))/(sum(n0) - 2)) # 共通の母標準偏差の推定値
es     <- m_diff/sigma; es # 標本効果サイズ
         1 
-0.2938698 

 各平均,分散,ケース数,標本平均の差といった基本情報を確認し,母標準偏差が共通であるとの仮定の下で標本効果サイズを求めた。

複数条件での検定力

 次に,各群でケース数が異なるので,小さい方のサンプルサイズ,大きい方のサンプルサイズ,算術平均,幾何平均の4通りのサンプルサイズをヴェクトルとして使用した。
 4通りのサンプルサイズそれぞれに対応した検定力が出力されている。
本文の結果ともほぼ一致している。

# 4通りの標本サイズ
n1 <- c(n0, floor(sum(n0)/2), floor(sqrt(prod(n0)))); n1
  1   2         
161 217 189 186 
power.t.test(n = n1, 
             sd = sigma, 
             delta = m_diff, 
             sig.level = .05, 
             power = NULL, 
             type = "two.sample")

     Two-sample t test power calculation 

              n = 161, 217, 189, 186
          delta = 0.5505911
             sd = 1.873589
      sig.level = 0.05
          power = 0.7481781, 0.8630816, 0.8131254, 0.8069331
    alternative = two.sided

NOTE: n is number in *each* group

パッケイジ”pwr”のpwr.t2n.test( )関数

 二群のサイズが異なる場合の検定力分析を正しく行う為の追加パッケイジ”pwr”とその関数pwr.t2n.test( )を使用してみよう。
これも本文とはスクリプトが少し異なるので注意してほしい。
PCにpwrをインストールした事が無い場合は,install.packages(“pwr”, repos=“https://cran.ism.ac.jp”)を一度だけ実行しておく。

# install.packages("pwr", repos="https://cran.ism.ac.jp")
library(pwr) # インストールした"pwr"パッケージを有効化する

pwr.t2n.test(n0[1], n0[2], 
             d = es, 
             sig.level = .05, 
             power = NULL, 
             alternative = "two.sided")

     t test power calculation 

             n1 = 161
             n2 = 217
              d = 0.2938698
      sig.level = 0.05
          power = 0.8045629
    alternative = two.sided

 以上の検定力や効果量の計算はいずれも二群の母分散が等しいと仮定した話であった。
次に,二群の母分散が等しいとは言えない場合について述べる。
 本文でも述べた様に,「平均的な標準偏差」を計算して利用する方法と,いずれか一方の標準偏差(分散)を「基準」として利用する方法がある。
男女別の世帯年収fincomeで例示しよう。

不等分散での検定力
m0 <- with(data01, by(fincome, sex, mean, na.rm = T))
v0 <- with(data01, by(fincome, sex, var,  na.rm = T))
n0 <- with(data01, by(fincome, sex, function(x) sum(!is.na(x))))
m0
sex: 1
[1] 697.4409
------------------------------------------------------------ 
sex: 2
[1] 737.0588
# 平均的な不偏分散を計算して,その平方根を平均的な標準偏差とする。
sd0 <- sqrt(sum(v0*(n0 - 1))/sum(n0 - 1))
sd1 <- c(sqrt(v0), sd0)
es1 <- diff(rev(m0))/sd1
sd1
       1        2          
377.5947 525.3713 467.9973 
pwr.t2n.test(n0[1], n0[2], 
             d = es1, 
             sig.level = .05, 
             power = NULL, 
             alternative = "two.sided")

     t test power calculation 

             n1 = 127
             n2 = 170
              d = 0.10492171, 0.07540930, 0.08465408
      sig.level = 0.05
          power = 0.14486618, 0.09821609, 0.11107465
    alternative = two.sided

 男性においてよりも女性において世帯年収の標準偏差が大きい。
大きな標準偏差を「基準」として用いるほど,それだけ偶然による変動のリスクが高まり,効果サイズは小さくなる(結果出力のd=の欄)。
効果サイズが小さくなれば,小さな効果を検出するのは難しくなるので,検定力も下がる(power=の欄)。
 しかしいずれにしても効果サイズは小さく,検定力は15%にも満たない。

検定可能な効果サイズ

検定力を.80(80%)として,このサンプルサイズで検定できる効果サイズを求めると,約0.33は必要である事が分かる。

pwr.t2n.test(n0[1], n0[2], 
             d = NULL, 
             sig.level = .05, 
             power = .80, 
             alternative = "two.sided")

     t test power calculation 

             n1 = 127
             n2 = 170
              d = 0.329665
      sig.level = 0.05
          power = 0.8
    alternative = two.sided
必要となる標本サイズ

 逆に,.085程度の小さな効果サイズを有意水準.05,検定力.80で検定しようと思うとどの程度の大きさの標本が必要かを計算してみよう。
 まずは,二群のサイズが等しい場合を標準のpower.t.test( )関数で計算する。
その次に,二群のサイズが異なる場合をpwr.t2n.test( )関数で計算する。

# 標準の検定力関数
power.t.test(n = NULL, 
             d = .085, 
             sig.level = .05, 
             power = .80, 
             alternative = "two.sided")

     Two-sample t test power calculation 

              n = 2173.661
          delta = 0.085
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group
# パッケイジpwr中のpwr.t2n.test( )関数
# 両方の標本サイズを同時に空欄には出来ないので,片方の標本サイズは3通りくらいで指定する
pwr.t2n.test(n1 = 1100, n2 = NULL, 
             d = .085, 
             sig.level = .05, power = .80,
             alternative = "two.sided")

     t test power calculation 

             n1 = 1100
             n2 = 87681.14
              d = 0.085
      sig.level = 0.05
          power = 0.8
    alternative = two.sided
pwr.t2n.test(n1 = 1500, n2 = NULL, 
             d = .085, 
             sig.level = .05,
             power = .80, 
             alternative = "two.sided")

     t test power calculation 

             n1 = 1500
             n2 = 3944.396
              d = 0.085
      sig.level = 0.05
          power = 0.8
    alternative = two.sided
pwr.t2n.test(n1 = 2100, n2 = NULL,
             d = .085,
             sig.level = .05, 
             power = .80,
             alternative = "two.sided")

     t test power calculation 

             n1 = 2100
             n2 = 2252.664
              d = 0.085
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

 二群のサイズが等しい場合には,それぞれの群のサイズが2174人,合計で4348人が必要である事が分かる。
 片方の群のサイズが小さい程,もう一つの群のサイズがかなり大きい事が必要になる。

“MBESS”パッケイジのconf.limits.nct( )関数

 最後に,標本効果サイズから母効果サイズの信頼区間を求める一つの方法を紹介しよう。
 MBESSパッケイジのconf.limits.nct( )関数は,第4章で解説した,非心t分布の非心パラメタの信頼区間を求める関数である。

# install.packages("MBESS", repos="https://cran.ism.ac.jp")
library("MBESS")

m0 <- with(data01, by(fincome, sex, mean, na.rm = T))
v0 <- with(data01, by(fincome, sex, var,  na.rm = T))
n0 <- with(data01, by(fincome, sex, function(x) sum(!is.na(x))))
m0
sex: 1
[1] 697.4409
------------------------------------------------------------ 
sex: 2
[1] 737.0588
sd0 <- sqrt(sum(v0*(n0 - 1))/sum(n0 - 1)) # 平均的な標準偏差
es0 <- diff(rev(m0))/sd0; es0 # 標本効果サイズ
          1 
-0.08465408 
t0  <- es0 / sqrt(sum(1/n0)) # t統計量
lambda <- conf.limits.nct(t.value = t0, 
                          df = (sum(n0) - 2), 
                          conf.level = .95, 
                          alpha.lower = NULL, 
                          alpha.upper = NULL)

lambda # 非心パラメタの下限と上限
$Lower.Limit
[1] -2.681983

$Prob.Less.Lower
    1 
0.025 

$Upper.Limit
[1] 1.239674

$Prob.Greater.Upper
    1 
0.025 
# 非心パラメタの上限と下限を効果サイズに換算する
lambda$Lower.Limit * sqrt(sum(1/n0)); lambda$Upper.Limit * sqrt(sum(1/n0))
[1] -0.3145632
[1] 0.1453983

発展1 正規性の検定(normality test)

 まずはシャピロ・ウィルクの検定とコルモゴロフ・スミルノフ検定の実例を示そう。
注意すべき点は,shapiro.test( )は変数が欠損値NAを含んでいても実行可能であるが,ks.test( )はNAを許容しないことである。

Shapiro-Wilk test
v1 <- data01$q1700b[data01$sex == 1] # 欠損値含むvector
shapiro.test(v1)

    Shapiro-Wilk normality test

data:  v1
W = 0.95646, p-value = 6.396e-05
Kolmogorov test
ks.test(v1, "pnorm", mean = mean(v1), sd = sd(v1))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  v1
D = NA, p-value = NA
alternative hypothesis: two-sided
v2 <- data01$q1700b[data01$sex == 1 & !is.na(data01$q1700b)] # 欠損値除外
ks.test(v2, "pnorm", mean = mean(v2), sd = sd(v2))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  v2
D = 0.15721, p-value = 0.0006992
alternative hypothesis: two-sided

 男女別の幸福度と世帯年収の正規性の検定は以下で実行できる。全て5%水準で有意となり,「母集団において正規分布している」とのゼロ仮説が棄却されてしまう。

v1 <- data01$q1700b[data01$sex == 1 & !is.na(data01$q1700b)]
shapiro.test(v1)

    Shapiro-Wilk normality test

data:  v1
W = 0.95646, p-value = 6.396e-05
ks.test(v1, "pnorm", mean = mean(v1), sd = sd(v1))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  v1
D = 0.15721, p-value = 0.0006992
alternative hypothesis: two-sided
v1 <- data01$q1700b[data01$sex == 2 & !is.na(data01$q1700b)]
shapiro.test(v1)

    Shapiro-Wilk normality test

data:  v1
W = 0.93864, p-value = 6.464e-08
ks.test(v1, "pnorm", mean = mean(v1), sd = sd(v1))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  v1
D = 0.18243, p-value = 1.068e-06
alternative hypothesis: two-sided
v1 <- data01$fincome[data01$sex == 1 & !is.na(data01$fincome)]
shapiro.test(v1)

    Shapiro-Wilk normality test

data:  v1
W = 0.93008, p-value = 5.623e-06
ks.test(v1, "pnorm", mean = mean(v1), sd = sd(v1))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  v1
D = 0.1215, p-value = 0.04704
alternative hypothesis: two-sided
v1 <- data01$fincome[data01$sex == 2 & !is.na(data01$fincome)]
shapiro.test(v1)

    Shapiro-Wilk normality test

data:  v1
W = 0.89074, p-value = 7.354e-10
ks.test(v1, "pnorm", mean = mean(v1), sd = sd(v1))

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  v1
D = 0.14576, p-value = 0.001458
alternative hypothesis: two-sided

 ヒストグラムとQ-Qプロットで視覚的に正規性を判断するスクリプトを紹介しよう。
本文の男性の幸福度のグラフの最低限のものは以下で描ける。
par (mfrow=c(1, 2))は,描画領域を1行2列に分割するというコマンドである。
例示したグラフではまず右側にhist(v1, freq = F)でヒストグラムを描き,そこにcurve(dnorm(x, mean=mean(v1), sd=sd(v1)), xlim = c(0, 10), add = T)で正規分布を重ね描きしている。
次に左側にqqnorm(v1)でQ-Qプロットを描き,qqline(v1)によって第1四分位と第3四分位を結んだ直線を重ねて描いている。

v1 <- data01$q1700b[data01$sex == 1 & !is.na(data01$q1700b)]
par(mfrow = c(1, 2))
hist(v1, freq = F)
curve(dnorm(x, mean = mean(v1), sd = sd(v1)), xlim = c(0, 10), add = T)
qqnorm(v1); qqline(v1)

女性の幸福度に例を変えて,もう少し装飾を施したグラフを描いてみよう。

v1 <- data01$q1700b[data01$sex == 2 & !is.na(data01$q1700b)]
par(mfrow = c(1, 2))

hist(v1, 
     main = "女性の幸福度ヒストグラム", 
     xlab="幸福度(女性)", 
     ylab="確率密度", 
     freq=F)

curve(dnorm(x, mean = mean(v1), sd = sd(v1)), 
      xlim = c(0, 10), add = T, col = "red")

qqnorm(v1, main = "女性の幸福度Q-Qプロット", sub = "一直線上に点が並ぶと良い")
qqline(v1, col = "red")

第6章の【練習問題】の解答

1) の答え

data01$q1700b <- c(0:10)[data01$q1700+1] # 欠損値処理

by(data01$q1700b, data01$edu1, var, na.rm=T)
data01$edu1: 0
[1] 3.767647
------------------------------------------------------------ 
data01$edu1: 1
[1] 2.630877
var.test(data01$q1700b ~ data01$edu1)

    F test to compare two variances

data:  data01$q1700b by data01$edu1
F = 1.4321, num df = 249, denom df = 121, p-value = 0.0265
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.043390 1.933145
sample estimates:
ratio of variances 
          1.432088 
t.test(data01$q1700b ~ data01$edu1)

    Welch Two Sample t-test

data:  data01$q1700b by data01$edu1
t = -1.5372, df = 282.24, p-value = 0.1254
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.67098830  0.08252928
sample estimates:
mean in group 0 mean in group 1 
        6.63200         6.92623 
t.test(data01$q1700b ~ data01$edu1, var.equal=T)

    Two Sample t-test

data:  data01$q1700b by data01$edu1
t = -1.4457, df = 370, p-value = 0.1491
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.6944225  0.1059635
sample estimates:
mean in group 0 mean in group 1 
        6.63200         6.92623 

2) の答え

平均値差は,大卒の幸福度-非大卒の幸福度 で計算する。

m0 <- by(data01$q1700b, data01$edu1, mean, na.rm=T); m0
data01$edu1: 0
[1] 6.632
------------------------------------------------------------ 
data01$edu1: 1
[1] 6.92623
v0 <- by(data01$q1700b, data01$edu1, var, na.rm=T)
n0 <- by(data01$q1700b, data01$edu1, function(x) sum(!is.na(x))); n0
data01$edu1: 0
[1] 250
------------------------------------------------------------ 
data01$edu1: 1
[1] 122
m_diff <- diff(m0); m_diff # 学歴別の幸福度の標本平均差
[1] 0.2942295
sigma <- sqrt(sum(v0*(n0 - 1))/(sum(n0)-2)) # 共通の母標準偏差の推定値
sd0 <- c(sqrt(v0), sigma); sd0
       0        1          
1.941043 1.621998 1.842795 
es <- m_diff/sd0; es # 標本効果サイズ
        0         1           
0.1515832 0.1813995 0.1596648 

3) の答え

 用いる標準偏差は3通りとする。検定力は30%以下であり,高くない。そもそも1)の検定で有意にならない程度の平均値差なのであるから,もしもそれを検出しようとするなら大きなサンプルサイズが必要である。

m0 <- by(data01$q1700b, data01$edu1, mean, na.rm=T); m0
data01$edu1: 0
[1] 6.632
------------------------------------------------------------ 
data01$edu1: 1
[1] 6.92623
v0 <- by(data01$q1700b, data01$edu1, var, na.rm=T)
n0 <- by(data01$q1700b, data01$edu1, function(x) sum(!is.na(x))); n0
data01$edu1: 0
[1] 250
------------------------------------------------------------ 
data01$edu1: 1
[1] 122
m_diff <- diff(m0); m_diff # 学歴別の幸福度の標本平均差
[1] 0.2942295
sigma <- sqrt(sum(v0*(n0 - 1))/(sum(n0)-2)) # 共通の母標準偏差の推定値
sd0 <- c(sqrt(v0), sigma); sd0
       0        1          
1.941043 1.621998 1.842795 
es <- m_diff/sd0; es # 標本効果サイズ
        0         1           
0.1515832 0.1813995 0.1596648 
pwr.t2n.test(n0[1], n0[2], d=es, sig.level=.05, power=NULL, alternative="two.sided")

     t test power calculation 

             n1 = 250
             n2 = 122
              d = 0.1515832, 0.1813995, 0.1596648
      sig.level = 0.05
          power = 0.2777058, 0.3740029, 0.3025693
    alternative = two.sided

4) の答え

スクリプトは本文にあるので,ここでは指示に従って2回繰り返す結果を示す(乱数の種seedを使用していないので実行する毎に結果が変わる)。母数の設定で明らかな様に,本当は母平均には差があるが,母標準偏差にもかなりの差がある。ケース数も異なる。この条件で,「母平均には差がない」と云うゼロ仮説を正しく棄却出来るか否かに関心がある。

n1 <- 20; n2 <- 80
mu1 <- 100; mu2 <- 115
sd1 <- 10; sd2 <- 40
# 1回目
x1 <-rnorm(n1, mean=mu1, sd=sd1)
x2 <-rnorm(n2, mean=mu2, sd=sd2)
t.test(x1, x2, var.equal=T)

    Two Sample t-test

data:  x1 and x2
t = -1.2201, df = 98, p-value = 0.2254
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -27.31703   6.51578
sample estimates:
mean of x mean of y 
  99.2602  109.6608 
t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -2.2286, df = 97.929, p-value = 0.02813
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -19.661927  -1.139322
sample estimates:
mean of x mean of y 
  99.2602  109.6608 
# 2回目
x1 <-rnorm(n1, mean=mu1, sd=sd1)
x2 <-rnorm(n2, mean=mu2, sd=sd2)
t.test(x1, x2, var.equal=T)

    Two Sample t-test

data:  x1 and x2
t = -1.9593, df = 98, p-value = 0.05291
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -37.7646232   0.2405548
sample estimates:
mean of x mean of y 
 102.3732  121.1352 
t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -3.6681, df = 96.134, p-value = 0.0004008
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -28.914856  -8.609213
sample estimates:
mean of x mean of y 
 102.3732  121.1352 

5) の答え

 ここではmu2は変えずに。n2とsd2を変えて検定結果がどの様に変わるかを示す。
 n2のサンプルサイズを半分にすると,ゼロ仮説を棄却出来ない事が多くなる(実行のたびに結果は変わる)。

n1 <- 20; n2 <- 40
mu1 <- 100; mu2 <- 115
sd1 <- 10; sd2 <- 40

x1 <-rnorm(n1, mean=mu1, sd=sd1)
x2 <-rnorm(n2, mean=mu2, sd=sd2)
t.test(x1, x2, var.equal=T)

    Two Sample t-test

data:  x1 and x2
t = -1.4037, df = 58, p-value = 0.1657
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -34.727004   6.097798
sample estimates:
mean of x mean of y 
 95.29327 109.60788 
t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -1.9349, df = 45.126, p-value = 0.05929
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -29.2141843   0.5849786
sample estimates:
mean of x mean of y 
 95.29327 109.60788 

 更にバラツキsd2を半分にすると,ゼロ仮説がかなり棄却出来る様になる。

n1 <- 20; n2 <- 40
mu1 <- 100; mu2 <- 115
sd1 <- 10; sd2 <- 20

x1 <-rnorm(n1, mean=mu1, sd=sd1)
x2 <-rnorm(n2, mean=mu2, sd=sd2)
t.test(x1, x2, var.equal=T)

    Two Sample t-test

data:  x1 and x2
t = -3.005, df = 58, p-value = 0.003919
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -20.417418  -4.091292
sample estimates:
mean of x mean of y 
 100.5294  112.7838 
t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -3.6145, df = 57.442, p-value = 0.0006347
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -19.042242  -5.466469
sample estimates:
mean of x mean of y 
 100.5294  112.7838 

同じ条件で何度も繰り返し実行してみると(シミュレイション),条件が同じであるにも関わらずかなり変動が大きい事が分かるだろう。これはサンプルサイズが比較的小さく設定してある事にもよる。逆に両方のサンプルサイズが大きければ,母平均の差がもっと小さくてもゼロ仮説をよく棄却する(大標本による検定は検定力が高い=感受性が高い)。試してみると良い。

ウェブ増補1 対応のある2群の母平均の差の検定

 第6章で解説したのは,「対応のない2群の母平均の差の検定」と呼ばれるものである。紙幅の都合上テクストでは割愛したが,2群の母平均の差の検定には,「対応のある2群」のものもある。ここで一応説明しておこう。

対応のある2群の母平均の差の検定

 以下のシミュレイションでは,まずサイズnの正規変量を整数化した変数xを生成し,そのxから平均1・分散1の正規変量を減じて整数化した変数yを生成する。一言で言えば,xとyはほぼ同じくらいの値であるが,yの方がやや小さくなる傾向がある。その事を,cor( ),plot( )で確認する。グラフから,yはxとほぼ同じ値であるが,xより僅かに小さい場合が多い事が確認出来る。

#### 対応のある二群の母平均の差の検定 ####
# 意図的に,xがyよりも僅かに大きくなる傾向のある2変数を作成する
n <- 20
x <- round(rnorm(n, mean=50, sd=20), 0)
y <- round(x - rnorm(n, mean=1, sd=1), 0)
cor(y, x); plot(y, x); abline(0, 1)
[1] 0.9985195

 まずはこの2変数xとyをそれぞれ別々に(あたかも対応が無いかの様に)要約し,箱ひげ図を描いてみよう。

summary(cbind(x, y)); boxplot(x, y) # それぞれの変数の個別の要約
       x               y        
 Min.   :19.00   Min.   :19.00  
 1st Qu.:32.75   1st Qu.:32.75  
 Median :44.50   Median :43.00  
 Mean   :46.20   Mean   :45.35  
 3rd Qu.:57.50   3rd Qu.:55.50  
 Max.   :88.00   Max.   :87.00  

t.test(x, y, var.equal=T) # 対応のない母平均差のスチューデントのt検定

    Two Sample t-test

data:  x and y
t = 0.1407, df = 38, p-value = 0.8888
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -11.37943  13.07943
sample estimates:
mean of x mean of y 
    46.20     45.35 
t.test(x, y) # 対応のない母平均差のウェルチのt検定

    Welch Two Sample t-test

data:  x and y
t = 0.1407, df = 38, p-value = 0.8888
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -11.37943  13.07943
sample estimates:
mean of x mean of y 
    46.20     45.35 

 箱ひげ図を見ても,x(1)とy(2)は殆ど同じ分布であり,要約統計量に見られる違いも僅かである。スチューデントのt検定とウェルチの検定の両方の結果を示しておいたがいずれにしても殆ど違わず(xとyがほぼ等分散である為),標本平均の小さな差は,統計的に全く有意にならない。つまり,「xとyの母平均は等しい」と云うゼロ仮説を棄却する証拠は全然得られなかったと云う結果になる。
 但し確認しておきたいのは,変数の作成の仕方や散布図からは,明らかにxはyより大きい傾向がある筈なのである。
 xとyはそもそも関係がある。こうした場合を「対応のある場合」と呼ぶ。実験データや社会調査データで言えば,同一個人から2回測定した場合,同一個人に二つの質問をした場合がその典型である。二つの変数を比較するのだが,一方の変数が明らかに他方の変数に基づいているか(実験における事前―事後の測定),或は同じ個人の性格や属性から2変数とも共通の影響を受けていると云った場合は,「対応のある場合」の検定を行うべきである。
 そうした場合は,各ケース(各ユニット)について,xとyの値の差を求めて,個人個人について計算されるその差について検討する。対応のない場合の母平均の差の検定は,2つの平均を求めてからその差について検討したが,対応のある場合の母平均の差の検定は,まず個々に差を求めてからその平均について検討するのである。
 よって以下ではまず,xとyの差x-yの要約統計量を求めて箱ひげ図を描く。これによって明らかにxはyよりも大きい傾向がある事が確認される。そしてその差についてt検定を行う。関数はt.test( )であり,これにオプションのpaired=TRUE を付けるだけで良い。

summary(x-y); boxplot(x-y) # 変数値の差の要約
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -1.00    0.00    1.00    0.85    1.00    3.00 

t.test(x, y, paired=T) # 対応のある母平均差の検定

    Paired t-test

data:  x and y
t = 3.6551, df = 19, p-value = 0.001683
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.3632698 1.3367302
sample estimates:
mean difference 
           0.85 

 検定結果は,1%水準で有意となり,平均における真の差(母平均差)は0ではないとの対抗仮説が採択される。差の95%信頼区間も求められているが,検定結果と平仄があっており,区間に0は含んでいない。本来僅かな差がある筈のデータについて,変数値に関連があると云う情報を考慮する事によって,その差を確かに検出する事が出来る様になったのである。

半分は自力でこの検定結果を再現する

 上では,t.test( )関数にpaired=Tオプションを指定するだけで済んだので簡単に実行出来たが,その分,何を計算しているのかなどは分からなくなってしまっている。よって以下で,上のt.test(x, y, paired=T) の結果を再現し,それを通して理解を深める事を目指す。
 上の検定は,「対応のある2群の母平均の差が0である」=「値の差x-yの平均が0である」と云うゼロ仮説を設定し,そこからのズレが偶然の範囲内と言えるかどうかを検討している。要するに,「或1つの変数(x-y)の母平均が0である」と云うゼロ仮説を検定しているだけなので,単純な「母平均の検定」に帰着する。ケース数はn,標準偏差はsd(x-y)で求められるので,標本平均の標準誤差(の推定値)は\(sd(x-y)/\sqrt{n}\) となり,ここから検定統計量\(t=\frac{mean(x-y)}{sd(x-y)/\sqrt{n}}\)が自由度n-1のt分布に近似的に従う事が言え,t統計量の値,それに対応する有意確率を求める事が出来る。以下がそのスクリプトで,95%信頼区間も求めている。

# 対応のある二群の母平均の差のt検定: マニュアル
mean(x-y) # 差の平均
[1] 0.85
t0 <- mean(x-y)/(sd(x-y)/sqrt(n)); t0 # t統計量
[1] 3.655147
2*(1-pt(abs(t0), df=n-1)) # 両側有意確率
[1] 0.001683286
# 95%信頼区間
mean(x-y)-qt(.975, df=n-1)*sd(x-y)/sqrt(n)
[1] 0.3632698
mean(x-y)+qt(.975, df=n-1)*sd(x-y)/sqrt(n)
[1] 1.33673

 先の”Paired t-test”の結果と見比べると,mean of the differences と mean(x-y)の値は0.7で一致し,t統計量の値(約3.0361),有意確率(約0.006796)も一致している。95%信頼区間の上限と下限も一致している事が確認出来る。「対応のある二群の母平均の差の検定」は,「引き算した値の母平均が0でないと言えるかどうか」を検定しているに過ぎない事がこれで確認出来た。

ウェブ増補2 対応のない2群の母比率の差の検定

 比率については,第3章の発展1で母比率の区間推定を紹介したのみだが,母平均の場合と同じ様に,二つの母比率に差があると言えるかどうかの検定も存在しており,やはり「対応のない場合」と「対応のある場合」に分けられる。この項目ではまず「対応のない場合」を説明する。

対応のない2群の母比率の差のカイ二乗検定

 例として,或TV番組についての,関東地区の世帯視聴率と関西地区の世帯視聴率を考えてみよう。調査世帯数は関東では900世帯,関西では600世帯であり,全く異なる別々の二つの集団から二値変数(0か1の値)を収集してその平均(つまり1の値を取る世帯の比率)を求めており,それぞれの個々の測定値には何の対応も関係も無い。検定の為のゼロ仮説は「関東の母比率と関西の母比率は等しい」である。
 関東での標本比率を.150,関西での標本視聴率を.130だったとしよう。以上の情報で,prop.test( )関数を用いて,ゼロ仮説の有意性検定を実行する事が出来る。この関数はデフォルトでは第5章発展1-2のイェーツの連続性修正同様の修正を行うので,それを行わない場合と行う(デフォルトの)場合の両方を示しておく。  

# 対応のない母比率の差の検定。関東(E)と関西(W)での視聴率
nE <- 900
nW <- 600
pE <- .150 # 計算の都合上小数第3位は0とする。
pW <- .130
prop.test(c(nE*pE, nW*pW), c(nE, nW), correct=F)
prop.test(c(nE*pE, nW*pW), c(nE, nW)) # デフォルトでは連続性修正

    2-sample test for equality of proportions without continuity correction

data:  c(nE * pE, nW * pW) out of c(nE, nW)
X-squared = 1.1819, df = 1, p-value = 0.277
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.0156135  0.0556135
sample estimates:
prop 1 prop 2 
  0.15   0.13 


    2-sample test for equality of proportions with continuity correction

data:  c(nE * pE, nW * pW) out of c(nE, nW)
X-squared = 1.0235, df = 1, p-value = 0.3117
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.01700239  0.05700239
sample estimates:
prop 1 prop 2 
  0.15   0.13 

 またしても,自動で検定してくれるのは便利ではあるが,正直何をどうしたのかは分からない。結果から,カイ二乗検定を行っているらしいと云う事は読み取れるが,その値がどの様に計算されたのかは分からない。因みに,上の結果はいずれも,関東と関西の母視聴率に差があるとは言えないと云う検定結果を示している。900世帯の標本比率15%と600世帯の標本比率13%では,母集団でも関東の方が視聴率が良かったと結論するには不十分なのである。比率の差の95%信頼区間も求められているが,0を挟んでマイナスからプラスの領域に広がっており,「関西の方が高いかも知れないし,関東の方が高いかも知れない」と云う結果を意味しているのである。

分割表を用いた比率の差のカイ二乗検定

 カイ二乗検定を用いているのだから,すぐに思いつくのは分割表(クロス表)の独立性の検定である。一方の変数は地域(関東と関西),他方の変数は視聴の有無なので,2行×2列の分割表で表現する事が出来,2変数に関連が無いかどうか(つまり2つの地域で比が等しいか否か)を検定する事が出来る筈である。以下のスクリプトはそれを実行するものであり,結果を上のprop.test( )と見比べると,連続性の修正をする場合としない場合でそれぞれ全く一致している事が分かる。
 第5章発展1-2で説明した通り,分割表の独立性のカイ二乗検定を行う代わりに,より正確なフィッシャーの正確検定(直接確率)を用いる事も出来る。

# 分割表を用いたカイ二乗検定
x1 <- matrix(c(nE*pE, nW*pW, nE*(1-pE), nW*(1-pW)), ncol=2); x1
chisq.test(x1, correct=F) # イェーツの連続性修正無し
chisq.test(x1) # イェーツの連続性修正有り
fisher.test(x1) # どうせなら正確検定
     [,1] [,2]
[1,]  135  765
[2,]   78  522

    Pearson's Chi-squared test

data:  x1
X-squared = 1.1819, df = 1, p-value = 0.277


    Pearson's Chi-squared test with Yates' continuity correction

data:  x1
X-squared = 1.0235, df = 1, p-value = 0.3117


    Fisher's Exact Test for Count Data

data:  x1
p-value = 0.2911
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.866734 1.616939
sample estimates:
odds ratio 
  1.180866 

対応のない2群の母比率の差のz検定

 分割表で検定出来てしまうのでこれ以上追究する意味も余り無いと感じられるが,他によく紹介される方法にも触れておく。
 ゼロ仮説が正しければそれぞれの母集団での視聴率は等しいと仮定されるので,その共通の母比率の推定値Pを,標本比率の加重平均として求める。そうすると母分散の推定値は\(P(1-P)\)となり,標本平均の一種とも言える標本比率の差の検定は,第6章基礎1-2のt統計量の計算と同様にして検定統計量を求める事が出来る。そこからするとこの検定統計量はt分布に従うと考えるべきであるが,比率の場合は母分散が小数になり,標本分散による母分散の過小推定の問題も小さいからか,標準正規分布に従うと考えて検定を行う。
 検定統計量が標準正規変量と見做せるならば,その二乗は,確率分布の定義から,自由度1のカイ二乗分布に従う事になる。よってカイ二乗分布を用いて検定を行っても全く同じである。これらをRで実行してみる。後半は連続性修正を行う方法を定義に従って実行している。

# 母比率は等しいとのゼロ仮説の下での,共通の母比率推定値P
P <- (nE*pE+nW*pW)/(nE+nW)
z <- abs(pE-pW)/sqrt(P*(1-P)*(1/nE+1/nW)); z # 標準正規変量
[1] 1.08716
2*(1-pnorm(z, 0, 1)) # 標準正規分布を用いた両側有意確率
[1] 0.2769661
2*(1-pt(z, df=nE+nW-2)) # t分布を用いた両側有意確率
[1] 0.277141
z^2 # 自由度1のカイ二乗統計量
[1] 1.181917
1-pchisq(z^2, df=1) # カイ二乗分布を用いた有意確率
[1] 0.2769661
# 比率の差の場合の連続性修正の公式
zc <- (abs(pE-pW)-.5*(1/nE+1/nW))/sqrt(P*(1-P)*(1/nE+1/nW)); zc # 標準正規変量
[1] 1.011663
2*(1-pnorm(zc, 0, 1)) # 標準正規分布を用いた両側有意確率
[1] 0.3116994
zc^2 # 自由度1のカイ二乗統計量
[1] 1.023461
1-pchisq(zc^2, df=1) # カイ二乗分布を用いた有意確率
[1] 0.3116994

 標準正規分布を用いた検定結果とカイ二乗分布を用いた検定結果は,有意確率が完全に一致しており,同じ事を行っているのだと云う事が納得出来るであろう。前半のt分布を用いた場合も有意確率の相違は極めて僅かである。
 prop.test( )の結果と見比べても,カイ二乗値や有意確率が一致している事が確認出来る。つまり,先のprop.test( )のカイ二乗検定は,母比率の差のz検定と等価であると言える。そして母比率の差のz検定は,母平均の差のt検定と同じ公式から求められるのである。対応のない二群の母比率の差の検定は,比率が0-1変数の平均である事に着目して,対応のない二群の母平均の差の検定を適用しても近似的に実行出来るのである。
 2行×2列の分割表に表現して独立性の検定を行った結果とも完全に一致している。検定するだけであれば分割表のカイ二乗検定が一番簡単にも思えるが,比率の差の区間推定を行いたい場合には,prop.test( )を用いるか,自分でz統計量を求めて信頼区間を求める事が必要になる。
 この様に,実質的には同じ事を検定する為に,一見色々と異なった方法が存在する事がある。それぞれがどの様に関係しているか,どの様に一致するかをこうしてきちんと確認してゆくと,単に公式を機械的に当てはめるのではなく,また検定や推定の中身をブラックボックス化させずに理解を深める事が出来るだろう。

ウェブ増補3 対応のある2群の母比率の差の検定

 対応のある2群と云うのはウェブ増補1で書いた通り,同一個人から2回測定した場合,同一個人に2つの質問をした場合などを指す。世帯視聴率で言えば,同じ関東地区の900世帯で,異なるドラマの視聴率(比率)に差があると言えるかどうかは,対応のある場合となる。
 対応のある2群の母比率の差の検定にも,3通りくらいの方法がある。一つは標準正規分布を用いた検定法,2つ目はマクネマー検定と呼ばれる有名な方法,3つ目は2項分布を用いた方法である。標準正規分布を用いる方法は標準誤差の計算が簡単には分からない。マクネマー検定は最もよく紹介されている様に思えるが,ユーザにとってはかなりブラックボックスな感覚になる。結局,2項分布を用いた方法が最もすっきりしているのではないかと思われる。
 関東地区の900世帯で世帯視聴率調査を行って,同じ週に放映された二つのドラマの視聴率を比較したとしよう。一方のドラマ1は16.0%,他方のドラマ2は13.0%であった。これらの視聴率は標本比率である。母集団は関東地区の約1,835万世帯であるが,母集団における視聴率(母比率)もドラマ1の方が高かったと言えるだろうか。(視聴率調査最大手のビデオリサーチについては,週間高世帯視聴率番組10や視聴率調査について(視聴率ハンドブック)を参照の事。)
 これがもし「対応のない2群」であれば,上の「ウェブ増補2」で紹介した様に,分割表にして独立性のカイ二乗検定を行えばよい。連続性修正ありの場合で言えば,\(χ^2_{(df=1)}=3.03\), 有意確率\(p=.082\)で5%で有意にならない。つまり,ドラマ1の方が母集団での視聴率が高かったと云う証拠としては不十分である。  

# 二つのTV番組の視聴率の例
n <- 900
p1 <- .160 # 便宜上,視聴率は小数第二位までとしておく
p2 <- .130
n1 <- n*p1; n2 <- n*p2; n1; n2
[1] 144
[1] 117
# 「対応がない場合」で計算してみる
x0 <- matrix(c( n1, n2, n-n1, n-n2), ncol=2)
addmargins(x0)
              Sum
    144  756  900
    117  783  900
Sum 261 1539 1800
chisq.test(x0)

    Pearson's Chi-squared test with Yates' continuity correction

data:  x0
X-squared = 3.0293, df = 1, p-value = 0.08177

 ここで重要なのは,同じ900世帯に対して,二つのTV番組の視聴率を調べていると云う点であり,世帯数は全部で900である。上の様に対応の無い場合のカイ二乗検定を行うとあたかも全部で1800世帯あるかの様になってしまう。
 論理的には,900世帯は,ドラマ1も2も見た世帯,1か2の一方のみ見た世帯,いずれも見なかった世帯の4通りに分類される。ドラマ1を見た/見ていない,ドラマ2を見た/見ていないの2行×2列の分割表で表現される筈である。これが対応のある場合の比率の差を考える出発点である。
 しかし,視聴率調査では,そこまで詳しい情報は一般には公開されていないので,我々にはこの分割表は作る事が出来ない。
 そこで以下では,3つの対照的な分割表xa, xb, xcについて分析して比較しよう。3つの分割表は全て周辺度数は(当然)同一であるが,分割表の中身は大きく異なる。

n <- 900
p1 <- .160; p2 <- .130
n1 <- n*p1; n2 <- n*p2

# ドラマ2を見た人は全てドラマ1も見たと云う極端なケースa
xa <- matrix(c((a11 <- min(n1, n2)), n2-a11, n1-a11, n-(n1+n2-a11)), ncol=2)
addmargins(xa)
            Sum
    117  27 144
      0 756 756
Sum 117 783 900
# ほぼ完全独立の極端なケースb
xb <- matrix(c((b11 <- 19), n2-b11, n1-b11, n-(n1+n2-b11)), ncol=2)
addmargins(xb)
            Sum
     19 125 144
     98 658 756
Sum 117 783 900
# 完全分離の極端なケースc
xc <- matrix(c((c11 <- 0), n2-c11, n1-c11, n-(n1+n2-c11)), ncol=2)
addmargins(xc)
            Sum
      0 144 144
    117 639 756
Sum 117 783 900

 対応のある2群の母比率の検定では,「ドラマ1もドラマ2も見た」人達,「両方とも見なかった」人達は考慮に入れない。この人達は母比率に差をもたらす事が定義上無いからである。注目されるのは,「ドラマ1は見たがドラマ2は見なかった」人達(1行2列)と「ドラマ1は見なかったがドラマ2は見た」人達(2行1列)だけである。例えば分割表xaでは1行2列は27人,2行1列は0人,合計で27人であるが,「この人達は本来,いずれのセルに入るかは五分五分であり,この標本ではたまたま1行2列に27人が入った」と云う確率が小さいか,そうでないかを検定する。分割表xbでは125人と98人で合計223人,分割表xcでは144人と117人で合計261人が対象である。

標準正規分布を用いた検定

 ドラマ1の視聴世帯数をn1,ドラマ2の視聴世帯数をn2,両方のドラマを見た世帯数をa11(またはb11, c11)とすると,ドラマ1だけを見た世帯数(1行2列)はn1-a11(またはn1-b11, n1-c11),ドラマ2だけを見た世帯数(2行1列)はn2-a11(またはn2-b11, n2-c11)であり,いずれか一方だけを見た世帯の総数は\(n1+n2-2*a11\)(または\(n1+n2-2*b11\), \(n1+n2-2*c11\))となる。
 標準正規分布による検定では,1行2列(または2行1列でも良い)に入る人数は本来は\(n1+n2-2*a11\)の半分であるが(=期待値),それが偶然にn1-a11になったと云う確率が小さいか否かを考える。この人数の標準誤差は\(\sqrt{(n1+n2-2*a11)*.5*(1-.5)}\)となる事が分かっているので2項分布の平均と分散,確率的に変動する「1行2列の人数」を標準化すると\(((n2-n11)-.5*(n1+n2-2*n11))/\sqrt{(n1+n2-2*n11)*.5*.5}\)となる。これが標準正規分布に従う事を利用して検定を行う。
 以下に3つの場合全てについてこの標準正規分布検定を示した。後でマクネマー検定との関係を示す為に,標準正規変量を二乗した値(自由度1のカイ二乗変量になる)も示してある。

# 標準正規変量によるz検定.標準誤差sqrt(n*p*(1-p))は二項分布の成功回数のもの
za <- ((n2-a11)-.5*(n1+n2-2*a11))/sqrt((n1+n2-2*a11)*.5*.5); za
[1] -5.196152
2*(1-pnorm(abs(za), 0, 1)) # 標準正規分布の両側有意確率
[1] 2.034555e-07
za^2 # 標準正規変量の二乗=自由度1のカイ二乗変量
[1] 27
zb <- ((n2-b11)-.5*(n1+n2-2*b11))/sqrt((n1+n2-2*b11)*.5*.5); zb
[1] -1.808054
2*(1-pnorm(abs(zb), 0, 1)) # 標準正規分布の両側有意確率
[1] 0.07059814
zb^2 # 標準正規変量の二乗=自由度1のカイ二乗変量
[1] 3.269058
zc <- ((n2-c11)-.5*(n1+n2-2*c11))/sqrt((n1+n2-2*c11)*.5*.5); zc
[1] -1.671258
2*(1-pnorm(abs(zc), 0, 1)) # 標準正規分布の両側有意確率
[1] 0.09467072
zc^2 # 標準正規変量の二乗=自由度1のカイ二乗変量
[1] 2.793103

 本来半々の筈が27対0になる事はおよそ偶然とは考えられないが,偶然に125対98(約1.28倍)や144対117(約1.23倍)になる事は20回に1回よりは多く起こり得るのでそれ程おかしくはないと云う結果になっている。

マクネマー検定(McNemar test)

 こうした検定を行う場合,よく紹介されているのはマクネマー検定と云う方法である。Rにもmcnemar( )関数が用意されているが,2行2列の分割表の独立性の検定のイェーツの連続性修正と同様,デフォルトでは連続性修正を行う。以下ではa,b,cのそれぞれについて連続性の修正を行わない場合と行う場合の両方の結果を示す。

# マクネマー検定
mcnemar.test(xa, correct=F)

    McNemar's Chi-squared test

data:  xa
McNemar's chi-squared = 27, df = 1, p-value = 2.035e-07
mcnemar.test(xa)

    McNemar's Chi-squared test with continuity correction

data:  xa
McNemar's chi-squared = 25.037, df = 1, p-value = 5.624e-07
mcnemar.test(xb, correct=F)

    McNemar's Chi-squared test

data:  xb
McNemar's chi-squared = 3.2691, df = 1, p-value = 0.0706
mcnemar.test(xb)

    McNemar's Chi-squared test with continuity correction

data:  xb
McNemar's chi-squared = 3.0314, df = 1, p-value = 0.08167
mcnemar.test(xc, correct=F)

    McNemar's Chi-squared test

data:  xc
McNemar's chi-squared = 2.7931, df = 1, p-value = 0.09467
mcnemar.test(xc)

    McNemar's Chi-squared test with continuity correction

data:  xc
McNemar's chi-squared = 2.59, df = 1, p-value = 0.1075

 これも自動的に検定を行ってくれるのでとても便利であるが,カイ二乗検定である事以外はここからはよく分からない。
 ただ,先の標準正規分布による検定のp値及びカイ二乗値(標準正規変量の二乗)と(修正なしの)マクネマー検定の結果を見比べると,3つとも一致している事が確認出来る。マクネマー検定の方は連続性の修正を行う事は出来るが,基本的には同じ様な事をしているのだと理解する事が出来る(検定結果については最後に纏める)。

2項検定(Binomial test)

 以上二つの考え方の基本になっているのは,いずれか一方だけを見た世帯の総数\(n1+n2-2*a11\)(または\(n1+n2-2*b11\), \(n1+n2-2*c11\))は成功確率.5の2項分布に従って1行2列のセルもしくは2行1列のセルに割り当てられると云う考えであり,実際の結果が偶然の範囲内だと言えるかどうかを調べるのが検定である。
 であれば,最初から素直に2項検定を行うのが分かり易い。Rには2項検定を簡単に実行するbinom.test( )関数がある。

# 2項検定
binom.test(n1-a11, n1+n2-2*a11, p=.5)

    Exact binomial test

data:  n1 - a11 and n1 + n2 - 2 * a11
number of successes = 27, number of trials = 27, p-value = 1.49e-08
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.8722971 1.0000000
sample estimates:
probability of success 
                     1 
binom.test(n1-b11, n1+n2-2*b11, p=.5)

    Exact binomial test

data:  n1 - b11 and n1 + n2 - 2 * b11
number of successes = 125, number of trials = 223, p-value = 0.08144
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4927172 0.6267280
sample estimates:
probability of success 
             0.5605381 
binom.test(n1-c11, n1+n2-2*c11, p=.5)

    Exact binomial test

data:  n1 - c11 and n1 + n2 - 2 * c11
number of successes = 144, number of trials = 261, p-value = 0.1074
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4891790 0.6130813
sample estimates:
probability of success 
             0.5517241 

 number of trials が問題の2つのセルの度数合計,number of successes が1行2列の実際の度数,そこに当てはまる確率が.5である。 p-value がその2項検定の結果の有意確率である。やはりxaでは高度に有意,xb,xcでは5%で有意にならない。
 最後のxcの場合を用いて,この2項検定をもう少し丁寧に確認してみよう。

# 2項検定の内容の確認
pbinom(143, size=261, prob=.5, lower.tail=F) # 144以上に大きくなる確率
[1] 0.05367957
pbinom(261-144, size=261, prob=.5) # 117以下に小さくなる確率
[1] 0.05367957
pbinom(261-144, size=261, prob=.5)*2 # 両側確率
[1] 0.1073591

 多い方のセルの度数以上になるか,もしくは少ない方のセルの度数以下になる確率の合計が両側有意確率であり,これを pbinom( )関数で求めると,上の binom.test( )の結果に一致している。全部自動で実行してくれる関数はとても便利であるが,少なくとも最初は,それが何をやっているのかをこの様に確認しながら理解していく様にしよう。
 最後に,以上の数種類の検定結果を簡単に纏めて比較しておく。

各種検定での有意確率・一覧
検定方法 分割表(xa) 分割表(xb) 分割表(xc)
標準正規検定 2.034555e-07 .07059814 .09467072
マクネマー検定(修正なし) 2.035e-07 .0706 .09467
マクネマー検定(修正あり) 5.624e-07 .08167 .1075
2項検定 1.49e-08 .08144 .1074

 xaの様な極端な場合を除いて,xb,xcでは標準正規検定や連続性修正無しのマクネマー検定では2項検定よりも有意確率が低く,ゼロ仮説が棄却されやすい傾向が見られる。実際より有意確率が低くなると云う事は,第1種の過誤(αエラー)を犯しやすくなると云う事である。但し,この傾向は絶対ではなく,上の例でも,1行1列のセルの度数を色々と変えて行ってみると,p(マクネマー修正なし)>p(2項検定)となる場合も無い訳では無い(一方のセルの度数がかなり小さくて有意確率が極めて小さくなる条件下)。
 連続性修正を行ったマクネマー検定は2項検定の結果にかなり近い。Rでは連続性修正を行うのがデフォルトでもあり,簡単に「対応のある2群の母比率の差の検定」を行う為に適してはいるが,2項検定も同じくらい簡単に実行出来るので,敢えてマクネマー検定を行う意味があるかどうかは不明である。