How to transform the result of stan in survival analysis

[yanmin]

i perform both a classic survival analysis with the function survreg of 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 the output similar to the result of frequentist approach.
the data is from the survival package. the codes 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 \sim Weibull(\rho,\mu_i) \\ log(\mu_i) = X * \beta \\ \beta \sim normal(0,100^2) \\ \rho \sim cauchy(1,10)

Any help is appriciated!

[Bob_Carpenter]

Hi, @yanmin and welcome to the Stan forums! Sorry it’s taken some time to get you an answer.

@yanmin

i believe that there is no significant difference between the two approaches. so, how to transform the result of rstan to obtain the output similar to the result of frequentist approach.

Can you say more specifically what it is you want?

I’m not exactly sure what you want to transfer and I don’t know the survreg package in R.

@yanmin

to_vector(beta) ~ normal(0,100^2);

You don’t need to apply to_vector because beta is already a vector here. Also, Stan uses a scale parameterization of the normal, so you probably want 100 rather than 100^2 unless you meant the standard deviation to be 10,000.

And you should only use a Cauchy for rho if you think it might take extreme values in the 1000s or bigger.

The Weibull in Stan takes a shape and a scale parameter in that order. I’d check those are in the right order.

Then, it’s more efficient to use the sampling versions, such as rho ~ cauchy(1, 10); (Stan doesn’t require the extra .0 on arguments.)

[yanmin]

Bob, thank you for your reply. i am performing a suvival analysis, and want to use the result from frequentist approach as robustness check. in general, the bayesian result is considered as robust if there is no significant difference between bayesian and frequentist approach. in the work of julia et al., there is no significant difference between bayesian from stran and frequentist approach from the function survreg of the survival package (Therneau, 2015). they transform the result from stan using a specific equation, as shown in their paper. i mailed to julia and was told that the transform equation should be different because they used jags. so, i posted and request for help.

[jroon]

As it happens I had this exact same issue yesterday. In the documentation for survreg you will find this:

So I believe to compare your Stan output to survreg output you get the exponential or log transform of the parameter estimates (depending on which parameter you talk about)

[yanmin]

bob, your answer is correct. thank you very much!