Weibull Model in JAGS: Adding in covariates

2014-02-12
2014-05-22
  • Daniel Eacker

    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))
    m <- jags.model("updated.txt", data=weib.data, inits=weib.inits,n.chains=1,n.adapt=20000)
    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

    Martyn Plummer - 2014-02-13

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

    http://www.openbugs.net/Examples/Mice.html
    http://www.openbugs.net/Examples/Kidney.html

    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

      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

    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

    Daniel Eacker - 2014-05-22

    Thank you atj409, I will check this book out.

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks