The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Introduction to coxmnar: Cox Regression with Missing Not At Random Failure Indicators

Shikhar Tyagi, Arvind Pandey, Bhupendra Singh, Vrijesh Tripathi

2026-07-23

1. Motivation

In survival analysis and epidemiological studies, cause-of-failure information (failure indicators) is often subject to missingness due to incomplete data collection, loss to follow-up, or diagnostic uncertainty. When the missingness mechanism depends on the unobserved failure indicator itself, the missingness is Missing Not At Random (MNAR). Standard estimation methods such as Complete-Case (CC) analysis or Missing At Random (MAR) imputation lead to severe estimation bias and invalid inference under MNAR.

The coxmnar package implements the two adjusted imputation-based estimating equations developed by Liu and Liu (2026):

  1. Adjusted Imputation (AI) estimator (\(\hat{\beta}_n^{AI}\))
  2. Adjusted Complete (AC) estimator (\(\hat{\beta}_n^{AC}\))

These estimators incorporate joint maximum likelihood estimation of the outcome and propensity models alongside Nadaraya-Watson kernel propensity smoothing to achieve consistency and asymptotic normality under MNAR.

2. Input Data Format

coxmnar() expects survival data containing: - Observed failure time (time). - Cause/failure indicator (cause), where 1 indicates an event of interest, 0 indicates right-censoring, and NA indicates an unobserved/missing failure indicator. - One or more covariates specified via a right-hand-side formula (e.g., cause ~ Z1 + Z2).

Here is a small 10-row example:

set.seed(123)
toy_df <- simulate_coxmnar_data(n = 10, beta0 = 0.2, mechanism = "mnar", seed = 123)
head(toy_df, 10)
#>          time cause          Z
#> 1  0.01862152     1  4.3136628
#> 2  0.07433200     1 11.8245770
#> 3  0.11412212     1  6.1346538
#> 4  0.03942608     1 13.2452611
#> 5  0.13533962     1 14.1070093
#> 6  0.09207125     1  0.6833475
#> 7  0.28754705     1  7.9215823
#> 8  0.21785067     1 13.3862857
#> 9  0.21321259     1  8.2715252
#> 10 0.01183403     1  6.8492210

3. Simulated Data Walkthrough

We simulate a dataset with \(n = 300\), true coefficient \(\beta_0 = 0.2\), censoring rate \(\approx 20\%\), and missingness rate \(\approx 20\%\) under an MNAR mechanism. We fit all five supported methods using coxmnar_compare():

set.seed(42)
sim_data <- simulate_coxmnar_data(n = 300, beta0 = 0.2, mechanism = "mnar", seed = 42)

comp <- coxmnar_compare(
  formula = cause ~ Z,
  time = "time",
  data = sim_data,
  methods = c("ai", "ac", "mar", "cc", "full"),
  B = 100L,
  seed = 42
)

summary(comp)
#> Call:
#> coxmnar_compare(formula = cause ~ Z, time = "time", data = sim_data, 
#>     methods = c("ai", "ac", "mar", "cc", "full"), B = 100L, seed = 42)
#> 
#> Comparative Coxmnar Fit Table:
#> 
#>  Method Covariate  Estimate         SE        z       pvalue
#>      AI         Z 0.1441008 0.02206463 6.530847 6.539880e-11
#>      AC         Z 0.1386539 0.02145286 6.463190 1.025182e-10
#>     MAR         Z 0.2042798 0.02122210 9.625806 0.000000e+00
#>      CC         Z 0.2040264 0.02474211 8.246120 2.220446e-16
#>    FULL         Z 0.1386539 0.02145286 6.463190 1.025182e-10

As demonstrated in the comparison table, the proposed AI and AC estimators closely recover the true coefficient (\(\beta_0 = 0.2\)), whereas the Complete-Case (CC) estimator exhibits significant negative bias.

4. Bootstrap vs. Asymptotic Standard Errors

Liu and Liu (2026) show that while asymptotic sandwich standard errors (ASE-T) are theoretically consistent, they can be unstable in finite samples. Nonparametric bootstrap standard errors (ASE) provide superior coverage probabilities:

fit_both <- coxmnar(
  formula = cause ~ Z,
  time = "time",
  data = sim_data,
  method = "ai",
  se_method = "both",
  B = 100L,
  seed = 42
)

cat("Bootstrap SE:", fit_boot_se <- sqrt(fit_both$vcov_boot[1, 1]), "\n")
#> Bootstrap SE: 0.02206463
cat("Asymptotic SE:", fit_asymp_se <- sqrt(fit_both$vcov_asymp[1, 1]), "\n")
#> Asymptotic SE: 0.01715322

5. Real-Data Illustration (survival::pbc)

We illustrate coxmnar() using the classic Primary Biliary Cirrhosis dataset (survival::pbc). To demonstrate MNAR missingness handling, we synthetically induce MNAR missingness on pbc$status:

data(pbc, package = "survival")
pbc_clean <- na.omit(pbc[, c("time", "status", "bili", "albumin")])
pbc_clean$cause <- ifelse(pbc_clean$status == 2, 1, 0)
pbc_clean$log_bili <- log(pbc_clean$bili)
pbc_clean$log_alb <- log(pbc_clean$albumin)

# Synthetically induce 25% MNAR missingness
set.seed(999)
prob_obs <- 1 / (1 + 0.3 * (pbc_clean$log_bili^(0.1 + 2 * pbc_clean$cause)))
xi <- rbinom(nrow(pbc_clean), 1, prob_obs)
#> Warning in rbinom(nrow(pbc_clean), 1, prob_obs): NAs produced
pbc_clean$cause[xi == 0] <- NA

fit_pbc_ai <- coxmnar(
  cause ~ log_bili + log_alb,
  time = "time",
  data = pbc_clean,
  method = "ai",
  B = 100L,
  seed = 999
)

summary(fit_pbc_ai)
#> Call:
#> coxmnar(formula = cause ~ log_bili + log_alb, time = "time", 
#>     data = pbc_clean, method = "ai", B = 100L, seed = 999)
#> 
#> Cox Regression with MNAR Failure Indicators (Method: AI)
#> n = 418, number of events = 108, missing indicators = 88 (21.1%)
#> 
#>                coef  exp(coef)  se(coef)         z     Pr(>|z|)   lower .95
#> log_bili  0.5519898 1.73670532 0.0942715  5.855320 4.760915e-09 1.443717063
#> log_alb  -3.2006060 0.04073751 0.8384715 -3.817191 1.349795e-04 0.007875647
#>          upper .95
#> log_bili 2.0891527
#> log_alb  0.2107185

We can plot the estimated survival curve alongside a baseline Kaplan-Meier curve:

plot(fit_pbc_ai, compare_km = TRUE)

6. References

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.