Liang KY. Odds ratio inference with dependent data. Biometrika 1985; 72(3): 678–682.
k x 2 x 2表の共通オッズ比の推定で、独立な積二項モデルの仮定のもとではCMLEもMH推定量も一致性があるが、層内で相関がある(ベータ二項モデルなど)とMH推定量しか一致性を持たない。
— Triad sou. (@triadsou) 2020年2月14日
と言ったしたものの確認したことがなかったのでシミュレーションしてみました。
マルコフ連鎖よりもベータ二項分布の方が簡単そうだったので、積ベータ二項分布モデルについて検討しています(原典と微妙な違いがあるかもしれません)。
シミュレーション用コード
require("VGAM") require("reshape2") require("foreach") require("doParallel") mhsim <- function(k, phi, rho, n = 5, m = 5, r = 5000) { p2k0 <- 0.2 cluster <- makeCluster(detectCores() - 1) registerDoParallel(cluster) res <- foreach(i = 1:r, .combine = "rbind") %dopar% { tab <- array(dim = c(2, 2, k)) p2 <- p2k0 for (l in 1:k) { p2 <- p2 + 0.6/k p1 <- (p2*phi) / (1 - p2 + p2*phi) n11 <- VGAM::rbetabinom(1, n, p2, rho) m11 <- VGAM::rbetabinom(1, m, p1, rho) tab[2, 1, l] <- n11 tab[1, 1, l] <- m11 tab[2, 2, l] <- n - n11 tab[1, 2, l] <- m - m11 } dat1 <- reshape2::melt(tab) dat1 <- data.frame(x=rep(dat1$Var1 - 1, dat1$value), y=rep(dat1$Var2 - 1, dat1$value), k=rep(dat1$Var3, dat1$value)) mh <- mantelhaen.test(tab) cmle <- mantelhaen.test(tab, exact = TRUE) mle <- glm(y ~ factor(x) + factor(k), family = binomial(link = "logit"), data = dat1) data.frame(mh = mh$estimate, cmle = cmle$estimate, mle = exp(mle$coefficients[2])) } stopCluster(cluster) res } summary(mhsim(k = 25, phi = 5, rho = 0.2)) summary(mhsim(k = 50, phi = 5, rho = 0.2)) summary(mhsim(k = 100, phi = 5, rho = 0.2)) summary(mhsim(k = 500, phi = 5, rho = 0.2, r = 100)) summary(mhsim(k = 100, phi = 5, rho = 0.0)) summary(mhsim(k = 500, phi = 5, rho = 0.0, n = 1, m = 1, r = 100))
相関の影響についての結果
オッズ比$\phi=5$、相関のパラメータ$\rho=0.2$、各層内のサンプルサイズが$N_k=M_k=5$、$K=25, 50, 100$のときを検討します。
mh:Mantel–Haenszel推定量、cmle:条件付き最尤推定量、mle:無条件の最尤推定量です。
Mantel–Haenszel推定量は$K$が大きければ$\phi=5$に近づいていきます。
条件付き最尤推定量は$K$が大きくなっても$\phi=5$から離れたままです。
無条件の最尤推定量はもともと一致性がありません。
> summary(mhsim(k = 25, phi = 5, rho = 0.2)) mh cmle mle Min. : 1.086 Min. : 1.098 Min. : 1.109 1st Qu.: 3.833 1st Qu.: 4.402 1st Qu.: 5.236 Median : 5.167 Median : 6.059 Median : 7.497 Mean : 6.018 Mean : 7.074 Mean : 9.313 3rd Qu.: 7.125 3rd Qu.: 8.459 3rd Qu.: 11.038 Max. :89.333 Max. :115.649 Max. :231.215 > summary(mhsim(k = 50, phi = 5, rho = 0.2)) mh cmle mle Min. : 2.042 Min. : 2.148 Min. : 2.340 1st Qu.: 4.130 1st Qu.: 4.786 1st Qu.: 5.764 Median : 5.064 Median : 5.968 Median : 7.395 Mean : 5.410 Mean : 6.381 Mean : 8.096 3rd Qu.: 6.261 3rd Qu.: 7.511 3rd Qu.: 9.606 Max. :35.167 Max. :32.358 Max. :54.375 > summary(mhsim(k = 100, phi = 5, rho = 0.2)) mh cmle mle Min. : 2.458 Min. : 2.770 Min. : 3.109 1st Qu.: 4.365 1st Qu.: 5.101 1st Qu.: 6.186 Median : 5.041 Median : 5.961 Median : 7.392 Mean : 5.218 Mean : 6.167 Mean : 7.730 3rd Qu.: 5.880 3rd Qu.: 7.012 3rd Qu.: 8.899 Max. :14.608 Max. :16.252 Max. :24.016 > summary(mhsim(k = 500, phi = 5, rho = 0.2, r = 100)) mh cmle mle Min. :3.928 Min. :4.636 Min. : 5.555 1st Qu.:4.757 1st Qu.:5.628 1st Qu.: 6.930 Median :5.061 Median :6.040 Median : 7.504 Mean :5.098 Mean :6.042 Mean : 7.514 3rd Qu.:5.433 3rd Qu.:6.390 3rd Qu.: 7.995 Max. :6.756 Max. :7.992 Max. :10.339
おまけ:無相関の積ベータ二項分布モデル
相関がない積ベータ二項分布モデルは積二項分布モデルに等しく、Mantel–Haenszel推定量も条件付き最尤推定量も一致性があります。
> summary(mhsim(k = 100, phi = 5, rho = 0.0)) mh cmle mle Min. :2.945 Min. :3.050 Min. : 3.465 1st Qu.:4.498 1st Qu.:4.495 1st Qu.: 5.363 Median :5.005 Median :5.007 Median : 6.059 Mean :5.090 Mean :5.075 Mean : 6.169 3rd Qu.:5.597 3rd Qu.:5.573 3rd Qu.: 6.842 Max. :9.277 Max. :9.877 Max. :13.188
おまけ:無条件のMLE
よく知られた結果ですが、無条件のMLEは1:1マッチングのときに$\phi^2$に収束します(Breslow, 1981)。
> summary(mhsim(k = 500, phi = 5, rho = 0.0, n = 1, m = 1, r = 100)) mh cmle mle Min. :3.423 Min. :3.423 Min. :11.72 1st Qu.:4.311 1st Qu.:4.310 1st Qu.:18.58 Median :4.872 Median :4.872 Median :23.73 Mean :5.013 Mean :5.012 Mean :25.99 3rd Qu.:5.457 3rd Qu.:5.457 3rd Qu.:29.78 Max. :7.966 Max. :7.965 Max. :63.45
- Breslow N. Odds ratio estimators when the data are sparse. Biometrika 1981; 68(1): 73–84.