if (instantiate::stan_cmdstan_exists()) {
if(requireNamespace("survival")){
library(survival)
data(E1684)
data(E1690)
## take subset for speed purposes
E1684 = E1684[1:100, ]
E1690 = E1690[1:50, ]
## replace 0 failure times with 0.50 days
E1684$failtime[E1684$failtime == 0] = 0.50/365.25
E1690$failtime[E1690$failtime == 0] = 0.50/365.25
E1684$cage = as.numeric(scale(E1684$age))
E1690$cage = as.numeric(scale(E1690$age))
data_list = list(currdata = E1690, histdata = E1684)
nbreaks = 3
probs = 1:nbreaks / nbreaks
breaks = as.numeric(
quantile(E1690[E1690$failcens==1, ]$failtime, probs = probs)
)
breaks = c(0, breaks)
breaks[length(breaks)] = max(10000, 1000 * breaks[length(breaks)])
d.bhm = pwe.bhm(
formula = survival::Surv(failtime, failcens) ~ treatment + sex + cage + node_bin,
data.list = data_list,
breaks = breaks,
chains = 1, iter_warmup = 500, iter_sampling = 1000
)
pwe.logml.map(
post.samples = d.bhm,
bridge.args = list(silent = TRUE),
chains = 1, iter_warmup = 1000, iter_sampling = 2000
)
}
}
Run the code above in your browser using DataLab