Triad sou.

層内に相関のあるK×2×2分割表の共通オッズ比の推測について

Liang KY. Odds ratio inference with dependent data. Biometrika 1985; 72(3): 678–682.


と言ったしたものの確認したことがなかったのでシミュレーションしてみました。
マルコフ連鎖よりもベータ二項分布の方が簡単そうだったので、積ベータ二項分布モデルについて検討しています(原典と微妙な違いがあるかもしれません)。

シミュレーション用コード

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.