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) $$
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
coxph)- 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 I(T \leq t, \epsilon=1 )/G_c(T_i-) - expit( X^T beta)) = 0 $$ for IPCW adjusted responses, with (default, type="II") an additional censoring augmentation term $$X \int E(Y(t)| T>s)/G_c(s) dM_c$$ with $$Y(t) = I(T \leq t, \epsilon=1 )$$
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.
Variance is based on $$ \sum w_i^2 $$ with IPCW adjustment, and naive.var is variance under known censoring model.
Censoring model may depend on strata.
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)
#>
#> n events
#> 408 160
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.180389 0.126750 -0.428815 0.068038 0.1547
#> tcell -0.419079 0.345487 -1.096221 0.258063 0.2251
#> platelet -0.436963 0.240975 -0.909265 0.035338 0.0698
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.83495 0.65128 1.0704
#> tcell 0.65765 0.33413 1.2944
#> platelet 0.64600 0.40282 1.0360
#>
#>
predict(out,data.frame(tcell=c(0,1),platelet=c(1,1)),se=TRUE)
#> pred se lower upper
#> 1 0.3503839 0.04848565 0.2553521 0.4454158
#> 2 0.2618392 0.06969476 0.1252374 0.3984409
##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)
#>
#> n events
#> 408 160
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.241503 0.131457 -0.499154 0.016148 0.0662
#> tcell -0.344892 0.368352 -1.066848 0.377064 0.3491
#> platelet -0.292631 0.262667 -0.807450 0.222187 0.2652
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.78545 0.60704 1.0163
#> tcell 0.70830 0.34409 1.4580
#> platelet 0.74630 0.44599 1.2488
#>
#>
##########################################
### 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)
#> strata: 1
#>
#> status 0 1 2
#> cause
#> 0 160 0 0
#> 1 0 161 0
#> 2 0 0 87
#> ------------------------------------------------------------
#> strata: 2
#>
#> status 0 1 2
#> cause
#> 0 160 0 0
#> 1 0 0 161
#> 2 0 87 0
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)
#>
#> n events
#> 408 85
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -1.321992 0.157775 -1.631225 -1.012759 0.0000
#> tcell 0.747597 0.352327 0.057049 1.438146 0.0338
#> platelet -0.019491 0.277007 -0.562415 0.523432 0.9439
#> age -0.072305 0.141684 -0.350001 0.205391 0.6098
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.26660 0.19569 0.3632
#> tcell 2.11192 1.05871 4.2129
#> platelet 0.98070 0.56983 1.6878
#> age 0.93025 0.70469 1.2280
#>
#>
summary(cifs2)
#>
#> n events
#> 408 85
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -1.321992 0.157775 -1.631225 -1.012759 0.0000
#> tcell 0.747597 0.352327 0.057049 1.438146 0.0338
#> platelet -0.019491 0.277007 -0.562415 0.523432 0.9439
#> age -0.072305 0.141684 -0.350001 0.205391 0.6098
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.26660 0.19569 0.3632
#> tcell 2.11192 1.05871 4.2129
#> platelet 0.98070 0.56983 1.6878
#> age 0.93025 0.70469 1.2280
#>
#>
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)
#>
#> n events
#> 816 245
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> factor(strata)1 -0.199015 0.130984 -0.455739 0.057709 0.1287
#> factor(strata)2 -1.321992 0.157775 -1.631225 -1.012759 0.0000
#> tcell -0.637921 0.356676 -1.336993 0.061151 0.0737
#> platelet -0.344154 0.246034 -0.826371 0.138063 0.1619
#> age 0.437363 0.107271 0.227116 0.647610 0.0000
#> factor(strata)2:tcell 1.385518 0.601030 0.207521 2.563516 0.0212
#> factor(strata)2:platelet 0.324662 0.432059 -0.522158 1.171483 0.4524
#> factor(strata)2:age -0.509668 0.208101 -0.917538 -0.101797 0.0143
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> factor(strata)1 0.81954 0.63398 1.0594
#> factor(strata)2 0.26660 0.19569 0.3632
#> tcell 0.52839 0.26263 1.0631
#> platelet 0.70882 0.43763 1.1480
#> age 1.54862 1.25498 1.9110
#> factor(strata)2:tcell 3.99690 1.23062 12.9814
#> factor(strata)2:platelet 1.38356 0.59324 3.2268
#> factor(strata)2:age 0.60070 0.39950 0.9032
#>
#>
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)
#> Estimate Std.Err 2.5% 97.5% P-value
#> p1 0.6602 0.2737 0.1237 1.197 0.01586