-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathREADME.Rmd
108 lines (83 loc) · 3.73 KB
/
README.Rmd
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
---
title: Targeted estimation of heterogeneous treatment effect in observational survival
analysis
output:
md_document:
variant: markdown_github
html_document:
df_print: paged
---
[](http://www.repostatus.org/#active)
[](http://opensource.org/licenses/MIT)
[](https://doi.org/10.1016/j.jbi.2020.103474)
This paper proposed a three-stage modular design for estimating the treatment effect heterogeneity in observational survival analysis. The method provides monotonic survival curves with adjustment for selection and censoring bias. We automate the identification of features contributing to the effect heterogeneity. We avoid the ad-hoc subgroup selection by non-parametrically estimating the conditional treatment effect. We provide evidence that oral anticoagulants confer protection against stroke and death on newly diagnosed non-valvular atrial fibrillation patients.
```{r message=FALSE, warning=FALSE, paged.print=FALSE}
require(dplyr)
require(MOSS) #devtools::install_github('wilsoncai1992/MOSS')
require(survival)
#require(simcausal) #if you want complicate simulation, please install from local directory install.packages("~/simcausal_0.5.5.tar", repos = NULL)
require(abind)
require(tidyverse)
```
Generate some sample
```{r}
n <- 100
W <- data.frame(W1 = runif(n), W2 = rbinom(n, 1, 0.5))
A <- rbinom(n, 1, 0.5)
EventTime <- rgeom(n,plogis(-4 + W$W1 * W$W2 - A)) + 1
CensorTime <- rgeom(n, plogis(-6 + W$W1)) + 1
T.tilde <- pmin(EventTime, CensorTime)
Delta <- as.numeric(T.tilde == EventTime)
df = data.frame(A = A, T.tilde = T.tilde, Delta = Delta, W1 = W$W1, W2=W$W2)
df$ID <- seq.int(nrow(df))
max_time = 30
head(df)
```
Run the method proposed in the step 1 and 2 of paper
```{r message=FALSE, warning=FALSE, paged.print=FALSE}
df <- df[df$T.tilde<= max_time & df$T.tilde>0,]
df <- df[complete.cases(df),]
adjustVars <- grep('W', colnames(df), value = T)
sl_lib_g <- c( "SL.earth","SL.gam") #choose your own esemble algorithm here
sl_lib_censor <- c( "SL.earth","SL.gam")
sl_lib_failure <- c( "SL.earth","SL.gam")
#df$T.tilde <- df$T.tilde + 1
k_grid <- 1:max(df$T.tilde)
#SL
sl_fit <- initial_sl_fit(
T_tilde = df$T.tilde,
Delta = df$Delta,
A = df$A,
W = data.frame(df[, adjustVars]),
#adjustVars = df[,c('W','W1')],
t_max = max(df$T.tilde),
sl_treatment = sl_lib_g,
sl_censoring = sl_lib_censor,
sl_failure = sl_lib_failure
)
sl_fit$density_failure_1$hazard_to_survival()
sl_fit$density_failure_0$hazard_to_survival()
sl_fit$density_failure_1$t <- k_grid
sl_fit$density_failure_0$t <- k_grid
sl_density_failure_1_marginal <- sl_fit$density_failure_1$clone(deep = TRUE)
sl_density_failure_0_marginal <- sl_fit$density_failure_0$clone(deep = TRUE)
sl_density_failure_1_marginal$survival <- matrix(colMeans(sl_density_failure_1_marginal$survival), nrow = 1)
sl_density_failure_0_marginal$survival <- matrix(colMeans(sl_density_failure_0_marginal$survival), nrow = 1)
out <- list(sl_fit_1 = sl_fit$density_failure_1$survival,
sl_fit_0 = sl_fit$density_failure_0$survival,
SL_diff = sl_fit$density_failure_1$survival-sl_fit$density_failure_0$survival)
```
Individual difference in survival probabilities
```{r}
head(out$SL_diff)
hist(out$SL_diff)
```
The average on the ITEs will be the TMLE adjusted ATE.
```{r}
#ATE
mean(out$SL_diff)
```
Please refer to the paper for BART variable importance measure and kernal grouping to get the CATE for subgroups.