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

\[z=\frac{(\bar{y}_{1}-\bar{y}_{2})-(\color{magenta}{\mu_{1}}-\color{magenta}{\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{magenta}{\mu_{1}}-\color{magenta}{\mu_{2}})}{\sigma\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}=\frac{(\bar{y}_{1}-\bar{y}_{2})-(\color{magenta}{\mu_{1}}-\color{magenta}{\mu_{2}})}{\sigma}\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}\;\sim\;N(0, 1^2)\]

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

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

data01 <- read.csv("practice.csv") # 先にworking directoryにデータファイルを置いておく
names(data01) # 読み込んだデータフレイムの中の変数名を確認

# まず最初に使用する変数の度数分布表を確認しておく。欠損値もあれば表示させる。
table(data01$sex, useNA="ifany")
table(data01$q1700, useNA="ifany")

 [1] "id"      "mode"    "sex"     "age"     "q0101"   "q0102"   "q0103"  
 [8] "q0200"   "q0301"   "q0302"   "q0401"   "q0402"   "q0403"   "edu1"   
[15] "edu2"    "edu"     "job"     "q0700"   "q0900"   "q1000"   "q1101"  
[22] "q1102"   "q1400"   "q1501"   "q1502"   "q1601"   "q1602"   "q1603"  
[29] "q1604"   "q1605"   "q1700"   "q1801"   "q1802"   "q1900"   "q2100"  
[36] "q2200"   "q2601"   "q2602"   "income"  "fincome" "q3000"  

  1   2 
165 219 

   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] # 変数のリコード法にはやや工夫をしている。
table(data01$q1700b, useNA="always") # 変換の結果を確認

   0    1    2    3    4    5    6    7    8    9   10 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    6 

# 男女別の幸福度の標本平均，不偏分散，有効ケース数
by(data01$q1700b, data01$sex, mean, na.rm=T)
by(data01$q1700b, data01$sex, var, na.rm=T)
by(data01$q1700b, data01$sex, function(x) sum(!is.na(x)))

data01$sex: 1
[1] 6.36646
-------------------------------------------------------- 
data01$sex: 2
[1] 6.917051
data01$sex: 1
[1] 4.083618
-------------------------------------------------------- 
data01$sex: 2
[1] 3.08568
data01$sex: 1
[1] 161
-------------------------------------------------------- 
data01$sex: 2
[1] 217

# tapply()関数でも同様にできる。
tapply(data01$q1700b, data01$sex, mean, na.rm=T)
tapply(data01$q1700b, data01$sex, var, na.rm=T)
tapply(data01$q1700b, data01$sex, function(x) sum(!is.na(x)))

       1        2 
6.366460 6.917051 
       1        2 
4.083618 3.085680 
  1   2 
161 217 

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

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

    F test to compare two variances

data:  data01$q1700b by data01$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 

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

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

var0 <- by(data01$q1700b, data01$sex, var, na.rm=T)
sum(var0*df0)/sum(df0)

[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{magenta}{\mu_{1}}-\color{magenta}{\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 <- t.test(data01$q1700b ~ data01$sex, var.equal=T)
test1

    Two Sample t-test

data:  data01$q1700b by data01$sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means 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 

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

boxplot(data01$q1700b ~ data01$sex, las=1, cex.main=1,
    main="男女別幸福度", sub="実線は中央値，赤い線は算術平均", varwidth=T,
    names=c("男性","女性"), par(family="serif"), horizontal=T)
m0 <- by(data01$q1700b, data01$sex, mean, na.rm=T)
segments(m0[1], 0.65, m0[1], 1.35, col="red")
segments(m0[2], 1.6, m0[2], 2.4, col="red")

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

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検定を行おう。

table(data01$fincome, useNA="always")
by(data01$fincome, data01$sex, mean, na.rm=T)
by(data01$fincome, data01$sex, var, na.rm=T)
df0 <- by(data01$fincome, data01$sex, function(x) sum(!is.na(x))); df0

var.test(data01$fincome ~ data01$sex)

   0   25   75  125  200  300  400  500  600  700  800  925 1125 1375 1750 
   5    1    7    8   17   19   28   36   32   27   23   33   31   13   11 
2500 <NA> 
   6   87 
data01$sex: 1
[1] 697.4409
-------------------------------------------------------- 
data01$sex: 2
[1] 737.0588
data01$sex: 1
[1] 142577.7
-------------------------------------------------------- 
data01$sex: 2
[1] 276015
data01$sex: 1
[1] 127
-------------------------------------------------------- 
data01$sex: 2
[1] 170

    F test to compare two variances

data:  data01$fincome by data01$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, df0[1]-1, df0[2]-1), type="l", bty="n",
    xlab="", ylab="", family="serif",
    main=sprintf("自由度%d, %dのF分布", df0[1]-1, df0[2]-1), 
    sub=sprintf("両側%d%%棄却域", p0*100))
segments(x1 <- c(seq(0, qf(p0/2, df0[1]-1, df0[2]-1), by=.01),
                        seq(qf(1-p0/2, df0[1]-1, df0[2]-1), 2, by=.01)), 0,
              x1, df(x1, df0[1]-1, df0[2]-1), col="red")
axis(side=1, las=2, family="serif", font=2, col.axis="red",
     at=round(c(qf(p0/2, df0[1]-1, df0[2]-1), qf(1-p0/2, df0[1]-1, df0[2]-1)), 2))

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

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

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

t.test(data01$fincome ~ data01$sex)

    Welch Two Sample t-test

data:  data01$fincome by data01$sex
t = -0.756, df = 294.6, p-value = 0.4503
alternative hypothesis: true difference in means 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検定と等価な分析を行う方法を紹介しています。

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

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

m0 <- by(data01$q1700b, data01$sex, mean, na.rm=T)
var0 <- by(data01$q1700b, data01$sex, var, na.rm=T)
n0 <- by(data01$q1700b, data01$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))

# 信頼区間上限値
# m_diff - qt(.025, df=sum(n0)-2)*sigma*sqrt(sum(1/n0))

# t検定の結果
test1 <- t.test(data01$q1700b ~ data01$sex, var.equal=T)
test1

[1] -0.9337912 -0.1673909

    Two Sample t-test

data:  data01$q1700b by data01$sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means 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 <- by(data01$fincome, data01$sex, mean, na.rm=T)
var0 <- by(data01$fincome, data01$sex, var, na.rm=T)
n0 <- by(data01$fincome, data01$sex, function(x) sum(!is.na(x)))

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

# ウェルチ検定の自由度計算
df1 <- prod(n0-1)*sum(var0/n0)^2/((n0[2]-1)*(var0[1]/n0[1])^2+(n0[1]-1)*(var0[2]/n0[2])^2)

m_diff - qt(c(.975, .025), df=df1)*se0 # 母平均差の95%信頼区間下限値と上限値
# m_diff - qt(.025, df=df1)*se0 # 信頼区間上限値

# Welch検定の結果と比較する
test2 <- t.test(data01$fincome ~ data01$sex); test2

[1] -142.7534   63.5176

    Welch Two Sample t-test

data:  data01$fincome by data01$sex
t = -0.756, df = 294.6, p-value = 0.4503
alternative hypothesis: true difference in means 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 <- by(data01$q1700b, data01$sex, mean, na.rm=T); m0
v0 <- by(data01$q1700b, data01$sex, var, na.rm=T); v0
n0 <- by(data01$q1700b, data01$sex, function(x) sum(!is.na(x))); n0

m_diff <- diff(rev(m0)); m_diff # 男女の幸福度の標本平均差
sigma <- sqrt(sum(v0*(n0 - 1))/(sum(n0)-2)) # 共通の母標準偏差の推定値
es <- m_diff/sigma; es # 標本効果サイズ

data01$sex: 1
[1] 6.36646
-------------------------------------------------------- 
data01$sex: 2
[1] 6.917051
data01$sex: 1
[1] 4.083618
-------------------------------------------------------- 
data01$sex: 2
[1] 3.08568
data01$sex: 1
[1] 161
-------------------------------------------------------- 
data01$sex: 2
[1] 217
         1 
-0.5505911 
         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

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

Warning: package 'pwr' was built under R version 3.5.3

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

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

# 平均的な不偏分散を計算して，その平方根を平均的な標準偏差とする。
sd0 <- sqrt(sum(v0*(n0-1))/sum(n0-1))
sd1 <- c(sqrt(v0), sd0)
es1 <- diff(rev(m0))/sd1
sd1

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

data01$sex: 1
[1] 697.4409
-------------------------------------------------------- 
data01$sex: 2
[1] 737.0588
       1        2          
377.5947 525.3713 467.9973 

     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

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

# 標準の検定力関数
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

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

Warning: package 'MBESS' was built under R version 3.5.3

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

data01$sex: 1
[1] 697.4409
-------------------------------------------------------- 
data01$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) の答え

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 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 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 = -0.91526, df = 98, p-value = 0.3623
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -23.968588   8.837886
sample estimates:
mean of x mean of y 
 103.5450  111.1104 

t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -1.7068, df = 96.585, p-value = 0.09108
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -16.36325   1.23255
sample estimates:
mean of x mean of y 
 103.5450  111.1104 

# 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 = -2.4476, df = 98, p-value = 0.01616
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -38.630349  -4.037095
sample estimates:
mean of x mean of y 
 97.27848 118.61220 

t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -4.6035, df = 95.522, p-value = 1.279e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -30.53321 -12.13424
sample estimates:
mean of x mean of y 
 97.27848 118.61220 

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.7247, df = 58, p-value = 0.0899
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -33.036404   2.455695
sample estimates:
mean of x mean of y 
 97.88225 113.17260 

t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -2.2834, df = 52.345, p-value = 0.0265
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -28.725432  -1.855277
sample estimates:
mean of x mean of y 
 97.88225 113.17260 

　更にバラツキ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.5706, df = 58, p-value = 0.0007237
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -27.547288  -7.755944
sample estimates:
mean of x mean of y 
 101.2986  118.9502 

t.test(x1, x2)

    Welch Two Sample t-test

data:  x1 and x2
t = -4.4494, df = 57.858, p-value = 3.973e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -25.593179  -9.710053
sample estimates:
mean of x mean of y 
 101.2986  118.9502 

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

対応のある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.9987412

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

summary(cbind(x, y)); boxplot(x, y) # それぞれの変数の個別の要約

       x               y        
 Min.   :29.00   Min.   :27.00  
 1st Qu.:41.25   1st Qu.:39.25  
 Median :54.00   Median :52.50  
 Mean   :53.20   Mean   :51.90  
 3rd Qu.:62.75   3rd Qu.:61.25  
 Max.   :86.00   Max.   :85.00  

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

    Two Sample t-test

data:  x and y
t = 0.23332, df = 38, p-value = 0.8168
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.97948 12.57948
sample estimates:
mean of x mean of y 
     53.2      51.9 

t.test(x, y) # 対応のない母平均差のウェルチのt検定

    Welch Two Sample t-test

data:  x and y
t = 0.23332, df = 37.991, p-value = 0.8168
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.979565 12.579565
sample estimates:
mean of x mean of y 
     53.2      51.9 

　箱ひげ図を見ても，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.0     1.0     1.0     1.3     2.0     3.0 

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

    Paired t-test

data:  x and y
t = 6.2962, df = 19, p-value = 4.815e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.8678446 1.7321554
sample estimates:
mean of the differences 
                    1.3 

　検定結果は，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] 1.3

t0 <- mean(x-y)/(sd(x-y)/sqrt(n)); t0 # t統計量

[1] 6.296187

2*(1-pt(abs(t0), df=n-1)) # 両側有意確率

[1] 4.815134e-06

#　95%信頼区間
mean(x-y)-qt(.975, df=n-1)*sd(x-y)/sqrt(n)

[1] 0.8678446

mean(x-y)+qt(.975, df=n-1)*sd(x-y)/sqrt(n)

[1] 1.732155

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

対応のない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統計量を求めて信頼区間を求める事が必要になる。
　この様に，実質的には同じ事を検定する為に，一見色々と異なった方法が存在する事がある。それぞれがどの様に関係しているか，どの様に一致するかをこうしてきちんと確認してゆくと，単に公式を機械的に当てはめるのではなく，また検定や推定の中身をブラックボックス化させずに理解を深める事が出来るだろう。

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

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