r - 将参数向量乘以 JAGS 中的自变量矩阵

Question

我正在使用 dirlichet 分布在 JAGS 中拟合多变量模型。我有一个包含 3 个物种比例丰度的矩阵 y。

#generate 3 columns of species proprotional abundance data
y    <- matrix(ncol = 3, nrow = 100)
y[,] <- abs(rnorm(length(y)))
for(i in 1:nrow(y)){
  y[i,] <- y[i,] / sum(y[i,])
}

我有一个预测值矩阵 x，第一个是截距。

#generate 2 columns of predictors and an intercept
x <- matrix(ncol = 2, nrow = 100)
x[,] <- rnorm(length(x), mean = 20, sd = 4)
x <- cbind(rep(1,nrow(x)),x)

我指定了一个多变量 jags 模型，jags.model:

jags.model = "
model {
    #setup parameter priors for each species * predictor combination.
    for(j in 1:N.spp){
      for(k in 1:N.preds){
        m[k,j] ~ dgamma(1.0E-3, 1.0E-3)
      }
    }

    #go ahead and fit means of species abundances as a linear combination of predictor and parameters.
    for(i in 1:N){
        for(j in 1:N.spp){
             log(a0[i,j]) <- m[,j] * x[i,]
           }
    y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp]) 
    }

} #close model loop.
"

我设置了 JAGS 数据对象，jags.data:

jags.data <- list(y = as.matrix(y), x = as.matrix(x),
                  N.spp = ncol(y), N.preds = ncol(x), N = nrow(y))

我使用 R 中的 runjags 包拟合 JAGS 模型。

jags.out <- runjags::run.jags(jags.model,
                              data=jags.data,
                              adapt = 100,
                              burnin = 200,
                              sample = 400,
                              n.chains=3,
                              monitor=c('m'))

我收到以下错误:

Error: The following error occured when compiling and adapting the model using rjags:
 Error in rjags::jags.model(model, data = dataenv, n.chains = length(runjags.object$end.state),  : 
  RUNTIME ERROR:
Invalid vector argument to exp

我在这里做错了什么？作为引用，通过预测组合拼出每个参数仍然很合适:

jags.model = "
model {
  #setup parameter priors for each species * predictor combination.
    for(j in 1:N.spp){
      for(k in 1:N.preds){
        m[k,j] ~ dgamma(1.0E-3, 1.0E-3)
      }
    }

  #go ahead and fit means of species abundances as a linear combination of predictor and parameters.
  for(i in 1:N){
    for(j in 1:N.spp){
      log(a0[i,j]) <- m[1,j] * x[i,1] + m[2,j] * x[i,2] + m[3,j] * x[i,3]
    }
    y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp]) 
  }

} #close model loop.
"

Answer 1

这个问题的解决方案是采用点积，或者 JAGS 中的内积。更改行:

log(a0[i,j]) <- m[,j] * x[i,]

到:

log(a0[i,j]) <- inprod(m[,j] , x[i,])

一切都应该工作正常。完整模型如下。

jags.model = "
model {
    #setup parameter priors for each species * predictor combination.
    for(j in 1:N.spp){
      for(k in 1:N.preds){
        m[k,j] ~ dgamma(1.0E-3, 1.0E-3)
      }
    }

    #go ahead and fit means of species abundances as a linear combination of predictor and parameters.
    for(i in 1:N){
        for(j in 1:N.spp){
             log(a0[i,j]) <- inprod(m[,j] , x[i,])
           }
    y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp]) 
    }

} #close model loop.
"