Dear list,
I am probably missing something obvious here, but I seem to have some inconsistent results using a Weibull regression model --- I'm trying to compare the results obtained using MLE (flexsurv), Bayesian analysis using INLA and Bayesian analysis using MCMC (via JAGS, of course).
I have a set of survival data for two treatment arms. Say these are data from a cancer trial, where the "event" is death but patients may also exit the observation because they progress to a worse health state. So my data look like this
Notice that patient no. 91 is censored because they progress (after 0.79 time units) and so we don't know whether they die (in this subset of the data).
My understanding is that INLA and JAGS have the same parameterisation for the Weibull model, in terms of a shape parameter and a function of the linear predictor, say lambda. So, in JAGS, my model is the following:
Here, I'm also rescaling the coefficients in the linear predictor (beta0 and beta1) to estimate the scale and the treatment effect so that these could be compared with the standard output from flexsurv. I can do similar re-scaling in INLA.
I manipulate the original data (which I can directly feed to both flexsurv and INLA) to create a suitable list for JAGS, like so:
dataJags <-list(t=ifelse(data$censored==0,data$time,NA),
censored=data$time,
is.censored=data$censored,
treat=data$treat,n=length(data$arm),# These are the same priors as in the default INLA
shape.inla=25,rate.inla=25,
prec.fe.inla=0.001,
prec.int.inla=0.000001)
The JAGS model runs, but (here's where I may be missing the obvious point!) I do get weird results for the coefficient beta1 (and hence its rescaled version effect).
Here's the output from flexsurv:
data mean est L95% U95% se exp(est) L95% U95%
shape NA 1.8164 1.6134 2.0449 0.1098 NA NA NA
scale NA 10.2210 9.1617 11.4026 0.5705 NA NA NA
treat 0.4850 0.3420 0.1744 0.5097 0.0855 1.4078 1.1905 1.6648
(NB: you could compute (Intercept)^(1/-alpha) to get a value that is comparable to the scale from flexsurv and similarly treat1^(1/-alpha) to get a similar value to the treatment effect from the MLE analysis)
As you can see, the estimate for the shape and beta0 (hence the scale) are kind of in line with the other methods. However, beta1 (and hence the effect) are way off the MLE and INLA analysis.
Is there an obvious reason why this should happen, that I just don't see? Or a blatant difference in the way that the censoring is taken into account (but even if so, I am not sure why this would affect the coefficient for the treatment effect only)?
Many thanks!
Gianluca
Last edit: Martyn Plummer 2016-01-04
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
There are two distinct parameterizations of the Weibull distribution and you cannot switch back and forth between them so easily. In fact you can only do this when the Weibull observations are independent and identically distributed. In a regression model, the two parameterizations have different interpretations: one is the accelerated life model and the other is the proportional hazards model.
R uses the shape/scale parameterization of the Weibull distribution. Hence the survreg function in the survival package and the flexsurvreg function from the flexsurv package use the same parameterization for regression models, which leads to the accelerated life model.
You can fit an accelerated life model in JAGS with the Weibull distribution but you must use the generalized gamma distribution dgen.gamma and not the Weibull distribution. A reproducible demonstration is given below
Note the long run for the JAGS model. Unfortunately the mixing is rather poor with this parameterization (This is fixable but will have to wait for a future release). If you run the chain long enough then you do get compatible results between JAGS and INLA, as shown below (and again removing some extraneous details).
Your code and explanations above has been incredibly helpful for building Weibull AFT models in JAGS. I'm hoping to run simulations based on this model. I've attached my model of interest below (it's a simpler case - I'm interested in whether people are treated or not - 'treat'). After hundreds of repeated runs, I've been running into "Error in node alpha - Current value is inconsistent with data". What am I missing here? Thanks!
This was happened around a seedNum, 57 and c=0.15 on my system.
numObs<-125# per armdat<-data.frame((matrix(NA,nrow=numObs*2,ncol=3)))colnames(dat)<-c("time","delta","trt")# Generate dataset.seed(seedNum)Y0<-rweibull(numObs,0.5,4)Y1<-rweibull(numObs,0.5,4+exp(0.5))# Cleandat$time<-c(Y0,Y1)dat$delta<-sample(c(0,1),size=numObs*2,prob=c(c,1-c),replace=T)dat$trt<-c(rep(0,numObs),rep(1,numObs))jags.data<-with(dat,list("is.censored"=(delta==0),"tcens"=time,"t"=ifelse(delta==1,time,NA),"trt"=trt+1))inits1<-with(jags.data,list("t"=ifelse(is.censored,tcens+0.1,NA_real_),".RNG.name"="base::Wichmann-Hill",.RNG.seed=1))inits2<-with(jags.data,list("t"=ifelse(is.censored,tcens+0.1,NA_real_),".RNG.name"="base::Wichmann-Hill",.RNG.seed=2))inits3<-with(jags.data,list("t"=ifelse(is.censored,tcens+0.1,NA_real_),".RNG.name"="base::Wichmann-Hill",.RNG.seed=3))params<-c("shape","scale","treat")cl<-makePSOCKcluster(3)tmp<-clusterEvalQ(cl,library(dclone))pm<-jags.parfit(cl,jags.data,params,modelString,inits=list(inits1,inits2,inits3))# I also saw this error before parallelizing the code. # m <- jags.model(textConnection(jags.model),# data=jags.data,# inits=list(inits1, inits2, inits3), n.chains=3)
Last edit: Arlene Jiang 2024-04-15
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The issue seems to be that you need to set the initial values for some of the parameters explicitly. In my modified version of your example below, I have explicitly set the parameters shape, scale, and treat and this seems to work.
There is another issue with your example. Your mechanism for removing observations is not censoring. They way you have set it up, the missing values are "missing completely at random", using the usual Little-Rubin terminology for missing data. This form of missingness is ignorable and you do not need to encode the missing data mechanism in your bugs code. In my version, the missing data are right-censored, i.e. any failure after the censoring time tcens is not observed: we know only that the failure did not occur before this time. To model this data you need the dinterval distribution in JAGS.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Dear list,
I am probably missing something obvious here, but I seem to have some inconsistent results using a Weibull regression model --- I'm trying to compare the results obtained using MLE (flexsurv), Bayesian analysis using INLA and Bayesian analysis using MCMC (via JAGS, of course).
I have a set of survival data for two treatment arms. Say these are data from a cancer trial, where the "event" is death but patients may also exit the observation because they progress to a worse health state. So my data look like this
Notice that patient no. 91 is censored because they progress (after 0.79 time units) and so we don't know whether they die (in this subset of the data).
My understanding is that INLA and JAGS have the same parameterisation for the Weibull model, in terms of a
shapeparameter and a function of the linear predictor, saylambda. So, in JAGS, my model is the following:Here, I'm also rescaling the coefficients in the linear predictor (
beta0andbeta1) to estimate thescaleand the treatmenteffectso that these could be compared with the standard output from flexsurv. I can do similar re-scaling in INLA.I manipulate the original data (which I can directly feed to both flexsurv and INLA) to create a suitable list for JAGS, like so:
and the initial values like so:
The JAGS model runs, but (here's where I may be missing the obvious point!) I do get weird results for the coefficient
beta1(and hence its rescaled version effect).Here's the output from flexsurv:
the one from INLA
Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld (Intercept) -3.9664 0.2016 -4.3405 -3.9642 -3.5959 -3.9544 0 treat1 -0.6005 0.1549 -0.9066 -0.5999 -0.2982 -0.5986 0 Model hyperparameters: mean sd 0.025quant 0.5quant 0.975quant mode alpha parameter for weibull 1.6968 0.0997 1.5070 1.6948 1.8989 1.6918(NB: you could compute
(Intercept)^(1/-alpha)to get a value that is comparable to the scale from flexsurv and similarlytreat1^(1/-alpha)to get a similar value to the treatment effect from the MLE analysis)and the one from JAGS
n.sims = 1000 iterations saved mu.vect sd.vect 2.5% 97.5% Rhat n.eff beta0 -4.013 0.226 -4.444 -3.568 1.004 380 beta1 -0.139 0.142 -0.413 0.127 1.007 220 effect 1.090 0.093 0.925 1.280 1.008 200 scale 11.027 0.720 9.690 12.510 1.006 370 shape 1.674 0.094 1.487 1.862 1.001 1000As you can see, the estimate for the shape and
beta0(hence the scale) are kind of in line with the other methods. However,beta1(and hence the effect) are way off the MLE and INLA analysis.Is there an obvious reason why this should happen, that I just don't see? Or a blatant difference in the way that the censoring is taken into account (but even if so, I am not sure why this would affect the coefficient for the treatment effect only)?
Many thanks!
Gianluca
Last edit: Martyn Plummer 2016-01-04
There are two distinct parameterizations of the Weibull distribution and you cannot switch back and forth between them so easily. In fact you can only do this when the Weibull observations are independent and identically distributed. In a regression model, the two parameterizations have different interpretations: one is the accelerated life model and the other is the proportional hazards model.
R uses the shape/scale parameterization of the Weibull distribution. Hence the
survregfunction in the survival package and theflexsurvregfunction from the flexsurv package use the same parameterization for regression models, which leads to the accelerated life model.You can fit an accelerated life model in JAGS with the Weibull distribution but you must use the generalized gamma distribution
dgen.gammaand not the Weibull distribution. A reproducible demonstration is given belowComparing results (and cutting out some details from the printed output) shows compatible results
INLA uses the alternative proportional hazards parameterization of the Weibull distribution, and so does the
dweibdistribution in BUGS.Note the long run for the JAGS model. Unfortunately the mixing is rather poor with this parameterization (This is fixable but will have to wait for a future release). If you run the chain long enough then you do get compatible results between JAGS and INLA, as shown below (and again removing some extraneous details).
Small differences between JAGS and INLA results may be due to the choice of prior: I have not made an effort to use the INLA prior in the JAGS model.
Last edit: Martyn Plummer 2016-01-05
Hi Martyn,
Your code and explanations above has been incredibly helpful for building Weibull AFT models in JAGS. I'm hoping to run simulations based on this model. I've attached my model of interest below (it's a simpler case - I'm interested in whether people are treated or not - 'treat'). After hundreds of repeated runs, I've been running into "Error in node alpha - Current value is inconsistent with data". What am I missing here? Thanks!
This was happened around a seedNum, 57 and c=0.15 on my system.
Last edit: Arlene Jiang 2024-04-15
The issue seems to be that you need to set the initial values for some of the parameters explicitly. In my modified version of your example below, I have explicitly set the parameters shape, scale, and treat and this seems to work.
There is another issue with your example. Your mechanism for removing observations is not censoring. They way you have set it up, the missing values are "missing completely at random", using the usual Little-Rubin terminology for missing data. This form of missingness is ignorable and you do not need to encode the missing data mechanism in your bugs code. In my version, the missing data are right-censored, i.e. any failure after the censoring time tcens is not observed: we know only that the failure did not occur before this time. To model this data you need the dinterval distribution in JAGS.
Thank you Martyn,
in your jags.model2:
jags.model2 <- "model {
for (i in 1:length(t)) {
is.censored[i] ~ dinterval(t[i], tcens[i])
t[i] ~ dweib(shape, lambda[i])
log(lambda[i]) <- alpha + beta[group[i]]
}
alpha ~ dnorm(0, 1.0E-6)
beta[1] <- 0
beta[2] ~ dnorm(0, 1.0E-6)
beta[3] ~ dnorm(0, 1.0E-6)
shape ~ dexp(0.01)
Intercept <- alpha
groupMedium <- beta[2]
groupPoor <- beta[3]
}"
how can we get scale like in jags. Model : scale <- exp(-alpha)?
thank you