如何在JAGS/BUGS中编写Tweedie发行版？

Question

我想使用JAGS到R在Tweedie分布变量上运行一个模型。我知道JAGS没有Tweedie分布作为标准，但可以指定一个作为复合Gamma/Poisson。不幸的是，我不知道如何在JAGS中编写它。我基于从各种来源收集的代码编写了以下代码，以简单地尝试并从Tweedie随机变量中恢复平均值、幂和phi参数。它目前不会运行，因为y上的父值无效，大概是因为yi出现在表达式的右侧和左侧。这是在源代码中写的，但我显然是误用了它。任何关于如何正确指定这个发行版的建议都会非常感谢，而且可能会有更广泛的用途，因为我还没有找到任何关于如何在JAGS中设置Tweedie模型的简单示例代码。

y = mgcv::rTweedie(mu=rep(2,100),p=1.33,phi=1)

jags_data = list(y=y,n=length(y))

jags_model = 
 "model{
    
    for (i in 1:n) {
      lambda[i] <- pow(mu,2-p)/(phi *(2-p))
      num[i] ~ dpois(lambda[i])
      shape[i,1] <- num[i]*((2-p)/(p-1))
      rate[i,1] <- 1/(phi*(p-1)*pow(mu,p-1))
      shape[i,2] <- 1
      rate[i,2] <- exp(-lambda[i])
      # Takes shape/rate parameter 1 if y > 0 and 2 if y = 0
      y[i] ~ dgamma(shape[i,1+equals(y[i],0)],rate[i,1+equals(y[i],0)]) 
    }
    
    mu    ~ dunif(0,100)
    p     ~ dunif(1,2)  ## Tweedie power parameter
    phi   ~ dunif(0,30) ## Dispersion parameter
    
  }
  
  "

model_file = tempfile(fileext = 'txt')
writeLines(jags_model,model_file)

jm = rjags::jags.model(
  file = model_file,
  data = jags_data,
  n.chains = 3,
  n.adapt = 1500
)

Answer 1

我能够让它发挥作用。它是否“有效”似乎更像是一个悬而未决的问题，但参数的后验均值非常接近真实值。我用runjags在R中做了这件事，但这里的所有部分都是在独立于R的JAGS中完成的。

首先，我生成了数据，并将形状矩阵放入具有一列NA值的数据中，以便它们可以在每次模拟迭代中被覆盖。我还添加了一个名为yind的变量，如果为y == 0，则为2，否则为1。这将用作索引shape和rate矩阵的值。在数据中这样做可能更有效，而不是让JAGS在每次迭代中对y的每个值进行几次求值和equals()函数，这些值在开始时是固定的。

y = mgcv::rTweedie(mu=rep(2,100),p=1.33,phi=1)
shape_mat <- matrix(NA, nrow=length(y), ncol=2)
shape_mat[,2] <- 1
jags_data = list(y=y,n=length(y), yind = 2-(y > 0), shape = shape_mat)

接下来，除了形状矩阵的处理方式和添加的yind之外，模型是相同的。

jags_model = 
  "model{
    
    for (i in 1:n) {
      lambda[i] <- pow(mu,2-p)/(phi *(2-p))
      num[i] ~ dpois(lambda[i])
      shape[i,1] <- num[i]*((2-p)/(p-1))
      rate[i,1] <- 1/(phi*(p-1)*pow(mu,p-1))
      ## moved to data
      # shape[i,2] <- 1
      rate[i,2] <- exp(-lambda[i])
      # Takes shape/rate parameter 1 if y > 0 and 2 if y = 0
      y[i] ~ dgamma(shape[i,yind[i]],rate[i,yind[i]]) 
    }
    
    mu    ~ dunif(0,100)
    p     ~ dunif(1,2)  ## Tweedie power parameter
    phi   ~ dunif(0,30) ## Dispersion parameter
    
  }
  "

无效的父值可能来自初始值。num泊松绘制必须与lambda保持一致，这是mu、p和phi的函数。因此，初始值对于确保所有这些值都是一致的非常重要。我将三个模型参数设置为下面的值，然后计算lambda。根据这三个参数的值，我将num设置为与lambda最接近的整数值。

inits <- list(mu = 5, p=1.5, phi=5,  
              num = rep(1, length(y)))

接下来，我运行模型并进行总结：

library(runjags)
out <- run.jags(model=jags_model, data = jags_data, monitor=c("mu", "p", "phi"), burnin=5000, sample=10000, 
                inits = list(inits, inits), n.chains = 2, keep.jags.files=TRUE)
summary(out)
# Lower95   Median Upper95     Mean         SD Mode       MCerr MC%ofSD SSeff        AC.10     psrf
# mu  1.606100 1.921855 2.24793 1.927986 0.16489395   NA 0.001534324     0.9 11550 -0.006514202 1.000242
# p   1.252230 1.377415 1.50807 1.379773 0.06563517   NA 0.002049511     3.1  1026  0.354243967 1.013941
# phi 0.871354 1.098075 1.36870 1.109978 0.12874006   NA 0.003880090     3.0  1101  0.322202621 1.004731