Randomised Complete Block Anova (JAGS)

This tutorial will focus on the use of Bayesian estimation to fit simple linear regression models. BUGS (Bayesian inference Using Gibbs Sampling) is an algorithm and supporting language (resembling R) dedicated to performing the Gibbs sampling implementation of Markov Chain Monte Carlo (MCMC) method. Dialects of the BUGS language are implemented within three main projects:

  1. OpenBUGS - written in component pascal.

  2. JAGS - (Just Another Gibbs Sampler) - written in C++.

  3. STAN - a dedicated Bayesian modelling framework written in C++ and implementing Hamiltonian MCMC samplers.

Whilst the above programs can be used stand-alone, they do offer the rich data pre-processing and graphical capabilities of R, and thus, they are best accessed from within R itself. As such there are multiple packages devoted to interfacing with the various software implementations:

  • R2OpenBUGS - interfaces with OpenBUGS

  • R2jags - interfaces with JAGS

  • rstan - interfaces with STAN

This tutorial will demonstrate how to fit models in JAGS (Plummer (2004)) using the package R2jags (Su et al. (2015)) as interface, which also requires to load some other packages.

Overview

Introduction

In the previous tutorial (nested ANOVA), we introduced the concept of employing sub-replicates that are nested within the main treatment levels as a means of absorbing some of the unexplained variability that would otherwise arise from designs in which sampling units are selected from amongst highly heterogeneous conditions. Such (nested) designs are useful in circumstances where the levels of the main treatment (such as burnt and un-burnt sites) occur at a much larger temporal or spatial scale than the experimental/sampling units (e.g. vegetation monitoring quadrats). For circumstances in which the main treatments can be applied (or naturally occur) at the same scale as the sampling units (such as whether a stream rock is enclosed by a fish proof fence or not), an alternative design is available. In this design (randomised complete block design), each of the levels of the main treatment factor are grouped (blocked) together (in space and/or time) and therefore, whilst the conditions between the groups (referred to as “blocks”) might vary substantially, the conditions under which each of the levels of the treatment are tested within any given block are far more homogeneous.

If any differences between blocks (due to the heterogeneity) can account for some of the total variability between the sampling units (thereby reducing the amount of variability that the main treatment(s) failed to explain), then the main test of treatment effects will be more powerful/sensitive. As an simple example of a randomised complete block (RCB) design, consider an investigation into the roles of different organism scales (microbial, macro invertebrate and vertebrate) on the breakdown of leaf debris packs within streams. An experiment could consist of four treatment levels - leaf packs protected by fish-proof mesh, leaf packs protected by fine macro invertebrate exclusion mesh, leaf packs protected by dissolving antibacterial tablets, and leaf packs relatively unprotected as controls. As an acknowledgement that there are many other unmeasured factors that could influence leaf pack breakdown (such as flow velocity, light levels, etc) and that these are likely to vary substantially throughout a stream, the treatments are to be arranged into groups or “blocks” (each containing a single control, microbial, macro invertebrate and fish protected leaf pack). Blocks of treatment sets are then secured in locations haphazardly selected throughout a particular reach of stream. Importantly, the arrangement of treatments in each block must be randomized to prevent the introduction of some systematic bias - such as light angle, current direction etc.

Blocking does however come at a cost. The blocks absorb both unexplained variability as well as degrees of freedom from the residuals. Consequently, if the amount of the total unexplained variation that is absorbed by the blocks is not sufficiently large enough to offset the reduction in degrees of freedom (which may result from either less than expected heterogeneity, or due to the scale at which the blocks are established being inappropriate to explain much of the variation), for a given number of sampling units (leaf packs), the tests of main treatment effects will suffer power reductions. Treatments can also be applied sequentially or repeatedly at the scale of the entire block, such that at any single time, only a single treatment level is being applied (see the lower two sub-figures above). Such designs are called repeated measures. A repeated measures ANOVA is to an single factor ANOVA as a paired t-test is to a independent samples t-test. One example of a repeated measures analysis might be an investigation into the effects of a five different diet drugs (four doses and a placebo) on the food intake of lab rats. Each of the rats (“subjects”) is subject to each of the four drugs (within subject effects) which are administered in a random order. In another example, temporal recovery responses of sharks to bi-catch entanglement stresses might be simulated by analyzing blood samples collected from captive sharks (subjects) every half hour for three hours following a stress inducing restraint. This repeated measures design allows the anticipated variability in stress tolerances between individual sharks to be accounted for in the analysis (so as to permit more powerful test of the main treatments). Furthermore, by performing repeated measures on the same subjects, repeated measures designs reduce the number of subjects required for the investigation. Essentially, this is a randomised complete block design except that the within subject (block) effect (e.g. time since stress exposure) cannot be randomised.

To suppress contamination effects resulting from the proximity of treatment sampling units within a block, units should be adequately spaced in time and space. For example, the leaf packs should not be so close to one another that the control packs are effected by the antibacterial tablets and there should be sufficient recovery time between subsequent drug administrations. In addition, the order or arrangement of treatments within the blocks must be randomized so as to prevent both confounding as well as computational complications. Whilst this is relatively straight forward for the classic randomized complete block design (such as the leaf packs in streams), it is logically not possible for repeated measures designs. Blocking factors are typically random factors that represent all the possible blocks that could be selected. As such, no individual block can truly be replicated. Randomised complete block and repeated measures designs can therefore also be thought of as un-replicated factorial designs in which there are two or more factors but that the interactions between the blocks and all the within block factors are not replicated.

Linear models

The linear models for two and three factor nested design are:

\[ y_{ij} = \mu + \beta_i + \alpha_j + \epsilon_{ij},\]

\[ y_{ijk} = \mu + \beta_i + \alpha_j + \gamma_k + (\beta\alpha)_{ij} + (\beta\gamma)_{ik} + (\alpha\gamma)_{jk} + (\alpha\beta\gamma)_{ijk} + \epsilon_{ijk}, \;\;\; \text{(Model 1)}\]

\[ y_{ijk} = \mu + \beta_i + \alpha_j + \gamma_k + (\alpha\gamma)_{jk} + \epsilon_{ijk}, \;\;\; \text{(Model 2)},\]

where \(\mu\) is the overall mean, \(\beta\) is the effect of the Blocking Factor B (\(\sum \beta=0\)), \(\alpha\) and \(\gamma\) are the effects of withing block Factor A and Factor C, respectively, and \(\epsilon \sim N(0,\sigma^2)\) is the random unexplained or residual component.

Tests for the effects of blocks as well as effects within blocks assume that there are no interactions between blocks and the within block effects. That is, it is assumed that any effects are of similar nature within each of the blocks. Whilst this assumption may well hold for experiments that are able to consciously set the scale over which the blocking units are arranged, when designs utilize arbitrary or naturally occurring blocking units, the magnitude and even polarity of the main effects are likely to vary substantially between the blocks. The preferred (non-additive or “Model 1”) approach to un-replicated factorial analysis of some bio-statisticians is to include the block by within subject effect interactions (e.g. \(\beta\alpha\)). Whilst these interaction effects cannot be formally tested, they can be used as the denominators in F-ratio calculations of their respective main effects tests. Proponents argue that since these blocking interactions cannot be formally tested, there is no sound inferential basis for using these error terms separately. Alternatively, models can be fitted additively (“Model 2”) whereby all the block by within subject effect interactions are pooled into a single residual term (\(\epsilon\)). Although the latter approach is simpler, each of the within subject effects tests do assume that there are no interactions involving the blocks and that perhaps even more restrictively, that sphericity holds across the entire design.

Assumptions

As with other ANOVA designs, the reliability of hypothesis tests is dependent on the residuals being:

  • normally distributed. Boxplots using the appropriate scale of replication (reflecting the appropriate residuals/F-ratio denominator should be used to explore normality. Scale transformations are often useful.

  • equally varied. Boxplots and plots of means against variance (using the appropriate scale of replication) should be used to explore the spread of values. Residual plots should reveal no patterns. Scale transformations are often useful.

  • independent of one another. Although the observations within a block may not strictly be independent, provided the treatments are applied or ordered randomly within each block or subject, within block proximity effects on the residuals should be random across all blocks and thus the residuals should still be independent of one another. Nevertheless, it is important that experimental units within blocks are adequately spaced in space and time so as to suppress contamination or carryover effects.

Simple RCB

Data generation

Imagine we has designed an experiment in which we intend to measure a response (y) to one of treatments (three levels; “a1,” “a2” and “a3”). Unfortunately, the system that we intend to sample is spatially heterogeneous and thus will add a great deal of noise to the data that will make it difficult to detect a signal (impact of treatment). Thus in an attempt to constrain this variability you decide to apply a design (RCB) in which each of the treatments within each of 35 blocks dispersed randomly throughout the landscape. As this section is mainly about the generation of artificial data (and not specifically about what to do with the data), understanding the actual details are optional and can be safely skipped.

> library(plyr)
> set.seed(123)
> nTreat <- 3
> nBlock <- 35
> sigma <- 5
> sigma.block <- 12
> n <- nBlock*nTreat
> Block <- gl(nBlock, k=1)
> A <- gl(nTreat,k=1)
> dt <- expand.grid(A=A,Block=Block)
> #Xmat <- model.matrix(~Block + A + Block:A, data=dt)
> Xmat <- model.matrix(~-1+Block + A, data=dt)
> block.effects <- rnorm(n = nBlock, mean = 40, sd = sigma.block)
> A.effects <- c(30,40)
> all.effects <- c(block.effects,A.effects)
> lin.pred <- Xmat %*% all.effects
> 
> # OR
> Xmat <- cbind(model.matrix(~-1+Block,data=dt),model.matrix(~-1+A,data=dt))
> ## Sum to zero block effects
> block.effects <- rnorm(n = nBlock, mean = 0, sd = sigma.block)
> A.effects <- c(40,70,80)
> all.effects <- c(block.effects,A.effects)
> lin.pred <- Xmat %*% all.effects
> 
> 
> 
> ## the quadrat observations (within sites) are drawn from
> ## normal distributions with means according to the site means
> ## and standard deviations of 5
> y <- rnorm(n,lin.pred,sigma)
> data.rcb <- data.frame(y=y, expand.grid(A=A, Block=Block))
> head(data.rcb)  #print out the first six rows of the data set
         y A Block
1 45.80853 1     1
2 66.71784 2     1
3 93.29238 3     1
4 43.10101 1     2
5 73.20697 2     2
6 91.77487 3     2

Exploratory data analysis

Normality and Homogeneity of variance

> boxplot(y~A, data.rcb)

Conclusions:

  • there is no evidence that the response variable is consistently non-normal across all populations - each boxplot is approximately symmetrical.

  • there is no evidence that variance (as estimated by the height of the boxplots) differs between the five populations. . More importantly, there is no evidence of a relationship between mean and variance - the height of boxplots does not increase with increasing position along the \(y\)-axis. Hence it there is no evidence of non-homogeneity

Obvious violations could be addressed either by:

  • transform the scale of the response variables (to address normality, etc). Note transformations should be applied to the entire response variable (not just those populations that are skewed).

Block by within-Block interaction

> library(car)
> with(data.rcb, interaction.plot(A,Block,y))
> 
> #OR with ggplot
> library(ggplot2)

> ggplot(data.rcb, aes(y=y, x=A, group=Block,color=Block)) + geom_line() +
+   guides(color=guide_legend(ncol=3))

> 
> residualPlots(lm(y~Block+A, data.rcb))

           Test stat Pr(>|Test stat|)
Block                                
A                                    
Tukey test   -1.4163           0.1567
> 
> # the Tukey's non-additivity test by itself can be obtained via an internal function
> # within the car package
> car:::tukeyNonaddTest(lm(y~Block+A, data.rcb))
      Test     Pvalue 
-1.4163343  0.1566776 
> 
> # alternatively, there is also a Tukey's non-additivity test within the
> # asbio package
> library(asbio)
> with(data.rcb,tukey.add.test(y,A,Block))

Tukey's one df test for additivity 
F = 2.0060029   Denom df = 67    p-value = 0.1613102

Conclusions:

  • there is no visual or inferential evidence of any major interactions between Block and the within-Block effect (A). Any trends appear to be reasonably consistent between Blocks.

Model fitting

Full parameterisation

\[ y_{ijk} \sim N(\mu_{ij}, \sigma^2), \;\;\; \mu_{ij}=\beta_0 + \beta_i + \gamma_{j(i)}, \]

where \(\gamma_{ij)} \sim N(0, \sigma^2_B)\), \(\beta_0, \beta_i \sim N(0, 1000000)\), and \(\sigma^2, \sigma^2_B \sim \text{Cauchy(0, 25)}\). The full parameterisation, shows the effects parameterisation in which there is an intercept (\(\beta_0\)) and two treatment effects (\(\beta_i\), where \(i\) is \(1,2\)).

Matrix parameterisation

\[ y_{ijk} \sim N(\mu_{ij}, \sigma^2), \;\;\; \mu_{ij}=\boldsymbol \beta \boldsymbol X + \gamma_{j(i)}, \]

where \(\gamma_{ij} \sim N(0, \sigma^2_B)\), \(\boldsymbol \beta \sim MVN(0, 1000000)\), and \(\sigma^2, \sigma^2_B \sim \text{Cauchy(0, 25)}\). The full parameterisation, shows the effects parameterisation in which there is an intercept (\(\alpha_0\)) and two treatment effects (\(\beta_i\), where \(i\) is \(1,2\)). The matrix parameterisation is a compressed notation, In this parameterisation, there are three alpha parameters (one representing the mean of treatment a1, and the other two representing the treatment effects (differences between a2 and a1 and a3 and a1). In generating priors for each of these three alpha parameters, we could loop through each and define a non-informative normal prior to each (as in the Full parameterisation version). However, it turns out that it is more efficient (in terms of mixing and thus the number of necessary iterations) to define the priors from a multivariate normal distribution. This has as many means as there are parameters to estimate (\(3\)) and a \(3\times3\) matrix of zeros and \(100\) in the diagonals.

\[ \boldsymbol \mu = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \;\;\; \sigma^2 \sim \begin{bmatrix} 1000000 & 0 & 0 \\ 0 & 1000000 & 0 \\ 0 & 0 & 1000000 \end{bmatrix}. \]

Hierarchical parameterisation

\[ y_{ijk} \sim N(\mu_{ij}, \sigma^2), \;\;\; \mu_{ij}= \beta_0 + \beta_i + \gamma_{j(i)}, \]

where \(\gamma_{ij} \sim N(0, \sigma^2_B)\), \(\beta_0, \beta_i \sim N(0, 1000000)\), and \(\sigma^2, \sigma^2_B \sim \text{Cauchy(0, 25)}\).

Rather than assume a specific variance-covariance structure, just like lme we can incorporate an appropriate structure to account for different dependency/correlation structures in our data. In RCB designs, it is prudent to capture the residuals to allow checks that there are no outstanding dependency issues following model fitting.

Full effect parameterisation

> modelString="
+ model {
+    #Likelihood
+    for (i in 1:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- beta0 + beta[A[i]] + gamma[Block[i]]
+       res[i] <- y[i]-mu[i]
+    }
+    
+    #Priors
+    beta0 ~ dnorm(0, 1.0E-6)
+    beta[1] <- 0
+    for (i in 2:nA) {
+      beta[i] ~ dnorm(0, 1.0E-6) #prior
+    }
+    for (i in 1:nBlock) {
+      gamma[i] ~ dnorm(0, tau.B) #prior
+    }
+    tau <- pow(sigma,-2)
+    sigma <- z/sqrt(chSq) 
+    z ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq ~ dgamma(0.5, 0.5)
+ 
+    tau.B <- pow(sigma.B,-2)
+    sigma.B <- z/sqrt(chSq.B) 
+    z.B ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq.B ~ dgamma(0.5, 0.5)
+  }
+ "
> 
> ## write the model to a text file
> writeLines(modelString, con = "fullModel.txt")

Arrange the data as a list (as required by JAGS). As input, JAGS will need to be supplied with: the response variable, the predictor matrix, the number of predictors, the total number of observed items. This all needs to be contained within a list object. We will create two data lists, one for each of the hypotheses.

> data.rcb.list <- with(data.rcb,
+         list(y=y,
+                  Block=as.numeric(Block),
+          A=as.numeric(A),
+          n=nrow(data.rcb),
+          nBlock=length(levels(Block)),
+                  nA = length(levels(A))
+          )
+ )

Define the nodes (parameters and derivatives) to monitor and the chain parameters.

> params <- c("beta0","beta",'gamma',"sigma","sigma.B","res")
> burnInSteps = 3000
> nChains = 2
> numSavedSteps = 3000
> thinSteps = 1
> nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

Start the JAGS model (check the model, load data into the model, specify the number of chains and compile the model). Load the R2jags package.

> library(R2jags)

Now run the JAGS code via the R2jags interface.

> data.rcb.r2jags.f <- jags(data = data.rcb.list, inits = NULL, parameters.to.save = params,
+     model.file = "fullModel.txt", n.chains = nChains, n.iter = nIter,
+     n.burnin = burnInSteps, n.thin = thinSteps)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 105
   Unobserved stochastic nodes: 42
   Total graph size: 582

Initializing model
> 
> print(data.rcb.r2jags.f)
Inference for Bugs model at "fullModel.txt", fit using jags,
 2 chains, each with 4500 iterations (first 3000 discarded)
 n.sims = 3000 iterations saved
          mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]     0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
beta[2]    27.923   1.236  25.587  27.086  27.918  28.736  30.363 1.001  3000
beta[3]    40.263   1.229  37.821  39.463  40.271  41.070  42.706 1.001  3000
beta0      41.834   2.154  37.519  40.406  41.833  43.338  46.002 1.001  3000
gamma[1]    3.731   3.394  -2.865   1.418   3.748   6.043  10.452 1.001  3000
gamma[2]    4.534   3.439  -2.033   2.182   4.508   6.865  11.317 1.003   690
gamma[3]   -3.951   3.464 -10.690  -6.324  -3.884  -1.600   2.829 1.001  3000
gamma[4]   -4.129   3.454 -10.758  -6.477  -4.200  -1.784   2.630 1.003   780
gamma[5]   -5.314   3.480 -12.143  -7.633  -5.325  -2.947   1.594 1.001  3000
gamma[6]   -6.050   3.377 -12.486  -8.331  -6.071  -3.700   0.524 1.001  3000
gamma[7]   -0.709   3.373  -7.083  -3.017  -0.728   1.585   6.032 1.001  3000
gamma[8]  -15.033   3.446 -21.689 -17.322 -15.065 -12.741  -7.939 1.001  3000
gamma[9]   27.856   3.444  20.996  25.498  27.927  30.226  34.525 1.001  3000
gamma[10]  12.830   3.591   5.798  10.453  12.809  15.249  20.065 1.001  3000
gamma[11] -14.936   3.427 -21.825 -17.228 -14.945 -12.635  -8.165 1.001  3000
gamma[12]  -7.878   3.427 -14.571 -10.161  -7.929  -5.551  -1.275 1.001  3000
gamma[13] -10.865   3.430 -17.569 -13.203 -10.852  -8.555  -4.076 1.001  2100
gamma[14]   9.153   3.466   2.557   6.679   9.161  11.473  15.941 1.002  1100
gamma[15]  -3.897   3.495 -10.811  -6.227  -3.882  -1.565   2.866 1.001  3000
gamma[16]   1.321   3.437  -5.486  -0.927   1.328   3.539   8.066 1.001  3000
gamma[17]  -4.137   3.445 -11.102  -6.451  -4.046  -1.802   2.384 1.001  2300
gamma[18]  -4.257   3.449 -10.970  -6.607  -4.280  -1.978   2.695 1.001  2900
gamma[19]  16.435   3.468   9.749  14.031  16.449  18.830  23.059 1.002  1600
gamma[20]  -5.108   3.439 -11.784  -7.392  -5.100  -2.851   1.708 1.001  3000
gamma[21]  18.935   3.517  12.023  16.632  18.857  21.210  25.944 1.001  3000
gamma[22] -20.654   3.459 -27.290 -23.041 -20.593 -18.341 -14.063 1.001  3000
gamma[23]   7.325   3.570   0.265   4.959   7.346   9.652  14.226 1.001  2500
gamma[24]  -1.300   3.500  -8.248  -3.680  -1.256   1.043   5.590 1.001  3000
gamma[25]  -6.114   3.419 -12.672  -8.442  -6.078  -3.725   0.554 1.001  2500
gamma[26]   1.038   3.451  -5.502  -1.351   1.018   3.415   7.820 1.001  3000
gamma[27]  -4.346   3.400 -10.942  -6.604  -4.351  -1.980   2.254 1.001  3000
gamma[28]  -4.721   3.368 -11.281  -6.939  -4.724  -2.572   2.057 1.001  3000
gamma[29] -12.328   3.513 -19.166 -14.660 -12.295 -10.096  -5.361 1.001  3000
gamma[30] -12.858   3.534 -19.927 -15.216 -12.782 -10.455  -5.999 1.001  3000
gamma[31]   0.272   3.457  -6.677  -2.020   0.279   2.562   7.061 1.002  1700
gamma[32]   8.682   3.389   1.905   6.358   8.769  10.986  15.304 1.002  1100
gamma[33]   0.315   3.433  -6.393  -2.013   0.292   2.624   7.101 1.001  3000
gamma[34]   9.586   3.491   2.775   7.232   9.621  11.954  16.438 1.002  1600
gamma[35]  23.906   3.429  17.230  21.620  23.976  26.208  30.709 1.004   660
res[1]      0.243   2.906  -5.276  -1.686   0.233   2.121   5.940 1.001  3000
res[2]     -6.771   2.869 -12.345  -8.647  -6.778  -4.904  -1.088 1.001  3000
res[3]      7.465   2.940   1.587   5.472   7.482   9.382  13.252 1.001  3000
res[4]     -3.267   2.927  -8.857  -5.249  -3.282  -1.245   2.325 1.003   710
res[5]     -1.085   2.946  -6.792  -3.130  -1.110   0.894   4.685 1.003   580
res[6]      5.144   2.913  -0.568   3.244   5.128   7.112  11.094 1.003   600
res[7]     -0.049   2.992  -5.806  -2.023  -0.081   1.917   5.915 1.002  1400
res[8]     -2.652   3.020  -8.617  -4.629  -2.690  -0.658   3.254 1.001  2100
res[9]      2.018   2.992  -3.965  -0.008   1.989   4.022   7.816 1.001  1900
res[10]    -2.071   2.894  -7.884  -4.016  -2.075  -0.143   3.626 1.003   790
res[11]     0.729   2.960  -5.122  -1.274   0.754   2.666   6.590 1.003   660
res[12]     0.287   2.974  -5.811  -1.700   0.396   2.232   5.977 1.003   700
res[13]    -2.939   2.956  -8.894  -4.938  -2.989  -0.959   2.831 1.001  3000
res[14]     4.213   2.943  -1.578   2.254   4.220   6.162  10.013 1.002  3000
res[15]    -2.451   2.949  -8.246  -4.449  -2.478  -0.503   3.535 1.001  3000
res[16]    -2.462   2.911  -8.169  -4.429  -2.422  -0.560   3.167 1.001  3000
res[17]     3.440   2.912  -2.164   1.559   3.394   5.390   9.051 1.001  3000
res[18]    -2.208   2.908  -7.897  -4.117  -2.183  -0.279   3.505 1.001  3000
res[19]    -5.249   2.902 -10.958  -7.149  -5.263  -3.289   0.456 1.001  2800
res[20]     4.201   2.886  -1.484   2.238   4.217   6.100   9.930 1.001  3000
res[21]     1.085   2.889  -4.579  -0.852   1.104   2.988   6.981 1.001  3000
res[22]     0.756   2.958  -5.385  -1.108   0.712   2.747   6.541 1.001  3000
res[23]     1.284   2.979  -4.882  -0.642   1.295   3.319   6.941 1.001  3000
res[24]    -5.388   2.941 -11.391  -7.320  -5.410  -3.425   0.336 1.001  3000
res[25]     3.141   2.979  -2.786   1.153   3.135   5.090   9.051 1.001  3000
res[26]    -4.587   2.978 -10.372  -6.597  -4.605  -2.589   1.268 1.001  3000
res[27]     7.011   2.963   1.359   4.983   6.994   8.965  12.803 1.001  3000
res[28]     7.495   3.018   1.549   5.443   7.508   9.538  13.331 1.002  1600
res[29]     0.730   3.013  -5.133  -1.311   0.723   2.759   6.603 1.001  2200
res[30]    -5.563   3.054 -11.597  -7.683  -5.577  -3.452   0.255 1.001  2100
res[31]    -3.927   2.962  -9.656  -5.903  -3.900  -1.965   1.925 1.001  3000
res[32]     2.986   2.928  -2.794   1.134   2.985   4.946   8.735 1.001  3000
res[33]    -1.871   2.951  -7.734  -3.824  -1.864   0.076   4.020 1.001  3000
res[34]    -0.528   2.951  -6.322  -2.505  -0.516   1.395   5.205 1.001  3000
res[35]    -1.472   2.949  -7.285  -3.350  -1.439   0.455   4.220 1.001  3000
res[36]     0.722   2.938  -4.922  -1.223   0.716   2.663   6.487 1.001  3000
res[37]    -0.493   2.928  -6.329  -2.503  -0.515   1.546   5.111 1.002   940
res[38]    -2.832   2.938  -8.610  -4.822  -2.827  -0.812   2.778 1.002  1300
res[39]     1.268   2.931  -4.339  -0.752   1.281   3.287   6.896 1.002  1200
res[40]     2.968   2.956  -2.952   0.972   2.987   4.925   8.768 1.003   540
res[41]    -2.427   2.942  -8.293  -4.383  -2.363  -0.534   3.629 1.003   660
res[42]     1.150   2.922  -4.453  -0.813   1.176   3.054   7.197 1.003   620
res[43]    -7.026   2.953 -12.859  -8.930  -7.097  -5.036  -1.170 1.001  3000
res[44]     2.862   2.915  -2.600   0.922   2.711   4.854   8.438 1.001  3000
res[45]     3.397   2.955  -2.273   1.404   3.364   5.361   9.127 1.001  3000
res[46]     1.390   2.972  -4.597  -0.607   1.321   3.366   7.167 1.001  2000
res[47]     2.490   2.991  -3.362   0.502   2.495   4.531   8.375 1.001  3000
res[48]    -3.582   2.963  -9.351  -5.573  -3.608  -1.606   2.286 1.001  2700
res[49]    -2.288   2.916  -7.729  -4.298  -2.366  -0.291   3.588 1.003  1000
res[50]    -1.083   2.906  -6.569  -3.067  -1.105   0.863   4.716 1.002  1300
res[51]     2.286   2.912  -3.216   0.331   2.244   4.222   8.075 1.002  1200
res[52]    -2.829   2.903  -8.540  -4.796  -2.804  -0.891   2.706 1.001  3000
res[53]     1.533   2.928  -4.178  -0.380   1.544   3.452   7.214 1.001  2700
res[54]     0.366   2.907  -5.326  -1.576   0.391   2.274   6.016 1.001  3000
res[55]     7.373   2.957   1.529   5.353   7.414   9.358  13.176 1.002  1300
res[56]    -3.029   2.984  -9.041  -4.985  -2.963  -1.047   2.888 1.002  1000
res[57]    -0.932   2.961  -6.889  -2.936  -0.824   0.995   4.741 1.002  1100
res[58]     0.955   2.941  -4.752  -1.089   0.981   2.907   6.650 1.001  3000
res[59]    -2.168   2.938  -8.052  -4.106  -2.164  -0.262   3.668 1.001  3000
res[60]    -0.054   2.957  -5.874  -2.011   0.024   1.879   5.755 1.001  3000
res[61]     4.652   3.013  -1.277   2.681   4.595   6.649  10.447 1.001  3000
res[62]     1.763   3.015  -4.198  -0.283   1.800   3.811   7.690 1.001  3000
res[63]    -2.627   2.975  -8.515  -4.603  -2.640  -0.650   3.124 1.001  3000
res[64]    -1.878   3.019  -7.923  -3.808  -1.879   0.207   3.832 1.001  3000
res[65]    -7.955   2.986 -14.124  -9.906  -7.908  -5.914  -2.185 1.001  3000
res[66]     5.629   3.022  -0.370   3.651   5.702   7.692  11.441 1.001  3000
res[67]    -9.447   2.997 -15.297 -11.427  -9.472  -7.515  -3.360 1.002  1100
res[68]     3.633   3.011  -2.241   1.639   3.558   5.638   9.647 1.002  1400
res[69]     7.139   3.005   1.189   5.081   7.118   9.082  13.136 1.002  1300
res[70]    -6.267   2.959 -11.956  -8.323  -6.249  -4.253  -0.493 1.001  3000
res[71]     6.538   2.978   0.563   4.494   6.525   8.512  12.391 1.001  3000
res[72]    -0.621   2.967  -6.432  -2.609  -0.683   1.364   5.125 1.001  3000
res[73]    -0.989   2.943  -6.875  -2.968  -0.975   0.949   4.951 1.001  3000
res[74]     1.375   2.946  -4.265  -0.614   1.336   3.269   7.187 1.001  2400
res[75]    -1.399   2.934  -7.135  -3.371  -1.401   0.618   4.478 1.001  2700
res[76]    -0.971   2.970  -6.912  -2.914  -1.004   1.027   4.841 1.001  3000
res[77]    -3.549   2.958  -9.315  -5.540  -3.567  -1.572   2.143 1.001  3000
res[78]     4.860   2.940  -0.988   2.887   4.805   6.835  10.606 1.001  3000
res[79]     6.984   2.964   1.461   4.967   6.920   9.004  12.824 1.001  3000
res[80]    -7.875   3.011 -13.687  -9.926  -7.936  -5.815  -1.954 1.001  3000
res[81]     0.160   2.947  -5.484  -1.824   0.127   2.155   6.047 1.001  3000
res[82]     2.734   2.915  -3.183   0.820   2.794   4.729   8.263 1.001  3000
res[83]     2.626   2.932  -3.161   0.748   2.624   4.515   8.266 1.001  3000
res[84]    -6.416   2.954 -12.297  -8.315  -6.382  -4.439  -0.826 1.001  3000
res[85]    -2.326   3.020  -8.242  -4.293  -2.308  -0.307   3.564 1.001  3000
res[86]    -1.054   3.030  -7.047  -3.116  -1.048   0.971   4.956 1.001  3000
res[87]     0.823   3.034  -5.036  -1.236   0.832   2.846   6.820 1.001  3000
res[88]    -3.700   3.021  -9.662  -5.732  -3.744  -1.656   2.141 1.001  3000
res[89]     5.124   3.000  -0.700   3.125   5.166   7.053  10.967 1.001  3000
res[90]    -3.973   3.026  -9.965  -5.967  -3.979  -1.965   2.004 1.001  3000
res[91]     6.800   2.936   1.114   4.888   6.726   8.792  12.683 1.001  2100
res[92]    -1.633   2.932  -7.382  -3.615  -1.645   0.303   4.342 1.002  1500
res[93]    -5.027   2.937 -10.801  -6.975  -5.090  -3.114   0.895 1.002  1600
res[94]    11.068   2.894   5.415   9.144  11.037  12.950  16.844 1.017   920
res[95]    -5.145   2.871 -10.758  -7.025  -5.209  -3.257   0.653 1.002   920
res[96]    -3.909   2.889  -9.670  -5.821  -3.975  -2.083   2.104 1.002  1000
res[97]     1.670   2.927  -4.070  -0.286   1.655   3.647   7.528 1.001  2600
res[98]    -1.855   2.902  -7.533  -3.814  -1.846   0.085   3.932 1.001  3000
res[99]     0.809   2.905  -4.984  -1.113   0.807   2.832   6.509 1.001  3000
res[100]    1.492   2.952  -4.294  -0.488   1.451   3.471   7.192 1.001  2000
res[101]    0.647   2.987  -5.287  -1.341   0.649   2.697   6.465 1.002  1500
res[102]   -0.289   2.974  -6.080  -2.262  -0.296   1.781   5.395 1.002  1600
res[103]   -1.309   2.976  -7.063  -3.273  -1.394   0.728   4.790 1.004   440
res[104]   11.580   2.959   5.725   9.638  11.567  13.515  17.434 1.003   770
res[105]   -5.108   2.939 -10.788  -7.110  -5.108  -3.165   0.707 1.006   490
sigma       5.090   0.453   4.294   4.775   5.059   5.360   6.091 1.002   980
sigma.B    11.494   1.491   8.926  10.487  11.365  12.348  14.912 1.002   920
deviance  637.702  11.556 618.449 629.190 636.668 645.140 663.417 1.001  3000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 66.8 and DIC = 704.5
DIC is an estimate of expected predictive error (lower deviance is better).

Matrix parameterisation

> modelString2="
+ model {
+    #Likelihood
+    for (i in 1:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- inprod(beta[],X[i,]) + gamma[Block[i]]
+     res[i] <- y[i]-mu[i]
+    } 
+    
+    #Priors
+    beta ~ dmnorm(a0,A0)
+    for (i in 1:nBlock) {
+      gamma[i] ~ dnorm(0, tau.B) #prior
+    }
+    tau <- pow(sigma,-2)
+    sigma <- z/sqrt(chSq) 
+    z ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq ~ dgamma(0.5, 0.5)
+ 
+    tau.B <- pow(sigma.B,-2)
+    sigma.B <- z/sqrt(chSq.B) 
+    z.B ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq.B ~ dgamma(0.5, 0.5)
+  }
+ "
> 
> ## write the model to a text file
> writeLines(modelString2, con = "matrixModel.txt")

Arrange the data as a list (as required by JAGS). As input, JAGS will need to be supplied with: the response variable, the predictor matrix, the number of predictors, the total number of observed items. This all needs to be contained within a list object. We will create two data lists, one for each of the hypotheses.

> A.Xmat <- model.matrix(~A,data.rcb)
> data.rcb.list <- with(data.rcb,
+         list(y=y,
+                  Block=as.numeric(Block),
+          X=A.Xmat,
+          n=nrow(data.rcb),
+          nBlock=length(levels(Block)),
+          a0=rep(0,3), A0=diag(3)
+          )
+ )

Define the nodes (parameters and derivatives) to monitor and the chain parameters.

> params <- c("beta",'gamma',"sigma","sigma.B","res")
> adaptSteps = 1000
> burnInSteps = 3000
> nChains = 2
> numSavedSteps = 3000
> thinSteps = 1
> nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

Now run the JAGS code via the R2jags interface.

> data.rcb.r2jags.m <- jags(data = data.rcb.list, inits = NULL, parameters.to.save = params,
+     model.file = "matrixModel.txt", n.chains = nChains, n.iter = nIter,
+     n.burnin = burnInSteps, n.thin = thinSteps)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 105
   Unobserved stochastic nodes: 40
   Total graph size: 910

Initializing model
> 
> print(data.rcb.r2jags.m)
Inference for Bugs model at "matrixModel.txt", fit using jags,
 2 chains, each with 4500 iterations (first 3000 discarded)
 n.sims = 3000 iterations saved
          mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]     0.550   1.010  -1.396  -0.163   0.543   1.212   2.618 1.001  3000
beta[2]     0.624   0.969  -1.267  -0.021   0.609   1.267   2.494 1.002   830
beta[3]     1.556   1.005  -0.458   0.886   1.564   2.248   3.495 1.001  3000
gamma[1]   64.957  12.115  41.593  56.722  64.856  72.937  89.019 1.001  3000
gamma[2]   65.491  12.061  40.343  57.573  65.450  73.443  89.230 1.001  3000
gamma[3]   57.107  11.937  34.004  49.204  56.969  65.191  79.840 1.001  3000
gamma[4]   56.660  11.705  32.676  48.922  56.796  64.418  79.244 1.002   980
gamma[5]   55.496  12.269  30.868  47.463  55.494  63.576  79.632 1.001  3000
gamma[6]   54.801  11.877  31.954  46.951  54.385  62.674  78.056 1.002  3000
gamma[7]   60.234  11.740  36.788  52.823  60.472  67.916  83.053 1.001  3000
gamma[8]   45.171  11.789  21.628  37.438  45.279  53.104  69.031 1.002  3000
gamma[9]   89.740  11.870  66.877  81.784  89.985  97.893 112.876 1.001  3000
gamma[10]  74.392  11.959  51.486  65.937  74.289  82.703  98.386 1.002  1500
gamma[11]  45.824  11.994  22.631  37.864  45.861  53.737  69.411 1.002  1100
gamma[12]  52.874  11.847  29.777  44.934  53.044  61.009  75.635 1.002  1800
gamma[13]  49.828  12.010  26.670  41.934  49.809  57.903  73.611 1.001  3000
gamma[14]  70.252  11.879  46.842  62.471  70.259  78.152  93.618 1.001  3000
gamma[15]  56.956  11.790  34.318  48.968  56.782  65.102  79.771 1.001  3000
gamma[16]  62.229  12.088  39.561  53.779  62.120  70.801  85.754 1.002  1500
gamma[17]  56.154  11.923  32.379  48.260  56.301  64.248  79.156 1.004   840
gamma[18]  56.302  11.789  32.839  48.585  56.477  64.399  79.221 1.001  3000
gamma[19]  78.011  12.127  53.246  69.932  78.041  86.393 100.917 1.001  3000
gamma[20]  55.445  11.822  32.414  47.480  55.477  63.566  78.456 1.002  3000
gamma[21]  79.975  11.935  56.573  72.184  80.069  87.728 103.091 1.001  3000
gamma[22]  39.667  11.800  16.288  31.836  39.723  47.341  62.882 1.001  3000
gamma[23]  68.605  11.860  45.831  60.540  68.721  76.702  91.376 1.001  3000
gamma[24]  59.858  12.057  36.038  51.799  59.900  67.940  83.312 1.001  2200
gamma[25]  54.974  11.970  31.161  47.135  55.153  62.882  78.327 1.002  1400
gamma[26]  62.397  11.915  38.745  54.113  62.260  70.437  86.239 1.002  3000
gamma[27]  56.526  11.968  32.943  48.610  56.627  64.491  80.067 1.003  1800
gamma[28]  56.062  12.002  33.315  48.254  56.158  64.016  79.264 1.001  3000
gamma[29]  47.976  11.787  24.940  39.984  47.840  55.752  71.140 1.001  3000
gamma[30]  47.866  11.894  24.408  40.161  47.877  55.955  70.787 1.001  3000
gamma[31]  61.528  12.021  37.815  53.617  61.620  69.800  84.596 1.001  3000
gamma[32]  70.047  11.805  46.872  62.230  70.224  78.030  93.180 1.001  3000
gamma[33]  61.830  11.934  38.654  53.718  61.828  69.563  85.450 1.002  1800
gamma[34]  70.909  12.075  47.408  62.788  70.805  78.794  95.338 1.002  3000
gamma[35]  85.532  12.138  61.716  77.422  85.479  93.434 109.336 1.001  3000
res[1]    -19.698  12.037 -43.684 -27.635 -19.684 -11.694   3.838 1.001  3000
res[2]      0.587  12.062 -23.374  -7.403   0.579   8.638  24.005 1.001  2500
res[3]     26.230  12.050   2.155  18.315  26.057  34.356  50.024 1.001  3000
res[4]    -22.940  11.997 -46.569 -30.787 -23.060 -15.111   2.226 1.001  3000
res[5]      6.542  11.989 -16.855  -1.311   6.542  14.381  31.699 1.001  3000
res[6]     24.179  11.977   0.862  16.176  24.040  31.902  49.149 1.001  3000
res[7]    -19.824  11.904 -42.484 -27.853 -19.648 -11.889   3.474 1.001  3000
res[8]      4.872  11.923 -17.801  -3.160   4.933  12.798  28.368 1.001  3000
res[9]     20.951  11.893  -1.904  12.827  21.045  28.785  44.008 1.001  3000
res[10]   -21.576  11.660 -44.079 -29.235 -21.708 -13.770   2.548 1.002   920
res[11]     8.523  11.658 -14.011   0.830   8.435  16.332  32.472 1.002  1100
res[12]    19.489  11.657  -2.829  11.774  19.266  27.343  43.385 1.002   910
res[13]   -22.465  12.167 -46.486 -30.516 -22.434 -14.565   2.114 1.001  3000
res[14]    11.986  12.183 -11.847   3.842  12.016  20.032  36.714 1.001  3000
res[15]    16.730  12.218  -7.281   8.522  16.737  24.706  41.406 1.001  3000
res[16]   -22.029  11.845 -45.221 -29.925 -21.621 -14.244   0.444 1.001  3000
res[17]    11.172  11.874 -12.338   3.326  11.433  18.964  33.977 1.001  3000
res[18]    16.933  11.886  -6.385   8.941  17.228  24.796  39.651 1.001  3000
res[19]   -24.908  11.723 -47.784 -32.552 -25.170 -17.491  -1.875 1.001  3000
res[20]    11.841  11.745 -11.265   4.268  11.662  19.369  35.065 1.001  3000
res[21]    20.133  11.744  -2.743  12.424  19.956  27.740  43.426 1.001  3000
res[22]   -18.164  11.729 -41.484 -26.023 -18.365 -10.358   5.425 1.001  3000
res[23]     9.664  11.740 -13.656   1.865   9.491  17.451  33.025 1.001  3000
res[24]    14.399  11.768  -8.920   6.550  14.224  22.160  38.381 1.001  3000
res[25]   -17.459  11.814 -40.240 -25.508 -17.731  -9.597   5.503 1.001  3000
res[26]     2.112  11.834 -20.695  -5.859   1.791   9.903  25.291 1.001  3000
res[27]    25.119  11.869   2.237  17.073  24.856  33.014  48.190 1.001  3000
res[28]   -12.783  11.927 -36.486 -21.081 -12.685  -4.484  10.101 1.002  1600
res[29]     7.751  11.951 -16.443  -0.606   7.820  16.063  30.397 1.001  2100
res[30]    12.866  11.971 -10.922   4.637  12.799  21.163  36.114 1.002  1600
res[31]   -23.404  11.959 -47.054 -31.381 -23.430 -15.471  -0.221 1.002  1200
res[32]    10.809  11.960 -12.765   2.757  10.883  18.657  33.990 1.002  1400
res[33]    17.359  11.972  -6.134   9.343  17.364  25.293  40.617 1.002  1100
res[34]   -19.997  11.800 -42.736 -28.106 -20.202 -12.124   3.343 1.002  2500
res[35]     6.359  11.812 -16.232  -1.893   6.081  14.115  29.804 1.002  1900
res[36]    19.960  11.807  -2.519  11.753  19.730  27.900  43.458 1.002  2600
res[37]   -19.902  11.980 -43.600 -27.871 -19.843 -12.058   3.587 1.001  3000
res[38]     5.059  12.005 -18.429  -3.025   5.164  12.756  28.837 1.001  3000
res[39]    20.566  11.996  -3.070  12.608  20.634  28.366  44.009 1.001  3000
res[40]   -16.847  11.831 -39.884 -24.766 -17.012  -9.239   6.596 1.001  3000
res[41]     5.057  11.855 -18.074  -2.883   4.927  12.687  28.394 1.001  3000
res[42]    20.042  11.830  -3.241  12.206  20.003  27.718  43.766 1.001  3000
res[43]   -26.596  11.746 -49.680 -34.580 -26.349 -18.627  -3.862 1.001  3000
res[44]    10.592  11.735 -12.047   2.531  10.659  18.547  33.288 1.001  3000
res[45]    22.535  11.759  -0.427  14.411  22.577  30.481  45.274 1.001  3000
res[46]   -18.234  12.031 -41.669 -26.696 -18.076  -9.978   4.443 1.002  1200
res[47]    10.165  12.061 -13.282   1.662  10.365  18.474  33.054 1.002  1500
res[48]    15.501  12.060  -7.938   7.051  15.690  23.578  38.325 1.002  1200
res[49]   -21.296  11.891 -44.247 -29.499 -21.509 -13.441   2.004 1.002  1000
res[50]     7.208  11.902 -15.618  -0.939   6.981  15.152  30.441 1.002   860
res[51]    21.986  11.895  -0.592  13.867  21.799  29.853  45.640 1.002  1000
res[52]   -22.104  11.729 -44.941 -29.987 -22.312 -14.241   1.133 1.001  3000
res[53]     9.556  11.751 -13.024   1.751   9.257  17.457  33.011 1.001  3000
res[54]    19.797  11.734  -2.966  11.838  19.629  27.542  43.295 1.001  3000
res[55]   -12.918  12.097 -35.949 -21.363 -12.982  -4.866  11.578 1.001  3000
res[56]     3.979  12.106 -19.166  -4.374   3.889  12.063  28.610 1.001  3000
res[57]    17.484  12.122  -5.688   8.981  17.300  25.557  42.114 1.001  3000
res[58]   -18.315  11.797 -41.135 -26.620 -18.209 -10.520   4.471 1.001  3000
res[59]     5.862  11.815 -16.900  -2.341   6.146  13.803  29.023 1.001  3000
res[60]    19.383  11.783  -3.464  11.144  19.447  27.072  42.559 1.001  3000
res[61]   -15.105  11.900 -38.125 -22.888 -15.284  -7.276   8.233 1.001  3000
res[62]     9.306  11.923 -13.585   1.507   9.128  17.140  32.775 1.001  3000
res[63]    16.323  11.956  -6.322   8.357  16.116  24.316  39.779 1.001  3000
res[64]   -20.915  11.793 -44.085 -28.598 -21.048 -13.147   2.444 1.001  3000
res[65]     0.307  11.820 -22.750  -7.598   0.334   7.975  23.864 1.001  3000
res[66]    25.298  11.813   2.156  17.551  25.247  33.151  48.893 1.001  3000
res[67]   -29.443  11.830 -52.213 -37.430 -29.692 -21.434  -6.575 1.001  3000
res[68]    10.936  11.848 -11.706   2.810  10.811  18.885  33.881 1.001  3000
res[69]    25.850  11.873   2.835  17.722  25.645  33.821  48.845 1.001  3000
res[70]   -26.142  11.991 -49.697 -34.072 -26.243 -18.049  -2.711 1.001  2300
res[71]    13.963  11.984  -9.101   6.160  13.990  22.029  37.049 1.001  3000
res[72]    18.211  12.015  -4.983  10.314  18.104  26.258  41.685 1.001  2300
res[73]   -20.794  11.932 -43.963 -28.725 -20.938 -12.958   3.106 1.002  1600
res[74]     8.870  11.937 -14.459   0.895   8.787  16.905  32.822 1.001  2000
res[75]    17.504  11.929  -5.274   9.402  17.317  25.457  41.765 1.002  1600
res[76]   -21.046  11.851 -44.910 -28.922 -21.131 -12.843   2.003 1.002  3000
res[77]     3.675  11.857 -19.985  -4.059   3.613  11.706  27.196 1.003  3000
res[78]    23.492  11.881  -0.028  15.552  23.399  31.576  46.776 1.003  3000
res[79]   -12.603  11.937 -35.961 -20.361 -12.719  -4.579  10.948 1.002  2900
res[80]    -0.163  11.955 -23.303  -8.120  -0.343   7.789  23.405 1.002  2200
res[81]    19.279  11.963  -4.001  11.319  19.160  27.293  42.755 1.002  3000
res[82]   -16.766  11.961 -39.955 -24.721 -16.887  -9.030   6.016 1.001  3000
res[83]    10.426  11.958 -13.000   2.504  10.510  18.221  32.859 1.001  3000
res[84]    12.792  11.939 -10.192   4.873  12.720  20.537  35.661 1.001  3000
res[85]   -21.347  11.725 -44.125 -28.923 -21.164 -13.609   1.848 1.001  3000
res[86]     7.224  11.690 -15.758  -0.362   7.382  15.079  30.544 1.001  3000
res[87]    20.510  11.739  -2.412  12.708  20.647  28.327  43.283 1.001  3000
res[88]   -23.140  11.838 -46.005 -31.274 -23.099 -15.502  -0.066 1.001  3000
res[89]    12.983  11.858  -9.972   4.729  13.009  20.680  36.860 1.001  3000
res[90]    15.293  11.892  -7.122   6.928  15.405  23.032  38.696 1.001  3000
res[91]   -13.173  11.935 -36.126 -21.376 -13.267  -5.416  10.293 1.001  3000
res[92]     5.694  11.935 -17.082  -2.485   5.640  13.450  29.247 1.001  3000
res[93]    13.708  11.936  -9.020   5.488  13.624  21.589  37.402 1.001  3000
res[94]    -9.013  11.766 -31.685 -16.863  -9.233  -1.091  14.418 1.001  3000
res[95]     2.073  11.780 -20.569  -5.756   1.871   9.830  25.639 1.001  3000
res[96]    14.717  11.782  -8.181   6.737  14.422  22.602  37.767 1.001  3000
res[97]   -18.561  11.906 -41.914 -26.360 -18.779 -10.492   4.073 1.002  1300
res[98]     5.213  11.908 -18.272  -2.559   5.068  13.190  28.414 1.002  1100
res[99]    19.285  11.939  -4.469  11.483  19.203  27.469  42.441 1.002  1300
res[100]  -18.547  12.018 -42.274 -26.484 -18.430 -10.631   5.083 1.001  3000
res[101]    7.907  12.020 -16.041  -0.084   8.072  15.733  31.471 1.001  3000
res[102]   18.379  12.034  -5.679  10.333  18.514  26.429  41.930 1.001  3000
res[103]  -21.652  12.104 -45.310 -29.603 -21.605 -13.692   1.709 1.001  2500
res[104]   18.537  12.095  -5.105  10.797  18.645  26.404  42.117 1.001  1900
res[105]   13.256  12.137 -10.949   5.489  13.253  21.341  36.897 1.001  2600
sigma      20.838   1.809  17.597  19.556  20.736  21.936  24.751 1.002  1000
sigma.B    63.500   7.812  50.201  58.005  62.978  68.138  80.806 1.001  3000
deviance  934.350  11.459 914.767 926.442 933.465 941.520 959.367 1.004   460

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 65.5 and DIC = 999.9
DIC is an estimate of expected predictive error (lower deviance is better).

For a simple model with only two hierarchical levels, the model is the same as above. If you want to include finite-population standard deviations in the model you can use the following code.

> modelString3="
+ model {
+    #Likelihood (esimating site means (gamma.site)
+    for (i in 1:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- gamma[Block[i]] + inprod(beta[], X[i,]) 
+       y.err[i]<- mu[i]-y[i]
+    }
+    for (i in 1:nBlock) {
+       gamma[i] ~ dnorm(0, tau.block)
+    }
+    #Priors
+    for (i in 1:nX) {
+      beta[i] ~ dnorm(0, 1.0E-6) #prior
+    }
+    sigma ~ dunif(0, 100)
+    tau <- 1 / (sigma * sigma)
+    sigma.block ~ dunif(0, 100)
+    tau.block <- 1 / (sigma.block * sigma.block)
+ 
+    sd.y <- sd(y.err)
+    sd.block <- sd(gamma)
+  }
+ "
> 
> ## write the model to a text file
> writeLines(modelString3, con = "SDModel.txt")
> 
> #data list
> A.Xmat <- model.matrix(~A,ddply(data.rcb,~Block,catcolwise(unique)))
> data.rcb.list <- with(data.rcb,
+         list(y=y,
+                  Block=Block,
+          X= A.Xmat,
+          n=nrow(data.rcb),
+          nBlock=length(levels(Block)),
+                  nX = ncol(A.Xmat)
+          )
+ )
> 
> #parameters and chain details
> params <- c("beta","sigma","sd.y",'sd.block','sigma.block')
> burnInSteps = 3000
> nChains = 2
> numSavedSteps = 3000
> thinSteps = 1
> nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
> 
> data.rcb.r2jagsSD <- jags(data = data.rcb.list, inits = NULL, parameters.to.save = params,
+     model.file = "SDModel.txt", n.chains = nChains, n.iter = nIter,
+     n.burnin = burnInSteps, n.thin = thinSteps)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 105
   Unobserved stochastic nodes: 40
   Total graph size: 899

Initializing model
> 
> print(data.rcb.r2jagsSD)
Inference for Bugs model at "SDModel.txt", fit using jags,
 2 chains, each with 4500 iterations (first 3000 discarded)
 n.sims = 3000 iterations saved
            mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]      41.715   2.196  37.449  40.231  41.710  43.183  45.995 1.001  3000
beta[2]      27.928   1.209  25.537  27.146  27.918  28.713  30.317 1.002   980
beta[3]      40.272   1.210  37.832  39.461  40.267  41.096  42.658 1.001  2300
sd.block     11.358   0.519  10.353  11.029  11.345  11.706  12.370 1.001  3000
sd.y          5.014   0.260   4.592   4.827   4.986   5.172   5.609 1.002  1300
sigma         5.074   0.443   4.322   4.752   5.045   5.350   6.036 1.001  3000
sigma.block  11.692   1.546   9.114  10.589  11.586  12.612  15.118 1.002  3000
deviance    637.262  10.949 618.392 629.413 636.321 644.252 660.930 1.001  2200

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 59.9 and DIC = 697.2
DIC is an estimate of expected predictive error (lower deviance is better).

Calculate \(R^2\) from the posterior of the model.

> data.rcb.mcmc.listSD <- as.mcmc(data.rcb.r2jagsSD)
> 
> Xmat <- model.matrix(~A, data.rcb)
> coefs <- data.rcb.r2jagsSD$BUGSoutput$sims.list[['beta']]
> fitted <- coefs %*% t(Xmat)
> X.var <- aaply(fitted,1,function(x){var(x)})
> Z.var <- data.rcb.r2jagsSD$BUGSoutput$sims.list[['sd.block']]^2
> R.var <- data.rcb.r2jagsSD$BUGSoutput$sims.list[['sd.y']]^2
> R2.marginal <- (X.var)/(X.var+Z.var+R.var)
> R2.marginal <- data.frame(Mean=mean(R2.marginal), Median=median(R2.marginal), HPDinterval(as.mcmc(R2.marginal)))
> R2.conditional <- (X.var+Z.var)/(X.var+Z.var+R.var)
> R2.conditional <- data.frame(Mean=mean(R2.conditional),
+    Median=median(R2.conditional), HPDinterval(as.mcmc(R2.conditional)))
> R2.block <- (Z.var)/(X.var+Z.var+R.var)
> R2.block <- data.frame(Mean=mean(R2.block), Median=median(R2.block), HPDinterval(as.mcmc(R2.block)))
> R2.res<-(R.var)/(X.var+Z.var+R.var)
> R2.res <- data.frame(Mean=mean(R2.res), Median=median(R2.res), HPDinterval(as.mcmc(R2.res)))
> 
> rbind(R2.block=R2.block, R2.marginal=R2.marginal, R2.res=R2.res, R2.conditional=R2.conditional)
                    Mean     Median      lower      upper
R2.block       0.2927774 0.29248056 0.24902731 0.33605200
R2.marginal    0.6500204 0.65101312 0.60509352 0.68965593
R2.res         0.0572022 0.05628758 0.04596228 0.07055798
R2.conditional 0.9427978 0.94371242 0.92944202 0.95403772

Planned comparisonsand pairwise tests

Since there are no restrictions on the type and number of comparisons derived from the posteriors, Bayesian analyses provide a natural framework for exploring additional contrasts and comparisons. For example, to compare all possible levels:

> coefs <- data.rcb.r2jags.m$BUGSoutput$sims.list[[c('beta')]]
> head(coefs)
           [,1]        [,2]       [,3]
[1,] -1.0697767 -0.46647636  0.4808020
[2,]  0.6186153  1.46210386  2.3592529
[3,] -1.5100302  0.09180824  1.1835869
[4,] -0.3127107  0.66392714 -0.5681012
[5,]  1.5552936  1.06785499  2.6443403
[6,]  0.7282182  0.59829747  2.8548669
> 
> newdata <- data.frame(A=levels(data.rcb$A))
> # A Tukeys contrast matrix
> library(multcomp)
> tuk.mat <- contrMat(n=table(newdata$A), type="Tukey")
> Xmat <- model.matrix(~A, data=newdata)
> pairwise.mat <- tuk.mat %*% Xmat
> pairwise.mat
      (Intercept) A2 A3
2 - 1           0  1  0
3 - 1           0  0  1
3 - 2           0 -1  1
> 
> comps <- coefs %*% t(pairwise.mat)
> 
> MCMCsum <- function(x) {
+    data.frame(Median=median(x, na.rm=TRUE), t(quantile(x,na.rm=TRUE)),
+               HPDinterval(as.mcmc(x)),HPDinterval(as.mcmc(x),p=0.5))
+ }
> 
> (comps <-plyr:::adply(comps,2,MCMCsum))
     X1    Median       X0.         X25.      X50.     X75.    X100.      lower
1 2 - 1 0.6093838 -2.556240 -0.020575421 0.6093838 1.267051 4.786166 -1.2766747
2 3 - 1 1.5638199 -1.833977  0.886430287 1.5638199 2.248195 4.835948 -0.4024791
3 3 - 2 0.9310770 -4.672228 -0.003765539 0.9310770 1.864871 5.592247 -1.5970184
     upper     lower.1  upper.1
1 2.479999  0.03512204 1.316762
2 3.539364  0.92897200 2.273364
3 3.728297 -0.03345687 1.823124

RCB (repeated measures) - continuous within

Data generation

Imagine now that we has designed an experiment to investigate the effects of a continuous predictor (\(x\), for example time) on a response (\(y\)). Again, the system that we intend to sample is spatially heterogeneous and thus will add a great deal of noise to the data that will make it difficult to detect a signal (impact of treatment). Thus in an attempt to constrain this variability, we again decide to apply a design (RCB) in which each of the levels of \(X\) (such as time) treatments within each of \(35\) blocks dispersed randomly throughout the landscape. As this section is mainly about the generation of artificial data (and not specifically about what to do with the data), understanding the actual details are optional and can be safely skipped.

> set.seed(123)
> slope <- 30
> intercept <- 200
> nBlock <- 35
> nTime <- 10
> sigma <- 50
> sigma.block <- 30
> n <- nBlock*nTime
> Block <- gl(nBlock, k=1)
> Time <- 1:10
> rho <- 0.8
> dt <- expand.grid(Time=Time,Block=Block)
> Xmat <- model.matrix(~-1+Block + Time, data=dt)
> block.effects <- rnorm(n = nBlock, mean = intercept, sd = sigma.block)
> #A.effects <- c(30,40)
> all.effects <- c(block.effects,slope)
> lin.pred <- Xmat %*% all.effects
> 
> # OR
> Xmat <- cbind(model.matrix(~-1+Block,data=dt),model.matrix(~Time,data=dt))
> ## Sum to zero block effects
> ##block.effects <- rnorm(n = nBlock, mean = 0, sd = sigma.block)
> ###A.effects <- c(40,70,80)
> ##all.effects <- c(block.effects,intercept,slope)
> ##lin.pred <- Xmat %*% all.effects
> 
> ## the quadrat observations (within sites) are drawn from
> ## normal distributions with means according to the site means
> ## and standard deviations of 5
> eps <- NULL
> eps[1] <- 0
> for (j in 2:n) {
+   eps[j] <- rho*eps[j-1] #residuals
+ }
> y <- rnorm(n,lin.pred,sigma)+eps
> 
> #OR
> eps <- NULL
> # first value cant be autocorrelated
> eps[1] <- rnorm(1,0,sigma)
> for (j in 2:n) {
+   eps[j] <- rho*eps[j-1] + rnorm(1, mean = 0, sd = sigma)  #residuals
+ }
> y <- lin.pred + eps
> data.rm <- data.frame(y=y, dt)
> head(data.rm)  #print out the first six rows of the data set
         y Time Block
1 282.1142    1     1
2 321.1404    2     1
3 278.7700    3     1
4 285.8709    4     1
5 336.6390    5     1
6 333.5961    6     1
> 
> ggplot(data.rm, aes(y=y, x=Time)) + geom_smooth(method='lm') + geom_point() + facet_wrap(~Block)

Exploratory data analysis

Normality and Homogeneity of variance

> boxplot(y~Time, data.rm)

> 
> ggplot(data.rm, aes(y=y, x=factor(Time))) + geom_boxplot()

Conclusions:

  • there is no evidence that the response variable is consistently non-normal across all populations - each boxplot is approximately symmetrical.

  • there is no evidence that variance (as estimated by the height of the boxplots) differs between the five populations. More importantly, there is no evidence of a relationship between mean and variance - the height of boxplots does not increase with increasing position along the \(y\)-axis. Hence it there is no evidence of non-homogeneity

Obvious violations could be addressed either by:

  • transform the scale of the response variables (to address normality, etc). Note transformations should be applied to the entire response variable (not just those populations that are skewed).

Block by within-Block interaction

> with(data.rm, interaction.plot(Time,Block,y))

> 
> ggplot(data.rm, aes(y=y, x=Time, color=Block, group=Block)) + geom_line() +
+   guides(color=guide_legend(ncol=3))

> 
> residualPlots(lm(y~Block+Time, data.rm))

           Test stat Pr(>|Test stat|)
Block                                
Time         -0.7274           0.4675
Tukey test   -0.9809           0.3267
> 
> # the Tukey's non-additivity test by itself can be obtained via an internal function
> # within the car package
> car:::tukeyNonaddTest(lm(y~Block+Time, data.rm))
      Test     Pvalue 
-0.9808606  0.3266615 
> 
> # alternatively, there is also a Tukey's non-additivity test within the
> # asbio package
> with(data.rm,tukey.add.test(y,Time,Block))

Tukey's one df test for additivity 
F = 0.3997341   Denom df = 305    p-value = 0.5277003

Conclusions:

  • there is no visual or inferential evidence of any major interactions between Block and the within-Block effect (Time). Any trends appear to be reasonably consistent between Blocks.

Sphericity

Since the levels of Time cannot be randomly assigned, it is likely that sphericity is not met. We can explore whether there is an auto-correlation patterns in the residuals. Note, as there was only ten time periods, it does not make logical sense to explore lags above \(10\).

> library(nlme)
> data.rm.lme <- lme(y~Time, random=~1|Block, data=data.rm)
> acf(resid(data.rm.lme), lag=10)

Conclusions:

The autocorrelation factor (ACF) at a range of lags up to \(10\), indicate that there is a cyclical pattern of residual auto-correlation. We really should explore incorporating some form of correlation structure into our model.

Model fitting

Full effect parameterisation

> modelString="
+ model {
+    #Likelihood
+    for (i in 1:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- beta0 + beta*Time[i] + gamma[Block[i]]
+       res[i] <- y[i]-mu[i]
+    }
+    
+    #Priors
+    beta0 ~ dnorm(0, 1.0E-6)
+    beta ~ dnorm(0, 1.0E-6) #prior
+    
+    for (i in 1:nBlock) {
+      gamma[i] ~ dnorm(0, tau.B) #prior
+    }
+    tau <- pow(sigma,-2)
+    sigma <- z/sqrt(chSq) 
+    z ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq ~ dgamma(0.5, 0.5)
+ 
+    tau.B <- pow(sigma.B,-2)
+    sigma.B <- z/sqrt(chSq.B) 
+    z.B ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq.B ~ dgamma(0.5, 0.5)
+  }
+ "
> 
> ## write the model to a text file
> writeLines(modelString, con = "fullModel2.txt")
> 
> data.rm.list <- with(data.rm,
+         list(y=y,
+                  Block=as.numeric(Block),
+          Time=Time,
+          n=nrow(data.rm),
+          nBlock=length(levels(Block))
+              )
+ )
> 
> params <- c("beta0","beta",'gamma',"sigma","sigma.B","res")
> burnInSteps = 3000
> nChains = 2
> numSavedSteps = 3000
> thinSteps = 1
> nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
> 
> data.rm.r2jags.f <- jags(data = data.rm.list, inits = NULL, parameters.to.save = params,
+     model.file = "fullModel2.txt", n.chains = nChains, n.iter = nIter,
+     n.burnin = burnInSteps, n.thin = thinSteps)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 350
   Unobserved stochastic nodes: 41
   Total graph size: 1815

Initializing model
> 
> print(data.rm.r2jags.f)
Inference for Bugs model at "fullModel2.txt", fit using jags,
 2 chains, each with 4500 iterations (first 3000 discarded)
 n.sims = 3000 iterations saved
           mu.vect sd.vect     2.5%      25%      50%      75%    97.5%  Rhat
beta        30.689   1.047   28.609   30.017   30.687   31.401   32.705 1.001
beta0      189.009  12.648  164.589  180.318  189.054  197.523  213.976 1.001
gamma[1]   -35.015  20.219  -74.991  -47.706  -35.021  -21.382    3.826 1.002
gamma[2]   -52.026  20.114  -91.663  -65.012  -52.008  -38.575  -12.685 1.001
gamma[3]    20.417  19.878  -19.313    7.462   20.587   33.579   60.228 1.001
gamma[4]     0.671  20.295  -38.799  -12.839    0.882   14.490   39.807 1.001
gamma[5]    67.812  19.967   29.683   54.189   67.514   81.150  109.392 1.001
gamma[6]    36.338  19.760   -2.575   22.780   36.203   49.381   76.772 1.001
gamma[7]    24.072  20.155  -14.701   10.740   24.009   36.728   63.397 1.001
gamma[8]   -31.199  20.016  -70.564  -44.687  -31.149  -17.691    7.011 1.001
gamma[9]    73.971  20.034   35.053   60.309   73.846   87.726  113.132 1.003
gamma[10]   58.034  19.900   19.397   44.730   58.085   71.380   97.283 1.001
gamma[11]  141.644  20.387  101.956  127.950  141.240  154.897  181.751 1.001
gamma[12]    5.655  20.094  -32.833   -7.787    5.349   19.065   47.017 1.001
gamma[13]  -44.187  20.168  -84.778  -57.576  -44.641  -30.571   -5.599 1.001
gamma[14]  -23.866  19.908  -63.673  -37.311  -23.653  -10.578   14.435 1.001
gamma[15]   30.407  20.239   -8.109   16.928   30.379   43.830   70.587 1.001
gamma[16]  103.433  20.123   64.608   90.052  103.087  116.736  143.495 1.002
gamma[17]   91.556  20.060   53.115   77.561   91.725  104.473  131.814 1.002
gamma[18]  -63.563  20.127 -102.913  -77.195  -63.210  -50.190  -24.916 1.002
gamma[19]   16.404  19.820  -21.892    2.880   16.232   29.420   55.497 1.001
gamma[20]  -26.858  19.837  -66.283  -39.890  -26.676  -13.651   12.760 1.002
gamma[21] -104.771  19.743 -143.174 -117.701 -105.347  -92.448  -64.620 1.001
gamma[22]  -14.307  19.903  -54.617  -27.704  -14.041   -0.901   23.918 1.001
gamma[23]  -81.493  19.863 -121.932  -94.860  -81.367  -67.618  -43.350 1.001
gamma[24]  -86.520  20.067 -125.826 -100.133  -86.297  -73.003  -47.481 1.001
gamma[25]  -47.166  20.417  -86.549  -61.568  -47.155  -33.479   -6.020 1.001
gamma[26]  -92.375  19.497 -130.380 -105.540  -92.310  -79.043  -54.017 1.002
gamma[27]   20.875  20.031  -18.328    6.625   20.873   34.275   60.661 1.002
gamma[28]   74.464  19.909   36.179   61.091   74.369   87.433  114.480 1.001
gamma[29]   -0.792  19.771  -39.202  -14.402   -0.589   12.677   36.999 1.001
gamma[30]  -75.855  20.350 -116.077  -89.286  -76.179  -62.465  -35.911 1.001
gamma[31]  -58.457  20.104  -97.479  -72.394  -58.390  -44.695  -19.492 1.001
gamma[32]  -53.274  20.080  -91.482  -66.984  -53.201  -40.309  -13.057 1.001
gamma[33]   15.405  20.131  -22.770    1.695   15.048   28.831   55.869 1.001
gamma[34]   59.338  20.334   19.750   45.810   59.517   72.790   99.194 1.001
gamma[35]   56.933  20.301   15.197   43.906   57.388   70.588   94.583 1.001
res[1]      97.432  17.921   62.629   85.359   97.337  109.111  134.580 1.002
res[2]     105.769  17.698   71.182   93.805  105.826  117.408  142.144 1.002
res[3]      32.710  17.534   -1.148   20.909   32.522   44.178   68.460 1.002
res[4]       9.121  17.431  -23.865   -2.637    8.871   20.532   44.239 1.002
res[5]      29.201  17.391   -4.021   17.641   28.835   40.540   64.071 1.002
res[6]      -4.531  17.414  -38.094  -16.092   -4.863    6.832   30.205 1.002
res[7]    -107.185  17.500 -140.651 -118.950 -107.516  -95.741  -72.125 1.002
res[8]     -37.350  17.647  -71.015  -49.217  -37.556  -25.708   -1.746 1.002
res[9]     -67.307  17.855 -101.635  -79.246  -67.326  -55.699  -31.035 1.002
res[10]    -78.614  18.121 -113.985  -90.735  -78.654  -66.764  -42.063 1.002
res[11]     49.884  18.193   13.728   37.834   50.102   61.789   84.845 1.001
res[12]     11.595  17.930  -23.183   -0.446   11.692   23.269   46.321 1.001
res[13]     61.820  17.726   27.566   49.860   61.924   73.375   96.299 1.001
res[14]     -3.458  17.582  -37.401  -15.192   -3.465    8.046   30.969 1.001
res[15]    -10.511  17.499  -44.587  -21.767  -10.459    1.089   24.112 1.001
res[16]     -2.243  17.479  -36.572  -13.604   -2.049    9.337   32.255 1.001
res[17]    -50.520  17.522  -85.179  -61.867  -50.373  -38.782  -16.278 1.001
res[18]    -62.585  17.627  -97.953  -73.969  -62.501  -50.564  -28.305 1.001
res[19]    -42.079  17.792  -77.762  -53.482  -42.166  -30.142   -7.621 1.001
res[20]      9.165  18.018  -26.988   -2.455    9.147   21.271   43.778 1.001
res[21]    -77.926  17.501 -113.047  -89.469  -78.006  -66.052  -43.688 1.002
res[22]    -73.189  17.257 -108.057  -84.536  -73.411  -61.442  -39.850 1.002
res[23]    -14.228  17.075  -48.379  -25.553  -14.413   -2.773   18.665 1.002
res[24]    -31.958  16.955  -66.594  -43.263  -32.141  -20.655    0.554 1.002
res[25]     -7.975  16.900  -42.439  -19.028   -8.192    3.408   24.762 1.002
res[26]     24.320  16.909   -9.262   13.308   24.202   35.576   57.313 1.002
res[27]     38.799  16.983    5.216   27.563   38.571   50.247   72.421 1.002
res[28]     69.516  17.120   35.582   58.330   69.291   80.872  103.259 1.002
res[29]     55.153  17.321   20.895   43.887   54.848   66.615   89.335 1.002
res[30]     28.976  17.581   -5.726   17.342   28.565   40.577   63.635 1.002
res[31]   -121.587  17.872 -156.152 -133.567 -121.526 -109.634  -87.382 1.001
res[32]   -100.256  17.618 -133.884 -111.966 -100.063  -88.659  -66.409 1.001
res[33]    -57.168  17.423  -90.542  -68.902  -57.110  -45.627  -22.964 1.001
res[34]    -17.580  17.290  -51.232  -29.462  -17.603   -6.075   16.449 1.001
res[35]    -40.581  17.219  -74.430  -52.498  -40.599  -29.045   -6.738 1.001
res[36]     57.619  17.212   23.794   45.812   57.457   69.114   91.210 1.001
res[37]     61.388  17.269   27.008   49.782   61.277   72.884   95.367 1.004
res[38]     56.259  17.389   21.473   44.753   56.294   67.702   90.416 1.001
res[39]    109.316  17.570   73.992   97.577  109.390  120.643  143.528 1.001
res[40]     52.088  17.811   16.206   39.944   52.135   63.655   86.709 1.001
res[41]    -38.914  17.643  -73.656  -50.658  -38.634  -26.860   -5.046 1.001
res[42]     77.326  17.420   42.616   65.720   77.575   89.215  111.120 1.001
res[43]     50.865  17.258   16.915   39.497   51.222   62.617   84.564 1.001
res[44]    110.679  17.159   76.469   99.239  111.186  122.203  144.204 1.001
res[45]      4.790  17.122  -29.711   -6.616    5.343   16.301   38.330 1.001
res[46]    -17.661  17.150  -52.491  -29.227  -17.113   -5.957   16.017 1.001
res[47]     -7.311  17.242  -42.494  -18.933   -6.757    4.462   26.166 1.001
res[48]     -3.089  17.396  -38.865  -14.616   -2.473    8.846   31.113 1.001
res[49]    -65.133  17.611 -101.697  -76.777  -64.519  -53.174  -30.398 1.001
res[50]    -63.661  17.886 -100.707  -75.693  -63.048  -51.501  -28.264 1.001
res[51]    -31.518  17.144  -65.465  -43.209  -31.355  -19.712    1.672 1.001
res[52]     14.815  16.893  -18.897    3.470   14.776   26.305   47.890 1.001
res[53]     70.347  16.704   36.921   58.939   70.092   81.766  102.877 1.001
res[54]    -50.853  16.579  -83.834  -62.055  -51.073  -39.619  -18.942 1.001
res[55]     25.083  16.520   -8.057   14.027   24.935   36.357   56.754 1.001
res[56]    -38.143  16.527  -71.222  -49.423  -38.357  -26.836   -6.450 1.001
res[57]     -4.070  16.600  -36.763  -15.317   -4.236    7.197   27.724 1.001
res[58]     33.306  16.738    0.513   22.191   33.191   44.522   65.598 1.001
res[59]     20.080  16.940  -12.937    8.828   19.836   31.550   52.764 1.001
res[60]    -12.900  17.204  -46.161  -24.276  -13.088   -1.463   20.472 1.001
res[61]    -17.368  17.885  -51.348  -29.241  -17.425   -5.291   18.072 1.001
res[62]      7.369  17.652  -26.193   -4.253    7.234   19.236   42.448 1.001
res[63]     49.245  17.479   15.639   37.647   49.046   60.880   83.698 1.001
res[64]    -64.174  17.368  -97.742  -75.711  -64.357  -52.534  -30.120 1.001
res[65]   -134.249  17.320 -167.266 -146.107 -134.437 -122.586 -100.604 1.001
res[66]    -37.107  17.334  -70.680  -48.991  -37.270  -25.543   -3.298 1.001
res[67]     21.279  17.412  -12.492    9.532   21.267   33.096   55.105 1.001
res[68]     37.284  17.552    3.245   25.564   37.265   49.184   71.420 1.001
res[69]     63.944  17.753   29.095   52.068   63.864   75.938   98.781 1.001
res[70]     95.234  18.013   60.203   83.123   95.270  107.455  130.632 1.001
res[71]    -48.396  17.620  -83.102  -60.471  -48.230  -36.426  -13.894 1.001
res[72]     16.818  17.403  -17.131    4.932   16.936   28.621   50.566 1.001
res[73]    -10.912  17.247  -44.871  -22.609  -10.795    0.606   22.307 1.001
res[74]      2.547  17.154  -31.372   -8.982    2.664   14.012   35.509 1.001
res[75]    -13.113  17.125  -47.212  -24.603  -12.951   -1.578   20.054 1.001
res[76]     32.577  17.159   -1.807   20.989   32.596   44.219   65.805 1.001
res[77]      7.970  17.257  -27.008   -3.791    7.911   19.682   41.186 1.001
res[78]     31.495  17.418   -3.602   19.914   31.400   43.395   65.506 1.001
res[79]      4.718  17.639  -30.550   -6.860    4.806   16.702   39.109 1.001
res[80]    -51.946  17.919  -87.119  -63.726  -51.922  -39.897  -17.025 1.001
res[81]    -63.210  17.647  -97.791  -75.690  -63.212  -51.190  -29.440 1.002
res[82]    -31.067  17.390  -65.386  -43.236  -31.079  -19.273    2.374 1.002
res[83]     43.678  17.193    9.558   31.772   43.540   55.495   76.981 1.002
res[84]     20.380  17.058  -13.497    8.805   20.360   32.176   53.934 1.002
res[85]     54.597  16.986   21.096   43.276   54.749   66.249   87.971 1.002
res[86]    124.353  16.979   90.938  113.063  124.583  135.902  157.927 1.002
res[87]     67.176  17.037   33.903   55.809   67.439   78.691  101.104 1.002
res[88]    -30.778  17.158  -64.046  -42.370  -30.407  -19.292    2.913 1.002
res[89]    -55.098  17.342  -88.751  -66.851  -54.716  -43.392  -21.208 1.002
res[90]    -73.427  17.586 -107.590  -85.203  -73.116  -61.615  -38.693 1.002
res[91]    -39.879  17.759  -75.389  -51.772  -39.985  -28.068   -4.896 1.001
res[92]     40.803  17.486    5.829   29.249   40.623   52.320   75.462 1.001
res[93]      3.288  17.272  -30.880   -8.346    3.017   14.813   37.091 1.001
res[94]      8.096  17.120  -26.071   -3.487    7.840   19.496   41.421 1.001
res[95]    -18.230  17.031  -51.683  -29.609  -18.372   -6.678   14.819 1.001
res[96]    -27.022  17.006  -60.187  -38.426  -27.134  -15.863    6.296 1.001
res[97]    -19.513  17.045  -52.825  -30.880  -19.632   -8.423   13.668 1.001
res[98]     37.064  17.149    3.831   25.570   37.124   48.265   70.659 1.001
res[99]     21.843  17.315  -12.203   10.180   21.878   33.093   55.675 1.001
res[100]    39.105  17.542    5.025   27.487   38.948   50.788   73.056 1.001
res[101]    29.449  17.991   -5.205   17.155   29.781   41.807   63.442 1.001
res[102]    49.685  17.766   15.218   37.673   50.072   61.814   82.819 1.001
res[103]    -8.723  17.601  -43.029  -20.739   -8.303    3.244   24.590 1.001
res[104]    54.477  17.498   20.833   42.653   54.759   66.407   87.773 1.001
res[105]     4.508  17.456  -29.390   -7.159    4.817   16.421   37.676 1.001
res[106]   -21.847  17.477  -55.720  -33.395  -21.535   -9.860   11.772 1.001
res[107]    32.423  17.561   -2.084   20.791   32.811   44.448   66.578 1.001
res[108]    70.203  17.706   35.619   58.406   70.311   82.156  104.571 1.001
res[109]   -18.915  17.912  -54.616  -30.740  -18.911   -6.868   15.481 1.001
res[110]   -79.501  18.176 -115.834  -91.501  -79.578  -67.080  -44.590 1.001
res[111]   -35.407  17.785  -70.692  -46.976  -35.584  -23.160   -0.930 1.001
res[112]   -16.834  17.533  -51.394  -28.356  -17.011   -4.647   17.012 1.001
res[113]    -2.964  17.342  -37.311  -14.283   -3.253    9.068   30.829 1.001
res[114]    17.958  17.211  -16.244    6.888   17.697   29.740   51.566 1.001
res[115]    43.960  17.144    9.743   32.909   43.779   55.749   77.773 1.002
res[116]     6.922  17.141  -27.302   -4.089    6.746   18.522   41.245 1.002
res[117]   -42.437  17.202  -76.634  -53.389  -42.592  -31.097   -7.471 1.002
res[118]    18.962  17.325  -15.305    7.723   18.882   30.432   53.655 1.003
res[119]    54.157  17.511   19.505   42.669   54.007   65.841   88.924 1.003
res[120]   -30.836  17.757  -66.091  -42.472  -31.031  -18.943    4.597 1.003
res[121]    29.694  17.913   -5.583   17.603   29.789   41.163   65.302 1.001
res[122]    -8.428  17.668  -42.894  -20.446   -8.260    2.790   26.838 1.001
res[123]   -97.805  17.483 -132.040 -109.658  -97.722  -86.688  -63.003 1.001
res[124]   -58.400  17.360  -92.551  -70.268  -58.013  -47.184  -24.118 1.001
res[125]   -38.480  17.298  -72.858  -50.159  -38.072  -27.272   -4.162 1.001
res[126]   -23.590  17.301  -58.378  -35.298  -23.103  -12.266   10.906 1.001
res[127]     3.860  17.366  -31.068   -7.655    4.299   15.220   38.206 1.001
res[128]    58.998  17.494   23.758   47.395   59.472   70.391   93.062 1.001
res[129]    69.127  17.683   33.354   57.321   69.673   80.505  103.147 1.001
res[130]    35.991  17.931   -0.747   24.279   36.423   47.825   70.397 1.001
res[131]   -18.707  17.753  -53.599  -30.549  -18.607   -6.714   16.819 1.003
res[132]   -14.747  17.516  -48.954  -26.537  -14.555   -2.845   20.061 1.003
res[133]    10.374  17.340  -23.749   -1.370   10.597   22.133   44.664 1.003
res[134]   -37.152  17.225  -70.611  -48.839  -37.026  -25.712   -3.136 1.002
res[135]   -32.541  17.174  -65.860  -44.108  -32.536  -21.077    1.265 1.002
res[136]    28.588  17.186   -4.848   17.066   28.625   39.992   62.044 1.002
res[137]    23.576  17.262  -10.078   11.987   23.698   35.092   57.117 1.002
res[138]   -10.078  17.401  -43.670  -21.709   -9.920    1.610   23.933 1.001
res[139]   -16.016  17.601  -49.961  -27.799  -15.770   -4.358   18.774 1.001
res[140]    48.626  17.861   14.303   36.410   48.905   60.610   83.764 1.001
res[141]     4.593  17.780  -29.625   -7.191    4.424   16.505   38.788 1.001
res[142]    57.463  17.544   23.749   45.968   57.368   69.300   91.358 1.001
res[143]    44.743  17.368   11.543   33.225   44.780   56.414   78.495 1.001
res[144]    47.987  17.253   14.659   36.564   48.140   59.625   81.776 1.001
res[145]     2.009  17.202  -31.053   -9.596    2.094   13.679   35.963 1.001
res[146]    23.279  17.214   -9.893   11.605   23.253   34.820   57.175 1.001
res[147]   -15.427  17.290  -49.433  -27.118  -15.449   -4.007   19.070 1.001
res[148]   -92.242  17.428 -126.419 -104.123  -92.223  -80.822  -57.405 1.001
res[149]   -76.403  17.628 -110.319  -88.293  -76.440  -64.957  -41.280 1.001
res[150]    27.022  17.887   -7.788   14.991   26.934   38.901   62.553 1.001
res[151]    56.524  17.435   23.497   44.650   56.376   68.008   91.937 1.002
res[152]    94.888  17.193   62.473   83.265   94.748  106.107  129.551 1.001
res[153]    85.122  17.012   53.069   73.707   84.983   96.084  119.032 1.001
res[154]    28.800  16.893   -3.079   17.754   28.668   39.844   62.380 1.001
res[155]     3.921  16.839  -28.102   -7.248    3.904   15.053   37.301 1.001
res[156]   -19.671  16.850  -51.711  -30.814  -19.649   -8.511   13.270 1.001
res[157]   -48.454  16.926  -81.115  -59.563  -48.405  -37.301  -15.678 1.001
res[158]   -12.976  17.067  -46.088  -24.257  -12.925   -1.706   20.012 1.001
res[159]   -79.807  17.269 -113.216  -91.035  -79.682  -68.569  -46.865 1.001
res[160]   -30.224  17.532  -64.395  -41.789  -29.847  -18.816    2.893 1.001
res[161]   -10.710  17.984  -46.042  -22.979  -10.276    1.572   24.298 1.002
res[162]   -82.452  17.740 -118.119  -94.601  -82.059  -70.317  -47.916 1.002
res[163]   -48.077  17.555  -83.627  -59.974  -47.795  -35.994  -13.951 1.002
res[164]    68.821  17.431   33.752   56.976   68.995   80.691  102.233 1.005
res[165]    12.830  17.369  -22.123    1.019   12.897   24.719   45.628 1.002
res[166]    38.006  17.370    3.280   26.085   37.983   49.937   71.003 1.002
res[167]   -23.347  17.435  -58.240  -35.377  -23.237  -11.414    9.865 1.002
res[168]    22.078  17.561  -13.256    9.887   22.185   34.157   55.218 1.001
res[169]    15.237  17.749  -20.525    3.105   15.325   27.296   48.552 1.001
res[170]    79.730  17.996   43.542   67.295   79.674   91.820  113.635 1.002
res[171]    63.717  17.877   28.767   51.218   63.625   75.992   98.346 1.002
res[172]    50.693  17.611   16.356   38.398   50.624   62.696   84.475 1.001
res[173]    16.355  17.405  -17.386    4.388   16.412   28.087   49.899 1.001
res[174]     5.229  17.259  -28.246   -6.690    5.281   16.848   38.937 1.001
res[175]   -25.424  17.176  -58.537  -37.419  -25.426  -13.955    7.970 1.001
res[176]   -60.299  17.157  -93.240  -72.218  -60.299  -48.754  -26.651 1.001
res[177]   -17.707  17.202  -50.497  -29.605  -17.777   -6.222   15.477 1.001
res[178]   -67.087  17.310  -99.859  -79.145  -67.067  -55.435  -33.398 1.001
res[179]    21.851  17.480  -10.981    9.770   21.871   33.713   56.070 1.001
res[180]   -40.647  17.711  -73.957  -52.771  -40.513  -28.613   -5.841 1.001
res[181]   -19.458  17.542  -54.333  -31.085  -19.422   -7.722   15.640 1.001
res[182]    13.382  17.297  -20.657    2.031   13.446   24.888   47.985 1.001
res[183]    42.203  17.113    8.488   30.978   42.207   53.606   75.709 1.001
res[184]    20.699  16.991  -12.648    9.734   20.572   32.070   53.897 1.001
res[185]    22.419  16.933  -10.796   11.505   22.206   33.885   55.335 1.001
res[186]    67.746  16.941   34.410   57.004   67.445   79.392  100.423 1.001
res[187]   -17.017  17.012  -50.239  -27.920  -17.298   -5.427   15.763 1.001
res[188]   -51.228  17.148  -84.441  -62.328  -51.548  -39.688  -18.069 1.001
res[189]   -23.628  17.345  -57.805  -34.767  -23.760  -11.997   10.684 1.001
res[190]   -39.945  17.603  -74.634  -51.215  -40.154  -28.208   -5.214 1.001
res[191]    52.530  17.512   18.626   41.139   52.206   64.088   87.150 1.006
res[192]    79.593  17.246   46.248   68.448   79.388   91.100  113.049 1.017
res[193]    71.051  17.040   37.868   59.980   70.793   82.571  104.153 1.021
res[194]   -14.917  16.897  -48.034  -25.902  -15.296   -3.462   18.041 1.005
res[195]    -7.135  16.817  -39.943  -18.047   -7.544    4.215   25.897 1.005
res[196]   -18.174  16.803  -51.133  -29.182  -18.623   -6.895   14.827 1.005
res[197]   -16.439  16.854  -48.969  -27.463  -16.762   -5.147   16.557 1.004
res[198]   -69.150  16.970 -102.173  -80.328  -69.361  -57.572  -35.980 1.004
res[199]   -27.445  17.149  -60.854  -38.683  -27.690  -15.800    6.299 1.003
res[200]   -71.100  17.389 -104.684  -82.727  -71.394  -59.189  -37.691 1.003
res[201]     1.426  17.309  -34.407   -9.571    1.332   12.962   35.011 1.001
res[202]    36.131  17.046    0.771   25.266   36.227   47.347   69.464 1.001
res[203]     5.510  16.844  -28.967   -5.374    5.685   16.535   38.753 1.001
res[204]    49.200  16.705   15.453   38.505   49.435   60.187   82.079 1.001
res[205]   -10.960  16.631  -44.546  -21.713  -10.492   -0.147   21.703 1.001
res[206]  -133.889  16.623 -167.539 -144.693 -133.301 -123.036 -101.431 1.001
res[207]   -68.634  16.681 -102.616  -79.683  -67.975  -57.769  -36.182 1.001
res[208]     2.212  16.804  -31.897   -8.838    2.830   13.283   34.752 1.001
res[209]     2.431  16.991  -31.627   -8.767    3.000   13.550   35.269 1.001
res[210]    41.968  17.240    6.811   30.539   42.440   53.467   75.061 1.001
res[211]   -67.622  17.508 -102.075  -79.236  -67.541  -56.136  -32.937 1.001
res[212]   -57.530  17.266  -91.092  -68.962  -57.500  -46.025  -22.873 1.001
res[213]  -140.313  17.084 -173.718 -151.700 -140.155 -128.926 -105.859 1.002
res[214]   -50.542  16.964  -83.840  -61.841  -50.386  -39.139  -16.241 1.002
res[215]    55.074  16.909   22.122   43.726   55.352   66.374   89.195 1.002
res[216]   100.133  16.919   67.118   88.860  100.592  111.422  133.899 1.002
res[217]    80.975  16.993   48.141   69.704   81.338   92.108  115.080 1.002
res[218]    65.212  17.132   32.285   53.430   65.470   76.536  100.152 1.001
res[219]   -21.674  17.332  -55.118  -33.640  -21.456  -10.022   13.856 1.002
res[220]    24.003  17.593   -9.530   11.811   24.101   35.907   59.852 1.002
res[221]    60.185  17.783   25.820   48.457   59.813   72.436   96.171 1.001
res[222]    26.825  17.539   -6.907   15.267   26.510   38.780   62.320 1.001
res[223]   -37.765  17.355  -71.270  -49.250  -38.173  -25.946   -2.142 1.001
res[224]   -33.962  17.232  -67.496  -45.390  -34.499  -22.046    1.130 1.001
res[225]   -58.522  17.173  -91.829  -70.044  -59.083  -46.732  -23.657 1.001
res[226]   -55.706  17.177  -88.390  -67.397  -56.263  -44.042  -21.182 1.001
res[227]   -94.620  17.245 -127.457 -106.411  -95.119  -82.969  -60.749 1.001
res[228]    19.371  17.375  -13.948    7.607   19.044   31.010   53.459 1.001
res[229]    25.246  17.568   -8.529   13.626   24.828   37.075   59.410 1.001
res[230]    84.355  17.820   50.467   72.405   83.982   96.570  118.968 1.001
res[231]   -31.497  18.030  -66.696  -43.447  -31.785  -19.738    3.853 1.001
res[232]   -33.555  17.750  -68.255  -45.064  -33.695  -22.067    1.221 1.001
res[233]   -46.454  17.529  -80.525  -58.085  -46.585  -35.121  -11.905 1.001
res[234]   -84.283  17.368 -118.128  -95.729  -84.398  -73.094  -50.480 1.001
res[235]    23.792  17.270   -9.429   12.004   23.713   34.911   57.625 1.001
res[236]   -37.979  17.234  -71.525  -49.693  -38.039  -26.826   -3.726 1.001
res[237]    -0.851  17.262  -34.217  -12.666   -0.823   10.179   33.380 1.001
res[238]    55.116  17.354   21.132   43.276   55.278   66.313   89.932 1.001
res[239]    66.340  17.507   31.871   54.397   66.435   77.816  102.154 1.001
res[240]    22.509  17.721  -12.331   10.444   22.467   34.188   58.690 1.001
res[241]    48.779  17.985   13.144   36.603   49.101   61.358   82.744 1.001
res[242]    54.607  17.728   19.723   42.634   54.853   67.014   88.097 1.001
res[243]    -2.573  17.531  -36.920  -14.310   -2.250    9.522   31.153 1.001
res[244]   -64.682  17.394  -98.704  -76.398  -64.353  -52.693  -31.136 1.001
res[245]    59.232  17.320   25.143   47.518   59.597   70.983   92.873 1.001
res[246]    19.963  17.309  -14.316    8.145   20.310   31.886   53.371 1.001
res[247]   -70.443  17.361 -105.373  -82.320  -70.049  -58.397  -36.717 1.001
res[248]   -23.463  17.476  -58.679  -35.360  -23.248  -11.587   10.598 1.001
res[249]     2.831  17.652  -32.482   -8.951    3.087   14.786   37.287 1.001
res[250]   -59.475  17.888  -95.014  -71.176  -59.315  -47.543  -24.802 1.001
res[251]  -118.664  17.326 -151.645 -130.208 -118.894 -107.309  -83.446 1.002
res[252]   -91.020  17.066 -123.512 -102.427  -91.199  -79.836  -57.008 1.001
res[253]    -6.258  16.868  -38.555  -17.830   -6.464    4.709   27.751 1.001
res[254]    36.251  16.734    4.248   24.639   35.981   47.226   70.217 1.001
res[255]    13.667  16.664  -18.023    2.118   13.400   24.523   47.428 1.001
res[256]   -21.601  16.659  -53.042  -33.062  -21.878  -10.646   12.267 1.001
res[257]     5.310  16.721  -26.340   -6.021    5.033   16.411   39.015 1.001
res[258]    21.015  16.847  -10.852    9.662   20.790   32.193   55.145 1.001
res[259]    57.059  17.037   24.518   45.611   56.831   68.241   92.442 1.002
res[260]    34.481  17.288    1.254   23.094   34.411   45.908   70.086 1.001
res[261]    50.413  17.715   15.780   38.279   50.907   62.488   84.294 1.002
res[262]     1.013  17.462  -32.559  -11.003    1.463   13.068   33.840 1.002
res[263]   -13.632  17.268  -46.927  -25.466  -13.251   -1.826   18.917 1.002
res[264]    28.083  17.137   -5.178   16.694   28.340   39.623   60.566 1.003
res[265]    73.774  17.069   40.522   62.477   73.840   85.174  106.109 1.003
res[266]   -36.234  17.065  -69.304  -47.505  -35.963  -24.747   -3.793 1.003
res[267]   -22.093  17.125  -55.763  -33.330  -21.828  -10.616   10.070 1.003
res[268]    14.160  17.249  -19.827    3.041   14.507   25.776   46.182 1.003
res[269]   -59.954  17.435  -93.778  -71.235  -59.514  -48.101  -27.315 1.003
res[270]   -22.812  17.681  -57.121  -34.293  -22.324  -10.851   10.918 1.003
res[271]  -125.907  17.458 -160.933 -137.188 -125.802 -114.461  -92.053 1.001
res[272]   -62.314  17.229  -97.496  -73.416  -62.473  -50.835  -29.085 1.001
res[273]   -35.666  17.061  -70.145  -46.545  -35.800  -24.427   -3.232 1.001
res[274]    -2.957  16.956  -37.392  -13.883   -2.813    8.075   29.213 1.001
res[275]    -9.344  16.916  -43.467  -20.207   -9.343    1.681   22.885 1.001
res[276]    22.554  16.940  -11.628   11.868   22.593   33.645   54.629 1.001
res[277]    73.779  17.029   39.473   63.039   73.780   84.850  106.184 1.001
res[278]   143.908  17.181  109.484  132.820  143.915  155.171  176.723 1.001
res[279]   100.141  17.395   65.573   88.965  100.141  111.576  133.506 1.001
res[280]   -46.044  17.669  -80.667  -57.400  -45.955  -34.590  -12.407 1.001
res[281]    -5.699  17.595  -39.468  -17.381   -5.770    5.983   28.301 1.001
res[282]     0.419  17.343  -33.174  -10.923    0.288   11.974   34.262 1.001
res[283]   -12.869  17.152  -46.211  -24.317  -13.049   -1.401   20.398 1.001
res[284]    12.545  17.023  -20.331    1.056   12.458   24.113   45.686 1.001
res[285]    54.857  16.958   22.320   43.310   54.741   66.049   88.406 1.001
res[286]    12.136  16.958  -20.207    0.644   11.947   23.530   45.348 1.001
res[287]   -10.984  17.022  -43.192  -22.616  -10.977    0.294   22.275 1.001
res[288]     4.979  17.150  -27.724   -6.762    5.048   16.454   38.007 1.001
res[289]   -29.791  17.340  -62.952  -41.651  -29.779  -18.299    3.958 1.001
res[290]   -25.671  17.591  -59.797  -37.708  -25.457  -14.066    8.558 1.001
res[291]    28.550  17.963   -7.697   16.728   28.845   40.812   62.212 1.001
res[292]    -9.089  17.712  -45.024  -20.553   -8.933    2.924   24.333 1.001
res[293]   -49.732  17.521  -85.124  -61.052  -49.636  -38.060  -16.508 1.001
res[294]   -58.677  17.390  -93.259  -69.572  -58.567  -47.134  -25.799 1.001
res[295]   -57.955  17.322  -92.707  -68.933  -57.842  -46.344  -25.049 1.001
res[296]    -3.734  17.317  -38.668  -14.858   -3.653    7.873   29.033 1.001
res[297]    69.494  17.375   34.053   58.293   69.643   81.139  102.719 1.001
res[298]    42.465  17.496    6.356   31.206   42.578   54.056   75.774 1.001
res[299]     7.231  17.678  -29.387   -4.288    7.395   19.167   40.819 1.001
res[300]   -23.337  17.919  -60.558  -34.994  -23.228  -11.302   10.882 1.001
res[301]   -51.903  17.982  -88.248  -63.876  -52.309  -39.714  -17.303 1.001
res[302]   -37.852  17.747  -73.574  -49.599  -38.337  -25.637   -3.466 1.001
res[303]     9.383  17.571  -25.880   -2.168    9.060   21.337   43.272 1.001
res[304]     6.786  17.457  -27.982   -4.692    6.466   18.794   40.534 1.001
res[305]   -83.287  17.404 -117.576  -94.688  -83.511  -71.384  -50.132 1.001
res[306]   -56.131  17.415  -90.605  -67.604  -56.382  -44.092  -22.768 1.001
res[307]    29.387  17.488   -5.255   17.630   29.168   41.485   62.770 1.001
res[308]    97.884  17.624   63.110   86.054   97.739  110.020  131.499 1.001
res[309]    53.516  17.820   18.505   41.582   53.468   65.807   87.676 1.001
res[310]   -20.057  18.074  -55.492  -32.169  -20.032   -7.574   14.442 1.001
res[311]   101.300  17.549   67.220   89.168  101.137  113.302  135.233 1.001
res[312]    83.175  17.312   49.506   71.360   82.808   95.098  116.774 1.001
res[313]    71.785  17.135   37.980   60.064   71.677   83.600  105.227 1.001
res[314]    88.437  17.022   55.339   76.901   88.506  100.014  122.116 1.001
res[315]    -0.110  16.972  -33.336  -11.341    0.071   11.200   33.192 1.001
res[316]   -26.794  16.987  -60.411  -38.159  -26.771  -15.599    6.418 1.001
res[317]   -88.891  17.067 -122.585 -100.389  -88.824  -77.593  -55.404 1.001
res[318]   -96.339  17.209 -130.866 -107.673  -96.328  -84.910  -62.716 1.001
res[319]   -61.837  17.414  -96.679  -73.288  -61.840  -50.019  -28.051 1.001
res[320]  -108.552  17.679 -143.474 -120.010 -108.694  -96.508  -74.344 1.001
res[321]   -74.410  17.625 -109.575  -85.816  -74.485  -62.663  -39.753 1.001
res[322]   -41.400  17.380  -75.871  -52.760  -41.519  -29.797   -7.178 1.001
res[323]   -74.807  17.196 -108.901  -86.261  -74.760  -63.308  -41.163 1.001
res[324]   -45.144  17.074  -79.209  -56.502  -45.026  -33.729  -12.103 1.001
res[325]     4.537  17.016  -29.468   -6.710    4.542   15.856   37.618 1.001
res[326]    59.794  17.022   26.116   48.375   59.758   71.282   93.240 1.001
res[327]    40.165  17.092    6.814   28.767   40.355   51.655   73.479 1.001
res[328]    30.285  17.226   -3.464   18.692   30.679   41.872   63.673 1.001
res[329]    22.589  17.422  -10.860   10.756   23.053   34.110   56.454 1.001
res[330]    92.703  17.678   58.616   80.766   93.284  104.245  126.965 1.001
res[331]    95.283  17.977   59.650   83.167   95.192  107.412  130.221 1.001
res[332]   112.719  17.720   77.424  100.669  112.724  124.693  147.222 1.001
res[333]    70.443  17.522   35.436   58.422   70.438   82.385  104.268 1.001
res[334]    69.514  17.384   34.563   57.587   69.353   81.410  102.556 1.001
res[335]    70.134  17.309   35.643   58.219   69.855   81.817  103.521 1.001
res[336]    -1.755  17.297  -35.726  -13.635   -1.976   10.018   31.941 1.001
res[337]   -93.736  17.349 -127.604 -105.548  -93.768  -81.835  -60.000 1.001
res[338]   -48.951  17.463  -82.779  -60.869  -49.029  -36.839  -14.994 1.001
res[339]  -121.819  17.638 -156.071 -133.910 -122.015 -109.591  -87.286 1.001
res[340]  -103.701  17.874 -138.279 -115.672 -103.713  -91.305  -68.978 1.001
res[341]   -69.317  18.059 -103.485  -81.639  -69.415  -57.406  -33.457 1.001
res[342]   -32.346  17.788  -65.785  -44.511  -32.362  -20.467    2.839 1.001
res[343]   -21.629  17.575  -54.850  -33.597  -21.808   -9.882   13.449 1.001
res[344]   -59.307  17.422  -92.307  -71.080  -59.415  -47.483  -24.425 1.001
res[345]    -3.627  17.331  -36.621  -15.370   -3.984    8.092   31.060 1.001
res[346]    78.394  17.304   46.202   66.537   78.021   90.261  112.620 1.001
res[347]   100.999  17.339   68.667   89.254  100.770  113.032  135.290 1.001
res[348]   -22.296  17.438  -55.281  -34.095  -22.507  -10.328   12.264 1.001
res[349]    46.092  17.598   12.847   34.241   46.001   58.157   80.422 1.001
res[350]    27.883  17.819   -5.915   15.648   27.735   40.121   63.077 1.001
sigma       55.917   2.244   51.705   54.351   55.829   57.468   60.369 1.003
sigma.B     64.474   8.406   50.251   58.466   63.695   69.675   83.144 1.006
deviance  3809.753   9.145 3794.047 3803.114 3809.077 3815.608 3829.461 1.003
          n.eff
beta       3000
beta0      3000
gamma[1]   1800
gamma[2]   3000
gamma[3]   2100
gamma[4]   3000
gamma[5]   3000
gamma[6]   3000
gamma[7]   3000
gamma[8]   3000
gamma[9]   2800
gamma[10]  3000
gamma[11]  2100
gamma[12]  3000
gamma[13]  3000
gamma[14]  3000
gamma[15]  2500
gamma[16]  1700
gamma[17]  1700
gamma[18]  1800
gamma[19]  3000
gamma[20]  1500
gamma[21]  3000
gamma[22]  3000
gamma[23]  3000
gamma[24]  3000
gamma[25]  3000
gamma[26]  1700
gamma[27]  1500
gamma[28]  3000
gamma[29]  2300
gamma[30]  3000
gamma[31]  3000
gamma[32]  3000
gamma[33]  2500
gamma[34]  3000
gamma[35]  3000
res[1]     1100
res[2]     1000
res[3]      990
res[4]      940
res[5]      910
res[6]      880
res[7]      860
res[8]      850
res[9]      840
res[10]     840
res[11]    3000
res[12]    3000
res[13]    3000
res[14]    3000
res[15]    3000
res[16]    3000
res[17]    3000
res[18]    3000
res[19]    3000
res[20]    3000
res[21]    1200
res[22]    1200
res[23]    1100
res[24]    1000
res[25]    1000
res[26]     970
res[27]     950
res[28]    1200
res[29]     920
res[30]     920
res[31]    3000
res[32]    3000
res[33]    3000
res[34]    3000
res[35]    3000
res[36]    3000
res[37]    3000
res[38]    3000
res[39]    3000
res[40]    3000
res[41]    3000
res[42]    3000
res[43]    3000
res[44]    3000
res[45]    3000
res[46]    3000
res[47]    3000
res[48]    3000
res[49]    3000
res[50]    3000
res[51]    3000
res[52]    3000
res[53]    3000
res[54]    3000
res[55]    3000
res[56]    3000
res[57]    3000
res[58]    3000
res[59]    3000
res[60]    3000
res[61]    3000
res[62]    3000
res[63]    3000
res[64]    3000
res[65]    3000
res[66]    3000
res[67]    3000
res[68]    3000
res[69]    3000
res[70]    3000
res[71]    3000
res[72]    3000
res[73]    3000
res[74]    3000
res[75]    3000
res[76]    3000
res[77]    3000
res[78]    3000
res[79]    3000
res[80]    3000
res[81]    1100
res[82]    1100
res[83]    1000
res[84]     970
res[85]     930
res[86]     930
res[87]    1100
res[88]     860
res[89]     850
res[90]     850
res[91]    3000
res[92]    3000
res[93]    3000
res[94]    3000
res[95]    3000
res[96]    3000
res[97]    3000
res[98]    3000
res[99]    3000
res[100]   3000
res[101]   2700
res[102]   2700
res[103]   2800
res[104]   3000
res[105]   3000
res[106]   3000
res[107]   3000
res[108]   3000
res[109]   3000
res[110]   3000
res[111]   3000
res[112]   3000
res[113]   3000
res[114]   3000
res[115]   3000
res[116]   3000
res[117]   3000
res[118]   3000
res[119]   3000
res[120]   3000
res[121]   3000
res[122]   3000
res[123]   3000
res[124]   3000
res[125]   3000
res[126]   3000
res[127]   3000
res[128]   3000
res[129]   3000
res[130]   3000
res[131]   3000
res[132]   3000
res[133]   3000
res[134]   3000
res[135]   3000
res[136]   3000
res[137]   3000
res[138]   3000
res[139]   3000
res[140]   3000
res[141]   2700
res[142]   2700
res[143]   2800
res[144]   3000
res[145]   3000
res[146]   3000
res[147]   3000
res[148]   3000
res[149]   3000
res[150]   3000
res[151]   1700
res[152]   2200
res[153]   2300
res[154]   1900
res[155]   1900
res[156]   2000
res[157]   2200
res[158]   2300
res[159]   2500
res[160]   2700
res[161]   1500
res[162]   1500
res[163]   1500
res[164]   1000
res[165]   1600
res[166]   1700
res[167]   1800
res[168]   1900
res[169]   2000
res[170]   1600
res[171]   3000
res[172]   1900
res[173]   1900
res[174]   2000
res[175]   2000
res[176]   2100
res[177]   2300
res[178]   2400
res[179]   2600
res[180]   2800
res[181]   3000
res[182]   3000
res[183]   3000
res[184]   3000
res[185]   3000
res[186]   3000
res[187]   3000
res[188]   3000
res[189]   3000
res[190]   3000
res[191]   1500
res[192]    630
res[193]    570
res[194]   1600
res[195]   1700
res[196]   1800
res[197]   1900
res[198]   2000
res[199]   2100
res[200]   2300
res[201]   3000
res[202]   3000
res[203]   3000
res[204]   3000
res[205]   3000
res[206]   3000
res[207]   3000
res[208]   3000
res[209]   3000
res[210]   3000
res[211]   2000
res[212]   1900
res[213]   1700
res[214]   1600
res[215]   1600
res[216]   1700
res[217]   1700
res[218]   1900
res[219]   1400
res[220]   1400
res[221]   3000
res[222]   3000
res[223]   3000
res[224]   3000
res[225]   3000
res[226]   3000
res[227]   3000
res[228]   3000
res[229]   3000
res[230]   3000
res[231]   3000
res[232]   3000
res[233]   3000
res[234]   3000
res[235]   3000
res[236]   3000
res[237]   3000
res[238]   3000
res[239]   3000
res[240]   3000
res[241]   3000
res[242]   3000
res[243]   3000
res[244]   3000
res[245]   3000
res[246]   3000
res[247]   3000
res[248]   3000
res[249]   3000
res[250]   3000
res[251]   1800
res[252]   1800
res[253]   1900
res[254]   1900
res[255]   2000
res[256]   2100
res[257]   2200
res[258]   2400
res[259]   3000
res[260]   2800
res[261]    930
res[262]    880
res[263]    830
res[264]    790
res[265]    750
res[266]    740
res[267]    720
res[268]    710
res[269]    710
res[270]    710
res[271]   3000
res[272]   3000
res[273]   3000
res[274]   3000
res[275]   3000
res[276]   3000
res[277]   3000
res[278]   3000
res[279]   3000
res[280]   3000
res[281]   2600
res[282]   2600
res[283]   2700
res[284]   2800
res[285]   3000
res[286]   3000
res[287]   3000
res[288]   3000
res[289]   3000
res[290]   3000
res[291]   3000
res[292]   3000
res[293]   3000
res[294]   3000
res[295]   3000
res[296]   3000
res[297]   3000
res[298]   3000
res[299]   3000
res[300]   3000
res[301]   3000
res[302]   3000
res[303]   3000
res[304]   3000
res[305]   3000
res[306]   3000
res[307]   3000
res[308]   3000
res[309]   3000
res[310]   3000
res[311]   3000
res[312]   3000
res[313]   3000
res[314]   3000
res[315]   3000
res[316]   3000
res[317]   3000
res[318]   3000
res[319]   3000
res[320]   3000
res[321]   2700
res[322]   2800
res[323]   2900
res[324]   3000
res[325]   3000
res[326]   3000
res[327]   3000
res[328]   3000
res[329]   3000
res[330]   3000
res[331]   3000
res[332]   3000
res[333]   3000
res[334]   3000
res[335]   3000
res[336]   3000
res[337]   3000
res[338]   3000
res[339]   3000
res[340]   3000
res[341]   3000
res[342]   3000
res[343]   3000
res[344]   3000
res[345]   3000
res[346]   3000
res[347]   3000
res[348]   3000
res[349]   3000
res[350]   3000
sigma       700
sigma.B    1000
deviance    720

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 41.8 and DIC = 3851.5
DIC is an estimate of expected predictive error (lower deviance is better).

Matrix parameterisation

> modelString2="
+ model {
+    #Likelihood
+    for (i in 1:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- inprod(beta[],X[i,]) + gamma[Block[i]]
+     res[i] <- y[i]-mu[i]
+    } 
+    
+    #Priors
+    beta ~ dmnorm(a0,A0)
+    for (i in 1:nBlock) {
+      gamma[i] ~ dnorm(0, tau.B) #prior
+    }
+    tau <- pow(sigma,-2)
+    sigma <- z/sqrt(chSq) 
+    z ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq ~ dgamma(0.5, 0.5)
+ 
+    tau.B <- pow(sigma.B,-2)
+    sigma.B <- z/sqrt(chSq.B) 
+    z.B ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq.B ~ dgamma(0.5, 0.5)
+  }
+ "
> 
> ## write the model to a text file
> writeLines(modelString2, con = "matrixModel2.txt")
> 
> Xmat <- model.matrix(~Time,data.rm)
> data.rm.list <- with(data.rm,
+         list(y=y,
+                  Block=as.numeric(Block),
+          X=Xmat,
+          n=nrow(data.rm),
+          nBlock=length(levels(Block)),
+          a0=rep(0,ncol(Xmat)), A0=diag(ncol(Xmat))
+          )
+ )
> 
> params <- c("beta",'gamma',"sigma","sigma.B","res")
> adaptSteps = 1000
> burnInSteps = 3000
> nChains = 2
> numSavedSteps = 3000
> thinSteps = 1
> nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
> 
> data.rm.r2jags.m <- jags(data = data.rm.list, inits = NULL, parameters.to.save = params,
+     model.file = "matrixModel2.txt", n.chains = nChains, n.iter = nIter,
+     n.burnin = burnInSteps, n.thin = thinSteps)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 350
   Unobserved stochastic nodes: 40
   Total graph size: 2521

Initializing model
> 
> print(data.rm.r2jags.m)
Inference for Bugs model at "matrixModel2.txt", fit using jags,
 2 chains, each with 4500 iterations (first 3000 discarded)
 n.sims = 3000 iterations saved
           mu.vect sd.vect     2.5%      25%      50%      75%    97.5%  Rhat
beta[1]      0.118   0.993   -1.825   -0.560    0.127    0.772    2.136 1.002
beta[2]      9.067   1.041    6.989    8.384    9.085    9.754   11.096 1.009
gamma[1]   268.496  28.003  214.165  249.963  268.515  287.365  322.199 1.002
gamma[2]   249.357  27.207  196.523  230.541  249.374  268.093  303.009 1.001
gamma[3]   326.337  28.472  271.718  306.769  325.694  345.668  383.423 1.002
gamma[4]   305.015  28.782  249.411  285.199  305.040  324.243  361.356 1.001
gamma[5]   377.709  27.831  324.335  358.940  377.037  396.616  432.150 1.001
gamma[6]   343.430  27.578  288.864  325.843  342.939  361.444  396.819 1.001
gamma[7]   332.166  27.783  278.248  313.377  331.664  351.388  385.777 1.002
gamma[8]   271.616  27.163  217.184  253.588  271.956  289.628  322.963 1.001
gamma[9]   384.589  27.386  332.007  366.136  384.502  402.093  440.209 1.001
gamma[10]  367.756  27.587  314.137  349.123  368.026  386.501  421.707 1.001
gamma[11]  457.715  27.438  405.652  439.504  457.268  476.098  513.423 1.001
gamma[12]  312.867  27.216  262.216  293.798  313.154  331.284  366.377 1.001
gamma[13]  258.804  27.818  201.324  240.475  259.644  276.264  313.330 1.002
gamma[14]  279.449  28.042  224.175  260.014  279.067  298.411  335.404 1.001
gamma[15]  337.363  28.099  284.866  317.715  336.998  356.109  394.348 1.001
gamma[16]  416.383  27.354  363.925  397.826  416.595  434.011  473.215 1.002
gamma[17]  403.080  27.501  349.003  384.921  403.079  421.627  456.974 1.002
gamma[18]  237.418  27.946  183.713  218.161  237.435  256.648  290.420 1.001
gamma[19]  323.886  27.728  270.777  305.233  323.114  342.641  377.533 1.001
gamma[20]  275.462  28.013  220.684  255.987  275.838  294.317  330.957 1.001
gamma[21]  194.059  27.398  141.539  175.756  194.167  211.807  248.181 1.002
gamma[22]  289.821  27.693  236.092  271.210  289.707  307.332  344.403 1.001
gamma[23]  217.663  28.050  164.441  197.887  217.932  235.812  274.184 1.002
gamma[24]  212.748  27.923  159.816  193.937  211.567  231.708  268.892 1.001
gamma[25]  255.843  26.839  204.539  237.181  255.916  274.542  309.265 1.003
gamma[26]  206.396  28.055  152.547  187.160  206.227  225.611  261.170 1.001
gamma[27]  327.559  27.486  273.350  309.633  327.509  345.364  380.375 1.002
gamma[28]  384.714  27.705  330.300  365.884  384.552  403.155  437.847 1.003
gamma[29]  304.783  28.077  251.111  284.961  304.255  323.437  360.045 1.001
gamma[30]  225.251  28.017  168.742  206.532  226.000  243.944  279.897 1.001
gamma[31]  242.036  28.513  184.632  223.534  242.114  260.552  297.438 1.001
gamma[32]  250.059  27.603  197.965  231.432  249.373  268.506  305.275 1.003
gamma[33]  321.386  27.863  266.714  302.647  321.462  340.165  377.552 1.003
gamma[34]  368.704  27.907  314.047  350.229  368.434  387.452  423.293 1.001
gamma[35]  365.968  27.738  311.338  347.444  365.930  385.018  420.415 1.005
res[1]       4.433  27.784  -48.877  -14.717    4.562   22.705   58.869 1.002
res[2]      34.393  27.628  -18.119   15.474   34.478   52.583   88.788 1.002
res[3]     -17.044  27.510  -69.264  -35.707  -16.883    0.972   37.088 1.002
res[4]     -19.010  27.431  -70.907  -37.537  -18.756   -0.634   35.109 1.002
res[5]      22.692  27.391  -29.376    4.107   22.942   41.160   77.308 1.001
res[6]      10.582  27.391  -41.038   -7.968   10.837   28.907   65.778 1.001
res[7]     -70.449  27.431 -121.778  -88.991  -69.960  -52.214  -14.834 1.001
res[8]      21.008  27.509  -30.615    2.408   21.241   39.197   76.881 1.001
res[9]      12.673  27.627  -39.023   -6.003   12.972   31.145   68.777 1.001
res[10]     22.988  27.783  -29.250    4.168   22.937   41.683   79.296 1.001
res[11]    -40.986  26.942  -94.515  -59.717  -41.262  -22.608   11.356 1.001
res[12]    -57.652  26.765 -110.392  -76.201  -57.806  -39.438   -5.452 1.001
res[13]     14.195  26.627  -37.606   -4.203   14.145   32.218   66.504 1.001
res[14]    -29.460  26.529  -81.001  -47.601  -29.683  -11.482   22.755 1.001
res[15]    -14.892  26.472  -66.711  -32.886  -15.051    3.149   37.784 1.001
res[16]     15.000  26.456  -36.965   -3.255   14.703   33.011   68.048 1.001
res[17]    -11.656  26.481  -63.530  -29.872  -12.202    6.352   41.041 1.001
res[18]     -2.098  26.546  -53.836  -20.736   -2.884   15.886   51.209 1.001
res[19]     40.030  26.652  -11.158   21.366   39.102   58.171   93.729 1.001
res[20]    112.896  26.799   61.236   94.306  112.162  131.227  167.261 1.001
res[21]   -173.333  28.266 -230.343 -192.502 -172.496 -153.934 -119.008 1.002
res[22]   -146.973  28.095 -203.369 -166.084 -146.394 -127.937  -93.044 1.002
res[23]    -66.390  27.961 -122.215  -85.409  -65.713  -47.492  -12.298 1.002
res[24]    -62.498  27.865 -117.833  -81.428  -61.625  -43.630   -7.786 1.001
res[25]    -16.892  27.808  -71.818  -35.704  -16.206    2.104   37.863 1.001
res[26]     37.026  27.790  -17.773   18.140   37.505   55.996   91.643 1.001
res[27]     73.127  27.811   18.616   54.532   73.671   91.725  127.823 1.001
res[28]    125.466  27.871   71.518  106.855  125.737  144.036  180.266 1.001
res[29]    132.725  27.969   78.584  113.924  133.004  151.381  187.753 1.001
res[30]    128.171  28.106   73.520  109.247  128.243  147.168  184.155 1.001
res[31]   -215.417  28.519 -271.907 -234.515 -215.477 -196.031 -160.014 1.001
res[32]   -172.464  28.331 -228.398 -191.559 -172.626 -153.223 -117.470 1.001
res[33]   -107.754  28.180 -163.519 -126.542 -107.988  -88.710  -53.225 1.001
res[34]    -46.543  28.067 -102.122  -65.407  -46.791  -27.844    8.278 1.001
res[35]    -47.922  27.992 -103.460  -66.793  -48.157  -29.112    7.029 1.001
res[36]     71.901  27.955   16.630   52.816   71.633   90.451  126.471 1.001
res[37]     97.292  27.958   42.344   78.098   97.282  115.785  151.917 1.001
res[38]    113.786  27.999   59.110   94.749  113.857  132.397  169.334 1.001
res[39]    188.465  28.078  134.278  169.268  188.455  207.026  244.767 1.001
res[40]    152.859  28.196   98.827  133.783  152.725  171.438  209.588 1.001
res[41]   -138.298  27.589 -192.214 -157.018 -137.668 -119.884  -85.511 1.001
res[42]     -0.435  27.414  -54.092  -19.066    0.070   17.947   51.755 1.001
res[43]     -5.275  27.278  -58.296  -23.835   -4.661   12.955   46.693 1.001
res[44]     76.163  27.182   23.390   57.931   76.757   94.384  128.160 1.001
res[45]     -8.105  27.124  -60.470  -26.335   -7.526   10.317   43.968 1.001
res[46]     -8.933  27.107  -61.592  -27.005   -8.421    9.419   43.130 1.001
res[47]     23.039  27.130  -29.532    5.000   23.698   41.485   75.548 1.001
res[48]     48.884  27.192   -3.759   30.945   49.511   67.086  101.297 1.001
res[49]      8.462  27.295  -44.134   -9.630    8.812   26.700   61.215 1.001
res[50]     31.557  27.436  -20.973   13.381   31.746   50.020   84.424 1.001
res[51]   -128.097  27.415 -180.770 -146.327 -127.910 -110.576  -73.568 1.001
res[52]    -60.141  27.254 -112.252  -78.239  -59.732  -42.653   -5.561 1.001
res[53]     17.013  27.132  -35.210   -1.266   17.578   34.265   71.525 1.001
res[54]    -82.564  27.049 -134.784 -100.658  -82.293  -65.282  -27.984 1.002
res[55]     14.994  27.006  -36.578   -3.173   15.282   32.302   69.851 1.002
res[56]    -26.609  27.004  -78.326  -44.704  -26.580   -9.314   28.128 1.002
res[57]     29.085  27.041  -22.987   10.998   29.221   46.243   83.856 1.002
res[58]     88.084  27.118   35.657   69.858   88.230  105.263  142.463 1.004
res[59]     96.481  27.235   43.718   78.058   96.354  113.557  150.554 1.004
res[60]     85.123  27.391   32.058   66.554   85.013  102.661  139.651 1.005
res[61]   -114.949  27.627 -168.618 -134.296 -114.394  -96.450  -60.949 1.001
res[62]    -68.590  27.512 -121.492  -87.723  -68.056  -50.063  -14.649 1.001
res[63]     -5.091  27.435  -57.658  -24.187   -4.767   13.171   48.423 1.001
res[64]    -96.888  27.399 -149.185 -115.932  -96.721  -78.893  -43.557 1.001
res[65]   -145.340  27.401 -197.553 -164.388 -145.424 -127.432  -91.411 1.001
res[66]    -26.577  27.444  -78.399  -45.694  -26.594   -8.455   27.555 1.001
res[67]     53.432  27.525    1.678   34.253   53.095   71.854  107.546 1.001
res[68]     91.059  27.646   39.251   71.830   90.780  109.558  145.538 1.003
res[69]    139.341  27.805   87.022  120.095  138.966  157.961  193.921 1.001
res[70]    192.254  28.001  139.613  173.011  191.889  211.174  247.685 1.001
res[71]   -140.697  26.945 -192.131 -158.462 -140.921 -122.682  -86.498 1.001
res[72]    -53.861  26.803 -104.692  -71.382  -53.867  -36.045   -0.026 1.001
res[73]    -59.969  26.701 -110.888  -77.531  -60.116  -42.419   -5.779 1.001
res[74]    -24.888  26.639  -75.800  -42.339  -25.220   -7.759   28.792 1.001
res[75]    -18.925  26.617  -70.569  -36.419  -19.228   -1.820   34.501 1.001
res[76]     48.388  26.637   -3.281   30.722   48.168   65.796  102.013 1.001
res[77]     45.403  26.697   -6.795   27.656   45.293   62.830   99.324 1.001
res[78]     90.550  26.797   38.062   72.851   90.464  108.076  144.636 1.001
res[79]     85.396  26.937   32.309   67.374   84.992  103.181  139.994 1.001
res[80]     50.353  27.116   -2.874   32.360   50.071   68.275  105.635 1.001
res[81]   -163.315  27.183 -218.512 -180.651 -163.268 -145.149 -111.213 1.001
res[82]   -109.549  27.050 -164.522 -126.836 -109.407  -91.446  -57.185 1.001
res[83]    -13.182  26.957  -68.058  -30.309  -13.185    4.672   38.759 1.001
res[84]    -14.857  26.904  -69.375  -32.207  -15.130    2.965   37.149 1.001
res[85]     40.982  26.891  -12.987   23.736   40.580   58.715   93.178 1.001
res[86]    132.360  26.918   78.198  115.020  131.926  150.358  184.466 1.001
res[87]     96.805  26.986   42.362   79.456   96.215  114.807  149.347 1.001
res[88]     20.474  27.093  -33.942    3.067   19.924   38.470   72.341 1.001
res[89]     17.776  27.240  -36.170    0.154   17.205   35.757   69.869 1.001
res[90]     21.069  27.426  -33.512    3.375   20.566   38.921   73.507 1.001
res[91]   -139.087  27.401 -192.251 -157.857 -139.429 -120.731  -85.751 1.001
res[92]    -36.783  27.249  -89.549  -55.363  -37.012  -18.606   16.371 1.001
res[93]    -52.676  27.136 -105.357  -71.142  -52.862  -34.505    0.149 1.001
res[94]    -26.246  27.062  -78.812  -44.487  -26.694   -8.200   27.216 1.001
res[95]    -30.949  27.028  -83.569  -49.064  -31.432  -12.883   22.822 1.001
res[96]    -18.119  27.035  -70.366  -36.181  -18.630   -0.071   35.799 1.001
res[97]     11.013  27.081  -42.062   -7.308   10.371   29.090   64.960 1.001
res[98]     89.212  27.167   36.723   70.825   88.487  107.474  143.357 1.001
res[99]     95.614  27.293   43.369   77.242   94.962  114.096  150.199 1.001
res[100]   134.498  27.457   80.966  116.052  134.051  153.002  189.417 1.001
res[101]   -76.108  27.204 -130.082  -94.332  -75.847  -57.878  -24.353 1.001
res[102]   -34.250  27.032  -87.681  -52.048  -34.152  -16.151   17.341 1.001
res[103]   -71.035  26.900 -124.029  -88.589  -71.012  -53.026  -19.557 1.001
res[104]    13.787  26.807  -38.887   -3.638   13.685   31.994   65.003 1.001
res[105]   -14.560  26.754  -67.146  -32.253  -14.734    3.512   36.383 1.001
res[106]   -19.293  26.742  -72.025  -36.855  -19.308   -1.263   32.398 1.001
res[107]    56.600  26.770    3.888   38.981   56.413   74.590  108.690 1.001
res[108]   116.002  26.839   62.639   98.370  115.911  133.987  168.211 1.001
res[109]    48.507  26.948   -4.467   30.803   48.406   66.425  100.281 1.001
res[110]     9.543  27.096  -43.482   -8.442    9.395   27.710   61.257 1.001
res[111]  -132.106  27.068 -185.510 -150.360 -132.361 -113.290  -81.768 1.002
res[112]   -91.911  26.941 -144.896 -110.112  -92.368  -73.497  -41.675 1.002
res[113]   -56.418  26.853 -110.121  -74.546  -56.821  -38.144   -6.066 1.002
res[114]   -13.873  26.805  -66.754  -32.098  -14.064    4.250   36.754 1.002
res[115]    33.751  26.798  -19.527   15.438   33.439   51.710   84.504 1.002
res[116]    18.335  26.831  -34.945   -0.042   17.997   36.508   69.231 1.003
res[117]    -9.402  26.904  -62.947  -27.827   -9.801    8.739   41.996 1.003
res[118]    73.620  27.018   20.102   55.043   73.349   91.735  125.589 1.003
res[119]   130.437  27.170   77.165  111.769  130.243  148.448  183.288 1.003
res[120]    67.067  27.362   14.395   48.375   66.985   85.315  120.963 1.004
res[121]   -62.784  27.602 -117.124  -80.547  -63.554  -44.541   -5.745 1.002
res[122]   -79.284  27.432 -132.901  -96.959  -80.152  -61.213  -23.216 1.002
res[123]  -147.038  27.300 -200.594 -164.755 -147.987 -129.337  -91.122 1.002
res[124]   -86.011  27.208 -139.363 -103.583  -87.079  -68.166  -30.561 1.001
res[125]   -44.469  27.155  -96.901  -61.934  -45.519  -26.882   10.835 1.001
res[126]    -7.957  27.143  -60.233  -25.501   -9.007    9.627   47.109 1.001
res[127]    41.116  27.170  -10.769   23.449   40.140   58.470   96.951 1.001
res[128]   117.876  27.237   66.278   99.824  117.103  135.377  173.458 1.003
res[129]   149.628  27.343   98.060  131.402  148.903  167.283  205.474 1.002
res[130]   138.114  27.489   85.817  119.560  137.512  155.791  194.856 1.002
res[131]  -111.510  27.838 -167.229 -130.194 -111.277  -92.014  -56.583 1.002
res[132]   -85.928  27.671 -141.664 -104.693  -85.568  -66.741  -30.867 1.002
res[133]   -39.184  27.543  -95.068  -57.826  -39.100  -20.158   15.452 1.002
res[134]   -65.087  27.454 -120.528  -83.757  -65.031  -46.398  -10.806 1.002
res[135]   -38.854  27.404  -94.235  -57.313  -38.905  -20.233   15.932 1.002
res[136]    43.897  27.394  -11.085   25.502   43.662   62.725   99.335 1.002
res[137]    60.508  27.423    5.661   42.259   60.002   79.149  116.033 1.002
res[138]    48.476  27.491   -6.511   30.204   47.961   67.171  103.838 1.002
res[139]    64.161  27.599    9.521   45.658   63.389   82.927  119.918 1.002
res[140]   150.425  27.745   95.799  131.748  149.642  169.252  206.380 1.003
res[141]   -91.850  27.881 -148.268 -110.602  -91.492  -72.233  -39.029 1.001
res[142]   -17.357  27.708  -73.303  -36.036  -17.162    1.868   35.260 1.001
res[143]    -8.456  27.572  -63.972  -27.298   -8.090   10.649   43.948 1.001
res[144]    16.411  27.476  -38.675   -2.242   16.566   35.463   68.828 1.001
res[145]    -7.945  27.418  -62.620  -26.453   -7.775   11.028   44.009 1.001
res[146]    34.948  27.400  -19.537   16.194   35.119   53.767   87.148 1.001
res[147]    17.864  27.422  -36.519   -0.769   17.959   36.562   70.397 1.001
res[148]   -37.328  27.483  -92.124  -55.967  -37.389  -18.410   15.393 1.001
res[149]     0.133  27.583  -54.930  -18.289    0.029   18.855   53.369 1.001
res[150]   125.181  27.722   69.744  106.772  124.931  144.185  179.294 1.001
res[151]   -45.913  27.112 -101.600  -63.384  -46.504  -27.406    6.233 1.002
res[152]    14.074  26.949  -41.422   -3.190   13.364   32.310   66.645 1.001
res[153]    25.930  26.825  -28.692    8.617   25.258   43.842   78.335 1.001
res[154]    -8.769  26.741  -62.803  -25.915   -9.134    8.810   43.512 1.001
res[155]   -12.026  26.698  -66.536  -29.146  -12.358    5.538   40.695 1.001
res[156]   -13.996  26.695  -68.442  -31.190  -14.374    3.511   39.057 1.001
res[157]   -21.157  26.732  -75.444  -38.415  -21.645   -3.442   32.176 1.001
res[158]    35.944  26.810  -17.938   18.462   35.409   53.573   89.820 1.001
res[159]    -9.265  26.928  -63.044  -26.958   -9.879    8.560   44.932 1.001
res[160]    61.941  27.086    8.272   44.148   61.146   79.841  116.112 1.001
res[161]  -111.720  27.299 -164.219 -130.281 -111.941  -93.352  -57.961 1.002
res[162]  -161.840  27.142 -214.803 -180.270 -162.077 -143.321 -108.335 1.002
res[163]  -105.843  27.025 -158.096 -124.374 -106.136  -87.584  -52.496 1.002
res[164]    32.678  26.947  -18.932   14.279   32.115   50.832   85.368 1.002
res[165]    -1.691  26.909  -53.352  -20.188   -2.310   16.517   51.266 1.002
res[166]    45.107  26.912   -6.211   26.808   44.518   63.279   97.999 1.002
res[167]     5.376  26.954  -46.126  -12.677    4.882   23.392   58.513 1.002
res[168]    72.424  27.037   20.741   54.308   72.204   90.456  125.956 1.002
res[169]    87.205  27.159   35.982   68.965   87.205  105.288  141.275 1.002
res[170]   173.321  27.321  121.880  155.039  173.059  191.695  228.030 1.003
res[171]   -26.750  27.741  -79.583  -45.635  -26.886   -7.432   26.470 1.001
res[172]   -18.152  27.612  -70.442  -37.010  -18.270    0.993   35.150 1.001
res[173]   -30.868  27.522  -82.411  -49.842  -30.897  -11.752   22.359 1.001
res[174]   -20.371  27.471  -72.112  -39.216  -20.399   -1.090   32.975 1.001
res[175]   -29.402  27.460  -81.284  -48.575  -29.553   -9.855   23.536 1.001
res[176]   -42.654  27.488  -94.488  -61.958  -42.932  -23.151   10.118 1.001
res[177]    21.559  27.555  -30.491    2.158   21.213   40.881   74.590 1.001
res[178]    -6.198  27.662  -58.988  -25.815   -6.244   13.036   47.593 1.001
res[179]   104.363  27.807   51.627   84.332  104.497  123.625  158.553 1.001
res[180]    63.487  27.989   10.897   43.415   63.657   82.927  118.256 1.002
res[181]  -116.426  27.477 -170.058 -135.254 -115.973  -98.040  -63.859 1.001
res[182]   -61.964  27.315 -115.028  -80.596  -61.367  -43.741   -8.748 1.001
res[183]   -11.520  27.191  -65.196  -29.915  -10.904    6.736   41.954 1.001
res[184]   -11.403  27.107  -64.157  -29.928  -10.938    7.033   42.438 1.001
res[185]    11.940  27.063  -40.738   -6.618   12.599   30.143   65.828 1.001
res[186]    78.890  27.059   26.180   60.481   79.578   97.004  132.471 1.001
res[187]    15.749  27.094  -36.961   -2.786   16.169   33.923   68.912 1.001
res[188]     3.160  27.170  -49.311  -15.573    3.325   21.357   56.992 1.001
res[189]    52.382  27.285   -0.331   33.687   52.699   70.616  106.503 1.001
res[190]    57.688  27.439    4.937   38.848   57.982   75.882  112.335 1.001
res[191]   -39.277  27.823  -93.946  -58.037  -39.479  -19.984   14.965 1.001
res[192]     9.409  27.662  -44.962   -9.063    9.312   28.464   63.569 1.001
res[193]    22.489  27.540  -32.351    3.890   22.479   41.387   76.748 1.001
res[194]   -41.856  27.457  -96.305  -60.387  -41.784  -23.180   12.866 1.001
res[195]   -12.452  27.413  -66.510  -31.155  -12.405    6.113   42.359 1.001
res[196]    -1.869  27.408  -55.985  -20.561   -1.913   16.685   53.271 1.001
res[197]    21.489  27.443  -32.127    2.711   21.389   39.871   76.814 1.001
res[198]    -9.600  27.518  -62.931  -28.632   -9.617    8.338   45.590 1.002
res[199]    53.727  27.631   -0.059   34.427   53.676   71.736  108.974 1.002
res[200]    31.695  27.783  -22.231   12.143   31.743   49.979   87.523 1.002
res[201]   -86.890  27.181 -140.177 -104.770  -87.118  -68.603  -35.106 1.001
res[202]   -30.563  27.031  -83.642  -48.350  -31.011  -12.345   20.910 1.002
res[203]   -39.562  26.920  -92.523  -57.243  -39.987  -21.390   11.887 1.002
res[204]    25.751  26.850  -27.313    8.074   25.467   43.881   77.877 1.002
res[205]   -12.787  26.819  -65.884  -30.405  -13.246    5.237   39.562 1.002
res[206]  -114.094  26.829 -166.864 -131.854 -114.195  -96.423  -61.166 1.002
res[207]   -27.215  26.880  -80.040  -44.971  -27.431   -9.544   26.079 1.002
res[208]    65.253  26.970   12.034   47.409   65.333   82.874  119.333 1.002
res[209]    87.094  27.100   33.577   69.194   87.081  104.945  141.467 1.003
res[210]   148.253  27.269   94.421  130.305  148.113  166.284  203.466 1.002
res[211]  -161.236  27.503 -215.889 -178.736 -161.122 -142.886 -108.375 1.001
res[212]  -129.522  27.348 -183.477 -146.819 -129.294 -111.282  -76.871 1.001
res[213]  -190.682  27.231 -243.930 -208.062 -190.524 -172.562 -138.450 1.001
res[214]   -79.289  27.154 -132.322  -96.647  -79.209  -61.143  -27.250 1.001
res[215]    47.949  27.116   -4.985   30.474   48.166   66.107  100.695 1.001
res[216]   114.631  27.119   61.912   96.825  114.839  132.864  167.496 1.001
res[217]   117.095  27.161   63.624   99.215  117.216  135.276  169.783 1.001
res[218]   122.954  27.243   68.976  105.183  123.040  140.999  175.833 1.001
res[219]    57.691  27.364    2.985   40.012   57.620   75.929  110.855 1.001
res[220]   124.990  27.525   70.453  107.235  124.767  143.105  178.606 1.001
res[221]   -28.457  27.809  -85.147  -46.709  -28.877   -8.749   24.203 1.001
res[222]   -40.195  27.642  -96.634  -58.154  -40.465  -20.895   12.714 1.001
res[223]   -83.163  27.513 -139.662 -101.235  -83.213  -63.910  -30.005 1.001
res[224]   -57.737  27.424 -113.878  -75.793  -57.651  -38.683   -5.179 1.001
res[225]   -60.675  27.373 -116.422  -78.560  -60.376  -41.718   -8.273 1.001
res[226]   -36.237  27.363  -91.706  -54.143  -35.975  -17.495   16.397 1.001
res[227]   -53.528  27.391 -108.591  -71.787  -53.313  -34.530   -0.302 1.001
res[228]    82.085  27.460   27.176   63.546   82.505  100.972  135.212 1.001
res[229]   109.582  27.567   54.642   90.994  110.106  128.598  162.713 1.001
res[230]   190.314  27.713  135.634  171.529  190.765  209.465  244.223 1.001
res[231]  -120.252  27.694 -175.928 -138.875 -119.474 -101.313  -68.003 1.001
res[232]  -100.688  27.530 -155.801 -119.226  -99.978  -81.979  -48.238 1.001
res[233]   -91.964  27.405 -146.450 -110.666  -91.280  -73.148  -39.619 1.001
res[234]  -108.171  27.319 -162.489 -126.710 -107.330  -89.436  -55.884 1.001
res[235]    21.527  27.273  -32.755    3.074   22.292   40.459   73.543 1.001
res[236]   -18.622  27.266  -73.068  -36.988  -18.139    0.348   33.406 1.001
res[237]    40.129  27.299  -14.569   21.865   40.523   59.257   92.027 1.001
res[238]   117.718  27.372   62.996   99.591  118.175  136.830  169.568 1.001
res[239]   150.564  27.484   95.949  132.308  150.955  169.833  203.577 1.001
res[240]   128.355  27.634   73.287  109.633  128.714  147.821  181.492 1.001
res[241]   -43.717  26.592  -96.442  -61.945  -43.778  -25.366    7.706 1.002
res[242]   -16.266  26.421  -68.915  -34.345  -16.394    1.657   34.738 1.002
res[243]   -51.824  26.290 -104.250  -69.666  -51.904  -34.191   -0.773 1.002
res[244]   -92.310  26.200 -144.512 -109.770  -92.356  -74.742  -41.007 1.002
res[245]    53.225  26.151    1.165   35.773   53.155   70.545  104.484 1.002
res[246]    35.579  26.143  -16.717   18.099   35.633   52.862   86.869 1.001
res[247]   -33.204  26.177  -85.575  -50.833  -33.098  -15.939   18.146 1.001
res[248]    35.397  26.252  -16.792   17.809   35.174   52.691   87.176 1.001
res[249]    83.314  26.368   30.491   65.761   83.032  100.631  135.542 1.001
res[250]    42.630  26.524   -9.672   24.976   42.309   59.916   95.055 1.001
res[251]  -206.923  27.800 -260.529 -225.880 -206.782 -187.968 -153.986 1.001
res[252]  -157.656  27.629 -211.249 -176.383 -157.567 -138.906 -104.887 1.001
res[253]   -51.271  27.496 -104.867  -70.030  -51.271  -32.412    1.659 1.001
res[254]    12.860  27.402  -40.341   -5.856   13.103   31.625   65.496 1.001
res[255]    11.898  27.348  -41.371   -7.218   11.887   30.504   64.701 1.001
res[256]    -1.747  27.333  -54.813  -20.847   -1.807   16.587   50.997 1.001
res[257]    46.786  27.357   -6.081   27.648   46.704   65.069  100.383 1.001
res[258]    84.114  27.422   31.253   65.110   83.938  102.405  137.813 1.001
res[259]   141.780  27.525   88.996  122.717  141.448  159.933  195.838 1.001
res[260]   140.825  27.667   88.214  121.349  140.766  158.831  195.301 1.001
res[261]   -45.758  27.262  -98.677  -63.456  -45.803  -28.120    8.017 1.002
res[262]   -73.536  27.103 -126.119  -91.086  -73.677  -55.880  -19.998 1.001
res[263]   -66.558  26.984 -118.944  -84.155  -66.732  -49.089  -12.830 1.001
res[264]    -3.221  26.904  -55.491  -20.926   -3.268   14.217   50.110 1.001
res[265]    64.093  26.864   12.275   46.220   64.035   81.487  117.920 1.001
res[266]   -24.292  26.865  -75.837  -42.281  -24.457   -6.832   29.548 1.001
res[267]    11.470  26.906  -40.398   -6.459   11.267   29.011   65.223 1.001
res[268]    69.347  26.987   17.413   51.396   68.896   87.025  123.112 1.001
res[269]    16.854  27.108  -35.179   -1.244   16.533   34.478   70.967 1.001
res[270]    75.619  27.268   23.067   57.629   75.329   93.227  129.655 1.001
res[271]  -225.643  27.480 -278.886 -243.823 -225.694 -207.386 -171.764 1.003
res[272]  -140.428  27.320 -193.206 -158.567 -140.609 -122.309  -86.874 1.003
res[273]   -92.158  27.198 -145.009 -110.346  -92.252  -74.095  -38.573 1.003
res[274]   -37.826  27.116  -90.730  -55.873  -38.139  -19.735   15.494 1.002
res[275]   -22.591  27.074  -75.243  -40.456  -22.982   -4.368   30.621 1.002
res[276]    30.930  27.072  -22.094   13.172   30.419   49.127   84.548 1.002
res[277]   103.777  27.109   50.599   85.820  102.988  121.948  157.742 1.002
res[278]   195.528  27.187  142.165  177.513  194.832  213.701  249.967 1.002
res[279]   173.383  27.304  119.414  155.128  172.627  191.796  228.591 1.001
res[280]    48.822  27.460   -5.443   30.590   47.952   67.271  104.943 1.001
res[281]  -100.761  27.924 -155.664 -119.211 -100.098  -81.104  -47.734 1.001
res[282]   -73.020  27.777 -127.369  -91.266  -72.520  -53.694  -20.563 1.001
res[283]   -64.686  27.669 -119.132  -82.823  -64.075  -45.417  -12.494 1.001
res[284]   -17.650  27.599  -72.206  -35.666  -17.169    1.377   35.020 1.001
res[285]    46.285  27.569   -7.652   28.252   46.815   65.328   99.260 1.001
res[286]    25.187  27.577  -28.572    7.101   25.553   44.083   78.569 1.001
res[287]    23.689  27.625  -29.512    5.427   23.865   42.895   77.192 1.001
res[288]    61.274  27.713    8.220   42.691   61.404   80.495  114.966 1.001
res[289]    48.126  27.838   -5.127   29.210   48.155   67.369  102.003 1.001
res[290]    73.868  28.002   20.444   54.993   73.670   93.109  128.079 1.001
res[291]   -62.043  27.792 -115.822  -80.494  -62.779  -43.673   -4.435 1.001
res[292]   -78.060  27.657 -131.681  -96.527  -78.957  -60.027  -21.083 1.001
res[293]   -97.080  27.560 -150.130 -115.697  -97.882  -79.080  -40.450 1.001
res[294]   -84.403  27.502 -137.305 -103.193  -85.127  -66.419  -27.873 1.001
res[295]   -62.058  27.484 -115.409  -80.974  -62.634  -43.885   -5.395 1.001
res[296]    13.785  27.505  -39.891   -5.375   13.338   32.209   69.885 1.001
res[297]   108.636  27.565   55.243   89.343  108.404  127.130  164.689 1.001
res[298]   103.229  27.665   49.562   83.789  102.942  122.033  159.137 1.001
res[299]    89.617  27.803   35.620   70.171   89.485  108.312  146.507 1.001
res[300]    80.672  27.979   26.414   60.959   80.808   99.562  137.744 1.001
res[301]  -141.883  28.305 -196.509 -160.420 -141.920 -123.442  -84.721 1.001
res[302]  -106.210  28.157 -160.120 -124.453 -106.309  -87.823  -49.394 1.001
res[303]   -37.352  28.047  -91.475  -55.335  -37.285  -19.229   19.895 1.001
res[304]   -18.327  27.975  -72.662  -36.286  -18.237   -0.397   38.713 1.001
res[305]   -86.778  27.942 -140.890 -105.114  -86.791  -68.788  -29.209 1.001
res[306]   -37.999  27.947  -91.849  -56.327  -37.867  -20.131   19.461 1.001
res[307]    69.141  27.991   14.641   50.665   69.178   86.859  126.516 1.001
res[308]   159.261  28.074  104.599  140.704  159.387  177.180  216.359 1.001
res[309]   136.516  28.195   80.873  118.244  136.554  154.819  193.892 1.001
res[310]    84.565  28.354   28.496   66.131   84.693  102.951  142.719 1.002
res[311]     8.480  27.336  -47.030   -9.792    8.978   27.073   59.498 1.003
res[312]    11.977  27.128  -43.062   -5.921   12.427   30.119   62.594 1.003
res[313]    22.210  26.959  -32.858    4.243   22.751   40.219   72.580 1.003
res[314]    60.484  26.829    5.199   42.716   60.942   78.598  110.330 1.002
res[315]    -6.440  26.739  -61.507  -24.246   -6.091   11.508   43.556 1.002
res[316]   -11.502  26.689  -66.443  -29.364  -11.033    6.418   38.397 1.002
res[317]   -51.977  26.680 -107.165  -69.981  -51.390  -33.782   -1.318 1.002
res[318]   -37.802  26.711  -93.392  -55.747  -37.309  -19.574   12.624 1.002
res[319]    18.322  26.783  -37.480    0.309   19.012   36.680   68.620 1.001
res[320]    -6.770  26.895  -62.754  -24.603   -5.931   11.437   44.809 1.001
res[321]  -169.878  27.641 -225.412 -188.557 -170.087 -150.812 -115.846 1.003
res[322]  -115.245  27.475 -170.155 -133.481 -115.610  -96.509  -61.338 1.002
res[323]  -127.030  27.348 -181.991 -145.028 -127.405 -108.460  -73.235 1.002
res[324]   -75.745  27.260 -130.493  -94.150  -76.056  -57.347  -21.877 1.002
res[325]    -4.442  27.211  -59.066  -22.893   -4.731   14.091   49.373 1.002
res[326]    72.438  27.202   17.923   54.115   72.083   90.962  125.624 1.002
res[327]    74.431  27.233   20.320   56.096   74.054   93.130  128.264 1.001
res[328]    86.173  27.303   32.254   68.056   86.103  104.893  140.473 1.001
res[329]   100.100  27.413   46.110   82.027  100.161  118.754  154.311 1.001
res[330]   191.836  27.562  137.087  173.612  191.900  210.527  246.650 1.001
res[331]    -3.568  27.677  -58.111  -22.218   -3.240   14.872   51.330 1.001
res[332]    35.490  27.505  -18.504   16.700   35.467   53.896   89.971 1.001
res[333]    14.836  27.372  -38.705   -3.945   14.883   33.326   69.479 1.001
res[334]    35.529  27.278  -17.795   16.742   35.624   53.845   89.880 1.001
res[335]    57.772  27.224    4.745   39.406   57.712   76.032  112.258 1.001
res[336]     7.505  27.209  -45.201  -11.277    7.241   25.550   61.787 1.001
res[337]   -62.853  27.234 -115.231  -81.656  -63.102  -44.830   -8.470 1.001
res[338]     3.554  27.299  -48.771  -15.196    3.261   21.590   57.887 1.001
res[339]   -47.691  27.403 -100.181  -66.462  -48.175  -29.624    6.870 1.001
res[340]    -7.951  27.546  -60.486  -26.990   -8.417   10.133   47.240 1.001
res[341]  -167.838  27.534 -221.260 -186.767 -167.997 -149.518 -113.530 1.005
res[342]  -109.245  27.409 -162.241 -128.030 -109.426  -91.138  -55.319 1.005
res[343]   -76.905  27.323 -129.335  -95.553  -76.780  -58.727  -23.086 1.004
res[344]   -92.961  27.276 -144.845 -111.585  -92.688  -74.742  -39.601 1.004
res[345]   -15.658  27.269  -67.414  -34.298  -15.410    2.716   37.316 1.004
res[346]    87.984  27.302   36.319   69.290   88.171  106.278  140.902 1.003
res[347]   132.212  27.374   80.204  113.503  132.167  150.380  185.832 1.003
res[348]    30.540  27.486  -21.770   11.869   30.461   48.883   84.590 1.003
res[349]   120.550  27.636   67.565  101.960  120.450  138.866  174.815 1.002
res[350]   123.963  27.825   70.604  105.055  123.691  142.633  178.437 1.002
sigma       86.053   4.095   78.239   83.300   85.994   88.737   94.268 1.013
sigma.B    317.817  38.958  252.469  291.832  313.873  339.888  405.125 1.002
deviance  4111.820  21.476 4069.687 4097.266 4111.765 4126.251 4152.352 1.008
          n.eff
beta[1]    1400
beta[2]     230
gamma[1]   1000
gamma[2]   3000
gamma[3]   1300
gamma[4]   2800
gamma[5]   3000
gamma[6]   3000
gamma[7]   1800
gamma[8]   3000
gamma[9]   1800
gamma[10]  2900
gamma[11]  3000
gamma[12]  1800
gamma[13]  1900
gamma[14]  3000
gamma[15]  3000
gamma[16]  1300
gamma[17]  1100
gamma[18]  3000
gamma[19]  3000
gamma[20]  3000
gamma[21]  3000
gamma[22]  3000
gamma[23]  1800
gamma[24]  3000
gamma[25]   810
gamma[26]  3000
gamma[27]  1500
gamma[28]   570
gamma[29]  3000
gamma[30]  3000
gamma[31]  3000
gamma[32]   660
gamma[33]   730
gamma[34]  3000
gamma[35]   340
res[1]     1100
res[2]     1200
res[3]     1400
res[4]     1700
res[5]     2100
res[6]     2700
res[7]     3000
res[8]     3000
res[9]     3000
res[10]    3000
res[11]    3000
res[12]    3000
res[13]    3000
res[14]    3000
res[15]    3000
res[16]    3000
res[17]    3000
res[18]    3000
res[19]    3000
res[20]    3000
res[21]    1200
res[22]    1400
res[23]    1700
res[24]    2000
res[25]    2500
res[26]    3000
res[27]    3000
res[28]    3000
res[29]    3000
res[30]    3000
res[31]    3000
res[32]    3000
res[33]    3000
res[34]    3000
res[35]    3000
res[36]    3000
res[37]    3000
res[38]    3000
res[39]    3000
res[40]    3000
res[41]    3000
res[42]    3000
res[43]    3000
res[44]    3000
res[45]    3000
res[46]    3000
res[47]    3000
res[48]    3000
res[49]    3000
res[50]    3000
res[51]    3000
res[52]    2800
res[53]    2200
res[54]    1700
res[55]    1400
res[56]    1200
res[57]    1000
res[58]     710
res[59]     630
res[60]     550
res[61]    2000
res[62]    2400
res[63]    3000
res[64]    3000
res[65]    3000
res[66]    3000
res[67]    3000
res[68]    3000
res[69]    3000
res[70]    3000
res[71]    3000
res[72]    3000
res[73]    3000
res[74]    3000
res[75]    3000
res[76]    3000
res[77]    3000
res[78]    3000
res[79]    3000
res[80]    3000
res[81]    2100
res[82]    2600
res[83]    3000
res[84]    3000
res[85]    3000
res[86]    3000
res[87]    3000
res[88]    3000
res[89]    3000
res[90]    3000
res[91]    3000
res[92]    3000
res[93]    3000
res[94]    3000
res[95]    3000
res[96]    3000
res[97]    3000
res[98]    3000
res[99]    3000
res[100]   3000
res[101]   3000
res[102]   3000
res[103]   3000
res[104]   3000
res[105]   3000
res[106]   3000
res[107]   3000
res[108]   3000
res[109]   3000
res[110]   3000
res[111]   1700
res[112]   1400
res[113]   1100
res[114]    980
res[115]    840
res[116]    740
res[117]    650
res[118]    580
res[119]    600
res[120]    480
res[121]   1600
res[122]   2000
res[123]   2400
res[124]   3000
res[125]   3000
res[126]   3000
res[127]   3000
res[128]   3000
res[129]   3000
res[130]   3000
res[131]   3000
res[132]   3000
res[133]   3000
res[134]   3000
res[135]   3000
res[136]   3000
res[137]   3000
res[138]   3000
res[139]   2700
res[140]   3000
res[141]   3000
res[142]   3000
res[143]   3000
res[144]   3000
res[145]   3000
res[146]   3000
res[147]   3000
res[148]   3000
res[149]   3000
res[150]   3000
res[151]   1600
res[152]   1900
res[153]   2300
res[154]   3000
res[155]   3000
res[156]   3000
res[157]   3000
res[158]   3000
res[159]   3000
res[160]   3000
res[161]   1000
res[162]   1200
res[163]   1400
res[164]   1700
res[165]   2100
res[166]   2600
res[167]   3000
res[168]   3000
res[169]   3000
res[170]   3000
res[171]   3000
res[172]   3000
res[173]   3000
res[174]   3000
res[175]   3000
res[176]   3000
res[177]   2700
res[178]   2200
res[179]   2000
res[180]   1500
res[181]   3000
res[182]   3000
res[183]   3000
res[184]   3000
res[185]   3000
res[186]   3000
res[187]   3000
res[188]   3000
res[189]   3000
res[190]   3000
res[191]   3000
res[192]   3000
res[193]   3000
res[194]   3000
res[195]   3000
res[196]   2500
res[197]   2000
res[198]   1700
res[199]   1400
res[200]   1200
res[201]   3000
res[202]   3000
res[203]   2700
res[204]   2100
res[205]   1700
res[206]   1400
res[207]   1200
res[208]   1000
res[209]    890
res[210]   1100
res[211]   3000
res[212]   3000
res[213]   3000
res[214]   3000
res[215]   3000
res[216]   3000
res[217]   3000
res[218]   3000
res[219]   3000
res[220]   3000
res[221]   1800
res[222]   2200
res[223]   2800
res[224]   3000
res[225]   3000
res[226]   3000
res[227]   3000
res[228]   3000
res[229]   3000
res[230]   3000
res[231]   3000
res[232]   3000
res[233]   3000
res[234]   3000
res[235]   3000
res[236]   3000
res[237]   3000
res[238]   3000
res[239]   3000
res[240]   3000
res[241]    900
res[242]   1000
res[243]   1200
res[244]   1400
res[245]   1800
res[246]   2200
res[247]   2800
res[248]   3000
res[249]   3000
res[250]   3000
res[251]   3000
res[252]   3000
res[253]   3000
res[254]   3000
res[255]   3000
res[256]   3000
res[257]   3000
res[258]   3000
res[259]   2600
res[260]   2100
res[261]   1600
res[262]   1900
res[263]   2400
res[264]   3000
res[265]   3000
res[266]   3000
res[267]   3000
res[268]   3000
res[269]   3000
res[270]   3000
res[271]    590
res[272]    650
res[273]    730
res[274]    830
res[275]    960
res[276]   1100
res[277]   1200
res[278]   1500
res[279]   1800
res[280]   2500
res[281]   3000
res[282]   3000
res[283]   3000
res[284]   3000
res[285]   3000
res[286]   3000
res[287]   3000
res[288]   3000
res[289]   3000
res[290]   2800
res[291]   3000
res[292]   3000
res[293]   3000
res[294]   3000
res[295]   3000
res[296]   3000
res[297]   3000
res[298]   3000
res[299]   3000
res[300]   3000
res[301]   3000
res[302]   3000
res[303]   3000
res[304]   3000
res[305]   3000
res[306]   3000
res[307]   3000
res[308]   2700
res[309]   2300
res[310]   1700
res[311]    610
res[312]    670
res[313]    760
res[314]    860
res[315]    990
res[316]   1200
res[317]   1400
res[318]   1700
res[319]   2100
res[320]   2600
res[321]    770
res[322]    870
res[323]   1000
res[324]   1200
res[325]   1400
res[326]   1600
res[327]   2000
res[328]   2600
res[329]   2700
res[330]   3000
res[331]   3000
res[332]   3000
res[333]   3000
res[334]   3000
res[335]   3000
res[336]   3000
res[337]   3000
res[338]   3000
res[339]   3000
res[340]   3000
res[341]    350
res[342]    380
res[343]    420
res[344]    460
res[345]    510
res[346]    570
res[347]    810
res[348]    740
res[349]   1200
res[350]   1500
sigma       130
sigma.B    3000
deviance    200

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 229.5 and DIC = 4341.4
DIC is an estimate of expected predictive error (lower deviance is better).

Given that Time cannot be randomized, there is likely to be a temporal dependency structure to the data. The above analyses assume no temporal dependency - actually, they assume that the variance-covariance matrix demonstrates a structure known as sphericity. Lets specifically model in a first order autoregressive correlation structure in an attempt to accommodate the expected temporal autocorrelation.

> modelString3="
+ model {
+    #Likelihood
+    y[1]~dnorm(mu[1],tau)
+    mu[1] <- eta1[1]
+    eta1[1] ~ dnorm(eta[1], taueps)
+    eta[1] <- inprod(beta[],X[1,]) + gamma[Block[1]]
+    res[1] <- y[1]-mu[1]
+    for (i in 2:n) {
+       y[i]~dnorm(mu[i],tau)
+       mu[i] <- eta1[i]
+       eta1[i] ~ dnorm(temp[i], taueps)
+       temp[i] <- eta[i] + -rho*(mu[i-1]-y[i-1])
+       eta[i] <- inprod(beta[],X[i,]) + gamma[Block[i]]
+     res[i] <- y[i]-mu[i]
+    } 
+    beta ~ dmnorm(a0,A0)
+    for (i in 1:nBlock) {
+      gamma[i] ~ dnorm(0, tau.B) #prior
+    }
+    rho ~ dunif(-1,1)
+    tau <- pow(sigma,-2)
+    sigma <- z/sqrt(chSq) 
+    z ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq ~ dgamma(0.5, 0.5)
+    taueps <- pow(sigma.eps,-2)
+    sigma.eps <- z/sqrt(chSq.eps) 
+    z.eps ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq.eps ~ dgamma(0.5, 0.5)
+    tau.B <- pow(sigma.B,-2)
+    sigma.B <- z/sqrt(chSq.B) 
+    z.B ~ dnorm(0, 0.0016)I(0,)  #1/25^2 = 0.0016
+    chSq.B ~ dgamma(0.5, 0.5)
+    sd.y <- sd(res)
+    sd.block <- sd(gamma)
+  }
+ "
> 
> ## write the model to a text file
> writeLines(modelString3, con = "matrixModel3.txt")
> 
> Xmat <- model.matrix(~Time,data.rm)
> data.rm.list <- with(data.rm,
+         list(y=y,
+                  Block=as.numeric(Block),
+          X=Xmat,
+          n=nrow(data.rm),
+          nBlock=length(levels(Block)),
+          a0=rep(0,ncol(Xmat)), A0=diag(ncol(Xmat))
+          )
+ )
> 
> params <- c("beta",'gamma',"sigma","sigma.B","res",'sigma.eps','rho','sd.y','sd.block')
> adaptSteps = 1000
> burnInSteps = 3000
> nChains = 2
> numSavedSteps = 3000
> thinSteps = 1
> nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
> 
> data.rm.r2jags.mt <- jags(data = data.rm.list, inits = NULL, parameters.to.save = params,
+     model.file = "matrixModel3.txt", n.chains = nChains, n.iter = nIter,
+     n.burnin = burnInSteps, n.thin = thinSteps)
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 350
   Unobserved stochastic nodes: 393
   Total graph size: 3931

Initializing model
> 
> data.rm.mt.mcmc <- data.rm.r2jags.mt$BUGSoutput$sims.matrix
> summary(as.mcmc(data.rm.mt.mcmc[,grep('beta|sigma|rho',colnames(data.rm.mt.mcmc))]))

Iterations = 1:3000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 3000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

              Mean      SD Naive SE Time-series SE
beta[1]     0.1093  1.0090 0.018422       0.018422
beta[2]     9.1917  1.0421 0.019026       0.019026
rho         0.6991  0.1635 0.002985       0.002985
sigma      63.1901  7.6813 0.140241       0.140241
sigma.B   314.5779 38.4731 0.702419       0.702419
sigma.eps  32.8263  9.4551 0.172625       0.172625

2. Quantiles for each variable:

              2.5%      25%       50%      75%    97.5%
beta[1]    -1.8402  -0.5698   0.09792   0.7920   2.0576
beta[2]     7.2315   8.4966   9.19938   9.8961  11.2232
rho         0.4549   0.5508   0.68313   0.8453   0.9746
sigma      49.8811  56.9814  62.70134  69.5262  77.4328
sigma.B   251.2601 286.7395 310.40901 338.2187 398.3167
sigma.eps  14.8598  25.5907  35.65555  40.2488  46.3558
> 
> #head(data.rm.r2jags.mt$BUGSoutput$sims.list[[c('beta','rho','sigma')]]) 
> #print(data.rm.r2jags.mt)
> data.rm.mcmc.list.mt <- as.mcmc(data.rm.r2jags.mt)
> Data.Rm.mcmc.list.mt <- data.rm.mcmc.list.mt
> 
> # R2 calculations
> Xmat <- model.matrix(~Time, data.rm)
> coefs <- data.rm.r2jags.mt$BUGSoutput$sims.list[['beta']]
> fitted <- coefs %*% t(Xmat)
> X.var <- aaply(fitted,1,function(x){var(x)})
> X.var[1:10]
       1        2        3        4        5        6        7        8 
814.3418 591.0181 634.1438 685.5362 900.1883 740.2397 864.5962 435.7952 
       9       10 
672.4743 584.2064 
> 
> Z.var <- data.rm.r2jags.mt$BUGSoutput$sims.list[['sd.block']]^2
> R.var <- data.rm.r2jags.mt$BUGSoutput$sims.list[['sd.y']]^2
> R2.marginal <- (X.var)/(X.var+Z.var+R.var)
> R2.marginal <- data.frame(Mean=mean(R2.marginal), Median=median(R2.marginal), HPDinterval(as.mcmc(R2.marginal)))
> R2.conditional <- (X.var+Z.var)/(X.var+Z.var+R.var)
> R2.conditional <- data.frame(Mean=mean(R2.conditional),
+    Median=median(R2.conditional), HPDinterval(as.mcmc(R2.conditional)))
> R2.block <- (Z.var)/(X.var+Z.var+R.var)
> R2.block <- data.frame(Mean=mean(R2.block), Median=median(R2.block), HPDinterval(as.mcmc(R2.block)))
> R2.res<-(R.var)/(X.var+Z.var+R.var)
> R2.res <- data.frame(Mean=mean(R2.res), Median=median(R2.res), HPDinterval(as.mcmc(R2.res)))
> 
> (r2 <- rbind(R2.block=R2.block, R2.marginal=R2.marginal, R2.res=R2.res, R2.conditional=R2.conditional))
                     Mean     Median      lower     upper
R2.block       0.52595295 0.52197768 0.40289527 0.6376026
R2.marginal    0.07190426 0.06983026 0.03887103 0.1087763
R2.res         0.40214279 0.40351594 0.28261806 0.5200679
R2.conditional 0.59785721 0.59648406 0.47993214 0.7173819

It would appear that the incorporation of a first order autocorrelation structure is indeed appropriate. The degree of correlation between successive points is \(0.733\). Let’s have a look at a summary figure.

> coefs <- data.rm.r2jags.mt$BUGSoutput$sims.list[['beta']]
> newdata <- with(data.rm, data.frame(Time=seq(min(Time, na.rm=TRUE), max(Time, na.rm=TRUE), len=100)))
> Xmat <- model.matrix(~Time, newdata)
> pred <- (coefs %*% t(Xmat))
> pred <- adply(pred, 2, function(x) {
+    data.frame(Mean=mean(x), Median=median(x, na.rm=TRUE), t(quantile(x,na.rm=TRUE)),
+               HPDinterval(as.mcmc(x)),HPDinterval(as.mcmc(x),p=0.5))
+ })
> newdata <- cbind(newdata, pred)
> #Also calculate the partial observations
> Xmat <- model.matrix(~Time, data.rm)
> pred <- colMeans(as.vector(coefs %*% t(Xmat))+data.rm.r2jags.mt$BUGSoutput$sims.list[['res']])
> part.obs <- cbind(data.rm,Median=pred)
> 
> ggplot(newdata, aes(y=Median, x=Time)) +
+   geom_point(data=part.obs, aes(y=Median))+
+   geom_ribbon(aes(ymin=lower, ymax=upper), fill='blue',alpha=0.2) +
+   geom_line()+
+   scale_x_continuous('Time') +
+   scale_y_continuous('Y') +
+   theme_classic() +
+   theme(axis.title.y = element_text(vjust=2, size=rel(1.2)),
+         axis.title.x = element_text(vjust=-2, size=rel(1.2)),
+         plot.margin=unit(c(0.5,0.5,2,2), 'lines'))

References

Plummer, Martyn. 2004. “JAGS: Just Another Gibbs Sampler.”
Su, Yu-Sung, Masanao Yajima, Maintainer Yu-Sung Su, and JAGS SystemRequirements. 2015. “Package ‘R2jags’.” R Package Version 0.03-08, URL Http://CRAN. R-Project. Org/Package= R2jags.
Assistant Professor in Statistics

Related