blendR-paper/Scripts/CODE_case_study.R at main · StatisticsHealthEconomics/blendR-paper · GitHub

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
## ----install packages------------------------------------------------------------
if (.Platform$OS.type == "windows") {
  if ("survHE" %in% rownames(installed.packages()) == FALSE) {
    install.packages("survHE",
 	                   repos = c("http://www.statistica.it/gianluca/R",
 		                         "https://cran.rstudio.org",
                              "https://inla.r-inla-download.org/R/stable"),
 	                   dependencies = TRUE)
   }
}

if (.Platform$OS.type == "windows") {
  if ("INLA" %in% rownames(installed.packages()) == FALSE) {
    install.packages("INLA",
                     repos = c(getOption("repos"),
                     INLA = "https://inla.r-inla-download.org/R/stable"),
                     dep = TRUE)
   }
}


## ----load packages, results='hide', warning=FALSE, message=FALSE----------------
library(survHE)
library(INLA)


## ----load data------------------------------------------------------------------
load("../Data/ta174_FCR.Rdata")
head(dat_FCR)


## ----source function, include=TRUE----------------------------------------------
source("../Scripts/Functions.R", local = knitr::knit_global())


## ----observed estimate----------------------------------------------------------
obs_Surv <- surv_est_inla(data = dat_FCR,
                          cutpoints = seq(0,180, by =5))

names(obs_Surv)


## ----observedplot, echo=TRUE, fig.cap="Survival estimates of the FCR arm based on the piecewise exponential model fitted to the observed data"----
ggplot() + geom_line(aes(obs_Surv$KM$time, obs_Surv$KM$surv, colour = "Kaplan-Meier"), size = 1.25, linetype = "dashed")+
  xlim(0, 180) + ylim(0,1) +
  geom_line(aes(obs_Surv$time,rowMeans(obs_Surv$S_obs), colour = "Data fitting"), size = 1)+
  geom_ribbon(aes(x=obs_Surv$time, y = rowMeans(obs_Surv$S_obs),
                  ymin = apply(obs_Surv$S_obs, 1, quantile, probs = 0.025),ymax = apply(obs_Surv$S_obs, 1, quantile, probs = 0.975)), alpha = 0.1)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_line(colour = "black"), text = element_text(size = 8))+
  xlab("Time (months)") + ylab("Survival")+
  scale_colour_manual(name = "model", values = c("Data fitting"="#7CAE00", "Kaplan-Meier"="brown"))+
  annotate("text", x = 150, y = 0.9, label = "FCR arm", colour = "orange", fontface = 2)


## ----external estimate, message=FALSE, warning=FALSE, results='hide'------------
ext_Surv <- ext_surv_est(t_info = 144, S_info = 0.05,
                         T_max = 180,
                         times_est = seq(0, 180), distr = "gom")

names(ext_Surv)


## ----external_object, echo=FALSE------------------------------------------------
names(ext_Surv)


## ----weight function------------------------------------------------------------
blending_interval <- list(min = 48, max = 150)
beta_params <- list(alpha = 3, beta = 3)

tp <- seq(0, 180)

# parameters for the weight function
wt_par <- list(a = blending_interval$min, b = blending_interval$max,
               shape1 = beta_params$alpha, shape2 = beta_params$beta)

weight <- with(wt_par,
               pbeta((tp - a)/(b - a), shape1, shape2))


## ----blended estimate-----------------------------------------------------------
n_sim <- ncol(ext_Surv$S_ext)

ble_Surv <- matrix(NA, nrow = 180 + 1, ncol = n_sim)

for (i in seq_len(n_sim)) {
  ble_Surv[, i] <-
    obs_Surv$S_obs[, i]^(1 - weight) * ext_Surv$S_ext[, i]^weight
}


## ----blendedplot, echo=TRUE, fig.cap="Blended survival curve based on short-term data and external information for FCR arm"----
ggplot() + geom_line(aes(obs_Surv$KM$time, obs_Surv$KM$surv, colour = "Kaplan-Meier"), size = 1.25, linetype = "dashed")+
  xlim(0, 180) + ylim(0,1) +
  geom_line(aes(tp,rowMeans(obs_Surv$S_obs), colour = "Data fitting"), size = 1, linetype = "twodash")+
  geom_ribbon(aes(x=tp, y = rowMeans(obs_Surv$S_obs),
                  ymin = apply(obs_Surv$S_obs, 1, quantile, probs = 0.025),
                  ymax = apply(obs_Surv$S_obs, 1, quantile, probs = 0.975)), alpha = 0.1)+
  geom_line(aes(tp, rowMeans(ext_Surv$S_ext), colour = "External info"), size = 1, linetype = "longdash")+
  geom_ribbon(aes(x=tp, y = rowMeans(ext_Surv$S_ext), ymin = apply(ext_Surv$S_ext, 1, quantile, probs = 0.025),
                  ymax = apply(ext_Surv$S_ext, 1, quantile, probs = 0.975)), alpha = 0.1)+
  geom_line(aes(tp, rowMeans(ble_Surv), colour = "Blended curve"), size = 1.25)+
  geom_ribbon(aes(x=tp, y = rowMeans(ble_Surv), ymin = apply(ble_Surv, 1, quantile, probs = 0.025),
                  ymax = apply(ble_Surv, 1, quantile, probs = 0.975)), alpha = 0.1)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_line(colour = "black"), text = element_text(size = 8))+
  xlab("Time (months)") + ylab("Survival")+
  scale_colour_manual(name = "model", values = c("Data fitting"="#7CAE00", "External info"="#00BFC4", "Blended curve"="#F8766D", "Kaplan-Meier"="brown"))+
  annotate("text", x = 150, y = 0.9, label = "FCR arm", colour = "orange", fontface = 2)