| Title: | Weighted NPMLE for Recurrent Events with a Competing Terminal Event |
|---|---|
| Description: | Provides regression modeling and prediction for the marginal mean of recurrent events in the presence of a competing terminal event based on the weighted nonparametric maximum likelihood estimator (wNPMLE) of Bellach and Kosorok (2026) <https://arxiv.org/abs/2605.25934>. Two large classes of transformation models are included: the Box-Cox transformation models and the logarithmic transformation models, extending the proportional means model of Ghosh and Lin (2002) <doi:10.17615/pt0g-y207> and the model of Zeng and Lin (2006) <doi:10.1093/biomet/93.3.627>. Parameter estimates are derived via automatic differentiation using the Template Model Builder (TMB) framework. Standard errors are computed via a sandwich variance estimator with correction for the estimated weights, following Bellach, Kosorok, Ruschendorf and Fine (2019) <doi:10.1080/01621459.2017.1401540>. |
| Authors: | Anna Bellach [aut, cre] |
| Maintainer: | Anna Bellach <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.1 |
| Built: | 2026-06-19 09:27:42 UTC |
| Source: | https://github.com/abellach/wnpmle |
AIC for wnpmle objects
## S3 method for class 'wnpmle' AIC(object, ..., k = 2)## S3 method for class 'wnpmle' AIC(object, ..., k = 2)
object |
A |
... |
Ignored. |
k |
Penalty per parameter (default 2 for AIC). |
Returns a data frame with the estimated cumulative baseline mean function Lambda(t) and its increments lambda(t), together with standard errors and pointwise 95% confidence intervals.
baseline(object, conf_level = 0.95, ...)baseline(object, conf_level = 0.95, ...)
object |
A |
conf_level |
Confidence level for the pointwise intervals (default 0.95). |
... |
Ignored. |
A data frame with columns:
time |
Recurrent event times. |
lambda |
Estimated baseline increments. |
Lambda |
Estimated cumulative baseline mean. |
se_Lambda |
Standard error of Lambda (if SE was estimated). |
lower |
Lower confidence band for Lambda. |
upper |
Upper confidence band for Lambda. |
## Not run: bdata <- bladder_prep() bdata_clean <- bdata[, c("id", "time", "status", "treat", "num", "size")] fit <- wnpmle_fit(Surv(time, status) ~ treat + num + size, data = bdata_clean, id = "id", model = "log", rho = 1) bl <- baseline(fit) head(bl) ## End(Not run)## Not run: bdata <- bladder_prep() bdata_clean <- bdata[, c("id", "time", "status", "treat", "num", "size")] fit <- wnpmle_fit(Surv(time, status) ~ treat + num + size, data = bdata_clean, id = "id", model = "log", rho = 1) bl <- baseline(fit) head(bl) ## End(Not run)
BIC for wnpmle objects
## S3 method for class 'wnpmle' BIC(object, ...)## S3 method for class 'wnpmle' BIC(object, ...)
object |
A |
... |
Ignored. |
Prepares the bladder1 dataset from the survival package for
use with wnpmle_fit. Subjects randomised to pyridoxine are
excluded (only placebo and thiotepa arms are retained). Status codes are
recoded to 0 (censored), 1 (recurrent tumour), 2 (terminal/dropout), and
event times beyond tau are truncated to tau.
bladder_prep(tau = 59)bladder_prep(tau = 59)
tau |
Follow-up truncation time (default: 59 months). Event times
beyond |
The bladder cancer trial (Veterans Administration Cooperative Urological
Research Group) randomised patients to placebo, pyridoxine, or thiotepa.
Only the placebo and thiotepa arms are used here (pyridoxine arm excluded).
Subjects with event times equal to tau and status 0 have their
times truncated to tau, consistent with administrative censoring.
A data frame with columns:
id |
Subject identifier (integer). |
status |
Event status: 0 = censored, 1 = recurrent event, 2 = terminal event. |
status0 |
Indicator: 1 if censored. |
status1 |
Indicator: 1 if recurrent event. |
status2 |
Indicator: 1 if terminal event. |
time |
Event time (months). |
treat |
Treatment indicator: 1 = thiotepa, 0 = placebo. |
num |
Initial number of tumours. |
size |
Initial tumour size (cm). |
Byar, D.P. (1980). The Veterans Administration study of chemoprophylaxis for recurrent stage I bladder tumors. In Bladder Tumors and Other Topics in Urological Oncology, 363-370. Plenum, New York.
bdata <- bladder_prep(tau = 59) head(bdata) table(bdata$status)bdata <- bladder_prep(tau = 59) head(bdata) table(bdata$status)
Extract coefficients from a wnpmle object
## S3 method for class 'wnpmle' coef(object, ...)## S3 method for class 'wnpmle' coef(object, ...)
object |
A |
... |
Ignored. |
Log-likelihood for wnpmle objects
## S3 method for class 'wnpmle' logLik(object, ...)## S3 method for class 'wnpmle' logLik(object, ...)
object |
A |
... |
Ignored. |
Fits the model over a fine grid of transformation parameter values and plots the profile log-likelihood for both the Box-Cox and logarithmic transformation models on a single plot. The log transformation parameter r is shown on the left (negative axis) and the Box-Cox parameter rho on the right (positive axis), meeting at zero where both models coincide.
plot_loglik( formula, data, id = "id", rho_grid = seq(0.01, 1.2, by = 0.01), r_grid = seq(0.01, 1.2, by = 0.01), tau = NULL, mark_points = TRUE, file = NULL, verbose = TRUE, ... )plot_loglik( formula, data, id = "id", rho_grid = seq(0.01, 1.2, by = 0.01), r_grid = seq(0.01, 1.2, by = 0.01), tau = NULL, mark_points = TRUE, file = NULL, verbose = TRUE, ... )
formula |
A formula as passed to |
data |
A data frame. |
id |
Name of the subject identifier column. |
rho_grid |
A numeric vector of rho values for the Box-Cox model
(default: |
r_grid |
A numeric vector of r values for the log model
(default: |
tau |
Optional truncation time. If |
mark_points |
Logical; if |
file |
Optional path to save the plot as a PDF (e.g.
|
verbose |
Logical; print progress (default |
... |
Additional arguments passed to |
A data frame with columns model, param and
loglik, invisibly.
## Not run: bdata <- bladder_prep() bdata_clean <- bdata[, c("id", "time", "status", "treat", "num", "size")] plot_loglik(Surv(time, status) ~ treat + num + size, data = bdata_clean, id = "id") ## End(Not run)## Not run: bdata <- bladder_prep() bdata_clean <- bdata[, c("id", "time", "status", "treat", "num", "size")] plot_loglik(Surv(time, status) ~ treat + num + size, data = bdata_clean, id = "id") ## End(Not run)
Plots the estimated cumulative baseline mean function Lambda(t) with optional pointwise confidence bands.
## S3 method for class 'wnpmle' plot(x, conf_bands = TRUE, ...)## S3 method for class 'wnpmle' plot(x, conf_bands = TRUE, ...)
x |
A |
conf_bands |
Logical; if |
... |
Additional graphical parameters passed to |
No return value, called for side effects (produces a plot).
Predict marginal mean for new covariate values
## S3 method for class 'wnpmle' predict(object, newdata = NULL, times = NULL, ...)## S3 method for class 'wnpmle' predict(object, newdata = NULL, times = NULL, ...)
object |
A |
newdata |
A data frame with the same covariates used in fitting.
If |
times |
Time points at which to evaluate the marginal mean.
If |
... |
Ignored. |
A data frame with columns time and one column per row of
newdata (or a single column mu for baseline).
Print method for wnpmle objects
## S3 method for class 'wnpmle' print(x, ...)## S3 method for class 'wnpmle' print(x, ...)
x |
A |
... |
Ignored. |
Summary method for wnpmle objects
## S3 method for class 'wnpmle' summary(object, tau_grid = TRUE, ...)## S3 method for class 'wnpmle' summary(object, tau_grid = TRUE, ...)
object |
A |
tau_grid |
Logical; if |
... |
Ignored. |
Returns the full variance-covariance matrix for (beta, Lambda).
To get only the beta part, use vcov(fit)[1:p, 1:p].
## S3 method for class 'wnpmle' vcov(object, ...)## S3 method for class 'wnpmle' vcov(object, ...)
object |
A |
... |
Ignored. |
Estimates the weighted nonparametric maximum likelihood estimator (wNPMLE) for the marginal mean function of recurrent events in the presence of a competing terminal event. Two transformation models are supported: the Box-Cox model and the logarithmic model.
wnpmle_fit( formula, data, id = "id", type = c("recurrent", "cmprsk", "cmprsk_ltrc"), model = c("boxcox", "log"), rho = 1, tau = NULL, se = c("sandwich_adj", "sandwich", "fisher", "none"), init_beta = NULL, control = list(), silent = TRUE )wnpmle_fit( formula, data, id = "id", type = c("recurrent", "cmprsk", "cmprsk_ltrc"), model = c("boxcox", "log"), rho = 1, tau = NULL, se = c("sandwich_adj", "sandwich", "fisher", "none"), init_beta = NULL, control = list(), silent = TRUE )
formula |
A formula of the form |
data |
A data frame containing the variables in |
id |
Name of the subject identifier column in |
type |
Type of analysis. Currently only |
model |
Transformation model: |
rho |
Transformation parameter. For |
tau |
Follow-up truncation time. Kaplan-Meier censoring weights are
truncated at |
se |
Variance estimation method: |
init_beta |
Initial values for regression coefficients (default: all zeros). |
control |
A list of control parameters passed to |
silent |
Suppress TMB output (default: |
Analysis types: Currently only type = "recurrent" is implemented,
for the recurrent events with competing terminal event setting. Support for
"cmprsk" (Bellach, Kosorok, Ruschendorf and Fine, 2019) and
"cmprsk_ltrc" (left truncation and right censoring) will be
added in future versions.
Transformation models: The two models differ in the link function :
Box-Cox: , with
as .
Logarithmic: , with
as .
Both are estimated via automatic differentiation using TMB. The sandwich
variance estimator accounts for the estimation of the censoring weights.
The censoring-corrected sandwich (sandwich_adj) adds an influence
function correction for the Kaplan-Meier censoring weight estimation.
An object of class "wnpmle" with components:
coefficients |
Estimated regression coefficients (beta). |
se |
Standard errors for the regression coefficients. |
Lambda |
Estimated cumulative baseline mean function, evaluated at the recurrent event times. |
lambda |
Estimated baseline increments. |
event_times |
Recurrent event times at which Lambda is estimated. |
loglik |
Log-likelihood at the optimum. |
vcov |
Full variance-covariance matrix for (beta, Lambda). |
model |
Transformation model used. |
rho |
Transformation parameter used. |
tau |
Truncation time used. |
n |
Number of subjects. |
n_events |
Named vector: recurrent events, terminal events, censored. |
convergence |
Convergence message from |
call |
The matched call. |
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.
Bellach, A., Kosorok, M.R., Ruschendorf, 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.
library(survival) data("bladder2", package = "survival") bladder2_prepped <- bladder_prep() fit <- wnpmle_fit(Surv(time, status) ~ treat + num + size, data = bladder2_prepped, id = "id", model = "log", rho = 1) summary(fit)library(survival) data("bladder2", package = "survival") bladder2_prepped <- bladder_prep() fit <- wnpmle_fit(Surv(time, status) ~ treat + num + size, data = bladder2_prepped, id = "id", model = "log", rho = 1) summary(fit)