I performed both a classic survival analysis with the function survreg() from the survival package, and a bayesian analysis with rstan. I believe that there is no significant difference between the two approaches. So, how to transform the result of rstan to obtain an output similar to the result of the frequentist approach?
The data is from the survival package. The code is as follows:
#rm(list=ls())
library(survival)
data_for_Stan <- list(n_cen = with(ovarian,sum(fustat)),
n_not_cen = with(ovarian,sum(1 - fustat)),
t_cen = with(ovarian,futime[fustat == 1]),
t_not_cen = with(ovarian,futime[fustat == 0]),
x_cen = ovarian %>% mutate(intercept = 1) %>%
filter(fustat == 1) %>%
dplyr::select(intercept,ecog.ps,rx),
x_not_cen = ovarian %>% mutate(intercept = 1) %>%
filter(fustat == 0) %>%
dplyr::select(intercept,ecog.ps,rx),
K= 3)
##Model with Bayessian framework
library(rstan)
Bayesian_model <-
"
data {
int < lower = 1 > n_cen;
int < lower = 1 > n_not_cen;
int < lower = 1 > K;
vector[n_cen] t_cen;
vector[n_not_cen] t_not_cen;
matrix[n_cen,K] x_cen;
matrix[n_not_cen,K] x_not_cen;
}
parameters {
vector[K] beta;
real < lower = 0 > rho;
}
model {
to_vector(beta) ~ normal(0,100^2);
target += cauchy_lpdf(rho|1.0, 10.0);
target += weibull_lpdf(t_not_cen | rho, exp(x_not_cen * beta));
target += weibull_lccdf(t_cen | rho, exp(x_cen * beta));
}
"
write(Bayesian_model,file = "Bayesian_model.stan")
fit_predict = stan(file = "Bayesian_model.stan",
data = data_for_Stan, warmup = 1000,
iter = 2000, chains = 4,
cores = 2, thin = 1)
fit <- survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist='weibull')
fit_predict
summary(fit)$table
my motivation is from the work of Julia Stander et al. (A Bayesian Survival Analysis of a Historical Dataset: How Long Do Popes Live?The American Statistician,2017)
My Bayesian model is:
T_i ~ Weibull(rho, mu_i)
log(mu_i) = X * beta
beta ~ normal(0,100^2)
rho ~ cauchy(1,10)