Menu

Weibull model: problem with slow mixing and effective sample size

2015-02-06
2015-02-06
  • Daniel Eacker

    Daniel Eacker - 2015-02-06

    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:
    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(i in 1:M){
        for(j in 1:N){  
             one[i, j] <- 1  ##Define indicator variable for right censoring
    }
    }
    }
    
    model{   
      for(i in 1 : M) {
         for(j in 1 : 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

    alternate text

    Here is the read out from codamenu() on the cross-correlations from one chain:
                deviance median[1] median[2] median[3] median[4] pos.control       r
    deviance      1.0000    0.2841   0.47575   0.45853   0.23974     0.05988 -0.3685
    median[1]     0.2841    1.0000   0.01663   0.03986   0.03732     0.68181  0.1627
    median[2]     0.4758    0.0166   1.00000   0.02804   0.00684    -0.00287  0.1571
    median[3]     0.4585    0.0399   0.02804   1.00000   0.00624     0.01447  0.1348
    median[4]     0.2397    0.0373   0.00684   0.00624   1.00000    -0.69373  0.1628
    pos.control   0.0599    0.6818  -0.00287   0.01447  -0.69373     1.00000 -0.0685
    r            -0.3685    0.1627   0.15710   0.13478   0.16284    -0.06854  1.0000
    test.sub     -0.1049    0.6499  -0.02454  -0.71847   0.00480     0.47715 -0.0868
    veh.control  -0.1472    0.6806  -0.71231   0.00370   0.02021     0.48215 -0.0135
            test.sub veh.control
    deviance     -0.1049     -0.1472
    median[1]     0.6499      0.6806
    median[2]    -0.0245     -0.7123
    median[3]    -0.7185      0.0037
    median[4]     0.0048      0.0202
    pos.control   0.4771      0.4821
    r            -0.0868     -0.0135
    test.sub      1.0000      0.4780
    veh.control   0.4780      1.0000
    

    Any help or suggestions would be greatly appreciated!

     

    Last edit: Daniel Eacker 2015-02-06
  • Martyn Plummer

    Martyn Plummer - 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

             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:

         beta[i] <- -r * log.scale[i])
         log.scale[i] ~ dnorm(0, 1.0E-6)
         scale[i] <- exp(log.scale[i])
    

    You will find this works much better.

     

    Last edit: Martyn Plummer 2015-02-06
    • Daniel Eacker

      Daniel Eacker - 2015-02-06

      Hi Martyn,

      I tried reparameterizing the model as you suggested in terms of the scale parameter:

      data{
      for(i in 1:M){
          for(j in 1:N){  
              one[i, j] <- 1  ##Define indicator variable for right censoring
      
         }
      }
      }
      
      model{   
        for(i in 1 : M) {
           for(j in 1 : 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.

       

Log in to post a comment.