---
title: "Introduction to coxmnar: Cox Regression with Missing Not At Random Failure Indicators"
author: "Shikhar Tyagi, Arvind Pandey, Bhupendra Singh, Vrijesh Tripathi"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to coxmnar}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4.5
)
library(coxmnar)
library(survival)
```

## 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:

```{r toy-data}
set.seed(123)
toy_df <- simulate_coxmnar_data(n = 10, beta0 = 0.2, mechanism = "mnar", seed = 123)
head(toy_df, 10)
```

---

## 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()`:

```{r sim-walkthrough}
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)
```

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:

```{r se-comparison}
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")
cat("Asymptotic SE:", fit_asymp_se <- sqrt(fit_both$vcov_asymp[1, 1]), "\n")
```

---

## 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`:

```{r pbc-illustration}
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)
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)
```

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

```{r pbc-plot, fig.width = 6, fig.height = 4.5}
plot(fit_pbc_ai, compare_km = TRUE)
```

---

## 6. References

- Liu, Y., & Liu, K. (2026). Estimation in the Cox proportional hazards model with missing not at random failure indicators. *Statistics and Computing*, 36, 112. <doi:10.1007/s11222-026-10857-1>