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.

Author

Thomas Scheike

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