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

问题描述

我想通过R使用JAGS在Tweedie分布变量上运行模型。我知道JAGS没有标准的Tweedie分布，但是可以将其指定为复合Gamma / Poisson。不幸的是，我无法弄清楚如何用JAGS编写代码。我基于从各种来源清除的代码编写了以下内容，以尝试从Tweedie随机变量中恢复均值，幂和phi参数。由于y上的父值无效，因此目前无法运行，大概是因为y [i]出现在表达式的右侧和左侧。这就像在源代码中写的一样，但是我显然在滥用它。我们将不胜感激任何有关如何正确指定此分布的指针，并且可能会更广泛地使用它，因为我无法找到有关如何在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,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
)

解决方法

我能够使它起作用。是否“有效”似乎更是一个悬而未决的问题，但参数的后验方法非常接近真实值。我使用runjags在R中进行了此操作，但是所有这些步骤都将在JAGS中独立于R进行。

首先，我生成了数据，并将形状矩阵放入具有一列NA值的数据中，以便可以在每次模拟迭代中将其覆盖。我还添加了一个名为yind的变量，如果y == 0为2，否则为1。这将用作索引shape和rate矩阵的值。在数据中执行此操作可能比让JAGS在每次迭代中针对开头固定的equals()的每个值在每次迭代中多次执行y函数更有效。

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

jags r r