binreg: Binomial Regression for Censored Competing Risks Data

Description

Fits a binomial regression model for a specific time point in the presence of right-censored data and competing risks. This function implements the Inverse Probability of Censoring Weighted (IPCW) estimating equation approach.

Usage

binreg(
  formula,
  data,
  cause = 1,
  time = NULL,
  beta = NULL,
  type = c("II", "I"),
  offset = NULL,
  weights = NULL,
  cens.weights = NULL,
  cens.model = ~+1,
  se = TRUE,
  kaplan.meier = TRUE,
  cens.code = 0,
  no.opt = FALSE,
  method = "nr",
  augmentation = NULL,
  outcome = c("cif", "rmst", "rmtl"),
  model = c("default", "logit", "exp", "lin"),
  Ydirect = NULL,
  ...
)

Value

An object of class "binreg" containing coefficients, variance-covariance matrix, influence functions (iid), and model details.

Arguments

formula: A formula object specifying the outcome and covariates. The outcome must be an Event object (Event(time, cause)).
data: A data frame containing the variables in the formula.
cause: Numeric vector or scalar indicating the cause of interest for the competing risks.
time: Numeric scalar indicating the time point of interest for the cumulative incidence.
beta: Optional numeric vector of starting values for the coefficients. Defaults to zeros.
type: Character string. Either "I" (classic binomial regression) or "II" (adds augmentation term).
offset: Optional numeric vector of offsets for the linear predictor.
weights: Optional numeric vector of weights for the score equations.
cens.weights: Optional numeric vector of pre-calculated censoring weights. If NULL, they are estimated internally.
cens.model: A formula specifying the censoring model. Defaults to ~+1. Can include strata (e.g., ~strata(group)).
se: Logical. If TRUE, computes standard errors based on IPCW influence functions.
kaplan.meier: Logical. If TRUE, uses Kaplan-Meier estimator for IPCW weights; otherwise uses exponential baseline.
cens.code: Numeric code representing censored observations in the status variable. Defaults to 0.
no.opt: Logical. If TRUE, optimization is skipped and starting values are used.
method: Character string. Optimization method: "nr" (Newton-Raphson) or "nlm".
augmentation: Optional numeric vector for additional augmentation terms.
outcome: Character string. Outcome type: "cif" (Cumulative Incidence Function), "rmst" (Restricted Mean Survival Time), or "rmtl" (Restricted Mean Time Lost).
model: Character string. Link function: "default" (auto-selects based on outcome), "logit", "exp", or "lin" (identity).
Ydirect: Optional numeric vector. If provided, this outcome is used instead of constructing one from outcome. Useful for custom IPCW adjustments.
...: Additional arguments passed to lower-level optimization functions.

Outcome definition

The observed outcome $Y_i$ constructed from outcome depends on whether competing risks are present. Competing risks are detected automatically: the data are considered to have competing risks when causes other than those specified in cause are observed (i.e. any(!(Causes %in% cause)), where Causes are all non-censoring event codes found in the data).

No competing risks (all observed events belong to cause):

"cif": $Y_i = I(T_i \leq t, \epsilon_i = 1)$. Cumulative incidence.
"rmst": $Y_i = \min(T_i, t)$. Restricted mean survival time.
"rmtl": $Y_i = t - \min(T_i, t)$. Restricted mean time lost.

Competing risks (causes beyond cause are observed in the data):

"cif": $Y_i = I(T_i \leq t, \epsilon_i \in \code{cause})$. Cumulative incidence for the cause(s) of interest.
"rmst" or "rmtl": $Y_i = I(\epsilon_i \in \code{cause})(t - \min(T_i, t))$. Cause-specific years lost, accumulated only for subjects who experience an event in cause.

The default link function (model = "default") is "logit" for "cif" and "exp" for "rmst" and "rmtl". If Ydirect is supplied, outcome construction is bypassed entirely.

Author

Thomas Scheike

Details

The model estimates the probability: $$P(T \leq t, \epsilon=1 | X ) = \text{expit}( X^T \beta) $$

Based on binomial regresion IPCW response estimating equation: $$ X ( \Delta^{ipcw}(t) I(T \leq t, \epsilon=1 ) - expit( X^T beta)) = 0 $$ where $$\Delta^{ipcw}(t) = I((min(t,T)< C)/G_c(min(t,T)-)$$ is IPCW adjustment of the response $$Y(t)= I(T \leq t, \epsilon=1 )$$. Two types of estimators are available:

type="I": Solves the standard IPCW estimating equation.
type="II": Adds a censoring augmentation term for efficiency gains, solving $$X \int E(Y(t)| T>s)/G_c(s) d \hat M_c$$.

Alternatively, logitIPCW performs a standard logistic regression with IPCW weights applied directly to the likelihood. Thus solving $$ X I(min(T_i,t) < G_i)/G_c(min(T_i ,t)) ( I(T \leq t, \epsilon=1 ) - expit( X^T beta)) = 0 $$ a standard logistic regression with IPCW weights.

The variance estimation is based on squared influence functions, with options for naive variance (assuming known censoring) and robust variance (accounting for censoring model estimation).

Censoring model may depend on strata (cens.model=~strata(gX)).

References

Blanche PF, Holt A, Scheike T (2022). "On logistic regression with right censored data, with or without competing risks, and its use for estimating treatment effects." Lifetime data analysis, 29, 441–482.
Scheike TH, Zhang MJ, Gerds TA (2008). "Predicting cumulative incidence probability by direct binomial regression." Biometrika, 95(1), 205–220.

Examples

Run this code

data(bmt); bmt$time <- bmt$time+runif(408)*0.001
# logistic regresion with IPCW binomial regression 
out <- binreg(Event(time,cause)~tcell+platelet,bmt,time=50)
summary(out)
head(iid(out))

predict(out,data.frame(tcell=c(0,1),platelet=c(1,1)),se=TRUE)

outs <- binreg(Event(time,cause)~tcell+platelet,bmt,time=50,cens.model=~strata(tcell,platelet))
summary(outs)

## glm with IPCW weights 
outl <- logitIPCW(Event(time,cause)~tcell+platelet,bmt,time=50)
summary(outl)

##########################################
### risk-ratio of different causes #######
##########################################
data(bmt)
bmt$id <- 1:nrow(bmt)
bmt$status <- bmt$cause
bmt$strata <- 1
bmtdob <- bmt
bmtdob$strata <-2
bmtdob <- dtransform(bmtdob,status=1,cause==2)
bmtdob <- dtransform(bmtdob,status=2,cause==1)
###
bmtdob <- rbind(bmt,bmtdob)
dtable(bmtdob,cause+status~strata)

cif1 <- cif(Event(time,cause)~+1,bmt,cause=1)
cif2 <- cif(Event(time,cause)~+1,bmt,cause=2)
plot(cif1)
plot(cif2,add=TRUE,col=2)

cifs1 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=1,time=50)
cifs2 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=2,time=50)
summary(cifs1)
summary(cifs2)

cifdob <- binreg(Event(time,status)~-1+factor(strata)+
	 tcell*factor(strata)+platelet*factor(strata)+age*factor(strata)
	 +cluster(id),bmtdob,cause=1,time=50,cens.model=~strata(strata))
summary(cifdob)
head(iid(cifdob)) 

newdata <- data.frame(tcell=1,platelet=1,age=0,strata=1:2,id=1)
riskratio <- function(p) {
  cifdob$coef <- p
  p <- predict(cifdob,newdata,se=0)
  return(p[1]/p[2])
}
lava::estimate(cifdob,f=riskratio)

predict(cifdob,newdata)
(p1 <- predict(cifs1,newdata))
(p2 <- predict(cifs2,newdata))
p1[1,1]/p2[1,1]

Run the code above in your browser using DataLab