| Title: | Cox Regression with Missing not at Random Failure Indicators |
| Version: | 0.1.0 |
| Description: | Implements estimation for the Cox (1972, 1975) <doi:10.1111/j.2517-6161.1972.tb00899.x> <doi:10.1093/biomet/62.2.269> proportional hazards model when the failure indicator (cause of failure) is missing not at random (MNAR), following the two adjusted imputation-based estimating equations of Liu and Liu (2026) <doi:10.1007/s11222-026-10857-1>. Also provided for comparison are the full-data partial-likelihood estimator of Andersen and Gill (1982) <doi:10.1214/aos/1176345976>, the complete-case estimator, and the missing-at-random imputation estimator of Liu and Wang (2010, Statistica Sinica, 20, 1125-1142). The probability models for the failure indicator and for the missingness mechanism are estimated jointly by maximum likelihood following Sun, Xie, and Liang (2013) <doi:10.1007/s11425-012-4492-x>, and a Nadaraya-Watson kernel-smoothed estimator of the missingness propensity is constructed following Qiu, Chen, and Zhou (2015) <doi:10.1016/j.spl.2014.12.006>. Both an asymptotic (sandwich-type) variance estimator and a nonparametric bootstrap variance estimator are provided. When failure indicators are fully observed the estimators reduce algebraically to the classical Cox partial-likelihood estimator. |
| License: | GPL (≥ 3) |
| Encoding: | UTF-8 |
| Language: | en-US |
| Depends: | R (≥ 4.1.0) |
| Imports: | survival (≥ 3.5-0), stats, Rdpack |
| Suggests: | testthat (≥ 3.0.0), knitr, rmarkdown, covr |
| RoxygenNote: | 7.3.3 |
| VignetteBuilder: | knitr |
| RdMacros: | Rdpack |
| Config/testthat/edition: | 3 |
| NeedsCompilation: | no |
| Packaged: | 2026-07-22 23:23:20 UTC; shikhar tyagi |
| Author: | Shikhar Tyagi |
| Maintainer: | Shikhar Tyagi <shikhar1093tyagi@gmail.com> |
| Repository: | CRAN |
| Date/Publication: | 2026-07-31 15:10:02 UTC |
coxmnar: Cox Regression with Missing Not at Random Failure Indicators
Description
Implements estimation for the Cox proportional hazards model when failure indicators are missing not at random (MNAR) following the methodology of Liu and Liu (2026).
Author(s)
Maintainer: Shikhar Tyagi shikhar1093tyagi@gmail.com (ORCID)
Authors:
Arvind Pandey arvindmzu@gmail.com
Bhupendra Singh bhupendra.rana@gmail.com
Vrijesh Tripathi vrijesh.tripathi@uwi.edu
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.
Asymptotic and Nonparametric Bootstrap Variance Estimation
Description
Asymptotic and Nonparametric Bootstrap Variance Estimation
Usage
asymptotic_variance(
val_data,
beta_hat,
theta_hat,
alpha_hat,
m_formula,
pi_formula,
w,
method,
control = coxmnar_control()
)
Breslow Baseline Cumulative Hazard and Survival Functions
Description
Breslow Baseline Cumulative Hazard and Survival Functions
Usage
breslow_baseline_hazard(Z, time, w, beta)
Fit probability models m(w, theta) and pi(w, alpha) jointly by maximum likelihood
Description
Fit probability models m(w, theta) and pi(w, alpha) jointly by maximum likelihood
Usage
build_design_matrix(formula, data)
References
Byrd RH, Lu P, Nocedal J, Zhu C (1995). “A limited memory algorithm for bound constrained optimization.” SIAM Journal on Scientific Computing, 16(5), 1190–1208. doi:10.1137/0916069.
Sun Z, Xie T, Liang H (2013). “Statistical inference for right-censored data with nonignorable missing censoring indicators.” Science China Mathematics, 56(6), 1263–1278. doi:10.1007/s11425-012-4492-x.
Reverse cumulative risk set calculations, Cox score, Hessian, and Newton-Raphson solver
Description
Reverse cumulative risk set calculations, Cox score, Hessian, and Newton-Raphson solver
Usage
compute_cox_risk_sums(Z, time, beta)
References
Cox DR (1972). “Regression models and life-tables.” Journal of the Royal Statistical Society, Series B, 34(2), 187–220. doi:10.1111/j.2517-6161.1972.tb00899.x.
Cox DR (1975). “Partial likelihood.” Biometrika, 62(2), 269–276. doi:10.1093/biomet/62.2.269.
Cox Regression with Missing Not At Random Failure Indicators
Description
Fits the Cox proportional hazards model when failure indicators (causes of failure) are missing not at random (MNAR), using the adjusted imputation estimators of Liu and Liu (2026). Also supports complete-case, missing-at-random (MAR), and benchmark full-data estimation.
Usage
coxmnar(
formula,
time,
cause = NULL,
observed = NULL,
data,
method = c("ai", "ac", "mar", "cc", "full"),
m_formula = NULL,
pi_formula = NULL,
kernel = c("gaussian", "triangular", "epanechnikov", "quartic", "triweight"),
bandwidth = NULL,
se_method = c("bootstrap", "asymptotic", "both", "none"),
B = 500L,
control = coxmnar_control(),
ties = c("efron", "breslow"),
x = FALSE,
seed = NULL
)
Arguments
formula |
Formula specifying covariates, e.g. |
time |
Name of time column in |
cause |
Optional name of failure indicator column in |
observed |
Optional name of missingness indicator column |
data |
Data frame containing the variables. |
method |
Estimation method: '"ai"' (Adjusted-Imputation MNAR), '"ac"' (Adjusted-Complete MNAR), '"mar"' (Missing-at-Random), '"cc"' (Complete-Case), '"full"' (Full-Data benchmark). Default is '"ai"'. |
m_formula |
One-sided formula for outcome model |
pi_formula |
One-sided formula for propensity model |
kernel |
Kernel function for Nadaraya-Watson estimator: '"gaussian"', '"triangular"', '"epanechnikov"', '"quartic"', '"triweight"'. Default is '"gaussian"'. |
bandwidth |
Optional numeric scalar bandwidth for Nadaraya-Watson estimator. Calculated automatically if NULL. |
se_method |
Standard error estimation method: '"bootstrap"', '"asymptotic"', '"both"', '"none"'. Default is '"bootstrap"'. |
B |
Number of bootstrap replications. Default is 500L. |
control |
List of control parameters from |
ties |
Tie-handling method: '"efron"' or '"breslow"'. Default is '"efron"'. |
x |
Logical indicating whether to return the covariate matrix. Default is FALSE. |
seed |
Optional random seed for bootstrap reproducibility. |
Value
An object of S3 class "coxmnar", containing estimated coefficients, variance-covariance matrix, standard errors, z-scores, p-values, fitted probability models, baseline hazard, and model diagnostics.
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.
Andersen PK, Gill RD (1982). “Cox's regression model for counting processes: a large sample study.” The Annals of Statistics, 10(4), 1100–1120. doi:10.1214/aos/1176345976.
Liu C, Wang Q (2010). “Semiparametric estimation for regression coefficients in the Cox model with failure indicators missing at random.” Statistica Sinica, 20(3), 1125–1142.
Sun Z, Xie T, Liang H (2013). “Statistical inference for right-censored data with nonignorable missing censoring indicators.” Science China Mathematics, 56(6), 1263–1278. doi:10.1007/s11425-012-4492-x.
Qiu Z, Chen X, Zhou Y (2015). “A kernel-assisted imputation estimating method for the additive hazards model with missing censoring indicator.” Statistics & Probability Letters, 98, 89–97. doi:10.1016/j.spl.2014.12.006.
Examples
set.seed(123)
sim_data <- simulate_coxmnar_data(n = 100, beta0 = 0.2, mechanism = "mnar", seed = 123)
fit <- coxmnar(cause ~ Z, time = "time", data = sim_data, method = "ai", B = 10L, seed = 123)
print(fit)
summary(fit)
S3 Methods for coxmnar Objects
Description
Print, summary, prediction, confidence intervals, survival curves, and plotting for coxmnar objects.
Usage
## S3 method for class 'coxmnar'
print(x, ...)
## S3 method for class 'coxmnar'
summary(object, conf.level = 0.95, ...)
## S3 method for class 'summary.coxmnar'
print(x, ...)
## S3 method for class 'coxmnar'
coef(object, ...)
## S3 method for class 'coxmnar'
vcov(object, ...)
## S3 method for class 'coxmnar'
confint(object, parm, level = 0.95, ...)
## S3 method for class 'coxmnar'
predict(
object,
newdata = NULL,
type = c("lp", "risk", "survival", "cumhaz"),
times = NULL,
...
)
## S3 method for class 'coxmnar'
survfit(formula, newdata = NULL, ...)
## S3 method for class 'coxmnar'
plot(x, compare_km = TRUE, ...)
## S3 method for class 'coxmnar_compare'
summary(object, ...)
## S3 method for class 'summary.coxmnar_compare'
print(x, ...)
Arguments
x |
An object of class |
... |
Additional arguments passed to methods. |
object |
An object of class |
conf.level |
Confidence level for |
parm |
Specification of parameters to include in |
level |
Confidence level for |
newdata |
Optional data frame of new covariate values for |
type |
Prediction type: '"lp"' (linear predictor), '"risk"' (relative risk exp(lp)), '"survival"', '"cumhaz"'. Default is '"lp"'. |
times |
Numeric vector of times at which to evaluate survival/cumulative hazard predictions. |
formula |
Object of class |
compare_km |
Logical indicating whether to overlay Kaplan-Meier curve in |
Value
Various objects depending on the method called.
Compare Multiple Coxmnar Estimation Methods
Description
Convenience wrapper to fit and compare multiple estimation methods on the same dataset.
Usage
coxmnar_compare(
formula,
time,
data,
methods = c("ai", "ac", "mar", "cc", "full"),
...
)
Arguments
formula |
Formula specifying covariates. |
time |
Name of time column or numeric vector. |
data |
Data frame. |
methods |
Vector of method names to fit. Default is |
... |
Additional arguments passed to |
Value
An object of class "coxmnar_compare", containing a list of fits and a comparison table.
Examples
set.seed(123)
sim_data <- simulate_coxmnar_data(n = 100, beta0 = 0.2, mechanism = "mnar", seed = 123)
comp <- coxmnar_compare(cause ~ Z, time = "time", data = sim_data, B = 10L, seed = 123)
summary(comp)
Control Parameters for coxmnar
Description
Specifies numerical optimization, tolerance, and control parameters for coxmnar.
Usage
coxmnar_control(
maxit = 200L,
tol = 1e-08,
eps = 1e-06,
nn_fallback_k = NULL,
boot_max_attempts_factor = 10L,
optim_method = "L-BFGS-B",
step_halving_max = 20L,
on_failure = c("error", "warn"),
verbose = FALSE
)
Arguments
maxit |
Maximum number of iterations for optimization/Newton-Raphson. Default is 200L. |
tol |
Convergence tolerance for Newton-Raphson score norm. Default is 1e-8. |
eps |
Small positive probability clamp parameter. Default is 1e-6. |
nn_fallback_k |
Number of nearest neighbors for Nadaraya-Watson fallback. Default is 'max(10, ceiling(0.1*n))'. |
boot_max_attempts_factor |
Maximum multiplier for bootstrap replicate attempts. Default is 10L. |
optim_method |
Optimization method passed to |
step_halving_max |
Maximum number of step-halvings in Newton-Raphson solver. Default is 20L. |
on_failure |
Character string specifying action on convergence failure: '"error"' (default) or '"warn"'. |
verbose |
Logical indicating whether to print diagnostic messages. Default is FALSE. |
Value
A validated control list.
Examples
ctrl <- coxmnar_control(maxit = 100L, tol = 1e-6)
Evaluate kernel functions
Description
Computes the kernel function value for 1D inputs. Supported kernels: '"gaussian"', '"triangular"', '"epanechnikov"', '"quartic"' (biweight), '"triweight"'.
Usage
eval_kernel(
u,
kernel = c("gaussian", "triangular", "epanechnikov", "quartic", "triweight")
)
Arguments
u |
Numeric vector or matrix. |
kernel |
Character string specifying kernel type. |
Value
Numeric array of same dimension as 'u'.
Nadaraya-Watson kernel propensity estimator u_hat_n(w)
Description
Nadaraya-Watson kernel propensity estimator u_hat_n(w)
Usage
nadaraya_watson(
W_df,
xi,
eval_W_df = W_df,
kernel = "gaussian",
bandwidth = NULL,
control = coxmnar_control()
)
References
Qiu Z, Chen X, Zhou Y (2015). “A kernel-assisted imputation estimating method for the additive hazards model with missing censoring indicator.” Statistics & Probability Letters, 98, 89–97. doi:10.1016/j.spl.2014.12.006.
Simulate Survival Data with MNAR or MAR Failure Indicators
Description
Generates synthetic survival data following the Cox proportional hazards model and missingness data-generating processes (DGPs) described in Sections 3.1-3.3 of Liu and Liu (2026).
Usage
simulate_coxmnar_data(
n,
beta0 = 0.2,
mechanism = c("mnar", "mar"),
covariate_dist = function(n) stats::runif(n, 0, 15),
censoring_rate = 0.2,
missing_rate = 0.2,
conditionally_independent_censoring = FALSE,
extended_missingness = FALSE,
seed = NULL
)
Arguments
n |
Sample size. |
beta0 |
True regression coefficient scalar or vector. Default is 0.2. |
mechanism |
Missingness mechanism: '"mnar"' (Missing Not At Random) or '"mar"' (Missing At Random). Default is '"mnar"'. |
covariate_dist |
Function generating covariate values for sample size 'n'. Default generates Uniform(0, 15). |
censoring_rate |
Target censoring rate in (0, 1). Default is 0.2. |
missing_rate |
Target missing indicator rate in (0, 1). Default is 0.2. |
conditionally_independent_censoring |
Logical indicating whether censoring depends on covariate Z. Default is FALSE. |
extended_missingness |
Logical indicating whether missingness depends on observed time X (Section 3.3). Default is FALSE. |
seed |
Optional random seed. |
Value
A data frame with columns 'time', 'cause' (with 'NA' where 'xi=0'), 'Z' (or 'Z1'), and attributes recording ground truth parameters.
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.
Examples
set.seed(123)
sim_df <- simulate_coxmnar_data(n = 200, beta0 = 0.2, mechanism = "mnar", seed = 123)
head(sim_df)