Гамма-распределение в JAGS - Ошибка в узле

Question

Я пытаюсь параметризовать гамма-распределение в JAGS - кусочно-линейным предиктором, но моя модель не запускается со следующим сообщением об ошибке:

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

Модель работает, когда timber.recovery извлекается из нормального распределения, но прогнозы нижнего квантиля меньше нуля, что биологически невозможно. Я попробовал несколько приемов, таких как добавление 0,001 к параметру "mu" на случай, если он рисует ноль, устанавливая начальные значения на основе выходных данных glm; но ни один не разрешает сообщение об ошибке. Любая идея будет принята с благодарностью [я использую R2jags]. Моя модель:

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")