カテゴリカルデータ解析読書メモ(第5章)

シンプソンのパラドックス

層別をしないで解析すると関連が見られるが、層別をして解析すると関連が見られなくなるような現象。層別に用いるカテゴリカル変数が、他のカテゴリカル変数全てに影響を与えているような場合に発生する、見かけ上の相関。連続変数で言うところの擬似相関と同様のもの。2つの変数の関係を調べる際には両方の変数に影響する他の変数が無いかどうかに注意する必要がある。

層別2 x 2表の解析

カリフォルニア州立大学バークレー校の1973年の入試結果データであるUCBAdmissionsについて、各クラスの性別ごとの合否結果を合算して集計すると、男性の方が合格率が高いという結果が得られる。

> GAdata <- apply(UCBAdmissions, c(2, 1), sum)
> GAdata
        Admit
Gender   Admitted Rejected
  Male       1198     1493
  Female      557     1278
> chisq.test(GAdata)

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

data:  GAdata
X-squared = 91.61, df = 1, p-value < 2.2e-16

> mosaicplot(GAdata)

f:id:Rion778:20160726215732p:plain
しかし、層別にプロットすると、A学科を除いて合格者と不合格者の男女比はあまり変わらない。

old <- par(mfrow = c(2, 3), mar = c(3, 3, 3, 3), oma = c(2, 2, 2, 2))
for(i in 1:6){
  mosaicplot(~ Admit + Gender, data = UCBAdmissions[,,i],
             main = paste("Department", LETTERS[i]))
}
mtext("Student admissions at UC Barkeley", outer = TRUE)
par(old)

f:id:Rion778:20160726215922p:plain

条件付き独立

2つの変数間の関係について、3つめの変数で層別した場合に関係が見られなくなるような場合、最初の2つの変数は条件付き独立であるという。
3つめの変数で層別された2 x 2表について、条件付き独立であることを帰無仮説とする検定（マンテル・ヘンセル検定）を、Rではmantelhaen.test()で実行できる。

> mantelhaen.test(UCBAdmissions)

	Mantel-Haenszel chi-squared test with continuity correction

data:  UCBAdmissions
Mantel-Haenszel X-squared = 1.4269, df = 1, p-value = 0.2323
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.7719074 1.0603298
sample estimates:
common odds ratio 
        0.9046968

ここでcommon odds ratioとして推定されているのは共通オッズ比と呼ばれるもので、層別した場合の各層におけるオッズ比は等しいと仮定した場合のオッズ比を指している。

オッズ比の均一性の検定

それぞれの層のオッズ比が等しいかどうかを検定するために、パッケージepiDisplayのmhor()関数が使用できる*1。

> library(epiDisplay)
> mhor(mhtable = UCBAdmissions)

Stratified analysis by  Dept 
                OR lower lim. upper lim.  P value
Dept A       0.350      0.197      0.592 1.67e-05
Dept B       0.803      0.294      2.004 6.77e-01
Dept C       1.133      0.845      1.516 3.87e-01
Dept D       0.921      0.679      1.250 5.99e-01
Dept E       1.221      0.806      1.839 3.60e-01
Dept F       0.828      0.433      1.576 5.46e-01
M-H combined 0.905      0.772      1.060 2.17e-01

M-H Chi2(1) = 1.52 , P value = 0.217 
Homogeneity test, chi-squared 5 d.f. = 17.96 , P value = 0.003

f:id:Rion778:20160726221638p:plain
また、テキストに記述のある「マンテル・ヘンセル統計量を用いた別の均一性の検定法(Breslow and Day, 1980)」は、現在はDescToolsパッケージのBreslowDayTest()関数で実行できる様子。

> library(DescTools)
> BreslowDayTest(UCBAdmissions)

	Breslow-Day test on Homogeneity of Odds Ratios

data:  UCBAdmissions
X-squared = 18.826, df = 5, p-value = 0.002071

層別したオッズ比を求める

層別した後のオッズ比はvcdパッケージのoddsratio()関数で計算でき、confint()による信頼区間の計算やplot()によるグラフ出力にも対応している。
当時よりバージョンが上がっているのか、テキストに記載のものに比べると出力が見やすくなっている印象。

> library(vcd)
> lor <- oddsratio(UCBAdmissions, log = FALSE)
> lor
 odds ratios for Admit and Gender by Dept 

        A         B         C         D         E         F 
0.3492120 0.8025007 1.1330596 0.9212838 1.2216312 0.8278727 
> confint(lor)
                                    2.5 %    97.5 %
Admitted:Rejected/Male:Female|A 0.2086756 0.5843954
Admitted:Rejected/Male:Female|B 0.3403815 1.8920166
Admitted:Rejected/Male:Female|C 0.8545328 1.5023696
Admitted:Rejected/Male:Female|D 0.6863345 1.2366620
Admitted:Rejected/Male:Female|E 0.8250748 1.8087848
Admitted:Rejected/Male:Female|F 0.4552059 1.5056335

f:id:Rion778:20160726222138p:plain

層別I x J表の解析

層別I x J表においても、条件付き独立性の検定にはmantelhaen.test()が、線形関連の検定にはcoinパッケージのlbl_test()関数が同様に使用できる。

> library(coin)
> mantelhaen.test(jobsatisfaction)

	Cochran-Mantel-Haenszel test

data:  jobsatisfaction
Cochran-Mantel-Haenszel M^2 = 10.2, df = 9, p-value = 0.3345

> lbl_test(jobsatisfaction)

	Asymptotic Linear-by-Linear Association Test

data:  Job.Satisfaction (ordered) by
	 Income (<5000 < 5000-15000 < 15000-25000 < >25000) 
	 stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided