在 R 的 JAGS 中指定层次模型

Question

我有一些相关变量 y 的数据，可以将其建模为协变量 x1 和 x2 的函数。 y 和 x1 在 "plot" 级别观察到，并且 x2 在 "site" 级别观察到。地块分层嵌套在站点内。这是 y 的 100 个观察结果以及相关的协变量数据。

#generate covariate data at plot and site scales.
x1 <- runif(100,0,1)  #100 plot level observations of x1
x2 <- runif(10,10,20) #10  site level observations of x2

#generate site values - in this case characters A:J
site_1 <- LETTERS[sort(rep(seq(1,10, by = 1),10))]
site_2 <- LETTERS[sort(seq(1,10, by = 1))]

#put together site level data - 10 observations for 10 sites.
site_data <- data.frame(site_2,x2)
colnames(site_data) <- c('site','x2')

#put together plot level data - 100 observations across 10 sites
plot_data <- data.frame(site_1,x1)
colnames(plot_data) <- c('site','x1')
plot_data <- merge(plot_data,site_data, all.x=T) #merge in site level data.

#pick parameter values.
b1 <- 10
b2 <- 0.2

#y is a function of the plot level covariate x1 and the site level covariate x2.
plot_data$y <- b1*plot_data$x1 + b2*plot_data$x2 + rnorm(100)

#check that the model fits. it does.
summary(lm(y ~ x1 + x2, data = plot_data))

我可以将 y 建模为 x1 和 x2 的函数，使用数据框 plot_data 在 jags 中没有问题，它基本上复制了 x2 每个站点10次。

然而，我真正想做的是分层适合模型，这样 y[i] ~ x1[i] + x2[j]，其中 [i] 表示地块级观察，[j] 索引站点。我如何修改下面的 JAGS 代码来执行此操作？

#fit a JAGS model
jags.model = "
model{
# priors
b1 ~ dnorm(0, .001)
b2 ~ dnorm(0, .001)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 100)

#normal model
for (i in 1:N){
y[i] ~ dnorm(y.hat[i], tau)
y.hat[i] <- b1*x1[i] + b2*x2[i]
}


} #end model
"

#setup jags data as a list
jd <- list(y=plot_data$y, x1=plot_data$x1, x2=plot_data$x2, N=length(plot_data$y))

library(runjags)
#run jags model
jags.out <- run.jags(jags.model,
                     data=jd,
                     adapt = 1000,
                     burnin = 1000,
                     sample = 2000,
                     n.chains=3,
                     monitor=c('b1', 'b2'))
summary(jags.out)

Answer 1

您只需要一个响应级别的索引向量，该级别对应于站点内的唯一级别（如果将其编码为一个因素，这将是最简单的）。以下模型与您已有的模型完全相同：

jags.model = "
model{
    # priors
    b1 ~ dnorm(0, .001)
    b2 ~ dnorm(0, .001)
    tau <- pow(sigma, -2)
    sigma ~ dunif(0, 100)

    # Response:
    for (i in 1:N){
        y[i] ~ dnorm(y.hat[i], tau)
        y.hat[i] <- b1*x1[i] + site_effect[plot_site[i]]
    }

    # Effect of site:
    for (s in 1:S){
        site_effect[s] <- b2 * x2_site[site_site[s]]
    }

}
"
# Ensure the site is coded as a factor with the same levels in both data frames:
plot_data$site <- factor(plot_data$site)
site_data$site <- factor(site_data$site, levels=levels(plot_data$site))

#setup jags data as a list
jd <- list(y=plot_data$y, x1=plot_data$x1, plot_site=plot_data$site, 
            site_site=site_data$site, x2_site=site_data$x2, 
            N=length(plot_data$y), S=nrow(site_data))

library(runjags)
#run jags model
jags.out <- run.jags(jags.model,
                     data=jd,
                     adapt = 1000,
                     burnin = 1000,
                     sample = 2000,
                     n.chains=3,
                     monitor=c('b1', 'b2'))
summary(jags.out)

分层方法的优点是网站的效果现在可以修改为例如结合随机效应或其他。

马特

编辑以添加随机效果示例

以下代码添加了站点的随机效果以及对应于 x2 的固定效果：

jags.model = "
model{
    # priors
    b1 ~ dnorm(0, .001)
    b2 ~ dnorm(0, .001)
    tau <- pow(sigma, -2)
    sigma ~ dunif(0, 100)
    tau.site <- pow(sigma.site, -2)
    sigma.site ~ dunif(0, 100)

    # Response:
    for (i in 1:N){
        y[i] ~ dnorm(y.hat[i], tau)
        y.hat[i] <- b1*x1[i] + site_effect[plot_site[i]]
    }

    # Effect of site (fixed and random effects):
    for (s in 1:S){
        site_effect[s] <- b2 * x2_site[site_site[s]] + random[site_site[s]]
        random[site_site[s]] ~ dnorm(0, tau.site)
    }

}
"
# Ensure the site is coded as a factor with the same levels in both data frames:
plot_data$site <- factor(plot_data$site)
site_data$site <- factor(site_data$site, levels=levels(plot_data$site))

#setup jags data as a list
jd <- list(y=plot_data$y, x1=plot_data$x1, plot_site=plot_data$site, 
            site_site=site_data$site, x2_site=site_data$x2, 
            N=length(plot_data$y), S=nrow(site_data))

library(runjags)
#run jags model
jags.out <- run.jags(jags.model,
                     data=jd,
                     adapt = 1000,
                     burnin = 1000,
                     sample = 2000,
                     n.chains=3,
                     monitor=c('b1', 'b2', 'sigma.site', 'sigma'))
summary(jags.out)

这可能是也可能不是适合您的应用程序的模型 - 这只是一个示例。在这种情况下，估计 sigma.site 非常小，因为它没有出现在数据模拟中。

在 R 的 JAGS 中指定层次模型

specifying a hierarchical model in JAGS for R

r

jags

runjags

rjags