library(iclust2prog)
library(predsurv)
library(magrittr)
library(dplyr)
library(forcats)
library(tidyr)
library(purrr)
library(modelr)
library(tidybayes)
library(ggplot2)
library(ggstance)
library(ggridges)
library(rstan)
library(rstanarm)
library(brms)
library(brmstools)
import::from(LaplacesDemon, invlogit)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
devtools::document()
data("ic2surv")
?ic2surv
ic2surv <- ic2surv %>%
dplyr::mutate(mastectomy = breast_surgery == "MASTECTOMY",
crtherapy = (chemotherapy == "YES" | radio_therapy == "YES"),
chemotherapy = chemotherapy == "YES",
radio_therapy = radio_therapy == "YES",
hormone_therapy = hormone_therapy == "YES")
#prepare long format dataset
long_ic2dat <- gen_stan_dat(ic2surv)
m1.map2stan <- post_surv(x = ic2surv, surv_form = c("hormone_therapy", "mastectomy") )
summary(m1.map2stan)
m1.gam <- mgcv::gam(status ~ 1 + offset(log_dtime) + s(time), long_ic2dat, family='poisson')
m1.stan_gam <- rstanarm::stan_gamm4(status~1+offset(log_dtime)+s(time) + chemotherapy, data = long_ic2dat, family='poisson')
prior <- get_prior(m1.stan_gam)
m1.brms_gam <- brms::brm(bf(status ~ 1 + offset(log_dtime) + s(time) + crtherapy),
data = long_ic2dat, family = poisson(), cores = 4, seed = 17,
iter = 4000, warmup = 1000, thin = 10, refresh = 0,
control = list(adapt_delta = 0.99))
the_stan_model <- make_stancode(bf(status ~ 1 + offset(log_dtime) + s(time) + crtherapy),
data = long_ic2dat, family = poisson())
(prior <- get_prior(brm(bf(status ~ 1 + offset(log_dtime) + s(time) + factor(chemotherapy)),
data = long_ic2dat, family = poisson())))
summary(m1.brms_gam)
#plot predictions
timepoints <- seq(0, max(ic2surv$time),
length.out = 100L)
newdata <- gen_new.frame(ic2surv, timepoints = timepoints)
chemdata <- rbind(newdata %>% mutate(chemotherapy = T) , newdata %>% mutate(chemotherapy = F,
sample_id = sample_id+nrow(ic2surv)) ) %>%
select(-status)
post <- posterior_linpred(m1.stan_gam, newdata = chemdata)
plot.matrix <- link.surv(post = post, longdata = chemdata)
surv.mean <- as.data.frame( t( apply(plot.matrix, 1, mean) ), colnames = surv.mean)
surv.pi <- as.data.frame( t( apply(plot.matrix, 1, quantile, probs = c(0.055, 0.945)) ), colnames = c( "lower", "upper") )
plot.frame <- get_plot.frame(post = plot.matrix, strata = c("chemotherapy"), obs = rbind(ic2surv, ic2surv))
obs.mortality <- get_km.frame(ic2surv, strata = c("chemotherapy") )
p <- ggplot2::ggplot(plot.frame ,
aes(time, postmean, group = chemotherapy, colour = chemotherapy)) +
geom_line(mapping=aes(group = chemotherapy, colour = chemotherapy), alpha = 0.8) +
geom_ribbon(aes( ymin = lower,
ymax = upper,
fill = chemotherapy), alpha = 0.5, colour=NA) +
geom_step(data=obs.mortality , aes(time,
surv,
colour = chemotherapy), alpha = 0.9)+
labs(title="Bayesian Modeling of survival outcome for iC-2",
subtitle="Hormone therapy observed survival and predictive posterior distribution with 89% credible interval",
caption="4000 posterior sample smooths from GAMM model shown") +
ylab("Survival probability") +
xlab("Time (months)") +
# guides() +
theme_bw() + theme(legend.position=c(0.75, 0.75))
p
pdf('Hormone.pdf')
p
dev.off()
#Conterfactual plot
ic2surv <- ic2surv %>%
mutate(Treatment = ifelse(chemotherapy & !hormone_therapy,
"only_chemo", ifelse(hormone_therapy & !chemotherapy, "only_hormone", ifelse(hormone_therapy & chemotherapy, "both", "none"))))
m2.stan_gam <- rstanarm::stan_gamm4(status~1+offset(log_dtime)+s(time) + mastectomy + hormone_therapy, data = long_ic2dat, family='poisson', seed = 17,iter = 4000, warmup = 1000, thin = 10, control = list(adapt_delta = 0.99))
summary(m2.stan_gam)
timepoints = seq(min(ic2surv$time), max(ic2surv$time), length.out = 100)
dtime = c(0, diff(timepoints))
mastectomy = c( rep(T, 100), rep(F, 100))
hormone_therapy = c( rep(F, 100), rep(T, 100))
newdata <- data.frame(
time = timepoints,
log_dtime = log(dtime),
mastectomy = mastectomy,
hormone_therapy = hormone_therapy
)
post <- posterior_linpred(m1.map2stan, newdata = newdata)
str(post)
plot.matrix <- link.surv(post = post, longdata = long_ic2dat)
surv.mean <- as.data.frame( t( apply(plot.matrix, 1, mean) ), colnames = surv.mean)
surv.pi <- as.data.frame( t( apply(plot.matrix, 1, quantile, probs = c(0.055, 0.945)) ), colnames = c( "lower", "upper") )
plot.frame <- get_plot.frame(post = plot.matrix, strata = c("hormone_therapy", "mastectomy") )
obs.mortality <- get_km.frame(ic2surv, strata = c("hormone_therapy", "mastectomy") )
pY <- ggplot2::ggplot(plot.frame %>%
dplyr::mutate(Treatment = ifelse(hormone_therapy & !mastectomy, "Breast conservation with\n hormone therapy", ifelse(!hormone_therapy & mastectomy, "Mastectomy without hormone therapy", "something else"))) %>%
filter(Treatment != "something else"),
aes(time, postmean, group = Treatment, colour = Treatment)) +
geom_line(mapping=aes(group = Treatment, colour = Treatment), alpha = 0.8) +
geom_ribbon(aes( ymin = lower,
ymax = upper,
fill = Treatment), alpha = 0.5, colour=NA) +
geom_step(data= obs.mortality %>%
dplyr::mutate(Treatment = ifelse(as.logical(hormone_therapy) & !as.logical(mastectomy), "Breast conservation with\n hormone therapy", ifelse(!as.logical(hormone_therapy) & as.logical(mastectomy), "Mastectomy without hormone therapy", "something else"))) %>%
filter(Treatment != "something else"), aes(time,
surv,
colour = Treatment), alpha = 0.9)+
labs(title="Bayesian Modeling of clinical outcome for iC-2",
subtitle="Observed vs predicted posterior survival distribution with 89% credible interval\n by type of surgery and hormone therapy") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Survival probability") +
xlab("Time (months)") +
theme_bw() + theme(legend.position=c(0.65, 0.85))
pdf('HormoneSurgeryInteraction.pdf')
pY
dev.off()
pN <- ggplot2::ggplot(plot.frame %>%
dplyr::filter(!mastectomy),
aes(time, postmean, group = hormone_therapy, colour = hormone_therapy)) +
geom_line(mapping=aes(group = hormone_therapy, colour = hormone_therapy), alpha = 0.8) +
geom_ribbon(aes( ymin = lower,
ymax = upper,
fill = hormone_therapy), alpha = 0.5, colour=NA) +
geom_step(data=obs.mortality %>%
dplyr::filter(!as.logical(mastectomy)), aes(time,
surv,
colour = hormone_therapy), alpha = 0.9)+
ggtitle("Practised breast conservation") +
theme(plot.title = element_text(hjust = 0.5)) +
ylab("Survival probability") +
xlab("Time (months)") +
theme_bw() + theme(legend.position=c(0.75, 0.75))
pN
require(gridExtra)
pdf('HormoneRadioInteraction.pdf')
grid.arrange(pN, pY, ncol=2)
dev.off()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.