title: WinBUGS Code for Bayesian hierarchical model
body: <p>I’m reading this <a href="http://www-bcf.usc.edu/~qianghua/papers/P18.pdf" rel="nofollow">paper</a> (which aims to model the nanowires (NW) growth using a Bayesian hierarchical approach). In page 7, the author proposed a model to describe the growth of nanowires. I’m trying to write a WinBUGS code (Please see below), but I’m not sure whether it’s correct. </p>
  
  <p>I quoted this from the paper as a summary: </p>
  
  <blockquote>
    <p>The case studies will demonstrate the procedure of hierarchical modeling and estimation of the NW growth process under uncertainties. ... The data therein were collected over time (t =15 s, 30 s, 180 s, 900 s) under six growth conditions (T = 365 ◦C, 380 ◦C, 400 ◦C, 420 ◦C, 430 ◦C, 440 ◦C). The first four conditions are used for model building. We do not consider the two high-temperature conditions because there are no observations at times 180 s and 900 s. “</p>
  </blockquote>
  
  <p>So we’re dealing with three variables, including Time (in seconds), Temperature (in C), and Length of nanowires (in cm). There are four levels of Temperature (365 ◦C, 380 ◦C, 400 ◦C, and 420 ◦C), and four time points (15 s, 30 s, 180 s, 900 s). I don’t have the exact values of length under each condition, so I generated some values by looking at Figure 6 in the paper.</p>
  
  <p>Here is the WinBUGS code. It runs perfectly, but I’m not sure if it’s correct. I appreciate any input.  Thanks!</p>
  
  <pre><code>model
  {
      for(j in 1 : 4) {
          for(i in 1 : 4) {
              X[i,j] ~ dnorm(mu[i,j], tau)
              mu[i,j] &lt;- (step(t0[j] - t[i]))  * ((a1[j]*(exp(-alpha2[j]/t[i])))) + (step(t[i] - t0[j])) * ((a3[j]*t[i]) + (a1[j]*(exp(-alpha2[j]/t0[j]))) -a3[j]*t0[j])
          }
  
          t0[j] ~ dunif(1,40)
          alpha2[j] ~ dunif(1,20)
  
          a1[j]~ dnorm(mua1[j], tau1)
          a3[j]~ dnorm(mua3[j], tau3)
  
          alpha1[j] &lt;- exp(a11)
          alpha3[j] &lt;- exp(a33)
  
          log(mua1[j]) &lt;- log(alpha1[j]) - (230/(T[j]*1.3806503*pow(10, -23)))
          log(mua3[j]) &lt;- log(alpha3[j]) - (230/(T[j]*1.3806503*pow(10, -23)))
      }
  
      a11~ dnorm(0,0.001)
      a33~ dnorm(0, 0.001)
  
      tau1 &lt;- pow(sigma1, -2)
      sigma1 ~ dunif(0,10)
  
      tau3 &lt;- pow(sigma3, -2)
      sigma3 ~ dunif(0,10)
  
      tau &lt;- pow(sigma, -2)
      sigma ~ dunif(0,20)
  }
  
  list(X= structure(.Data= c(0.00011E+00, 0.00022E+00, 2.00000E-02, 3.50000E-02, 7.00000E-03, 7.00000E-03, 3.00000E-01, 3.00000E-01, 3.50000E-02, 7.00000E-02, 2.50000E-01, 5.00000E-01, 2.00000E-01, 2.25000E-01, 5.50000E-01, 2.30000E+00), .Dim=c(4, 4)), T=c(3.65000E+02, 3.80000E+02, 4.00000E+02, 4.20000E+02), t=c(1.50000E+01, 3.00000E+01, 1.80000E+02, 9.00000E+02))
  </code></pre>