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,
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,
...
)Arguments
- formula
formula with outcome (see
coxph)- data
data frame
- cause
cause of interest (numeric variable)
- time
time of interest
- beta
starting values
- 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
- ...
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.
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 $$ also with IPCW adjustment, and naive.var is variance under known censoring model.
Censoring model may depend on strata.
Examples
data(bmt)
# 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.178145 0.127205 -0.427462 0.071173 0.1614
#> tcell -0.426717 0.348029 -1.108841 0.255407 0.2202
#> platelet -0.442046 0.242375 -0.917092 0.033001 0.0682
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.83682 0.65216 1.0738
#> tcell 0.65265 0.32994 1.2910
#> platelet 0.64272 0.39968 1.0336
#>
#>
predict(out,data.frame(tcell=c(0,1),platelet=c(1,1)),se=TRUE)
#> pred se lower upper
#> 1 0.3497382 0.04857251 0.2545361 0.4449403
#> 2 0.2598194 0.07038970 0.1218556 0.3977833
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.180709 0.127413 -0.430434 0.069016 0.1561
#> tcell -0.365862 0.350626 -1.053077 0.321353 0.2967
#> platelet -0.433605 0.240229 -0.904446 0.037236 0.0711
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.83468 0.65023 1.0715
#> tcell 0.69360 0.34886 1.3790
#> platelet 0.64817 0.40477 1.0379
#>
#>
## 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.241415 0.131460 -0.499072 0.016241 0.0663
#> tcell -0.344384 0.368385 -1.066404 0.377637 0.3499
#> platelet -0.292936 0.262664 -0.807748 0.221875 0.2647
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.78552 0.60709 1.0164
#> tcell 0.70866 0.34424 1.4588
#> platelet 0.74607 0.44586 1.2484
#>
#>
##########################################
### 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=1,time=50)
cifs2 <- binreg(Event(time,cause)~tcell+platelet+age,bmt,cause=2,time=50)
summary(cifs1)
#>
#> n events
#> 408 160
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -0.195401 0.131030 -0.452216 0.061413 0.1359
#> tcell -0.637058 0.359000 -1.340685 0.066569 0.0760
#> platelet -0.352264 0.247108 -0.836588 0.132059 0.1540
#> age 0.419772 0.106900 0.210253 0.629292 0.0001
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.82250 0.63622 1.0633
#> tcell 0.52885 0.26167 1.0688
#> platelet 0.70309 0.43319 1.1412
#> age 1.52161 1.23399 1.8763
#>
#>
summary(cifs2)
#>
#> n events
#> 408 85
#>
#> 408 clusters
#> coeffients:
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) -1.318339 0.158459 -1.628914 -1.007764 0.0000
#> tcell 0.742321 0.359124 0.038451 1.446191 0.0387
#> platelet -0.030930 0.279381 -0.578507 0.516647 0.9118
#> age -0.093271 0.142475 -0.372518 0.185976 0.5127
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> (Intercept) 0.26758 0.19614 0.3650
#> tcell 2.10081 1.03920 4.2469
#> platelet 0.96954 0.56073 1.6764
#> age 0.91095 0.68900 1.2044
#>
#>
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.195401 0.131030 -0.452216 0.061413 0.1359
#> factor(strata)2 -1.318339 0.158459 -1.628914 -1.007764 0.0000
#> tcell -0.637058 0.359000 -1.340685 0.066569 0.0760
#> platelet -0.352264 0.247108 -0.836588 0.132059 0.1540
#> age 0.419772 0.106900 0.210253 0.629292 0.0001
#> factor(strata)2:tcell 1.379379 0.602993 0.197534 2.561224 0.0222
#> factor(strata)2:platelet 0.321334 0.431926 -0.525225 1.167893 0.4569
#> factor(strata)2:age -0.513043 0.207027 -0.918808 -0.107278 0.0132
#>
#> exp(coeffients):
#> Estimate 2.5% 97.5%
#> factor(strata)1 0.82250 0.63622 1.0633
#> factor(strata)2 0.26758 0.19614 0.3650
#> tcell 0.52885 0.26167 1.0688
#> platelet 0.70309 0.43319 1.1412
#> age 1.52161 1.23399 1.8763
#> factor(strata)2:tcell 3.97244 1.21840 12.9517
#> factor(strata)2:platelet 1.37897 0.59142 3.2152
#> factor(strata)2:age 0.59867 0.39899 0.8983
#>
#>
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.6639 0.2764 0.1221 1.206 0.01631