JAGSの階層モデルをR

Question

に指定しています。依存関係のあるベクトルyには、共変量x1とx2の関数としてモデル化できるデータがあります。 yおよびx1が「プロット」レベルで観察され、x2が「サイト」レベルで観察される。プロットはサイト内で階層的にネストされます。関連する共変量データを伴う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)) 

Iは基本的部位当たりx2 10倍のサイトレベルの観測を複製データフレームplot_dataを用いジャグにx1とx2問題なしの関数としてyをモデル化することができます。私は本当にしかし、やりたい何

は階層的にな[i]サイトプロットレベルの観察と[j]インデックスを示し y[i] ~ x1[i] + x2[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はデータシミュレーションでは機能しなかったため、非常に小さいと見積もられます。