ate: Compute the average treatment effects via the g-formula

Description Usage Arguments Value Author(s) Examples

View source: R/ate.R

Description

Use the g-formula to estimate the average treatment effect based on Cox regression with or without competing risks

Usage

1
2
3
4
5
ate(object, data, formula, treatment, contrasts = NULL, times, cause,
  landmark, conf.level = 0.95, se = TRUE, band = FALSE, B = 0,
  bootci.method = "perc", nsim.band = ifelse(band, 1000, 0), seed,
  handler = "foreach", mc.cores = 1, verbose = TRUE, store.iid = "full",
  ...)

Arguments

object

outcome model which describes how event risk depends on treatment and covariates. The object carry its own call and have a predictRisk method. See examples.

data

data set in which to evaluate risk predictions based on the outcome model

formula

For analyses with time-dependent covariates, the response formula. See examples.

treatment

name of the treatment variable

contrasts

the levels of the treatment variable to be compared

times

time points at which to evaluate risks

cause

the cause of interest

landmark

for models with time-dependent covariates the landmark time(s) of evaluation. In this case, argument time may only be one value and for the prediction of risks it is assumed that that the covariates do not change between landmark and landmark+time.

conf.level

Numeric value between 0 and 1 (default is 0.05). Confidence level of the confidence intervals.

se

Logical. If TRUE compute standard errors and confidence intervals

band

Logical. If TRUE compute confidence bands across time points.

B

the number of bootstrap replications used to compute the confidence intervals. If it equals 0, then Wald-type confidence intervals are computed. They rely on the standard error estimated using the influence function of the estimator.

bootci.method

Character. Method for constructing bootstrap confidence intervals. Either "perc" (the default), "norm", "basic", "stud", or "bca". Argument passed to boot::boot.ci.

nsim.band

the number of simulations used to compute the quantiles for the confidence bands.

seed

An integer used to generate seeds for bootstrap and to achieve reproducible results.

handler

parallel handler for bootstrap. Either "mclapply" or "foreach". If "foreach" use doParallel to create a cluster.

mc.cores

Passed to parallel::mclapply or doParallel::registerDoParallel. The number of cores to use, i.e. at most how many child processes will be run simultaneously. The option is initialized from environment variable MC_CORES if set.

verbose

Logical. If TRUE inform about estimated run time.

store.iid

Implementation used to estimate the standard error. Can be "full" or "minimal". "minimal" requires less memory but can only estimate the standard for the difference between treatment effects (and not for the ratio).

...

passed to predictRisk

Value

An object of class ate containing:

Author(s)

Brice Ozenne [email protected] and Thomas Alexander Gerds [email protected]

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
library(survival)
library(rms)

set.seed(10)

#### Survival settings  ####
#### ATE with Cox model ####

## generate data
n <- 100
dtS <- sampleData(n, outcome="survival")
dtS$time <- round(dtS$time,1)
dtS$X1 <- factor(rbinom(n, prob = c(0.3,0.4) , size = 2), labels = paste0("T",0:2))

## estimate the Cox model
fit <- cph(formula = Surv(time,event)~ X1+X2,data=dtS,y=TRUE,x=TRUE)

## compute the ATE at times 5, 6, 7, and 8 using X1 as the treatment variable
## Not run: 
## only punctual estimate (argument se = FALSE)
ateFit1a <- ate(fit, data = dtS, treatment = "X1", times = 5:8,
               se = FALSE)

## standard error / confidence intervals computed using the influence function
## (argument se = TRUE and B = 0)
ateFit1b <- ate(fit, data = dtS, treatment = "X1", times = 5:8,
               se = TRUE, B = 0)

## bootstrap confidence intervals: studentized Wald type 
ateFit1c <- ate(fit, data = dtS, treatment = "X1", times = 5,
               seed=3,se = TRUE, B = 100)
## bootstrap confidence intervals: studentized Wald type 
ateFit1d <- ate(fit, data = dtS, treatment = "X1", times = 5,
                seed=3,bootci.method="quantile",se = TRUE, B = 100)

## same as before with in addition the confidence bands for the ATE
## (argument band = TRUE)
ateFit1c <- ate(fit, data = dtS, treatment = "X1", times = 5:8,
               se = TRUE, band = TRUE, B = 0)

## standard error / confidence intervals computed using 100 boostrap samples
## (argument se = TRUE and B = 100) 
ateFit1d <- ate(fit, data = dtS, treatment = "X1",
                times = 5:8, se = TRUE, B = 100)
## NOTE: for real applications 100 bootstrap samples is not enougth 

## same but using 2 cpus for generating and analyzing the boostrap samples
## (parallel computation, argument mc.cores = 2) 
ateFit1e <- ate(fit, data = dtS, treatment = "X1",
                times = 5:8, se = TRUE, B = 100, mc.cores = 2)

## End(Not run)

#### Survival settings without censoring ####
#### ATE with glm                        ####

## generate data
n <- 100
dtS <- sampleData(n, outcome="survival")
dtS[, event5 := eventtime<=5]
dtS[, X2 := as.numeric(X2)]

## estimate the Cox model
fit <- glm(formula = event5 ~ X1+X2, data=dtS, family = "binomial")

## compute the ATE at times 5 using X1 as the treatment variable
## only punctual estimate (argument se = FALSE)
ateFit1a <- ate(fit, data = dtS, treatment = "X1", times = 5,
               se = FALSE)
ateFit1a

## Not run: 
## standard error / confidence intervals computed using the influence function
ateFit1b <- ate(fit, data = dtS, treatment = "X1", times = 5,
               se = TRUE, B = 0)
ateFit1b

## standard error / confidence intervals computed using 100 boostrap samples
ateFit1d <- ate(fit, data = dtS, treatment = "X1",
                times = 5, se = TRUE, B = 100)
ateFit1d

## using lava
ateLava <- estimate(fit, function(p, data){
a <- p["(Intercept)"] ; b <- p["X11"] ; c <- p["X2"] ;
R.X11 <- expit(a + b + c * data[["X2"]])
R.X10 <- expit(a + c * data[["X2"]])
list(risk0=R.X10,risk1=R.X11,riskdiff=R.X10-R.X11)},
average=TRUE)
ateLava

## End(Not run)

#### Competing risks settings               ####
#### ATE with cause specific Cox regression ####

## Not run: 
## generate data
n <- 500
dt <- sampleData(n, outcome="competing.risks")
dt$time <- round(dt$time,1)
dt$X1 <- factor(rbinom(n, prob = c(0.2,0.3,0.2) , size = 3), labels = paste0("T",0:3))

## estimate cause specific Cox model
fitCR <-  CSC(Hist(time,event)~ X1+X8,data=dt,cause=1)

## compute the ATE at times 5, 6, 7, and 8 using X1 as the treatment variable
ateFit2a <- ate(fitCR, data = dt, treatment = "X1", times = 5:8, cause = 1,
               se = FALSE)

## standard error / confidence intervals computed using the influence function
## (argument se = TRUE and B = 0)
ateFit2b <- ate(fitCR, data = dt, treatment = "X1", times = 5:8, cause = 1,
               se = TRUE, B = 0)

## same as before with in addition the confidence bands for the ATE
## (argument band = TRUE)
ateFit2c <- ate(fitCR, data = dt, treatment = "X1", times = 5:8, cause = 1,
               se = TRUE, band = TRUE, B = 0)

## standard error / confidence intervals computed using 100 boostrap samples
## (argument se = TRUE and B = 100) 
ateFit2d <- ate(fitCR, data = dt, treatment = "X1", times = 5:8, cause = 1,
                se = TRUE, B = 100)
## NOTE: for real applications 100 bootstrap samples is not enougth 

## same but using 2 cpus for generating and analyzing the boostrap samples
## (parallel computation, argument mc.cores = 2) 
ateFit2e <- ate(fitCR, data = dt, treatment = "X1", times = 5:8, cause = 1,
                se = TRUE, B = 100, mc.cores = 2)

## End(Not run)

#### time-dependent covariates ###
## Not run: 
library(survival)
fit <- coxph(Surv(time, status) ~ celltype+karno + age + trt, veteran)
vet2 <- survSplit(Surv(time, status) ~., veteran,
                       cut=c(60, 120), episode ="timegroup")
fitTD <- coxph(Surv(tstart, time, status) ~ celltype+karno + age + trt,
               data= vet2,x=1)
set.seed(16)
resVet <- ate(fitTD,formula=Hist(entry=tstart,time=time,event=status)~1,
          data = vet2, treatment = "celltype", contrasts = NULL,
        times=5,verbose=1,
        landmark = c(0,30,60,90), cause = 1, B = 20, se = 1,
        band = FALSE, mc.cores=1)
resVet

## End(Not run)

## Not run: 
set.seed(137)
d=sampleDataTD(127)
library(survival)
d[,status:=1*(event==1)]
## ignore competing risks
cox1TD <- coxph(Surv(start,time, status,type="counting") ~ X3+X5+X6+X8, data=d)
resTD1 <- ate(cox1TD,formula=Hist(entry=start,time=time,event=status)~1,
        data = d, treatment = "X3", contrasts = NULL,
        times=.5,verbose=1,
        landmark = c(0,0.5,1), B = 20, se = 1,
        band = FALSE, mc.cores=1)
resTD1
## adjust for competing risks
cscTD <- CSC(Hist(time=time, event=event,entry=start) ~ X3+X5+X6+X8, data=d)
set.seed(16)
resTD <- ate(cscTD,formula=Hist(entry=start,time=time,event=event)~1,
        data = d, treatment = "X3", contrasts = NULL,
        times=.5,verbose=1,
        landmark = c(0,0.5,1), cause = 1, B = 20, se = 1,
        band = FALSE, mc.cores=1)
resTD

## End(Not run)

riskRegression documentation built on April 19, 2018, 5:04 p.m.