# Understanding Discrepancies Between Frequentist and Bayesian Parametric Weibull Models in Accelerated Failure Time Analysis

Currently in the process of identifying a Bayesian Weibull Survival AFT model that is equivalent to the survreg() Weibull AFT model results from the Survival package in R. I came across this publication and I am trying to match the survreg results with the use of low informative priors. I have not succeeded yet.

Are these two models equivalent (bayewei vs. freqwei)? If yes, why the results are different and how can I acquire comparable results?

 mc<- structure(list(V1 = c(5.3, 6.8, 5.7, 3.7, 6, 6.3, 7.9, 5.9, 7.6,
7.1, 6, 6.9, 6.1, 8.8, 5.6, 6.1, 5, 6.6, 5.6, 6.1),
V2 = structure(c(1L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L,
2L, 2L, 2L), levels = c("0", "1"), class = "factor"),
V3 = c(0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1),
X = c(1.65,
1.79, 1.77, 2.74, 4.62, -0.86, 2.7, 3.47, 1.95, 4.3, 2.52, 1.95,
2.78, 2.7, 4.25, -0.19, 3.76, 0.68, 1.48, -0.06),
Time = c(38.863,
24.2, 120.373, 16.5, 115.533, 8.613, 106.997, 105.71, 105.963,
35.42, 107.613, 14.113, 105.82, 14.377, 99.253, 34.793, 96.657,
94.347, 91.19, 8.613),
Status = c(0, 0, 0, 1, 0, 1, 0, 0, 0,
1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1)),
row.names = c(1L, 2L, 4L,
5L, 7L, 9L, 11L, 12L, 14L, 20L, 21L, 24L, 25L, 31L, 32L, 33L,
37L, 40L, 41L, 42L), class = "data.frame")

library(rjags)
library(survival)

BayesSurv<-function(form="X", data=mc, model="weibull",
prior_mu=0, prior_prec=1,var,
eventvariable="Status",
timevariable="Time",
niter=10000, nthin=10,nburn=1000){
.formula <- reformulate(form)
data$$event<-data[,eventvariable] data$$time<-data[,timevariable]
X <- model.matrix(object=.formula, data = data)

time_un <- data$$time[data$$event==1]
cen_un <- data$$time[data$$event==1]
is.censored_un <- rep(0, length(cen_un))
X_un <- X[data$event==1, ] cen_cen <- data$$time[data$$event==0] time_cen <- rep(NA, length(cen_cen)) is.censored_cen <- rep(1, length(cen_cen)) X_cen <- X[data$event==0, ]

time <- c(time_un, time_cen)
cen <- c(cen_un, cen_cen)
is.censored <- c(is.censored_un, is.censored_cen)
Xnew <- rbind(X_un, X_cen)

cat("model{
for(i in 1:n){
is.censored[i]~dinterval(time[i],cen[i])
time[i]~dweib(alpha,lambda[i])
lambda[i]<-exp(-mu[i]*alpha)
mu[i]<-inprod(beta[],X[i,])}

#Prior distributions
for(l in 1:(Nbetas)){
beta[l]~dnorm(mu.betaMain[l],prec.betaMain[l])}
alpha ~ dunif(0, 10)
}",file = "AFT.txt")

d.jags <- list(n = nrow(X), time = time, cen = cen, X = Xnew, is.censored = is.censored, Nbetas = ncol(X),mu.betaMain=prior_mu,prec.betaMain=prior_prec)
i.jags <- function() list(beta = rnorm(ncol(X),sd = 1), alpha = runif(1,min = 0.1,max = 10))
p.jags <- c("beta", "alpha")

m1 <- rjags::jags.model(data = d.jags, file = "AFT.txt", inits = i.jags, n.chains = 1)

# Burn-in and update
update(m1, nburn)
res <- coda.samples(m1, variable.names = p.jags, n.iter = niter, thin = nthin)
results <- as.mcmc(do.call(rbind, res))

return(results)
}

freqwei<-summary(survreg(Surv(Time,Status)~V1 + V2 + V3 + X, data=mc))

outn<-BayesSurv(form="V1 + V2 + V3 + X",prior_mu=c(0,rep(0,5)),prior_prec=c(rep(0.001,5)),
data = mc,niter = 20000)
bayewei=summary(outn)
bayewei

BOut<-cbind(freqwei$$coefficients, bayewei[2]$$quantiles[,3][-1])
BOut