Перевести отрицательную биномиальную дисперсионную модель типа I и II из glmmTMB в JAGS

Question

Я хочу перевести модель nbinom1 из glmmTMB в JAGS. Но я не знаю, как перевести формулу дисперсии в JAGS. В этом примере я использую nbinom2, потому что он проще, но на самом деле я бы предпочел реализовать модель дисперсии из семейства nbinom1 от glmmTMB до JAGS.

Обратите внимание, что пример регрессии здесь полностью вне контекста и не имеет смысл. Это только для того, чтобы предоставить код, поскольку фактическая модель, с которой я работаю, намного сложнее.

У нас есть отрицательная биномиальная модель, которая определена в R следующим образом:

library(glmmTMB)
library(rjags)

Owls <- transform(Owls,Nest=reorder(Nest,NegPerChick),NCalls=SiblingNegotiation,FT=FoodTreatment)

fit_nbinom2 <- glmmTMB(
  NCalls~ FT,
  data=Owls,
  dispformula = ~SexParent,
  family=nbinom2)


runjags::template.jags(formula = NCalls ~ FT + ArrivalTime
                       ,data = Owls
                       ,family = "negative binomial")

Модифицированная вывод из runjags :: template.jags выглядит следующим образом:

model{
  # In the BUGS/JAGS language we must use an explicit for loop:
  for (i in 1:N) {
    # These lines describe the response distribution and linear model terms:
    NCalls[i] ~ dpois(regression_fitted[i])
    regression_residual[i] <- NCalls[i] - regression_fitted[i]
    dispersion[i] ~ dgamma(k, k)
    regression_fitted[i] <- regression_mean[i] * dispersion[i]
    # Note: this formulation of a gamma-Poisson is exactly equivalent to a Negative Binomial
    log(regression_mean[i]) <-
      intercept + ArrivalTime_coefficient * ArrivalTime[i] + FT_effect[FT[i]]
  }

  # These lines give the prior distributions for the parameters to be estimated:
  k ~ dgamma(0.001, 0.001)
  intercept ~ dnorm(0, 10 ^ -6)
  ArrivalTime_coefficient ~ dnorm(0, 10 ^ -6)

  FT_effect[1] <- 0    # Factor level Deprived
  FT_effect[2] ~ dnorm(0, 10 ^ -6)    # Factor level Satiated
}

Как вы могли бы реализовать dispformula = ~SexParent часть из модели glmmTMB?

Справка ?glmmTMB::sigma.glmmTMB говорит это:

nbinom1

возвращает параметр чрезмерной дисперсии (обычно обозначается как альфа, как в Hardin and Hilbe (2007)): ​​так, что дисперсия равна mu (1 + альфа).

nbinom2

возвращает параметр избыточной дисперсии (обычно обозначаемый как тета или k); в отличие от большинства других семейств, большая тета соответствует более низкой дисперсии, которая равна mu (1 + mu / theta).

Из этого я попытался сделать параметр дисперсии зависимым от переменной SexParent. Вот так

model{
  # In the BUGS/JAGS language we must use an explicit for loop:
  for (i in 1:N) {
    # These lines describe the response distribution and linear model terms:
    NCalls[i] ~ dpois(regression_fitted[i])
    regression_residual[i] <- NCalls[i] - regression_fitted[i]
    dispersion[i] ~ dgamma(k[SexParent[i]], k[SexParent[i]])
    regression_fitted[i] <- regression_mean[i] * dispersion[i]
    # Note: this formulation of a gamma-Poisson is exactly equivalent to a Negative Binomial
    log(regression_mean[i]) <-
      intercept + ArrivalTime_coefficient * ArrivalTime[i] + FT_effect[FT[i]]
  }

  # These lines give the prior distributions for the parameters to be estimated:
  #k ~ dgamma(0.001, 0.001)
  intercept ~ dnorm(0, 10 ^ -6)
  ArrivalTime_coefficient ~ dnorm(0, 10 ^ -6)

  FT_effect[1] <- 0    # Factor level Deprived
  FT_effect[2] ~ dnorm(0, 10 ^ -6)    # Factor level Satiated
  for (i_sex in 1:N.parent){
    k[i_sex] ~ dgamma(0.001, 0.001)
  }  

}

Но иногда я сталкиваюсь с проблемами, потому что k должно быть> 0, или я получаю разные результаты для параметров, как в summary(fit_nbinom2). Я просто не знаю, как воспроизвести эти результаты.

Правильно ли я смоделировал дисперсию? У меня такое чувство, что я что-то упускаю. Было бы здорово, если бы кто-то мог помочь.