1

可以说我有一些关于 3 个物种的相对丰度的数据。物种 1 与年平均温度 ( mat) 呈负相关。相对丰度数据存储在矩阵中r.spp.y：

#Simulate some species count data.
y1 <- round(rnorm(100, 100, 5)) 
mat <- rnorm(length(y1), 10, 3)
y1 <- y1 + round(mat*-3)
y2 <- round(rnorm(100, 100, 5)) 
y3 <- round(rnorm(100, 100, 5)) 
spp.y <- as.matrix(data.frame(y1,y2,y3))
r.spp.y <- spp.y / rowSums(spp.y)

我可以在 JAGS 中一次对每个物种进行 beta 回归，以表明mat物种 1 的相对丰度之间存在负相关，但mat使用 JAGS 时，物种 2 和 3 的相对丰度之间存在正相关：

library(runjags)
beta.model = "
model{
  #priors
  m0 ~ dnorm(0, 1.0E-3) #intercept prior
  m1 ~ dnorm(0, 1.0E-3) #      mat prior
  t0 ~ dnorm(0, .01)
  tau <- exp(t0) 

  #drop beta
  for (i in 1:N){
  y[i] ~ dbeta(p[i], q[i])
  p[i] <- mu[i] * tau
  q[i] <- (1 - mu[i]) * tau
  logit(mu[i]) <- m0 + m1 * mat[i]
  }

}#close model loop.
"
output <- list()
for(i in 1:ncol(r.spp.y)){
  beta.data <- list(y = r.spp.y[,i], mat = mat, N = nrow(r.spp.y))
  jags.beta <- run.jags(beta.model,
                        data=beta.data,
                        adapt = 200,
                        burnin = 1000,
                        sample = 1000,
                        n.chains = 3,
                        monitor = c('m0','m1'))
  output[[i]] <- summary(jags.beta)
}
names(output) <- c('species 1','species 2','species 3')

这表明每个都有一个截距值和一些关系mat。具体来说，物种 1 与 物种 2 和物种 3 呈负相关，mat而物种 2 和 3 表现出正相关，mat这反映在m1参数值中。

$`species 1`
       Lower95      Median     Upper95        Mean          SD        Mode        MCerr MC%ofSD SSeff     AC.10     psrf
m0 -0.71012252 -0.66164521 -0.61324227 -0.66175203 0.025504290 -0.65828829 0.0025395302    10.0   101 0.4951512 1.005070
m1 -0.04608964 -0.04139645 -0.03671938 -0.04138209 0.002453942 -0.04148652 0.0002420001     9.9   103 0.5013018 1.006634

$`species 2`
       Lower95      Median     Upper95        Mean          SD        Mode        MCerr MC%ofSD SSeff     AC.10     psrf
m0 -0.74993189 -0.70794571 -0.67146076 -0.70773078 0.019593399 -0.70704999 0.0018992840     9.7   106 0.4913926 1.006231
m1  0.01479929  0.01849042  0.02212619  0.01844102 0.001845884  0.01825645 0.0001812411     9.8   104 0.4964613 1.006281

$`species 3`
       Lower95      Median     Upper95        Mean          SD        Mode       MCerr MC%ofSD SSeff     AC.10     psrf
m0 -0.71786559 -0.67713045 -0.63342213 -0.67726218 0.021507927 -0.67767494 0.001908073     8.9   127 0.4398105 1.024932
m1  0.01160666  0.01560139  0.01958145  0.01561193 0.002017519  0.01576153 0.000187649     9.3   116 0.4513887 1.026890

我想使用多元模型同时拟合所有丰度。这可以使用 dirlichet 分布来完成，它是 beta 分布的多元泛化。为此，我在 JAGS 中设置了一个 dirlichet 模型，类似于上面的 beta：

dirlichet.model = "
model {
  #setup priors for each species
  for(j in 1:N.spp){
    m0[j] ~ dnorm(0, 1.0E-3) #intercept prior
    m1[j] ~ dnorm(0, 1.0E-3) #      mat prior
  }

  #implement dirlichet
  for(i in 1:N){
    for(j in 1:N.spp){
      logit(a0[i,j]) <- m0[j] + m1[j] * mat[i]
  }
  y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp]) 
  }

} #close model loop.
"

然后我设置了一个稍微不同的数据对象，并运行模型：

jags.data <- list(y = r.spp.y,mat = mat,N = nrow(r.spp.y), N.spp = ncol(r.spp.y))
jags.out <- run.jags(dirlichet.model,
                     data=jags.data,
                     adapt = 200,
                     burnin = 2000,
                     sample = 2000,
                     n.chains=3,
                     monitor=c('m0','m1'))
summary(jags.out)

但是，返回的m0和m1参数与循环 beta 模型返回的参数完全不同。似乎所有参数都反映了提供的（平坦）先验分布，并且数据并未以任何方式限制参数。该模型无法捕捉mat物种 1 之间的负相关关系和相对丰度。它甚至似乎根本没有限制数据：

> summary(jags.out)
         Lower95    Median  Upper95      Mean       SD      Mode     MCerr MC%ofSD SSeff         AC.10     psrf
m0[1] -49.005761  6.181091 63.01256  6.393293 28.43734  6.628375 0.5359853     1.9  2815 -0.0018363868 1.000590
m0[2] -45.469227  6.723323 67.85340  6.978617 28.52249  5.924171 0.5417058     1.9  2772  0.0129357542 1.000718
m0[3] -49.227870  6.233448 61.60683  6.599494 28.43137  4.742497 0.5279870     1.9  2900 -0.0057574414 1.000546
m1[1]  -1.869337 24.268826 64.16340 27.391725 19.13417 20.322548 0.4572891     2.4  1751  0.0109859218 1.006568
m1[2]  -1.465451 23.725515 64.21979 27.071232 18.98391 15.823296 0.4294849     2.3  1954 -0.0005437107 1.002363
m1[3]  -1.907332 24.031360 62.18261 27.017599 18.57743 18.597688 0.4290431     2.3  1875 -0.0086604393 1.001171

值得注意的是，如果我拟合一个仅截距的 dirlichet 模型（没有mat协变量），它可以很好地捕获相对丰度，但前提是我使用 log-link，而不是 beta 案例中使用的 logit-link（例如在这里）。扩展 log-link 案例并添加mat协变量并不能解决问题。我在这里想念什么？

4

0 回答 0