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)) 

나는 기본적으로 사이트 당 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) 

계층 접근 방식의 장점은 사이트의 효과가 지금에 변경 될 수 있다는 점이다 예를 들면, 무작위 효과 또는 무엇이든 통합하십시오.

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) 

:

매트

---

편집 랜덤 효과

---

다음 코드 (X2)에 대응하는 고정 효과와 함께 사이트의 임의의 효과를 추가하는 예를 추가

이것은 응용 프로그램에 적합한 모델 일 수도 있고 그렇지 않을 수도 있습니다. 이는 단지 예일뿐입니다. 이 경우 sigma.site는 데이터 시뮬레이션에 포함되지 않았기 때문에 매우 작다고 평가됩니다.