Package {wnpmle}

wnpmle: Weighted NPMLE for Recurrent Events with a Competing Terminal Event

Description

Implements the weighted nonparametric maximum likelihood estimator (wNPMLE) for the marginal mean function of recurrent events in the presence of a competing terminal event. Two semiparametric transformation models are provided: the Box-Cox transformation model and the logarithmic transformation model. Both use automatic differentiation via TMB for efficient and exact gradient computation.

Details

The main user-facing function is wnpmle_fit.

Author(s)

Maintainer: Anna Bellach abellach.biostat@gmail.com

References

Bellach, A. and Kosorok, M.R. (2026). Weighted NPMLE for the marginal mean of recurrent events with a competing terminal event. JASA, to appear.

See Also

Useful links:

• https://github.com/abellach/wnpmle

• Report bugs at https://github.com/abellach/wnpmle/issues

AIC for wnpmle objects

Description

AIC for wnpmle objects

Usage

## S3 method for class 'wnpmle'
AIC(object, ..., k = 2)

Arguments

Value

A numeric scalar giving the AIC value.

BIC for wnpmle objects

Description

BIC for wnpmle objects

Usage

## S3 method for class 'wnpmle'
BIC(object, ...)

Arguments

Value

A numeric scalar giving the BIC value.

Extract the estimated baseline mean function

Description

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.

Usage

baseline(object, conf_level = 0.95, ...)

Arguments

Value

A data frame with columns:

Examples

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)

Prepare bladder cancer data for wnpmle analysis

Description

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.

Usage

bladder_prep(tau = 59)

Arguments

Details

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.

Value

A data frame with columns:

References

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.

Examples

bdata <- bladder_prep(tau = 59)
head(bdata)
table(bdata$status)

Extract coefficients from a wnpmle object

Description

Extract coefficients from a wnpmle object

Usage

## S3 method for class 'wnpmle'
coef(object, ...)

Arguments

Value

A named numeric vector of regression coefficients.

Log-likelihood for wnpmle objects

Description

Log-likelihood for wnpmle objects

Usage

## S3 method for class 'wnpmle'
logLik(object, ...)

Arguments

Value

An object of class "logLik" with the log-likelihood value, degrees of freedom (df) equal to the number of regression coefficients, and number of observations (nobs).

Plot method for wnpmle objects

Description

Plots the estimated cumulative baseline mean function Lambda(t) with optional pointwise confidence bands.

Usage

## S3 method for class 'wnpmle'
plot(x, conf_bands = TRUE, ...)

Arguments

Value

No return value, called for side effects (produces a plot).

Log-likelihood profile plot for the transformation parameter

Description

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.

Usage

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,
  ...
)

Arguments

Value

A data frame with columns model, param and loglik, invisibly.

Examples

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")

Predict marginal mean for new covariate values

Description

Predict marginal mean for new covariate values

Usage

## S3 method for class 'wnpmle'
predict(object, newdata = NULL, times = NULL, ...)

Arguments

Value

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

Description

Print method for wnpmle objects

Usage

## S3 method for class 'wnpmle'
print(x, ...)

Arguments

Value

Invisibly returns x.

Summary method for wnpmle objects

Description

Summary method for wnpmle objects

Usage

## S3 method for class 'wnpmle'
summary(object, tau_grid = TRUE, ...)

Arguments

Value

Invisibly returns object.

Extract variance-covariance matrix from a wnpmle object

Description

Returns the full variance-covariance matrix for (beta, Lambda). To get only the beta part, use vcov(fit)[1:p, 1:p].

Usage

## S3 method for class 'wnpmle'
vcov(object, ...)

Arguments

Value

A numeric matrix containing the variance-covariance matrix for the regression coefficients and cumulative baseline mean function. Returns an error if se = "none" was used in wnpmle_fit.

Fit Weighted NPMLE for Survival Data with Recurrent or Competing Events

Description

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.

Usage

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
)

Arguments

Details

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

• Box-Cox: G(x) = ((1 + x)^\rho - 1) / \rho, with G(x) \to \log(1+x) as \rho \to 0.

• Logarithmic: G(x) = \log(1 + r \cdot x) / r, with G(x) \to x as r \to 0.

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.

Value

An object of class "wnpmle" with components:

References

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.

Examples

 
  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)