README

An R package implementing the weighted nonparametric maximum likelihood estimator (wNPMLE) for the marginal mean of recurrent events when a competing terminal event (e.g. death) is present.

Background

In many clinical studies, subjects experience repeated events (e.g. COPD exacerbations, bladder tumor recurrences, recurrent hospitalizations) while they are also exposed to a competing terminal event. Treating death as independent censoring leads to overestimation of the expected number of recurrences. This package implements the weighted nonparametric maximum likelihood estimators by Bellach and Kosorok (2026), which correctly account for the terminal event through inverse probability of censoring weighting.

When no terminal events are present, the package automatically reduces to the unweighted NPMLE, which is equivalent to the Zeng-Lin model (Zeng and Lin, 2006).

Both models are estimated via automatic differentiation (TMB). Standard errors are available via the adjusted sandwich estimator with correction for the estimated weights, a sandwich estimator without correction as an approximation, and the inverse Fisher information.

Installation

# install.packages("remotes")
remotes::install_github("abellach/wnpmle")

Quick start

library(wnpmle)

# Prepare bladder cancer data
bdata       <- bladder_prep()
bdata_clean <- bdata[, c("id", "time", "status", "treat", "num", "size")]

# --- Ghosh-Lin model (Box-Cox, rho = 1) ---
fit_bc <- wnpmle_fit(
  Surv(time, status) ~ treat + num + size,
  data  = bdata_clean,
  id    = "id",
  model = "boxcox",
  rho   = 1,
  se    = "sandwich_adj"
)
summary(fit_bc)
plot(fit_bc)
baseline(fit_bc)
AIC(fit_bc)
BIC(fit_bc)

# --- Proportional odds model (Log, r = 1) ---
fit_log <- wnpmle_fit(
  Surv(time, status) ~ treat + num + size,
  data  = bdata_clean,
  id    = "id",
  model = "log",
  rho   = 1,
  se    = "sandwich_adj"
)
summary(fit_log)
plot(fit_log)
baseline(fit_log)
AIC(fit_log)
BIC(fit_log)

# --- Prediction: marginal mean at new covariate values ---
newdat <- data.frame(treat = c(0, 1), num = c(1, 1), size = c(1, 1))
pred <- predict(fit_bc, newdata = newdat, times = seq(1, 50, by = 1))
plot(pred$time, pred$mu_1, type = "s", lwd = 2,
     xlab = "Time (months)", ylab = "Marginal mean number of recurrences",
     ylim = range(pred[, -1]))
lines(pred$time, pred$mu_2, lwd = 2, lty = 2, col = "firebrick")
legend("topleft", legend = c("Placebo", "Thiotepa"),
       lty = c(1, 2), col = c("black", "firebrick"), bty = "n")

# --- Profile log-likelihood for transformation parameter ---
result <- plot_loglik(
  Surv(time, status) ~ treat + num + size,
  data     = bdata_clean,
  id       = "id",
  tau      = 59,
  rho_grid = seq(0.01, 1.2, by = 0.01),
  r_grid   = seq(0.01, 1.2, by = 0.01)
)

# Table of log-likelihood values
head(result)

# Optimal transformation parameters
result[result$model == "boxcox", ][which.max(result$loglik[result$model == "boxcox"]), ]
result[result$model == "log",    ][which.max(result$loglik[result$model == "log"   ]), ]

Status codes

Data format

Users can bring their own data. The required format is one row per event per subject with columns:

Citation

Status	Meaning
0	Censored
1	Recurrent event
2	Terminal event (e.g. death)

Bellach, A. and Kosorok, M.R. (2026). Weighted NPMLE for the marginal mean of recurrent events with a competing terminal event. arXiv preprint arXiv:2605.25934

BibTeX

@misc{bellach2026wnpmle,
  title         = {Weighted {NPMLE} for the marginal mean of recurrent events
                   with a competing terminal event},
  author        = {Bellach, Anna and Kosorok, Michael R.},
  year          = {2026},
  eprint        = {2605.25934},
  archivePrefix = {arXiv},
  url           = {https://arxiv.org/abs/2605.25934}
}

References

Bellach, A., Kosorok, M.R., Rüschendorf, L. and Fine, J.P. (2019). Weighted NPMLE for the subdistribution of a competing risk. Journal of the American Statistical Association, 114(525), 259-270. doi:10.1080/01621459.2017.1401540

Ghosh, D. and Lin, D.Y. (2002). Marginal regression models for recurrent and terminal events. Statistica Sinica, 12, 663-688. doi:10.17615/pt0g-y207

Zeng, D. and Lin, D.Y. (2006). Efficient estimation of semiparametric transformation models for counting processes. Biometrika, 93(3), 627-640. doi:10.1093/biomet/93.3.627

wnpmle