シンプソンのパラドックス
層別をしないで解析すると関連が見られるが、層別をして解析すると関連が見られなくなるような現象。層別に用いるカテゴリカル変数が、他のカテゴリカル変数全てに影響を与えているような場合に発生する、見かけ上の相関。連続変数で言うところの擬似相関と同様のもの。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)
しかし、層別にプロットすると、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)
条件付き独立
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
また、テキストに記述のある「マンテル・ヘンセル統計量を用いた別の均一性の検定法(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
層別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
*1:テキストに記載のepicalcパッケージの新しいもの。epicalcパッケージに入っていた関数はそのまま使える様子。cf:r - Why was package 'epicalc' removed from CRAN? - Stack Overflow cf2:整形した結果を表示する - 盆栽日記