| Type: | Package |
| Title: | Correlated Destructive Negative Binomial Cure Rate Model |
| Version: | 1.1.1 |
| Maintainer: | Rodrigo M. R. de Medeiros <rodrigo.matheus@ufrn.br> |
| Description: | Provides tools for modeling time-to-event data with a cure fraction under correlated destructive negative binomial cure rate models. The models assume multiple latent competing causes with possible dependence and allow for elimination (inactivation) of some initial causes. Estimation is performed via an Expectation-Maximization algorithm, and diagnostic tools based on Cox-Snell residuals are provided. |
| Encoding: | UTF-8 |
| LazyData: | true |
| URL: | https://github.com/rdmatheus/cdnbcr |
| BugReports: | https://github.com/rdmatheus/cdnbcr/issues |
| Suggests: | survival, testthat (≥ 3.0.0) |
| Config/testthat/edition: | 3 |
| RoxygenNote: | 7.3.2 |
| Imports: | Formula (≥ 1.2.4), pracma (≥ 2.3.3) |
| Depends: | R (≥ 3.5) |
| License: | GPL (≥ 3) |
| NeedsCompilation: | no |
| Packaged: | 2026-03-04 13:19:07 UTC; rodri |
| Author: | Rodrigo M. R. de Medeiros
|
| Repository: | CRAN |
| Date/Publication: | 2026-03-09 11:10:02 UTC |
Correlated Destructive Negative Binomial Cure Rate Model
Description
Fits the correlated destructive negative binomial cure rate (CDNBCR) model to right-censored survival data using maximum likelihood estimation via the expectation–maximization (EM) algorithm.
Usage
cdnbcr(
formula,
data,
subset,
na.action,
theta.link = "log",
p.link = "logit",
alpha = TRUE,
control = control_EM(...),
y = TRUE,
x = TRUE,
...
)
## S3 method for class 'cdnbcr'
print(x, digits = getOption("digits"), ...)
Arguments
formula |
A model formula following the syntax of the |
data |
A data frame containing the variables used in the model. |
subset |
An optional vector specifying a subset of observations to be used in the model. |
na.action |
A function indicating how to handle missing values. See |
theta.link, p.link |
Link functions for the regression submodels of |
alpha |
Logical; if |
control |
A list of control parameters for the EM algorithm, as returned by
|
y |
Logical; if |
x |
Logical; if |
... |
Additional arguments passed to |
digits |
a non-null value for digits specifies the minimum number of significant digits to be printed in values. |
Details
The CDNBCR model assumes that the number of initial competing causes follows a negative binomial
distribution with mean \theta > 0 and dispersion parameter \phi > 0. Each initial
cause may remain active with probability p \in (0, 1), and the activation indicators
may be correlated through the parameter \alpha \in [0, 1]. The event time is defined
as the minimum of the latent failure times associated with the remaining active causes.
Individuals with no active causes are considered cured and will never experience the event
of interest. Regression structures can be specified for both the expected number of initial
competing causes (\theta) and the activation probability (p).
The CDNBCR model reduces to the (uncorrelated) destructive negative binomial cure rate model
(Rodrigues et al., 2011) when \alpha = 0. In practice, this is obtained by setting
alpha = FALSE in cdnbcr().
The response must be specified using Surv(time, status), where time is the observed
follow-up time and status is the event indicator (1 = event, 0 = right-censored).
Regression structures can be specified separately for \theta and p using the
| operator. For example,
Surv(time, status) ~ x1 + x2 | z1 + z2
specifies a model in which \theta depends on covariates x1 and x2, while
p depends on covariates z1 and z2.
Value
An object of class "cdnbcr" with components:
- coefficients
A list with the estimated regression coefficients for the
thetaandpsubmodels.- phi, alpha, mu, sigma
Maximum likelihood estimates of the additional model parameters.
- fitted
Fitted values of
thetaandpfor each observation.- links
A named list with the link functions used in the
thetaandpregression submodels.- vcov
Estimated covariance matrix associated with the parameter estimates.
- logLik
Log-likelihood of the fitted model.
- nobs
Number of observations used in the fit.
- df.null
Residual degrees of freedom of the null model.
- df.residual
Residual degrees of freedom of the fitted model.
- convergence
Convergence code of the EM algorithm (0 = successful, 1 = maximum iterations reached).
- inits
Initial values used in the EM algorithm.
- control
Control parameters used in the EM algorithm.
- iterations
Number of EM iterations.
- latent
A list with components
M,D, andYcontaining posterior expectations of the latent variables used in the EM algorithm. Here,Mdenotes the initial number of competing causes,Dis the number of remaining active competing causes, andYis a latent Bernoulli variable that governs the activation regime of the destructive mechanism (and induces dependence throughalpha).- call
Matched function call.
- formula
Model formula.
- terms
Terms objects for the regression submodels.
- y
Response vector (if
y = TRUE).- x
Model matrices (if
x = TRUE).
Author(s)
Diego I. Gallardo diego.gallardo.mateluna@gmail.com
Rodrigo M. R. de Medeiros rodrigo.matheus@ufrn.br
References
Cox, D. R., and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society B, 30, 248–265.
De Medeiros, R. M. R., Bourguignon, M., Gómez, Y. M., and Gallardo, D. I. (2026). A correlated approach to cancer cell counting in cure rate models.
Rodrigues, J., De Castro, M., Balakrishnan, N., and Cancho, V. G. (2011). Destructive weighted Poisson cure rate models. Lifetime Data Analysis, 17, 333–346
See Also
summary.cdnbcr for detailed model summaries,
residuals.cdnbcr to extract Cox–Snell residuals (Cox and Snell, 1968),
plot.cdnbcr for diagnostic plots based on Cox–Snell residuals, and
predict.cdnbcr for prediction under the CDNBCR model.
Additional methods for "cdnbcr" objects are documented in cdnbcr-methods.
Examples
## Loading survival package
library(survival)
## Dataset: e1690 (see ?e1690)
head(e1690)
## Correlated destructive fit
fit <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690)
fit
## Uncorrelated destructive fit
fit0 <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690, alpha = FALSE)
fit0
Extract Information From a CDNBCR Fit
Description
Methods for "cdnbcr" objects.
Usage
## S3 method for class 'cdnbcr'
model.frame(formula, ...)
## S3 method for class 'cdnbcr'
model.matrix(object, parm = c("full", "theta", "p"), ...)
## S3 method for class 'cdnbcr'
coef(object, parm = c("full", "theta", "p"), ...)
## S3 method for class 'cdnbcr'
vcov(object, parm = c("full", "theta", "p"), ...)
## S3 method for class 'cdnbcr'
residuals(object, ...)
## S3 method for class 'cdnbcr'
logLik(object, ...)
## S3 method for class 'cdnbcr'
AIC(object, ..., k = 2)
Arguments
formula |
|
... |
Additional arguments passed to or from other methods. |
object |
An object of class |
parm |
A character indicating which regression structure should be used. It can be
|
k |
Numeric, the penalty per parameter to be used; the default
|
Value
-
model.framereturns adata.framecontaining the variables required byformulaand any additional arguments provided via.... -
model.matrixreturns the design matrix used in the regression structure, as specified by theparmargument. -
coefreturns a numeric vector of estimated regression coefficients, based on theparmargument. Ifparm = "full", it returns a list with the components"theta"and"p", each containing the corresponding coefficient estimates. -
vcovreturns the asymptotic covariance matrix of the regression coefficients, based on theparmargument. -
residualsreturnas a"Surv"object with the Cox-Snell residuals (Cox and Snell, 1968). If the model is well fitted to the data, the Cox-Snell residuals are expected to be distributed as a censored random sample from the exponential distribution with mean 1. -
logLikreturns the log-likelihood value of the fitted model. -
AICreturns a numeric value representing the Akaike Information Criterion (AIC), Bayesian Information Criterion, or another criterion, depending onk.
Author(s)
Rodrigo M. R. de Medeiros <rodrigo.matheus@ufrn.br>
References
Cox, D. R., and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society B, 30, 248–265.
De Medeiros, R. M. R., Bourguignon, M., Gómez, Y. M., and Gallardo, D. I. (2026). A correlated approach to cancer cell counting in cure rate models.
Rodrigues, J., De Castro, M., Balakrishnan, N., and Cancho, V. G. (2011). Destructive weighted Poisson cure rate models. Lifetime Data Analysis, 17, 333–346
Examples
## Loading survival package
library(survival)
## Dataset: e1690 (see ?e1690)
head(e1690)
## Correlated destructive fit
fit <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690)
# Model frame
mf <- model.frame(fit)
mf
# Model matrices
model.matrix(fit, parm = "theta")
model.matrix(fit, parm = "p")
# Coef
coef(fit)
coef(fit, parm = "theta")
coef(fit, parm = "p")
# vcov
vcov(fit)
vcov(fit, parm = "theta")
vcov(fit, parm = "p")
# residuals
residuals(fit)
# Log-likelihood value
logLik(fit)
# AIC and BIC
AIC(fit)
AIC(fit, k = log(fit$nobs))
The Correlated Destructive Negative Binomial Cure Rate Model
Description
Provides the main distributional components of the correlated destructive negative binomial cure rate (CDNBCR) model, including the probability density function, cumulative distribution function, cure fraction, and random number generation.
Usage
dcdnbcr(x, theta, phi, p, alpha, mu, sigma, log.p = FALSE)
pcdnbcr(q, theta, phi, p, alpha, mu, sigma, lower.tail = TRUE, log.p = FALSE)
rcdnbcr(n, theta, phi, p, alpha, mu, sigma)
cure_rate(theta, phi, p, alpha)
Arguments
x |
Vector of positive event times. |
theta, phi |
Positive parameters related to the expected initial number of competing causes and their dispersion. |
p, alpha |
Parameters controlling the probability that a cause remains active. The dependence
among active causes is governed by |
mu, sigma |
Strictly positive parameters of the Weibull distribution assumed for the time-to-event of active competing causes. |
log.p |
Logical; if |
q |
Vector of quantiles. |
lower.tail |
Logical; if |
n |
Number of random values to generate. |
Details
The CDNBCR model describes the time until the occurrence of an event while accounting for a cure fraction, assuming multiple initial competing causes with possible dependence. The model is termed destructive because some of the initial competing causes contributing to the occurrence of the event may be eliminated or inactivated. This formulation is particularly useful in applications where interventions (e.g., treatments) may remove or deactivate some competing causes.
The number of initial competing causes is modeled by a negative binomial random variable
with mean theta and variance theta * (1 + phi * theta). Each initial competing cause
remains active with probability p. Dependence in the activation mechanism is introduced
through the parameter alpha. Let Y_j denote the survival time associated with the
jth remaining competing cause, j = 1, \ldots, D, where D is a latent random variable.
We assume that Y_j follows a Weibull distribution with mean mu and variance
mu^2 * (Gamma(2/sigma + 1) / Gamma(1/sigma + 1)^2 - 1).
The observed uncensored survival time is given by T = \min(Y_1, \ldots, Y_D).
Value
dcdnbcr returns the probability density function, pcdnbcr returns the
distribution function, cure_rate returns the cure fraction of the model, and
rcdnbcr generates random observations.
Author(s)
Diego I. Gallardo <diego.gallardo.mateluna@gmail.com>
Rodrigo M. R. de Medeiros <rodrigo.matheus@ufrn.br>
See Also
cdnbcr for parameter estimation in the CDNBCR model.
Examples
# Parameters
theta <- 1.5
phi <- 1.2
p <- 0.6
alpha <- 0.3
mu <- 5
sigma <- 1.5
n <- 1000
tobs <- rcdnbcr(n, theta, phi, p, alpha, mu, sigma)
# Censoring
Cens <- 10
delta <- ifelse(tobs <= Cens, 1, 0)
tobs[delta == 0] <- Cens
op <- par()$mfrow
par(mfrow = c(1, 2))
# Histogram
hist(tobs, prob = TRUE, main = " ", xlab = expression(T[obs]))
curve(dcdnbcr(x, theta, phi, p, alpha, mu, sigma), add = TRUE, col = 2, lwd = 2)
points(Cens, mean(1 - delta), pch = 16, col = 2, lwd = 2)
# Survival function
plot(survival::survfit(survival::Surv(tobs, delta) ~ 1, se.fit = FALSE),
xlab = expression(T[obs]), ylab = "Survival function")
curve(1 - pcdnbcr(x, theta, phi, p, alpha, mu, sigma), add = TRUE, col = 2, lwd = 2)
# Cure rate
abline(h = cure_rate(theta, phi, p, alpha), col = 4)
legend("bottomleft", c("Empirical", "Theoretical", "Cure rate"),
col = c(1, 2, 4), bty = "n", lty = 1)
par(mfrow = op)
Auxiliary for Controlling a CDNBCR Fit via EM Algorithm
Description
Optimization parameters that control the fitting of the correlated destructive negative binomial cure rate model (CDNBCR) via Expectation-Maximization (EM) algorithm.
Usage
control_EM(method = "BFGS", maxit = 10000, start = NULL, prec = 5e-05, ...)
Arguments
method |
Optimization method for the "M" steps, specifying the |
maxit |
Maximum number of EM algorithm iterations. The default is |
start |
An optional vector of initial values for the algorithm. The expected order of the
parameters is |
prec |
Numeric tolerance for convergence in the EM iterations. |
... |
Additional arguments passed to |
Value
A list with components named as the arguments.
Author(s)
Rodrigo M. R. de Medeiros <rodrigo.matheus@ufrn.br>
See Also
cdnbcr for parameter estimation in the CDNBCR model via EM algorithm.
Examples
## Loading the survival package
library(survival)
## Dataset: (see ?e1690)
head(e1690)
## Correlated destructive fit (default fit)
fit1 <- cdnbcr(Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
trt + thickness, data = e1690)
## Changing the initial values
beta1 <- c(1, 1.5, 2)
phi <- 2
beta2 <- c(0.5, -0.5, 0.5)
alpha <- 0.5
mu <- 3
sigma <- 2
fit2 <- cdnbcr(Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
trt + thickness, data = e1690,
start = c(beta1, phi, beta2, alpha, mu, sigma))
## Compare the fits:
fit1
fit2
Phase III cutaneous melanoma clinical trial
Description
A well-known dataset from a Phase III cutaneous melanoma clinical trial (Ibrahim et al., 2001). The study was conducted by the Eastern Cooperative Oncology Group (ECOG) and aimed to evaluate the effectiveness of the Interferon alpha-2b chemotherapy in preventing recurrence after surgery. After eliminating missing data, the observations include 417 patients from 1991 to 1995, with follow-up until 1998.
Usage
e1690
Format
A data frame with 417 rows and 9 columns:
- trt
Indicates whether the patient was treated with Interferon alpha-2b chemotherapy (chemotherapy) or not (control).
- time
Post-surgery survival or censoring time, in years.
- status
Censoring status.
- age
Age, in years.
- sex
Sex of the patient.
- thickness
Tumor thickness, in mm.
- nodeII, nodeIII, nodeIV
Binary variables indicating the nodal category.
References
Ibrahim, J. G., Chen, M., Sinha, D. (2001). Bayesian Survival Analysis. Springer.
Examples
data(e1690)
plot(time ~ trt, e1690, xlab = "Treatment", ylab = "Time")
Diagnostic Plots for the CDNBCR Model
Description
Produces two diagnostic plots to assess the goodness-of-fit of a Correlated Destructive Negative Binomial Cure Rate model fit based on the Cox-Snell residuals. Available plots include the Kaplan-Meier estimate and the cumulative hazard plot.
Usage
## S3 method for class 'cdnbcr'
plot(
x,
which = 1:2,
ask = prod(graphics::par("mfcol")) < length(which) && grDevices::dev.interactive(),
col.lines = c("black", "#56B4E9"),
lwd = 2,
pch = 16,
cex = 0.8,
...
)
Arguments
x |
An object of class |
which |
Numeric; if a subset of the plots is required, specify a subset
of the numbers |
ask |
Logical; if |
col.lines |
A vector with dimension two with the color for empirical and expected lines. |
pch, cex, lwd, ... |
Graphical parameters (see |
Details
The function produces two diagnostic plots for assessing model fit:
The Kaplan-Meier estimate for Cox-Snell residuals, compared to the expected survival function of an exponential distribution with mean 1;
The cumulative hazard function against the Cox-Snell residuals, which should align approximately with the identity line if the model is well-fitted.
Value
plot method for "cdnbcr" objects returns two types of
diagnostic plots based on the Cox-Snell residuals.
Author(s)
Diego I. Gallardo <diego.gallardo.mateluna@gmail.com>
Rodrigo M. R. de Medeiros <rodrigo.matheus@ufrn.br>
References
Cox, D. R., and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society B, 30, 248–265.
De Medeiros, R. M. R., Bourguignon, M., Gómez, Y. M., and Gallardo, D. I. (2026). A correlated approach to cancer cell counting in cure rate models.
Rodrigues, J., De Castro, M., Balakrishnan, N., and Cancho, V. G. (2011). Destructive weighted Poisson cure rate models. Lifetime Data Analysis, 17, 333–346
Examples
## Loading survival package
library(survival)
## Dataset: e1690 (see ?e1690)
head(e1690)
## Correlated destructive fit
fit <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690)
## Plot of the Cox-Snell residuals
oldpar <- par(mfrow = c(1, 2))
plot(fit, ask = FALSE)
par(oldpar)
Prediction Method for CDNBCR Models
Description
Computes fitted values and predictions from a correlated destructive negative
binomial cure rate (CDNBCR) model fitted with cdnbcr. Predictions
may include the survival function, cure rate, expected number of initial
competing causes, or the activation probability of competing causes.
Usage
## S3 method for class 'cdnbcr'
predict(
object,
newdata = NULL,
type = c("survival", "cure", "theta", "p"),
time,
na.action = stats::na.pass,
...
)
Arguments
object |
An object of class |
newdata |
An optional data frame containing covariate values for which predictions are required. If omitted, predictions are computed for the data used in the model fit. |
type |
Character string indicating the type of prediction. Possible values are:
|
time |
Numeric vector of time points at which the survival function is
evaluated. Only used when |
na.action |
A function specifying how missing values in |
... |
Further arguments passed to or from other methods. |
Value
A numeric vector or matrix of predicted values. The structure of the output
depends on the selected type:
"survival"A matrix of survival probabilities with rows corresponding to observations in
newdata(or the original data) and columns corresponding to the values intime."cure"A numeric vector with the predicted cure fractions.
"theta"A numeric vector with the predicted expected numbers of initial competing causes.
"p"A numeric vector with the predicted activation probabilities.
Author(s)
Diego I. Gallardo <diego.gallardo.mateluna@gmail.com>
Rodrigo M. R. de Medeiros <rodrigo.matheus@ufrn.br>
References
Cox, D. R., and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society, Series B, 30, 248–265.
De Medeiros, R. M. R., Bourguignon, M., Gómez, Y. M., and Gallardo, D. I. (2026). A correlated approach to cancer cell counting in cure rate models.
Rodrigues, J., De Castro, M., Balakrishnan, N., and Cancho, V. G. (2011). Destructive weighted Poisson cure rate models. Lifetime Data Analysis, 17, 333–346.
Examples
## Loading survival package
library(survival)
## Dataset: e1690 (see ?e1690)
head(e1690)
## Correlated destructive fit
fit <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690)
## New data for predictions
newdata <- data.frame(trt = c("Control", "Control", "Chemotherapy", "Chemotherapy"),
age = median(e1690$age),
sex = c("Male", "Female", "Male", "Female"),
thickness = median(e1690$thickness),
nodeII = c(0, 0, 0, 0),
nodeIII = c(0, 0, 0, 0),
nodeIV = c(1, 1, 1, 1))
newdata
## Fitted survival curves
pred <- predict(fit, newdata)
plot(pred[1, ], type = "l", ylim = c(0, 1), xlab = "Time", ylab = "Survival")
lines(pred[2, ], col = 2, lty = 2)
lines(pred[3, ], col = 3, lty = 3)
lines(pred[4, ], col = 4, lty = 4)
legend("topright", legend = c("trt: Control, sex: Male",
"trt: Control, sex: Female",
"trt: Chemotherapy, sex: Male",
"trt: Chemotherapy, sex: Female"),
col = 1:4, lty = 1:4)
## Predicted cure rates
predict(fit, newdata, type = "cure")
## Predicted expected number of initial competing causes
predict(fit, newdata, type = "theta")
## Predicted activation probability of an initial competing cause
predict(fit, newdata, type = "p")
Summarizing a CDNBCR Fit
Description
summary method for class "cdnbcr".
Usage
## S3 method for class 'cdnbcr'
summary(object, ...)
## S3 method for class 'summary.cdnbcr'
print(x, digits = getOption("digits"), ...)
Arguments
object |
An object of class |
... |
Further arguments passed to or from other methods. |
x |
An object of class |
digits |
A non-null value for digits specifies the minimum number of significant digits to be printed in values. |
Value
The function summary.cdnbcr returns an object of class "summary.cdnbcr",
a list with the following components:
- call
The original function call.
- theta
Summary of the regression coefficients for the
thetasubmodel.- p
Summary of the regression coefficients for the
psubmodel.- par
Estimates and standard errors for the additional model parameters:
phi(dispersion of the initial number of competing causes),alpha(dependence parameter, when estimated), and the baseline parametersmuandsigma(reparameterized Weibull).- links
Named list with the link functions used in the
thetaandpregression submodels.- residuals
A
"Surv"object with the Cox-Snell residuals (Cox and Snell, 1968). If the model is well fitted to the data, the Cox-Snell residuals are expected to be distributed as a censored random sample from the exponential distribution with mean 1.- iterations
Number of EM iterations.
- logLik
Log-likelihood of the fitted model.
- AIC, BIC
Akaike and Bayesian information criteria.
Author(s)
Diego I. Gallardo diego.gallardo.mateluna@gmail.com
Rodrigo M. R. de Medeiros rodrigo.matheus@ufrn.br
References
Cox, D. R., and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society B, 30, 248–265.
De Medeiros, R. M. R., Bourguignon, M., Gómez, Y. M., and Gallardo, D. I. (2026). A correlated approach to cancer cell counting in cure rate models.
Rodrigues, J., De Castro, M., Balakrishnan, N., and Cancho, V. G. (2011). Destructive weighted Poisson cure rate models. Lifetime Data Analysis, 17, 333–346
Examples
## Loading survival package
library(survival)
## Dataset: e1690 (see ?e1690)
head(e1690)
## Correlated destructive fit
fit <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690)
out <- summary(fit)
out
### Summary table for the regression coefficients
out$theta
out$p
### Summary table for the additional parameters
out$par
### Cox-Snell residuals
out$residuals
class(out$residuals)
plot(out$residuals)
curve(pexp(x, 1, lower.tail = FALSE), add = TRUE, col = "blue")
## Uncorrelated destructive fit
fit0 <- cdnbcr(formula = Surv(time, status) ~ nodeII + nodeIII + nodeIV - 1 |
sex + trt + thickness + age, data = e1690, alpha = FALSE)
summary(fit0)