Package 'cdnbcr'

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 Formula package. The left-hand side must be a Surv object specifying the observed follow-up time and censoring indicator. The right-hand side may contain two regression components separated by the | operator: the first for theta (expected number of initial competing causes) and the second for p (activation probability of an initial competing cause). If only one component is provided, it is used for theta and p is modeled with an intercept only.
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 model.frame.
theta.link, p.link	Link functions for the regression submodels of theta and p. Defaults are "log" for theta.link and "logit" for p.link. Alternative links can be specified using make.link.
alpha	Logical; if TRUE (default), the correlation parameter alpha is estimated. If FALSE, alpha is fixed at zero, yielding the uncorrelated destructive model.
control	A list of control parameters for the EM algorithm, as returned by control_EM.
y	Logical; if TRUE (default), the response vector is returned.
x	Logical; if TRUE (default), the model matrices are returned. For print(), x is a fitted model object of class "cdnbcr".
...	Additional arguments passed to control_EM.
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 theta and p submodels.

phi, alpha, mu, sigma

Maximum likelihood estimates of the additional model parameters.

fitted

Fitted values of theta and p for each observation.

links

A named list with the link functions used in the theta and p regression 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, and Y containing posterior expectations of the latent variables used in the EM algorithm. Here, M denotes the initial number of competing causes, D is the number of remaining active competing causes, and Y is a latent Bernoulli variable that governs the activation regime of the destructive mechanism (and induces dependence through alpha).

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 [email protected]

Rodrigo M. R. de Medeiros [email protected]

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	A model Formula or terms object or an "cdnbcr" object.
...	Additional arguments passed to or from other methods.
object	An object of class "cdnbcr".
parm	A character indicating which regression structure should be used. It can be "theta" for the expected initial competing causes regression structure, "p" for the activation probability of an initial competing cause regression submodel, or "full" for both regression structures.
k	Numeric, the penalty per parameter to be used; the default k = 2 is the classical AIC. See AIC.

Value

• model.frame returns a data.frame containing the variables required by formula and any additional arguments provided via ....

• model.matrix returns the design matrix used in the regression structure, as specified by the parm argument.

• coef returns a numeric vector of estimated regression coefficients, based on the parm argument. If parm = "full", it returns a list with the components "theta" and "p", each containing the corresponding coefficient estimates.

• vcov returns the asymptotic covariance matrix of the regression coefficients, based on the parm argument.

• residuals returnas 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.

• logLik returns the log-likelihood value of the fitted model.

• AIC returns a numeric value representing the Akaike Information Criterion (AIC), Bayesian Information Criterion, or another criterion, depending on k.

Author(s)

Rodrigo M. R. de Medeiros <[email protected]>

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 alpha. Both p and alpha are restricted to the closed unit interval [0, 1]. When alpha = 0, activations are assumed to be independent.
mu, sigma	Strictly positive parameters of the Weibull distribution assumed for the time-to-event of active competing causes.
log.p	Logical; if TRUE, probabilities p are returned as log(p).
q	Vector of quantiles.
lower.tail	Logical; if TRUE (default), probabilities are P(X <= x); otherwise, P(X > x).
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 $j$th 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 <[email protected]>

Rodrigo M. R. de Medeiros <[email protected]>

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 method argument passed to optim.
maxit	Maximum number of EM algorithm iterations. The default is 10.000 iterations. If this limit is reached, a warning message is displayed.
start	An optional vector of initial values for the algorithm. The expected order of the parameters is (beta1, phi, beta2, alpha, mu, sigma), where beta1 and beta2 are vectors of coefficients associated with the regression structures of theta and p, respectively.
prec	Numeric tolerance for convergence in the EM iterations.
...	Additional arguments passed to optim.

Value

A list with components named as the arguments.

Author(s)

Rodrigo M. R. de Medeiros <[email protected]>

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 "cdnbcr", a result of a call to cdnbcr.
which	Numeric; if a subset of the plots is required, specify a subset of the numbers 1:2.
ask	Logical; if TRUE, the user is asked before each plot.
col.lines	A vector with dimension two with the color for empirical and expected lines.
pch, cex, lwd, ...	Graphical parameters (see par).

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 <[email protected]>

Rodrigo M. R. de Medeiros <[email protected]>

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 "cdnbcr", as returned by cdnbcr.
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: "survival" Predicted survival function evaluated at the values supplied in time. "cure" Predicted cure fraction (probability of being cured). "theta" Predicted expected number of initial competing causes. "p" Predicted activation probability of an initial competing cause.
time	Numeric vector of time points at which the survival function is evaluated. Only used when type = "survival".
na.action	A function specifying how missing values in newdata should be handled. The default is to return NA for predictions involving incomplete cases.
...	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 in time.

"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 <[email protected]>

Rodrigo M. R. de Medeiros <[email protected]>

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 "cdnbcr", a result of a call to cdnbcr.
...	Further arguments passed to or from other methods.
x	An object of class "summary.cdnbcr", a result of a call to summary.cdnbcr.
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 theta submodel.

p

Summary of the regression coefficients for the p submodel.

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 parameters mu and sigma (reparameterized Weibull).

links

Named list with the link functions used in the theta and p regression 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 [email protected]

Rodrigo M. R. de Medeiros [email protected]

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)

---