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

第5章　2変数の関連の推定と検定

0 全体の章構成

• 第 1章　1変量の記述統計の基礎
• 第 2章　2変数の関連の記述統計
• 第 3章　推測統計の基礎
• 第 4章　統計的検定の一般型
• 第 5章　2変数の関連の推定と検定
• 第 6章　2群の母平均の差のt検定
• 第 7章　平均値の差の分散分析(ANOVA)
• 第 8章　単回帰分析(SRA)
• 第 9章　重回帰分析（Ⅰ）(MRA)
• 第10章　重回帰分析（Ⅱ）(MRA)
• 第11章　主成分分析と因子分析
• 第12章　ロジスティック回帰分析
• 第13章　階層構造のあるデータの分析

1 積率相関係数についての推測統計

積率相関

\[r_{xy}=\frac{s_{xy}}{s_{x}s_{y}}=\frac{1}{n}\sum_{i=1}^{n}(\frac{x_{i}-\bar{x}}{s_{x}})(\frac{y_{i}-\bar{y}}{s_{y}})\]

　本文の標本相関係数の標本抽出分布（一部）を描くスクリプトを紹介しておこう。単回帰分析を深く理解していないと意味が分からないと思うので，解読できなくて良い。シミュレイションなので，実行するたびに少しずつ違う結果が得られる。NやR，n，timesを変えて実行してみると良い。

N <- 1000 # 母集団サイズ
R <- 0.3 # 欲しい母相関係数の値

y <- rnorm(N, mean=0, sd=1)
x <- rnorm(N, mean=0, sd=1)

slope_xy <- cor(x, y, method="pearson")*sd(y)/sd(x)
intercept_xy <- mean(y) - slope_xy*mean(x)
e <- y - (intercept_xy + slope_xy*x)

sd_x <- sqrt(sd(x)^2*(N-1)/N)
sd_e <- sqrt(sd(e)^2*(N-1)/N)

std_x <- (x-mean(x))/sd_x
std_e <- (e-mean(e))/sd_e

u <- R*std_x + sqrt(1-R^2)*std_e
sd_u <- sqrt(sd(u)^2*(N-1)/N)

## サンプリング
n <- 50 # サンプルサイズの指定
times <- 100000 # 抽出回数
corrs <- vector(length=times) # 標本相関係数を格納するヴェクトル

for (i in 1:times) {
     j <- sample(N, n)
     corrs[i] <- cor(std_x[j], u[j], method="pearson")
}

hist(corrs, xlim=c(-1,1), col="#99CC0080", border="#99CC00", 
     main=sprintf("標本相関係数%d個のヒストグラム\n（母相関係数ρ=%.2f, 標本相関係数の平均%.3f）", times, R, mean(corrs)),
     sub=sprintf("母集団サイズ%d, 標本サイズ%.d", N, n), cex.main=1.0,
     xlab="標本相関係数rの値", ylab="標本の個数", family="serif")
axis(side=1, R, font=2, family="serif")
text(c(min(corrs), max(corrs)), 1000, round(c(min(corrs), max(corrs)),3), family="serif")

1-1 積率相関係数のt検定(t-test)

積率相関係数の検定統計量 \[\frac{r}{\sqrt{1-r^2}}\sqrt{n-2}\;\sim\;t_{(n-2)}\]

　標本サイズと標本相関係数から，ゼロ仮説「母相関係数=0」のt検定を行う（有意確率を表示する）スクリプトを示す。片側検定と両側検定の両方を表示しているが，両側検定での有意確率は片側検定での有意確率の2倍である。よって，両側検定で有意になる場合，同じ有意水準の片側検定では必ず有意になるが，片側検定で有意になったとしても，同じ有意水準の両側検定で有意になるとは限らない。
nやrに自由な値を設定して試してみよう。rの絶対値はもちろん1を超えてはならない。

n <- 50 # 標本サイズ
r <- -.30 # 標本相関係数

t0 <- r/sqrt(1-r^2)*sqrt(n-2); t0 # 検定統計量t

[1] -2.178819

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

[1] 0.03428618

# 片側検定の有意確率
# 標本相関係数の符号と同じ側に棄却域を設定する
if (r < 0) prob <- pt(t0, n-2)
if (r >=0) prob <- 1-pt(t0, n-2)
prob

[1] 0.01714309

1-2 標本サイズと有意確率

　本文の表のように，n，rの幾つかの組合せに対して，95%限界値，有意確率を計算してみよう。本文よりも，限界値に対応する標本相関係数の値など少し追加してある。

r <- .20 # 標本相関係数
n <- 50 # サンプルサイズ
p0 <- .95 # 信頼水準（信頼係数）
t0 <- r/sqrt(1-r^2)*sqrt(n-2)
t0 # t検定統計量

[1] 1.414214

qt((1-p0)/2, n-2); qt((1+p0)/2, n-2) # 100p%限界値

[1] -2.010635

[1] 2.010635

pt(-1*t0, n-2)+(1-pt(t0, n-2)) # 両側検定での有意確率

[1] 0.1637531

# t統計量の100p%限界値tcに対応する標本相関係数rcの値（正の値のみ）
tc <- qt((1+p0)/2, n-2)
sqrt((tc^2/(n-2))/(1+tc^2/(n-2)))

[1] 0.2787106

　上の結果を解説しよう。n=50，r=.20の場合，検定統計量tは約1.414となり，95%範囲の両端の±2.010635よりずっと中心寄りであり，両端の限界値の外側に位置する5%棄却域には入らない。即ち，「母相関係数はゼロである」とのゼロ仮説を棄却出来ない。両側有意確率は約16.4%である。このサンプルサイズで5%棄却域に入る（＝ゼロ仮説が棄却できる）のは，r=.2787くらいになってからである。

1-3 母相関係数の区間推定

　サンプルサイズと標本相関係数が与えられた時に，信頼区間を計算する方法は，本文脚注にも書いてある通りである。フィッシャーのZ変換を行ったものはそのままZと表記する事が多いが[南風原 2002: 116]，ここでは他の多くの標準正規変量と区別し易いようにあえてLと表記しておく。

\[L=\frac{1}{2}ln\frac{1+r}{1-r}\sim\;N(\mu_{L}, \sigma_{L}^2)=N\left(\frac{1}{2}ln\frac{1+\rho}{1-\rho}, \left(\frac{1}{\sqrt{n-3}}\right)^2\right)\]　 \[\frac{L-\mu_{L}}{\sigma_{L}}\sim\;N(0, 1^2)\] 95%信頼区間は， \[L-1.96\sigma_{L}<\mu_{L}<L+1.96\sigma_{L}\] 母相関係数\(\rho\)の信頼区間は，下記の関係式（逆関数）から計算する。
\[\rho=\frac{e^{2\mu_{L}}-1}{e^{2\mu_{L}}+1}\]

　

r <- .30
n <- 50
lwr <- 1/2*log((1+r)/(1-r))-1.96/sqrt(n-3) # Lの信頼区間下限
upr <- 1/2*log((1+r)/(1-r))+1.96/sqrt(n-3) # Lの信頼区間上限
lwr; upr

[1] 0.02362422

[1] 0.595415

# 逆変換して相関係数の値に戻す
rhoL <- (exp(2*lwr)-1)/(exp(2*lwr)+1)
rhoU <- (exp(2*upr)-1)/(exp(2*upr)+1)
rhoL; rhoU

[1] 0.02361983

[1] 0.5337789

ここで，フィッシャーのZ変換を行ったものが実際に正規分布で近似できるかをシミュレイションで確かめてみよう。
上の1で使用した模擬母集団とそこからのランダムサンプリングの結果を利用しよう。

hist(corrs, xlim=c(-1,1),  col="#99CC0080", border="#99CC00", 
     main=sprintf("標本相関係数%d個のヒストグラム\n（母相関係数ρ=%.2f, 標本相関係数の平均%.3f）", times, R, mean(corrs)),
     sub=sprintf("母集団サイズ%d, 標本サイズ%.d", N, n), cex.main=1.0,
     xlab="標本相関係数rの値", ylab="標本の個数", family="serif")
axis(side=1, R, font=2, family="serif")
text(c(min(corrs), max(corrs)), 1000, round(c(min(corrs), max(corrs)),3), family="serif")

L <- 1/2*log((1+corrs)/(1-corrs))-1.96/sqrt(n-3) # Lの信頼区間下限
muL <- 1/2*log((1+R)/(1-R))-1.96/sqrt(n-3) # Lの期待値

hist(L, xlim=c(-1,1), breaks=seq(-1, 1, by=.05),  col="#ffcc3399", border="#ffcc33", 
     main=sprintf("フィッシャーZ変換後のL：%d個のヒストグラム\n（Lの期待値=%.2f, Lの平均%.3f）", times, muL, mean(L)),
     sub=sprintf("母集団サイズ%d, 標本サイズ%.d", N, n), cex.main=1.0,
     xlab="Lの値", ylab="標本の個数", family="serif")
text(c(min(L), max(L)), 1000, round(c(min(L), max(L)),3), family="serif")

確かに，かなり正規分布に近付いている様子が分かる。
　本文で「t検定の結果と95%信頼区間，標本相関係数の値は，Rのcor.test(x, y, method=”pearson”)によってすべて出力される。」と書いたので，その出力例を以下に挙げておく。

data01 <- read.csv("practice.csv") # working directoryにデータファイルを置いておく
# names(data01)

cor.test(data01$age, data01$fincome, method="pearson", conf=.95)

    Pearson's product-moment correlation

data:  data01$age and data01$fincome
t = 1.9832, df = 295, p-value = 0.04827
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.0009066566 0.2255742479
sample estimates:
     cor 
0.114707 

2 分割表の独立性についてのカイ二乗検定

分割表の例

カイ二乗統計量 \[\chi^2=\sum_{i=1}^{I}\sum_{j=1}^{J}\left(\frac{n_{ij}-f_{ij}}{\sqrt{f_{ij}}}\right)^2\;\sim\;\chi^2_{df=(I-1)(J-1)}\]

2-1 カイ二乗分布(chi-squared distribution)

　本文の，「自由度1,4,6,8,9,12のカイ二乗分布」のグラフを描くスクリプトは以下である。

x <- seq(0, end <- 20, by=0.1)

df0 <- 1
# 最初に自由度1で基本のグラフを描く
plot(x, dchisq(x, df=df0), type="l", ylim=c(0, .20),
     bty="n", xlab="カイ二乗変量", ylab="確率密度", las=1,
     main="自由度1,4,6,8,9,12のカイ二乗分布",
     family="serif")

for (df1 in c(4, 6, 8, 9, 12)) {
par(new=T) # グラフを重ね描きする為のコマンド
plot(x, dchisq(x, df=df1), type="l", ylim=c(0, .20),
     bty="n", axes=F, ylab="", xlab="", lwd=1+df1/5,
     lty=floor(df1/2))
}

　次の，棄却域を塗り潰したカイ二乗分布のグラフは下記。自由度と有意水準は変更する事が出来る。

df0 <- 6 # 自由度を指定する
p0 <- .95 # 有意水準を指定する
x <- seq(0, end <- 20, by=0.1) # endはx軸の右端
plot(x, dchisq(x, df=df0), type="l", ylim=c(0, .20),
     bty="n", xlab="カイ二乗変量", ylab="確率密度", las=1,
     main=paste("自由度", df0, "のカイ二乗分布の", (1-p0)*100, "%棄却域"),
     family="serif", cex.main=1)
abline(h=0, col="gray")
rej <- seq(qchisq(p0, df0), end, by=.1) # 棄却域を設定
segments(rej, 0, rej, dchisq(rej, df0), col="#CC330080")
segments(min(rej), 0, max(rej), 0, lwd=3, col="#CC3300")
axis(side=1, at=round(min(rej),1), las=2, font=2, family="serif")

　実際のデータ（演習用データ）から分割表のカイ二乗検定を行うには下記の様にする。欠損値処理が入っているので少し手間がかかっている。ついでにクラメールのVやファイ係数の計算もしておこう。

data01 <- read.csv("practice.csv") # working directoryにデータファイルを置いておく
names(data01)

 [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"  

data01$q1601b <- data01$q1601 # q1601をq1601bに複製して，
data01$q1601b[data01$q1601 ==6] <- NA # 欠損値処理をする

table1 <- table("性別"=data01$sex, "機会平等なら格差やむをえない"=data01$q1601b)
dimnames(table1)[[1]] <- c("男性", "女性")
dimnames(table1)[[2]] <- c("そう思う", "どち…思う", "どちらとも", "どち…思わない", "そう思わない")
table1 # 分割表

      機会平等なら格差やむをえない
性別   そう思う どち…思う どちらとも どち…思わない そう思わない
  男性       53         57         32             12            8
  女性       35         81         76             13           12

ch <- chisq.test(table1); ch # カイ二乗検定の結果

    Pearson's Chi-squared test

data:  table1
X-squared = 19.041, df = 4, p-value = 0.0007715

names(ch)

[1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed" 
[7] "expected"  "residuals" "stdres"   

ch$residual # 標準化残差

      機会平等なら格差やむをえない
性別     そう思う どち…思う どちらとも どち…思わない そう思わない
  男性  2.5085606 -0.2586892 -2.0846038      0.4019594   -0.1877030
  女性 -2.1674657  0.2235146  1.8011553     -0.3473041    0.1621806

ch$stdres # 調整済み標準化残差

      機会平等なら格差やむをえない
性別     そう思う どち…思う どちらとも どち…思わない そう思わない
  男性  3.7834451 -0.4287253 -3.2579785      0.5496547   -0.2548786
  女性 -3.7834451  0.4287253  3.2579785     -0.5496547    0.2548786

V <- sqrt(ch$statistic/min(dim(table1)-1)/sum(table1))
phi <- sqrt(ch$statistic/sum(table1))
names(V) <- "Cramer's V"; names(phi) <- "phi coefficient"
V; phi

Cramer's V 
 0.2241438 

phi coefficient 
      0.2241438 

自由dによって変化するカイ二乗分布

カイ二乗分布の定義： \[\chi_{(df=\nu)}^2=z_{1}^2+z_{2}^2+...+z_{\nu}^2=\left(\frac{x_{1}-\bar{x}_{1}}{s_{x_{1}}}\right)^2+\left(\frac{x_{2}-\bar{x}_{2}}{s_{x_{2}}}\right)^2+...+\left(\frac{x_{\nu}-\bar{x}_{\nu}}{s_{x_{\nu}}}\right)^2\]

2-2 カイ二乗検定の前提条件

ポワソン分布（計数分布）\(x_{ij}\;\sim\;Poisson(\lambda=f_{ij}),\;mean=\lambda,\;s.d.=\sqrt{\lambda}\)
\(\lambda\geq5\)ならば，\(x_{ij}\;\sim\;Poisson(\lambda=f_{ij})\;\sim\;N(\lambda,\lambda)\)
標準化すると， \[\frac{x_{ij}-f_{ij}}{\sqrt{f_{ij}}}\;\sim\;N(0.1^2)\]

　分割表の独立性についてのカイ二乗検定において，検定統計量がカイ二乗分布で近似できる条件は，各セルの度数がポワソン分布に従い，そのポワソン分布が正規分布で近似できるということである。本文中の，二つのポワソン分布を重ね描きするスクリプトは下記である。

x1 <- c(0:15)
lambda1 <- 2
lambda2 <- 5
plot(x1, dpois(x1, lambda1), type="h", bty="n",
     lwd=5, col="#00000060", xlab="", ylab="", 
     las=1, family="serif", cex.main=1.0, ylim=c(0, .3),
     main=sprintf("平均%dと%dのポワソン分布", lambda1, lambda2))
par(new=T)
plot(x1, dpois(x1, lambda2), type="b", bty="n", axes=F, 
     lwd=1, col="#000000", xlab="", ylab="",
     las=1, family="serif", ylim=c(0, .3))

　少し改変したスクリプトとその結果も例示しておこう。

lambda <- c(3, 6) # 二つのポワソン分布の平均をヴェクトルで指定
x1 <- c(0: (max(lambda)*2))
plot(x1-.1, dpois(x1, min(lambda)), type="h", bty="n", xlim=range(x1),
     lwd=5, col="#0000ff60", xlab="", ylab="", las=1, family="serif", 
     cex.main=1.0, ylim=c(0, Y <- dpois(min(lambda), min(lambda))),
     main=sprintf("平均%dと%dのポワソン分布", min(lambda), max(lambda)))
par(new=T)
plot(x1+.1, dpois(x1, max(lambda)), type="h", bty="n", axes=F, 
     lwd=5, col="#ff000060", xlab="", ylab="", xlim=range(x1), 
     las=1, family="serif", ylim=c(0, Y))

2-3 データに対するモデルの適合度検定(Goodness of Fit test)

　社会学では，「カイ二乗検定」と云うと「独立性の検定」と同一視する傾向があるが，これはデータに対して「2変数独立モデル」が適合するかどうかを検定していると考えれば，「適合度検定」の一種であるとも言える。カイ二乗統計量は社会統計学のより発展的な場面で色々と登場してくる事になるが，いずれも「データとモデルのズレ」を表す統計量として統一的に考えるのが良い。ここでは，chisq.test( )関数を使用した最も単純な適合度の検定を紹介しよう。
　少人数の無作為標本に対して，「夫は外で働き，妻は家庭を守るべきである」という考え方について，そう思うかそう思わないか（性別役割分業意識）を，「そう思う／どちらかと言えばそう思う／どちらかと言えばそう思わない／そう思わない」の4件法で尋ねた結果，それぞれ5人，8人，15人，12人であったとする。
　「この結果では偶然ばらついたが，本当はどの回答も等確率で出現する」，つまり同じ比率ずつ分布していると仮定した時，その等比率（もしくは等確率）モデルと実際のデータのズレが偶然の範囲内と言えるかどうかをカイ二乗検定する事が出来るが，これが最も単純な適合度検定である。

# 適合度のカイ二乗検定1： 等比率（等確率）モデルがデータに適合するか
x <- c(5, 8, 15, 12) # 回答の結果のヴェクトル
sum(x) # 一応総人数を確認しておく

[1] 40

ch01 <- chisq.test(x); ch01 # 適合度のカイ二乗検定

    Chi-squared test for given probabilities

data:  x
X-squared = 5.8, df = 3, p-value = 0.1218

ch01$expected # 等比率だとした場合の期待度数

[1] 10 10 10 10

　等確率だとしたら全ての選択肢が10人ずつに選ばれる筈であり，実際の結果はそこからカイ二乗統計量にして5.8ずれているが，これは自由度3(=カテゴリ数4 - 1)の\(χ^2\)分布に近似的に従うと考えられる。自由度3の\(χ^2\)分布で5.8以上の範囲の面積（確率）は約12.2%ある。つまり本当は（母集団では）等比率で分布しているのにたまたまこの様な無作為標本が得られることは10回に1回より多く起こると考えられ，特に不自然であるとは言えないので，「本当は（母集団では）等比率である」と云うゼロ仮説は棄却出来ない。
　これを，等比率モデルの適合度は悪くないと考えるのが適合度検定である。「モデルがデータに適合している」と云うゼロ仮説が棄却されなかった事をもって「モデルは正しい」と解釈する事は統計の誤用であるとする立場からするとこの適合度検定の考え方には疑問や批判もある。その立場からすると，「等比率モデルは間違っていると自信をもって言えるだけの証拠が無かっただけで，等比率モデルが正しいと言える訳ではなく，精々判断が保留されるだけである」となる。この辺りは研究者によって考え方が違ってくるかも知れない。
　上の例では，適合度検定と独立性の検定の違いがいまいちはっきりしないかも知れないので，少し違った例を示そう。
　上では等比率仮説は棄却されないが，流石に現在の日本ではこの質問に対して4つの選択肢で同じ比率ずつ意見が分布しているとは考えられない。本当の意見分布を2:3:4:5であると仮定して，データと齟齬を来すかどうかを検討しよう。Rスクリプトと結果は以下の通りである。chisq.test( )関数のオプションp=で，それぞれのカテゴリの生起確率を指定する点だけが異なる。

ch02 <- chisq.test(x, p=c(2, 3, 4, 5)/14)
ch02

    Chi-squared test for given probabilities

data:  x
X-squared = 1.6092, df = 3, p-value = 0.6573

ch02$expected

[1]  5.714286  8.571429 11.428571 14.285714

　結果，性別役割分業に否定的な人ほど数が多いと云うこのモデルも有意性検定で棄却されず，データとの適合は悪くないと云う結果になった。むしろズレの指標は等比率モデルよりも遥かに小さく，有意確率もずっと高い。いずれのモデルも（そして他にも多くの可能なモデルが）データに適合する（適合しないとは言い切れない）と云う結果になる。そのいずれかが正しいと証明される訳でも無く，カイ二乗統計量や有意確率の値からいずれかが「より正しい」と言えるかどうかにも議論がある。
　適合度検定と考えられるものには，第6章の基礎1-1や第7章基礎1-5の等分散性の検定（各群の母分散が等しいモデルが適合するかどうか），第6章発展1の正規性の検定（変数の分布が正規分布に適合するかどうか）がある。それ以外にも，第11章発展1-1の構造方程式モデリング(SEM; Structural Equation Modelling)では，設定したモデルの適合度がカイ二乗検定される（それ以外に，GFI[Goodness-of-Fit Index]を始めとして幾つかの適合度指標が算出されるが，GFIらについては統計的検定がされる訳では無い）。第12章基礎1-4では，2項ロジスティック回帰分析に対する適合度検定としてHosmer-Lemeshow検定を紹介した。第12章発展1のログリニア・モデル(log-linear model)も「モデルがデータに適合するか」を重視した分析手法であると言える。
　適合度検定の多くは，統計的検定一般がそうである様に，サンプルサイズが大きいと有意になり易い。大標本では「僅かなズレでも敏感に検出してしまう」のである。適合度検定で統計的に有意になると云う事は，「モデルがデータに適合している」と仮定すると（ゼロ仮説），ズレの指標が偶然とは思えない程に大きくなるので，ゼロ仮説は棄却されると云う事であり，つまり「このモデルはデータには適合しない」と云う結論になるので，モデルが当てはまる事を期待していた分析者にとっては嬉しい結果ではない事も多いだろう。

発展1 分割表のさらなる検定

発展1-1　分割表の残差分析(analysis of residual)

　上の2-1でデータフレイムdata01，その中の変数data01$sex，data01$q1601b，分割表table1などを既に準備済みであるので，それを利用して調整済み標準化残差（ハーバーマン残差）を表示したりモザイクプロットを表示したりしてみよう。

table1 # 分割表

      機会平等なら格差やむをえない
性別   そう思う どち…思う どちらとも どち…思わない そう思わない
  男性       53         57         32             12            8
  女性       35         81         76             13           12

ch <- chisq.test(table1); ch # カイ二乗検定の結果

    Pearson's Chi-squared test

data:  table1
X-squared = 19.041, df = 4, p-value = 0.0007715

names(ch)

[1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed" 
[7] "expected"  "residuals" "stdres"   

ch$"stdres"

      機会平等なら格差やむをえない
性別     そう思う どち…思う どちらとも どち…思わない そう思わない
  男性  3.7834451 -0.4287253 -3.2579785      0.5496547   -0.2548786
  女性 -3.7834451  0.4287253  3.2579785     -0.5496547    0.2548786

round(ch$"stdres", 3) # 小数第三位までに四捨五入して表示

      機会平等なら格差やむをえない
性別   そう思う どち…思う どちらとも どち…思わない そう思わない
  男性    3.783     -0.429     -3.258          0.550       -0.255
  女性   -3.783      0.429      3.258         -0.550        0.255

mosaicplot(table1, type="pearson", shade=T,
     main="モザイクプロット", las=2)

発展1-2　イェーツの連続性修正(Yates’ correction for continuity)とフィッシャーの正確検定(Fisher’s exact test)

　模擬的に2行2列の分割表を作成し，カイ二乗検定の関数を適用してみよう。correct=というオプションを省略した場合，TRUEにした場合，FALSEにした場合をよく比較して欲しい。2行2列分割表ではcorrect=を省略すると自動的にcorrect=TRUEとなる。行数か列数が2を超えると連続修正は存在しない。

t1 <- matrix(c(10, 30, 20, 40), ncol=2); t1

     [,1] [,2]
[1,]   10   20
[2,]   30   40

chisq.test(t1) # デフォルトでのカイ二乗検定

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

data:  t1
X-squared = 0.44643, df = 1, p-value = 0.504

chisq.test(t1, correct=T) # 連続性修正のオプションをTRUE

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

data:  t1
X-squared = 0.44643, df = 1, p-value = 0.504

chisq.test(t1, correct=F) # 連続性修正のオプションをFALSE

    Pearson's Chi-squared test

data:  t1
X-squared = 0.79365, df = 1, p-value = 0.373

同じ分割表に対してフィッシャーの正確検定を行ってみる。fisher.test( )関数の出力するodds ratioは注意が必要なので，通常のオッズ比を計算するコマンドも最後に付け加えてある。

t1

     [,1] [,2]
[1,]   10   20
[2,]   30   40

fisher.test(t1)

    Fisher's Exact Test for Count Data

data:  t1
p-value = 0.5045
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.2420054 1.7643388
sample estimates:
odds ratio 
 0.6693434 

(odds.ratio <- t1[1,1]*t1[2,2]/t1[1,2]/t1[2,1]) # 素朴なオッズ比

[1] 0.6666667

　演習用データの学歴（3行）×性別役割分業意識（4列）の分割表にカイ二乗検定や正確検定を行ってみよう。2行2列ではないのでオッズ比は計算されない。

table1 <- table("学歴3区分"=data01$edu2, "性別役割分業"=data01$q1900)
dimnames(table1)[[1]] <- c("1 初等", "2 中等", "3 高等")
dimnames(table1)[[2]] <- c("1 そう思う", "2 どち…思う", "3 どち…思わない", "4 思わない")
table1 # まず分割表（特に行変数と列変数）を確認する。

         性別役割分業
学歴3区分 1 そう思う 2 どち…思う 3 どち…思わない 4 思わない
   1 初等          6           17               26         25
   2 中等          3           15               21         22
   3 高等          3           19               15         31

chisq.test(table1) # カイ二乗検定

Warning in chisq.test(table1): Chi-squared approximation may be incorrect

    Pearson's Chi-squared test

data:  table1
X-squared = 5.1487, df = 6, p-value = 0.5249

fisher.test(table1) # 正確検定

    Fisher's Exact Test for Count Data

data:  table1
p-value = 0.5283
alternative hypothesis: two.sided

fisher.test(table1, simulate.p.value=T) # 正確検定（近似計算）

    Fisher's Exact Test for Count Data with simulated p-value (based on
    2000 replicates)

data:  table1
p-value = 0.5322
alternative hypothesis: two.sided

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

1) の答え

r <- .25 # 標本相関係数
n <-　 30 # サンプルサイズ
p0 <- .95 # 信頼水準（信頼係数）
t0 <- r/sqrt(1-r^2)*sqrt(n-2)
t0 # t検定統計量

[1] 1.36626

qt((1-p0)/2, n-2); qt((1+p0)/2, n-2) # 100p%限界値

[1] -2.048407

[1] 2.048407

pt(-1*t0, n-2)+(1-pt(t0, n-2)) # 両側検定での有意確率

[1] 0.1827305

　有意確率は18.3%程度であり，有意水準を1%にしようが5%にしようが10%にしようが有意にはならない。つまりゼロ仮説は棄却出来ない。片側検定として考えるなら有意水準を半分にすればよいが，それでも5%では有意にならない。

2) の答え

r <- .25
n <- 70
p0 <- .95
t0 <- r/sqrt(1-r^2)*sqrt(n-2)
t0

[1] 2.129163

qt((1-p0)/2, n-2); qt((1+p0)/2, n-2)

[1] -1.995469

[1] 1.995469

pt(-1*t0, n-2) + (1-pt(t0, n-2))

[1] 0.03686488

　両側5%検定で有意になり，ゼロ仮説は棄却される。1)とは標本相関係数の値は同じであるが，標本サイズが2倍以上になっており，弱い関連に対してより敏感に検出出来る様になっている。

3) の答え

t01 <- matrix(c(49,47,84,37,37,17,17,14,6,25,25,17), ncol=4)
t01

     [,1] [,2] [,3] [,4]
[1,]   49   37   17   25
[2,]   47   37   14   25
[3,]   84   17    6   17

ch1 <- chisq.test(t01); ch1

    Pearson's Chi-squared test

data:  t01
X-squared = 30.377, df = 6, p-value = 3.333e-05

高度に有意となり，母集団では独立であるとのゼロ仮説は棄却される。
独立性が否定されたので，行比率や列比率を計算して，どの様な偏りがあるのかを確認しておく。

addmargins(prop.table(t01, margin=1))[1:3,]

                                         Sum
 0.3828125 0.2890625 0.1328125 0.1953125   1
 0.3821138 0.3008130 0.1138211 0.2032520   1
 0.6774194 0.1370968 0.0483871 0.1370968   1

addmargins(prop.table(t01, margin=2))[,1:4]

                                           
    0.2722222 0.4065934 0.4594595 0.3731343
    0.2611111 0.4065934 0.3783784 0.3731343
    0.4666667 0.1868132 0.1621622 0.2537313
Sum 1.0000000 1.0000000 1.0000000 1.0000000

4) の答え

　カイ二乗検定のみならず正確検定も行っておこう。

t01s <- round(t01/3, 0); t01s

     [,1] [,2] [,3] [,4]
[1,]   16   12    6    8
[2,]   16   12    5    8
[3,]   28    6    2    6

chisq.test(t01s)

Warning in chisq.test(t01s): Chi-squared approximation may be incorrect

    Pearson's Chi-squared test

data:  t01s
X-squared = 9.5046, df = 6, p-value = 0.1471

fisher.test(t01s)

    Fisher's Exact Test for Count Data

data:  t01s
p-value = 0.1437
alternative hypothesis: two.sided

addmargins(prop.table(t01s, margin=1))[1:3,]

                                          Sum
 0.3809524 0.2857143 0.14285714 0.1904762   1
 0.3902439 0.2926829 0.12195122 0.1951220   1
 0.6666667 0.1428571 0.04761905 0.1428571   1

いずれにしても有意にはならず，独立モデルは棄却出来ない。3)と比較すると，分割表の行数・列数は等しく，したがって自由度も同じである。各セルの度数を一律に三分の一にしているので，行比率の表もほぼ同じである。
　しかしカイ二乗統計量の値は，30.377と9.5046と3倍余り異なり，その結果検定結果は1%有意と10%でも有意でないと云う様に対照的である。サンプルサイズが小さいと，同じ様な関連が「偶然に過ぎないかも知れない」と見做されてしまうのである。

5) の答え

　調整済み標準化残差（ハーバーマン残差）が標準正規分布に近似的に従う事を利用する。
調整済み標準化残差は，カイ二乗検定の結果を格納したオブジェクトのstdresで取り出せる。

ch1$stdres

          [,1]      [,2]       [,3]       [,4]
[1,] -2.711816  1.508699  1.5961802  0.6057539
[2,] -2.650739  1.835024  0.6875128  0.8682866
[3,]  5.378449 -3.351823 -2.2948443 -1.4770222

有意水準を両側5%とすると，絶対値が1.96を超えているセルが，有意に独立モデルから乖離しているセルであり，[3行, 1列]が期待値より有意に，かなり多い。[3, 2], [3, 3], [1, 1], [1, 2]は期待値より有意に少ない。つまり，3行1列のセルに著しく集中している事が分かる。
　実はこの分割表は，演習用データの学歴3区分と従業上の地位4区分のものであり，高学歴者ほど常時雇用者が多いと云う結果を示しているのである。

addmargins(t01)

                 Sum
     49 37 17 25 128
     47 37 14 25 123
     84 17  6 17 124
Sum 180 91 37 67 375

tbl <- table("学歴3区分"=data01$edu2, "従業上の地位"=data01$job)
dimnames(tbl)[[1]] <- c("1 初等", "2 中等", "3 高等")
dimnames(tbl)[[2]] <- c("1 常時雇用", "2 臨時雇用", "3 自営", "4 無職")
addmargins(tbl)

         従業上の地位
学歴3区分 1 常時雇用 2 臨時雇用 3 自営 4 無職 Sum
   1 初等         49         37     17     25 128
   2 中等         47         37     14     25 123
   3 高等         84         17      6     17 124
   Sum           180         91     37     67 375

round(prop.table(tbl, margin=1)*100, 1)

         従業上の地位
学歴3区分 1 常時雇用 2 臨時雇用 3 自営 4 無職
   1 初等       38.3       28.9   13.3   19.5
   2 中等       38.2       30.1   11.4   20.3
   3 高等       67.7       13.7    4.8   13.7

ウェブ増補1　xtabs( )関数による分割表の作成

　第2章と第5章を中心として，本書では分割表を作成する為に table( )関数を使用しているが，他によく使用される関数として，xtabs( )があるのでここで補っておく。

分割表の基本

　まずは7ケース分のデータ・フレイムdata01を作成してから，table( )とxtabs( )で三元分割表を作成してみる。その後，ftable( )を利用して三元分割表を見易くする。
　コマンドの文法と出力スタイルの両方を見比べてみよう。

data01 <- data.frame("sex"=c( 1, 1, 1, 1, 2, 2, 2),
    "edu"=c(12,14,16,16,12,14,16),
    "outcome"=c( 3, 4, 4, 5, 5, 4, 3))
data01 # 7ケース分のデータ・フレイム, NA無し

table(edu, outcome, sex, data=data01) # table( )ではこの書き方はNG

Error in table(edu, outcome, sex, data = data01):  オブジェクト 'edu' がありません 

table(data01$edu, data01$outcome, data01$sex)

, ,  = 1

    
     3 4 5
  12 1 0 0
  14 0 1 0
  16 0 1 1

, ,  = 2

    
     3 4 5
  12 0 0 1
  14 0 1 0
  16 1 0 0

# xtabs( )での書き方. table( )よりも若干親切な出力表示と，lm( )に似た構文.
xtabs(~ edu + outcome + sex, data=data01)

, , sex = 1

    outcome
edu  3 4 5
  12 1 0 0
  14 0 1 0
  16 0 1 1

, , sex = 2

    outcome
edu  3 4 5
  12 0 0 1
  14 0 1 0
  16 1 0 0

# 三元分割表をftableで表示するには，変数の順序に注意
ftable(table(data01$edu, data01$outcome, data01$sex))

      1 2
         
12 3  1 0
   4  0 0
   5  0 1
14 3  0 0
   4  1 1
   5  0 0
16 3  0 1
   4  1 0
   5  1 0

ftable(table(data01$sex, data01$edu, data01$outcome))

      3 4 5
           
1 12  1 0 0
  14  0 1 0
  16  0 1 1
2 12  0 0 1
  14  0 1 0
  16  1 0 0

ftable(xtabs(~ edu + outcome + sex, data=data01))

            sex 1 2
edu outcome        
12  3           1 0
    4           0 0
    5           0 1
14  3           0 0
    4           1 1
    5           0 0
16  3           0 1
    4           1 0
    5           1 0

　この限りでは，特に重回帰分析（一般線型モデル）などに慣れている人ほど，xtabs( )の方を好むだろう。出力の表記だけをとってもxtabs( )の方がやや親切である。
　但しtable( )の場合は，以下の様に with( ) と組合せる事も出来る。

with(data01, table(edu, outcome, sex))

, , sex = 1

    outcome
edu  3 4 5
  12 1 0 0
  14 0 1 0
  16 0 1 1

, , sex = 2

    outcome
edu  3 4 5
  12 0 0 1
  14 0 1 0
  16 1 0 0

NAの処理

　実際の調査データ分析では，欠損値（欠測missing）の処理が必須である。この点も違いがあるので比べておこう。

data02 <- data.frame("sex"=c( 1, 1, 1,NA, 2, 2, 2),
    "edu"=c(12,14,NA,16,12,14,16),
    "outcome"=c( 3, 4, 4, 5, 5,NA, 3))
data02 # 7ケース分のデータ・フレイム, NA有り

table(data02$edu, data02$outcome, data02$sex)

, ,  = 1

    
     3 4 5
  12 1 0 0
  14 0 1 0
  16 0 0 0

, ,  = 2

    
     3 4 5
  12 0 0 1
  14 0 0 0
  16 1 0 0

xtabs(~ edu + outcome + sex, data=data02)

, , sex = 1

    outcome
edu  3 4 5
  12 1 0 0
  14 0 1 0
  16 0 0 0

, , sex = 2

    outcome
edu  3 4 5
  12 0 0 1
  14 0 0 0
  16 1 0 0

table(data02$edu, data02$outcome, data02$sex, useNA="ifany")

, ,  = 1

      
       3 4 5 <NA>
  12   1 0 0    0
  14   0 1 0    0
  16   0 0 0    0
  <NA> 0 1 0    0

, ,  = 2

      
       3 4 5 <NA>
  12   0 0 1    0
  14   0 0 0    1
  16   1 0 0    0
  <NA> 0 0 0    0

, ,  = NA

      
       3 4 5 <NA>
  12   0 0 0    0
  14   0 0 0    0
  16   0 0 1    0
  <NA> 0 0 0    0

xtabs(~ edu + outcome + sex, data=data02, na.action=na.pass, exclude=NULL)

, , sex = 1

      outcome
edu    3 4 5 <NA>
  12   1 0 0    0
  14   0 1 0    0
  16   0 0 0    0
  <NA> 0 1 0    0

, , sex = 2

      outcome
edu    3 4 5 <NA>
  12   0 0 1    0
  14   0 0 0    1
  16   1 0 0    0
  <NA> 0 0 0    0

, , sex = NA

      outcome
edu    3 4 5 <NA>
  12   0 0 0    0
  14   0 0 0    0
  16   0 0 1    0
  <NA> 0 0 0    0

ftable(table(data02$edu, data02$outcome, data02$sex, useNA="ifany"))

       1 2 NA
             
12 3   1 0  0
   4   0 0  0
   5   0 1  0
   NA  0 0  0
14 3   0 0  0
   4   1 0  0
   5   0 0  0
   NA  0 1  0
16 3   0 1  0
   4   0 0  0
   5   0 0  1
   NA  0 0  0
NA 3   0 0  0
   4   1 0  0
   5   0 0  0
   NA  0 0  0

ftable(xtabs(~ edu + outcome + sex, data=data02, na.action=na.pass, exclude=NULL))

            sex 1 2 NA
edu outcome           
12  3           1 0  0
    4           0 0  0
    5           0 1  0
    NA          0 0  0
14  3           0 0  0
    4           1 0  0
    5           0 0  0
    NA          0 1  0
16  3           0 1  0
    4           0 0  0
    5           0 0  1
    NA          0 0  0
NA  3           0 0  0
    4           1 0  0
    5           0 0  0
    NA          0 0  0

　NAの処理は，xtabs( )の方が何だか面倒である。

データ・フレイムへの（逆）変換

　分割表からデータ・フレイムへの変換（逆変換？）を行う関数，as.data.frame.table( )があるが，これを使用する場合には，table( )の結果よりもxtabs( )の結果の方が便利である。上のdata01を使って実際に示してみるので，どう云う事か違いをよく見比べて読み取って欲しい。

t1 <- table(data01$edu, data01$outcome, data01$sex); t1

, ,  = 1

    
     3 4 5
  12 1 0 0
  14 0 1 0
  16 0 1 1

, ,  = 2

    
     3 4 5
  12 0 0 1
  14 0 1 0
  16 1 0 0

xt1 <- xtabs(~ edu + outcome + sex, data=data01); xt1

, , sex = 1

    outcome
edu  3 4 5
  12 1 0 0
  14 0 1 0
  16 0 1 1

, , sex = 2

    outcome
edu  3 4 5
  12 0 0 1
  14 0 1 0
  16 1 0 0

as.data.frame.table(t1) # データフレイムに逆変換

as.data.frame.table(xt1) # データフレイムに逆変換

ウェイト付け

　上のas.data.frame.table( )関数の結果をよく見ると，出力されたデータ・フレイムにはFreqと云う変数が新たに加わっている。この例では，それぞれのパタンに該当するのはせいぜい1ケースしかなかったのでFreqは0か1の値しかないが，全く同じ変数値プロファイルを持つケースが複数あれば，Freqもその2以上のケース数を表す。
　逆に，その様なFreq変数（変数名は別でも良い）を持つデータ・フレイムから，その度数によって重み付けを行った分割表を作成する事も xtabs( )によって可能になる。
　

# weightを付ける事が出来る。
data01$Freq <- c(1,2,3,4,5,6,7)
data01

xtabs(~ edu + outcome + sex, data=data01)

, , sex = 1

    outcome
edu  3 4 5
  12 1 0 0
  14 0 1 0
  16 0 1 1

, , sex = 2

    outcome
edu  3 4 5
  12 0 0 1
  14 0 1 0
  16 1 0 0

xtabs(Freq ~ edu + outcome + sex, data=data01)

, , sex = 1

    outcome
edu  3 4 5
  12 1 0 0
  14 0 2 0
  16 0 3 4

, , sex = 2

    outcome
edu  3 4 5
  12 0 0 5
  14 0 6 0
  16 7 0 0

　これを利用すると，単に度数で重みを付ける（頻度ウェイティング）だけでなく，サンプルの歪みの補正にしばしば使われる重み付け（確率ウェイティング）を行う事も出来る。原理的には，重回帰分析などのウェイト付けと同じである。ウェイト付けを行わない場合と行った場合では，そもそも分布が違ったものになるので，分析結果も異なり得る。

# 単純に度数ではなく，ケースへのweight付けを同様に行える。
data01$weight <- c(1.6, .4, 1.8 , 0.2, .4, .6, 2.0)
data01

xtabs(weight ~ edu + outcome + sex, data=data01)

, , sex = 1

    outcome
edu    3   4   5
  12 1.6 0.0 0.0
  14 0.0 0.4 0.0
  16 0.0 1.8 0.2

, , sex = 2

    outcome
edu    3   4   5
  12 0.0 0.0 0.4
  14 0.0 0.6 0.0
  16 2.0 0.0 0.0

# 線型モデル lm( )におけるweight付け
data01$fsex <- factor(data01$sex) # 性別は要因型（名義尺度）に
data01$cedu <- data01$edu - mean(data01$edu) # 教育年数は中心化
lm(outcome ~ fsex * cedu, data=data01) # weight付け無し

Call:
lm(formula = outcome ~ fsex * cedu, data = data01)

Coefficients:
(Intercept)        fsex2         cedu   fsex2:cedu  
    3.92208     -0.06494      0.36364     -0.86364  

lm(outcome ~ fsex * cedu, data=data01, weights=weight) # weight付け有り

Call:
lm(formula = outcome ~ fsex * cedu, data = data01, weights = weight)

Coefficients:
(Intercept)        fsex2         cedu   fsex2:cedu  
     3.6734       0.1838       0.2725      -0.7725