在JAGS中为R拟合多元dirlichet模型

问题描述

我正在尝试使用R中实现的JAGS将多元模型拟合到物种组成数据.我有3个物种相对丰度(界于[0,1]之间)的数据，其中两个具有相关性.这是生成相似数据的代码.

I am trying to fit a a multivariate model to species composition data using JAGS, implemented in R. I have data on 3 species relative abundances (bounded between [0,1]), two of which are correlated. Here is code to generate similar data.

#generate some correlated fractional composition data.
y1 <- runif(100,10,200)
y2 <- y1*1.5 + rnorm(100, sd = 5)
y3 <- runif(100,10,200)
total <- y1+y2+y3
y1 <- y1/(total)
y2 <- y2/(total)
y3 <- y3/(total)
y <- data.frame(y1,y2,y3)

y是我的三个因变量y1，y2和y3的data.frame.我想对这些数据使用仅截距模型，考虑使用dirlichet分布(β分布的多元扩展)来解释因变量之间的协方差.

y is a data.frame of my three dependent variables, y1,y2 and y3. I would like to fit an intercept only model to these data, accounting for the covariance among the dependent variables using a dirlichet disitribution, the multivariate extension of the beta distribution.

我有点卡住了.我可以使用R中的runjags包，使用beta分布精细法将其编码为单个因变量，如下所示:

I'm sort of stuck. I can code this up for a single dependent variable using a beta distribution fine using the runjags package in R as follows:

library(runjags)

#Write JAGS model, save as R object.
jags.model = "
model{
  # priors
  a0 ~ dnorm(0, .001)
  t0 ~ dnorm(0, .01)
  tau <- exp(t0)

  # likelihood for continuous component - predicted value on interval (0,1)
  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]) <- a0
  }
}
"

#generate JAGS data as list.
jags.data <- list(y = y1,
                  N = length(y1))

#Fit a JAGS model using run.jags
jags.out <- run.jags(jags.model,
                     data=jags.data,
                     adapt = 1000,
                     burnin = 1000,
                     sample = 2000,
                     n.chains=3,
                     monitor=c('a0'))

Mu问题是:如何使用R中实现的JAGS中的dirlichet分布将其扩展到多变量情况?如果我们可以考虑y矩阵中的协方差，那么将有好处.

Mu question is: how can I extend this to the multivariate case, using the dirlichet distribution in JAGS, implemented in R? Bonus if we can account for covariance in the y matrix.

推荐答案

使用原始问题中提供的伪数据，这是JAGS中的直接解决方案.

Here is a straightforward solution in JAGS, using the pseudo-data provided in the original question.

JAGS数据对象:

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

JAGS模型作为R对象jags.model:

JAGS model as the R object jags.model:

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

    #implement dirlichet
    for(i in 1:N){
    for(j in 1:N.spp){
         log(a0[i,j]) <- m0[j] ## eventually add linear model here
      }
     y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp])
    }

} #close model loop.
"

继续进行模型拟合，监视m0截距参数，该参数是一个向量，表示每个物种组的截距.

Go ahead and fit the model, monitoring the m0 intercept parameter, which is a vector, an intercept for each species group.

#Fit a JAGS model using run.jags
jags.out <- run.jags(jags.model,
                     data=jags.data,
                     adapt = 100,
                     burnin = 100,
                     sample = 200,
                     n.chains=3,
                     monitor=c('m0'))

这篇关于在JAGS中为R拟合多元dirlichet模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！