0

I'm trying to parameterise a gamma distribution in JAGS - with a piecewise linear predictor but my model fails to run with the following error message: 

Error: Error in node (ashape/(aexp(mu[59]))) Invalid parent values

The model works when timber.recovery is drawn from a normal distribution, but the lower quantile predictions is less than zero, which is not biologically possible. I've tried a few tricks like adding 0.001 to the "mu" parameter in case it was drawing a zero, setting initial values based on outputs from a glm; but neither resolves the error message. any insights would be greatly appreciated [i'm using R2jags]. My model:

cat (
"model {

    # UNINFORMATIVE PRIORS
    sd_plot ~ dunif(0, 100)
    tau_plot <- 1/(sd_plot * sd_plot) 
      # precision for plot level variance

    alpha ~ dnorm(0, 1e-06) 
      # normal prior for intercept term 
    shape ~ dunif(0, 100) 
      # shape parameter for gamma     

    log_intensity ~ dnorm(0, 1e-06) 
      # uninformative prior for logging intensity

    beta_1 ~ dnorm (0, 1e-06) 
      # uninformative prior; change in slope for first segment : <=3.6 years
    beta_2 ~ dnorm (0, 1e-06) 
      # uninformative prior; change in slope for first segment : >3.6 years
    InX_1 ~ dnorm (0, 1e-06) 
      # uniformative prior for interaction between tsl and log_intensity : <=3.6 years
    InX_2 ~ dnorm (0, 1e-06) 
      # uniformative prior for interaction between tsl and log_intensity : >3.6 years

    # PLOT LEVEL RANDOM EFFECTS
    for (i in 1:nplots) {

      plot_Eff[i] ~ dnorm(0,tau_plot)

    }


    for (i in 1:Nobs) {

      # PIECEWISE LINEAR PREDICTOR     
      mu[i] <- 
        alpha + 
        beta_1 * (time.since.logged[i] * tsl.DUM1[i]) + 
        log_intensity * log.volume [i] + 
        beta_2 * (time.since.logged[i] * tsl.DUM2[i] - 3.6) + 
        beta_1 * (time.since.logged[i] * tsl.DUM2[i]) + 
        plot_Eff[plot.id[i]] + 
        InX_1 * (time.since.logged[i] * tsl.DUM1[i]) * log.volume [i] + 
        InX_2 * (time.since.logged[i] * tsl.DUM2[i] - 3.6) * log.volume[i] + 
        InX_1 * (time.since.logged[i] * tsl.DUM2[i]) * log.volume[i]


      timber.recovery[i] ~ dgamma(shape,shape/exp(mu[i])) 
        # observed recovery


      pred_timber_recovery[i] ~  dgamma(shape,shape/exp(mu[i])) 
        # posterior predictive distribution

      pearson.residual[i] <- 
        (timber.recovery[i] - pred_timber_recovery[i]) / (sqrt(timber.recovery[i]))
    }


 }",
 fill = TRUE, 
 file = "outputs/piecewise_TIMBER_MODEL_FINAL_GAMMA.txt")
Geoff
  • 4,676
  • 3
  • 26
  • 38
Anand Roopsind
  • 581
  • 2
  • 8
  • 11
  • Seems to be failing when i=59. Have you checked that data values at i=59 and the resulting mu? – John K. Kruschke Dec 12 '16 at 18:14
  • I have - i even modified the data values for i=59, but it breaks at other points. It's curious as these are the highest value of one of the predictors [log.volume]; the reason i wanted to move to a gamma from a normal was that some of my predictions [simulated with variance] was going below zero. I'm considering moving to a truncated normal, if i can't figure out the gamma issue. Thankyou for taking the time to have a look at my code. – Anand Roopsind Dec 13 '16 at 18:48

0 Answers0