## Weibull Model in JAGS: Adding in covariates

2014-02-12
2014-05-22

• Daniel Eacker
2014-02-12

Hi all,

I have searched for help on the internet on running a Weibull model using Bayesian analysis in JAGS. While there are some good examples out there, none of them make it clear how to code in covariates. Also, it is a little unclear what the lambda is in the following JAGS code: dweib(v, lambda). Is it the scale parameter or does it enter into the following equation: scale<-pow(lambda,-1/shape)?

I have used the trick provided by martyn_plummer for dealing with the right-sensored observations in the survival analysis. I also get somewhat reasonable estimates of log(shape) and log(scale) when I run the following code (as compared to the phreg() function in the R package eha):

## Jags code for Weibull Model

data {
for(j in 1:277){
one[j] <- 1
}
}

model {
for(j in 1:277){
one[j] ~ dinterval(y[j], y.cens[j])
y[j]~dweib(shape,lambda)
}
scale<-pow(lambda,-1/shape)
mean<-scale*exp(loggam(1+1/shape))
logscale<-log(scale)
logshape<-log(shape)
# Priors
shape~dunif(0,5)
lambda~dunif(0,1)
}
library(rjags)
weib.data<-list(y=Y,y.cens=Y.cens)
weib.inits <- list(shape=0.8,lambda=0.5)
weib.inits\$y <- with(weib.data, ifelse(is.na(y), y.cens+1, NA))
update(m, 10000)
s <- coda.samples(m, c("lambda","mean","logscale","logshape"),n.iter=20000)

The data is for elk calf survival, and includes a vector of failure times, sensoring times, and sex (coded as 0=female, male =1).

While I understand the parameters in the Weibull distribution (i.e., shape, scale), I am confused about how to:
1) code in covariates into the model (sex.beta in the data below)?
2) whether the above code is the correct formulation for a proportional hazards weibull model?

Here is the data:
weib.data
\$y
[1] NA NA NA NA NA NA 68 NA NA 6 234 NA NA NA NA NA NA 62
[19] NA NA NA NA 54 NA NA 18 104 35 283 152 NA 258 NA NA 215 NA
[37] NA 237 NA 218 NA NA NA NA 292 225 83 NA 14 6 NA NA 15 NA
[55] NA NA 88 71 NA NA NA NA 61 NA NA NA 97 244 15 NA 19 NA
[73] NA 115 14 NA NA NA NA NA NA NA NA NA NA NA NA 262 NA NA
[91] NA NA NA NA NA NA NA 8 NA 12 NA NA NA NA NA NA 142 NA
[109] NA NA NA 121 NA NA NA 25 NA NA NA NA 24 NA 161 145 100 NA
[127] NA NA 173 83 NA NA 3 NA 22 NA NA NA NA 83 86 21 27 16
[145] 22 NA NA 97 25 14 29 NA NA 15 25 NA NA 40 NA 67 NA NA
[163] NA 224 NA 247 320 195 316 257 NA NA 193 203 NA NA 192 NA 11 NA
[181] NA NA NA 15 NA 156 88 NA 100 211 NA NA NA 21 NA NA 34 NA
[199] NA NA 84 NA NA NA 20 61 94 28 NA NA NA NA 144 136 13 NA
[217] 14 NA NA NA NA NA NA NA NA 273 NA NA NA 266 NA 13 4 NA
[235] NA NA 12 NA NA NA NA 9 13 71 NA NA NA NA NA 37 19 173
[253] NA NA 157 154 NA 110 NA NA 24 48 NA NA NA 233 143 NA NA 153
[271] 78 NA NA NA NA 128 NA

\$y.cens
[1] 136 126 45 91 67 199 0 143 91 0 0 172 143 152 77 172 163 0
[19] 199 143 199 102 0 219 193 0 0 0 0 0 363 0 363 363 0 363
[37] 363 0 363 0 363 363 363 363 0 0 0 70 0 0 198 232 0 218
[55] 360 168 0 0 113 346 163 57 0 42 77 14 0 0 0 330 0 337
[73] 51 0 0 11 296 42 360 360 360 360 360 360 360 360 360 0 360 360
[91] 360 190 360 300 76 263 168 0 263 0 263 263 158 263 263 263 0 102
[109] 263 263 67 0 263 263 52 0 263 263 263 263 0 65 0 0 0 263
[127] 263 263 0 0 61 3 0 88 0 219 132 112 363 0 0 0 0 0
[145] 0 84 81 0 0 0 0 258 221 0 0 95 63 0 258 0 143 132
[163] 363 0 363 0 0 0 0 0 363 363 0 0 248 363 0 363 0 37
[181] 232 118 300 0 47 0 0 220 0 0 38 198 146 0 360 54 0 132
[199] 39 48 0 118 53 302 0 0 0 0 20 27 289 78 0 0 0 294
[217] 0 190 360 360 360 360 360 360 360 0 356 360 360 0 263 0 0 263
[235] 263 263 0 263 263 263 263 0 0 0 263 263 108 263 263 0 0 0
[253] 263 263 0 0 120 0 263 263 0 0 263 263 263 0 0 93 263 0
[271] 0 263 263 263 263 0 263

\$sex.beta
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Levels: 0 1

And the output:
chain1

Iterations = 20001:40000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

```       Mean        SD  Naive SE Time-series SE
```

lambda 0.007494 2.587e-03 1.829e-05 1.951e-04
logscale 6.391826 1.543e-01 1.091e-03 7.226e-03
logshape -0.257124 8.137e-02 5.754e-04 7.550e-03
mean 711.618629 1.473e+02 1.041e+00 1.084e+01

2. Quantiles for each variable:

```       2.5%        25%        50%        75%      97.5%
```

lambda 0.003443 0.005573 0.007228 0.009047 0.01341
logscale 6.115110 6.283464 6.385314 6.490020 6.71774
logshape -0.420408 -0.310682 -0.258977 -0.201057 -0.10120
mean 489.420815 607.137902 689.058124 789.291926 1065.68333

Thank you for looking and any helpful comments would be greatly appreciated!

• Martyn Plummer
2014-02-13

See the KIDNEY and MICE examples from OpenBUGS for how to define a proportional hazards model:

These examples have been translated into the JAGS dialect of the BUGS language in the tarball that you can download here:

https://sourceforge.net/projects/mcmc-jags/files/Examples/3.x/

• Daniel Eacker
2014-02-13

Thank you very much for the response Martyn. I have been working with the Mice and Kidney examples, and did explore the other related posts on the topic here on the JAGS forum. I will give them another look, and check out the translation website above. I have been able to run the model in JAGS with one covariate (sex) and will work on incorporating some continuous covariates as well.

Thank you.

• atj409
2014-03-27

Also check out this website that has companion BUGS code for the book, "Bayesian Survival Analysis." There are some additional and highly useful examples.
http://www.stat.uconn.edu/~mhchen/survbook/

• Daniel Eacker
2014-05-22

Thank you atj409, I will check this book out.