binreg {mets} | R Documentation |
Binomial Regression for censored competing risks data
Description
Simple version of comp.risk function of timereg for just one time-point thus fitting the model
P(T \leq t, \epsilon=1 | X ) = expit( X^T beta)
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"),
model = "exp",
Ydirect = NULL,
...
)
Arguments
formula |
formula with outcome (see |
data |
data frame |
cause |
cause of interest (numeric variable) |
time |
time of interest |
beta |
starting values |
type |
"II" adds augmentation term, and "I" classic binomial regression |
offset |
offsets for partial likelihood |
weights |
for score equations |
cens.weights |
censoring weights |
cens.model |
only stratified cox model without covariates |
se |
to compute se's based on IPCW |
kaplan.meier |
uses Kaplan-Meier for IPCW in contrast to exp(-Baseline) |
cens.code |
gives censoring code |
no.opt |
to not optimize |
method |
for optimization |
augmentation |
to augment binomial regression |
outcome |
can do CIF regression "cif"=F(t|X), "rmst"=E( min(T, t) | X) , or "rmst-cause"=E( I(epsilon==cause) ( t - mint(T,t)) ) | X) |
model |
possible exp model for E( min(T, t) | X)=exp(X^t beta) , or E( I(epsilon==cause) ( t - mint(T,t)) ) | X)=exp(X^t beta) |
Ydirect |
use this Y instead of outcome constructed inside the program (e.g. I(T< t, epsilon=1)), then uses IPCW vesion of the Y, set outcome to "rmst" to fit using the model specified by model |
... |
Additional arguments to lower level funtions |
Details
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 )
.
(type="I") sovlves this estimating equation using a stratified Kaplan-Meier for the censoring distribution. For (type="II") the default an additional censoring augmentation term
X \int E(Y(t)| T>s)/G_c(s) d \hat M_c
is added.
logitIPCW instead considers
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 weights that adjust for IPCW.
The variance is based on the squared influence functions that are also returned as the iid component. naive.var is variance under known censoring model.
Censoring model may depend on strata (cens.model=~strata(gX)).
Author(s)
Thomas Scheike
Examples
library(mets)
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)
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)
bplot(cif1)
bplot(cif2,add=TRUE,col=2)
cifs1 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=2,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)+cluster(id))
summary(cifdob)
riskratio <- function(p) {
Z <- rbind(c(1,0,1,1,0,0,0,0), c(0,1,1,1,0,1,1,0))
lp <- c(Z %*% p)
p <- lava::expit(lp)
return(p[1]/p[2])
}
lava::estimate(cifdob,f=riskratio)