I have been experiencing some problems with mixing and effective sample size running Weibull models in JAGS, especially in regard to the shape parameter and any covariates estimated as a function of the Weibull hazard. Besides noticing this problem in my own data, I have also noticed it here in the mice data as well. I am not sure what could be causing this problem, but there seems to be high correlation among some parameters.
Here is the mice data:
sp.data <-list(t =structure(.Data =c(12,1,21,25,11,26,27,30,13,12,21,20,23,25,23,29,35,NA,31,36,32,27,23,12,18,NA,NA,38,29,30,NA,32,NA,NA,NA,NA,25,30,37,27,22,26,NA,28,19,15,12,35,35,10,22,18,NA,12,NA,NA,31,24,37,29,27,18,22,13,18,29,28,NA,16,22,26,19,NA,NA,17,28,26,12,17,26),.Dim =c(4,20)),
t.cen =structure(.Data =c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,40,0,0,0,0,0,0,0,40,40,0,0,0,40,0,40,40,40,40,0,0,0,0,0,0,10,0,0,0,0,0,0,0,0,0,24,0,40,40,0,0,0,0,0,0,0,0,0,0,0,20,0,0,0,0,29,10,0,0,0,0,0,0),.Dim =c(4,20)),
M =4, N =20)
And the model:
data{for(iin1:M){for(jin1:N){one[i,j]<-1##Define indicator variable for right censoring}}}model{for(iin1:M){for(jin1:N){one[i,j]~dinterval(t[i,j],t.cen[i,j])t[i,j]~dweib(r,mu[i])}mu[i]<-exp(beta[i])beta[i]~dnorm(0.0,0.001)median[i]<-pow(log(2)*exp(-beta[i]),1/r)}r~dexp(0.001)veh.control<-beta[2]-beta[1]test.sub<-beta[3]-beta[1]pos.control<-beta[4]-beta[1]}
And to run the model:
require(rjags)
sp.inits <-list(beta =c(-1,-1,-1,-1), r =1,
t=with(sp.data,ifelse(is.na(t), t.cen+1,NA)))
m <- jags.model("Mice_JAGS.txt", data=sp.data, inits=sp.inits)
s <- coda.samples(m,c("median","pos.control","r","test.sub","veh.control"), n.iter=10000, thin=10)
And the results from the MCMC:
Inference for Bugs model at "Mice_JAGS.txt", fit using jags,
3 chains, each with 10000 iterations (first 1000 discarded)
n.sims = 27000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
median[1] 24.630 2.423 20.346 22.964 24.454 26.091 29.901 1.002 1700
median[2] 25.106 2.541 20.554 23.363 24.960 26.677 30.569 1.002 2200
median[3] 28.794 2.983 23.623 26.730 28.540 30.596 35.422 1.004 760
median[4] 27.303 2.646 22.627 25.481 27.133 28.937 32.998 1.002 1700
pos.control -0.271 0.347 -0.945 -0.505 -0.272 -0.040 0.416 1.002 3200
r 2.629 0.260 2.136 2.447 2.620 2.811 3.141 1.025 87
test.sub -0.409 0.360 -1.120 -0.648 -0.409 -0.167 0.297 1.002 1900
veh.control -0.050 0.356 -0.747 -0.290 -0.050 0.186 0.655 1.001 12000
deviance 475.415 7.127 463.340 470.337 474.792 479.699 491.354 1.002 1400
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 25.4 and DIC = 500.8
DIC is an estimate of expected predictive error (lower deviance is better).
Here are the results from the mice example, which used a slightly different method for
dealing with the right censoring
Here is the read out from codamenu() on the cross-correlations from one chain:
Yes this is a known problem. Unfortunately the MICE example runs very slowly due to the poor mixing of the shape parameter of the Weibull distribution.
The underlying problem here is the choice of parameterization. The usual parameterization of the Weibull distribution (e.g. the one used by R and explained on Wikipedia) is the shape & scale parameterization. In survival analysis this corresponds to the accelerated life model where the second parameter (scale) corresponds to the scaling of time.
However, the accelerated life model is not popular in biostatistics, where the proportional hazards model is ubiquitous in survival analysis. Hence the BUGS authors chose the alternate parameterization that corresponds to the proportional hazards model (The Weibull distribution is unique in having both an accelerated life parameterization and a proportional hazards parameterization).
In the MICE example, with
t[i, j] ~ dweib(r, mu[i])
the ratio mu[i]/mu[1] is a hazard ratio (i.e. rate ratio) parameter. The variables veh.control, test.sub, and pos.control are all log rate ratios.
It turns out that the proportional hazards parameterization is not a good one for Gibbs sampling when we try to update r and each element of mu[i] separately.
You can use the accelerated life parameterization with a parameter transformation. Instead of putting a prior on beta[i] directly, define it in terms of the shape parameter r and the log of the scale parameter, given a vague normal prior:
I tried reparameterizing the model as you suggested in terms of the scale parameter:
data{for(iin1:M){for(jin1:N){one[i,j]<-1##Define indicator variable for right censoring}}}model{for(iin1:M){for(jin1:N){one[i,j]~dinterval(t[i,j],t.cen[i,j])t[i,j]~dweib(r,mu[i])}mu[i]<-pow(scale[i],-r)scale[i]<-exp(log.scale[i])log.scale[i]~dnorm(0,1.0E-6)beta[i]<-(-r*log.scale[i])median[i]<-pow(log(2)*exp(-beta[i]),1/r)}r~dexp(0.001)veh.control<-beta[2]-beta[1]test.sub<-beta[3]-beta[1]pos.control<-beta[4]-beta[1]}
This does improve the mixing for the shape parameter considerably, and I don't get a message about low effective sample size. I have a couple of factors in the Weibull model that I am fitting to my data that I am using to predict survival probability up to a certain time. I will see if I can incorporate these factors into this reparameterization to solve the mixing problem.
Thank you very much for the help.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Dear JAGS users,
I have been experiencing some problems with mixing and effective sample size running Weibull models in JAGS, especially in regard to the shape parameter and any covariates estimated as a function of the Weibull hazard. Besides noticing this problem in my own data, I have also noticed it here in the mice data as well. I am not sure what could be causing this problem, but there seems to be high correlation among some parameters.
Here is the mice data:
And the model:
And to run the model:
And the results from the MCMC:
Here are the results from the mice example, which used a slightly different method for
dealing with the right censoring
Here is the read out from codamenu() on the cross-correlations from one chain:
Any help or suggestions would be greatly appreciated!
Last edit: Daniel Eacker 2015-02-06
Yes this is a known problem. Unfortunately the MICE example runs very slowly due to the poor mixing of the shape parameter of the Weibull distribution.
The underlying problem here is the choice of parameterization. The usual parameterization of the Weibull distribution (e.g. the one used by R and explained on Wikipedia) is the shape & scale parameterization. In survival analysis this corresponds to the accelerated life model where the second parameter (scale) corresponds to the scaling of time.
However, the accelerated life model is not popular in biostatistics, where the proportional hazards model is ubiquitous in survival analysis. Hence the BUGS authors chose the alternate parameterization that corresponds to the proportional hazards model (The Weibull distribution is unique in having both an accelerated life parameterization and a proportional hazards parameterization).
In the MICE example, with
the ratio
mu[i]/mu[1]
is a hazard ratio (i.e. rate ratio) parameter. The variablesveh.control
,test.sub
, andpos.control
are all log rate ratios.It turns out that the proportional hazards parameterization is not a good one for Gibbs sampling when we try to update
r
and each element ofmu[i]
separately.You can use the accelerated life parameterization with a parameter transformation. Instead of putting a prior on
beta[i]
directly, define it in terms of the shape parameterr
and the log of the scale parameter, given a vague normal prior:You will find this works much better.
Last edit: Martyn Plummer 2015-02-06
Hi Martyn,
I tried reparameterizing the model as you suggested in terms of the scale parameter:
This does improve the mixing for the shape parameter considerably, and I don't get a message about low effective sample size. I have a couple of factors in the Weibull model that I am fitting to my data that I am using to predict survival probability up to a certain time. I will see if I can incorporate these factors into this reparameterization to solve the mixing problem.
Thank you very much for the help.