title: specifying a hierarchical model in JAGS for R
body: <p>I have data on some dependent varaible <code>y</code> that can be modeled as a function of covariates <code>x1</code> and <code>x2</code>. <code>y</code> and <code>x1</code> are observed at the "plot" level, and and <code>x2</code> is observed at the "site" level. Plot is nested within site, hierarchically. Here are 100 observations of <code>y</code> with associated covariate data.</p>
  
  <pre><code>#generate covariate data at plot and site scales.
  x1 &lt;- runif(100,0,1)  #100 plot level observations of x1
  x2 &lt;- runif(10,10,20) #10  site level observations of x2
  
  #generate site values - in this case characters A:J
  site_1 &lt;- LETTERS[sort(rep(seq(1,10, by = 1),10))]
  site_2 &lt;- LETTERS[sort(seq(1,10, by = 1))]
  
  #put together site level data - 10 observations for 10 sites.
  site_data &lt;- data.frame(site_2,x2)
  colnames(site_data) &lt;- c('site','x2')
  
  #put together plot level data - 100 observations across 10 sites
  plot_data &lt;- data.frame(site_1,x1)
  colnames(plot_data) &lt;- c('site','x1')
  plot_data &lt;- merge(plot_data,site_data, all.x=T) #merge in site level data.
  
  #pick parameter values.
  b1 &lt;- 10
  b2 &lt;- 0.2
  
  #y is a function of the plot level covariate x1 and the site level covariate x2.
  plot_data$y &lt;- 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))
  </code></pre>
  
  <p>I can model <code>y</code> as a function of <code>x1</code> and <code>x2</code> no problem in jags using the data frame <code>plot_data</code> which basically replicates the site level observations of  <code>x2</code> 10 times per site. </p>
  
  <p>What I would really like to do however is fit the model <em>hierarchically</em> such that
  <code>y[i] ~ x1[i] + x2[j]</code>, where <code>[i]</code> indicates the plot level observation and <code>[j]</code> indexes the site. How can I modify the below JAGS code to do this?</p>
  
  <pre><code>#fit a JAGS model
  jags.model = "
  model{
  # priors
  b1 ~ dnorm(0, .001)
  b2 ~ dnorm(0, .001)
  tau &lt;- pow(sigma, -2)
  sigma ~ dunif(0, 100)
  
  #normal model
  for (i in 1:N){
  y[i] ~ dnorm(y.hat[i], tau)
  y.hat[i] &lt;- b1*x1[i] + b2*x2[i]
  }
  
  
  } #end model
  "
  
  #setup jags data as a list
  jd &lt;- 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 &lt;- run.jags(jags.model,
                       data=jd,
                       adapt = 1000,
                       burnin = 1000,
                       sample = 2000,
                       n.chains=3,
                       monitor=c('b1', 'b2'))
  summary(jags.out)
  </code></pre>