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,
...
)formula with outcome (see coxph)
data frame
cause of interest (numeric variable)
time of interest
starting values
"II" adds augmentation term, and "I" classic binomial regression
offsets for partial likelihood
for score equations
censoring weights
only stratified cox model without covariates
to compute se's based on IPCW
uses Kaplan-Meier for IPCW in contrast to exp(-Baseline)
gives censoring code
to not optimize
for optimization
to augment binomial regression
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)
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)
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
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)).
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)
#>
#> n events
#> 408 160
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.180650 0.127412 -0.430373 0.069073 0.1562
#> tcell -0.366536 0.350628 -1.053754 0.320682 0.2959
#> platelet -0.432322 0.240356 -0.903411 0.038767 0.0721
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.83473 0.65027 1.0715
#> tcell 0.69313 0.34863 1.3781
#> platelet 0.64900 0.40519 1.0395
#>
#>
## 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)
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet
bplot(cif2,add=TRUE,col=2)
#> Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet
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