Jump to main navigation

# Tutorial 9.2b - Nested ANOVA (Bayesian)

23 Dec 2015

If you are completely ontop of the conceptual issues pertaining to Nested ANOVA, and just need to use this tutorial in order to learn about Nested ANOVA in R, you are invited to skip down to the section on Nested ANOVA in R.

## Overview

You are strongly advised to review the information on the nested design in tutorial 9.2a. I am not going to duplicate the overview here.

## Linear models

So called 'random effects' are modelled differently from 'fixed effects' in that rather than estimate their individual effects, we instead estimate the variability due to these 'random effects'. Since technically all variables in a Bayesian framework are random, some prefer to use the terms 'fixed effects' and 'varying effects'. As random factors typically represent 'random' selections of levels (such as a set of randomly selected sites), incorporated in order to account for the dependency structure (observations within sites are more likely to be correlated to one another - not strickly independent) to the data, we are not overly interested in the individual differences between levels of the 'varying' (random) factor. Instead (in addition to imposing a separate correlation structure within each nest), we want to know how much variability is attributed to this level of the design.

The linear models for two and three factor nested design are:
\begin{align} y_{ijk}&=\mu+\alpha_i + \beta_{j(i)} + \varepsilon_{ijk} &\hspace{2em} \varepsilon_{ijk}&\sim\mathcal{N}(0,\sigma^2), \hspace{1em}\beta_{j(i)}\sim\mathcal{N}(0,\sigma_B^2)\\ y_{ijkl}&=\mu+\alpha_i + \beta_{j(i)} + \gamma_{k(j(i))} + \varepsilon_{ijkl} &\hspace{2em} \varepsilon_{ijkl}&\sim\mathcal{N}(0,\sigma^2), \hspace{1em}\beta_{j(i)}\sim\mathcal{N}(0,\sigma_B^2), \hspace{1em}\gamma_{k(j(i))}\sim\mathcal{N}(0,\sigma_C^2)\\ \end{align} where $\mu$ is the overall mean, $\alpha$ is the effect of Factor A, $\beta$ is the variability of Factor B (nested within Factor A), $\gamma$ is the variability of Factor C (nested within Factor B) and $\varepsilon$ is the random unexplained or residual component that is assumed to be normally distributed with a mean of zero and a constant amount of standard deviation ($\sigma^2$).

The subscripts are iterators. For example, the $i$ represents the number of effects to be estimated for Factor A. Thus the first formula can be read as 'the $k$th observation of $y$ is drawn from a normal distribution (with a specific level of variability) and mean proposed to be determined by a base mean ($\mu$ - mean of the first treatment across all nests) plus the effect of the $i$th treatment effect plus the variability 'the model proposes that given a base mean ($\mu$) and knowing the effect of the $i$th treatment (factor A) and which of the $j$th nests within the treatment 'the $k$th observation from Block $j$ (factor B) within treatment effect

## Variance components

As previously alluded to, it can often be useful to determine the relative contribution (to explaining the unexplained variability) of each of the factors as this provides insights into the variability at each different scales.

These contributions are known as Variance components and are estimates of the added variances due to each of the factors. For consistency with leading texts on this topic, I have included estimated variance components for various balanced nested ANOVA designs in the above table. However, variance components based on a modified version of the maximum likelihood iterative model fitting procedure (REML) is generally recommended as this accommodates both balanced and unbalanced designs.

While there are no numerical differences in the calculations of variance components for fixed and random factors, fixed factors are interpreted very differently and arguably have little biological meaning (other to infer relative contribution). For fixed factors, variance components estimate the variance between the means of the specific populations that are represented by the selected levels of the factor and therefore represent somewhat arbitrary and artificial populations. For random factors, variance components estimate the variance between means of all possible populations that could have been selected and thus represents the true population variance.

## Assumptions

Consistent with the hierarchical nature of the design (in which the site means are the appropriate level of replication for the main treatments), it is important that diagnostics associated with each factor reflect this hierarchy. Hence, the replicates for the treatment effects are the site means, and the replicates of the sites are the quadrat means.

The assumptions are:
1. All of the observations are independent - this must be addressed at the design and collection stages. Importantly, to be considered independent replicates, the replicates must be made at the same scale at which the treatment is applied. For example, if the experiment involves subjecting organisms housed in tanks to different water temperatures, then the unit of replication is the individual tanks not the individual organisms in the tanks. The individuals in a tank are strictly not independent with respect to the treatment.
2. The response variable (and thus the residuals) should be normally distributed for each sampled populations (combination of factors). Boxplots of each treatment combination are useful for diagnosing major issues with normality.
3. The response variable should be equally varied (variance should not be related to mean as these are supposed to be estimated separately) for each combination of treatments. Again, boxplots are useful.

## Design balance

In a Bayesian framework, issues of design balance essentially evaporate.

## Nested ANOVA in R

### Scenario and Data

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'). The treatments occur at a spatial scale (over an area) that far exceeds the logistical scale of sampling units (it would take too long to sample at the scale at which the treatments were applied). The treatments occurred at the scale of hectares whereas it was only feasible to sample $y$ using 1m quadrats.

Given that the treatments were naturally occurring (such as soil type), it was not possible to have more than five sites of each treatment type, yet prior experience suggested that the sites in which you intended to sample were very uneven and patchy with respect to $y$.

In an attempt to account for this inter-site variability (and thus maximize the power of the test for the effect of treatment, you decided to employ a nested design in which 10 quadrats were randomly located within each of the five replicate sites per three treatments. 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. Consequently, I have folded (toggled) this section away.

Random data incorporating the following properties
• the number of treatments = 3
• the number of sites per treatment = 5
• the total number of sites = 15
• the sample size per treatment= 5
• the number of quadrats per site = 10
• the mean of the treatments = 40, 70 and 80 respectively
• the variability (standard deviation) between sites of the same treatment = 12
• the variability (standard deviation) between quadrats within sites = 5
library(plyr)
set.seed(1)
nTreat <- 3
nSites <- 15
nSitesPerTreat <- nSites/nTreat
nQuads <- 10
site.sigma <- 12
sigma <- 5
n <- nSites * nQuads
sites <- gl(n=nSites,k=nQuads, lab=paste0('S',1:nSites))
A <- gl(nTreat, nSitesPerTreat*nQuads, n, labels=c('a1','a2','a3'))
a.means <- c(40,70,80)
## the site means (treatment effects) are drawn from normal distributions
## with means of 40, 70 and 80 and standard deviations of 12
A.effects <- rnorm(nSites, rep(a.means,each=nSitesPerTreat),site.sigma)
#A.effects <- a.means %*% t(model.matrix(~A, data.frame(A=gl(nTreat,nSitesPerTreat,nSites))))+rnorm(nSites,0,site.sigma)
Xmat <- model.matrix(~sites -1)
lin.pred <- Xmat %*% c(A.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.nest <- data.frame(y=y, A=A, Sites=sites,Quads=1:length(y))
head(data.nest)  #print out the first six rows of the data set

         y  A Sites Quads
1 32.25789 a1    S1     1
2 32.40160 a1    S1     2
3 37.20174 a1    S1     3
4 36.58866 a1    S1     4
5 35.45206 a1    S1     5
6 37.07744 a1    S1     6

library(ggplot2)
ggplot(data.nest, aes(y=y, x=1)) + geom_boxplot() + facet_grid(.~Sites)


### Exploratory data analysis

#### Normality and Homogeneity of variance

#Effects of treatment
library(plyr)
boxplot(y~A, ddply(data.nest, ~A+Sites,numcolwise(mean, na.rm=T)))

#Site effects
boxplot(y~Sites, ddply(data.nest, ~A+Sites+Quads,numcolwise(mean, na.rm=T)))

## with ggplot2
ggplot(ddply(data.nest, ~A+Sites,numcolwise(mean, na.rm=T)), aes(y=y, x=A)) +
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).

### JAGS

For non-hierarchical linear models, uniform priors on variance (standard deviation) parameters seem to work reasonably well. Gelman (2006) warns that the use of the inverse-gamma family of non-informative priors are very sensitive to $\epsilon$ particularly when variance is close to zero and this may lead to unintentionally informative priors. When the number of groups (treatments or varying/random effects) is large (more than 5), Gelman (2006) advocated the use of either uniform or half-cauchy priors. Yet when the number of groups is low, Gelman (2006) indicates that uniform priors have a tendency to result in inflated variance estimates. Consequently, half-cauchy priors are generally recommended for variances.

Full parameterizationMatrix parameterizationHeirarchical parameterization
$$\begin{array}{rcl} y_{ijk}&\sim&\mathcal{N}(\mu_{ij},\sigma^2)\\ \mu_{ij} &=& \alpha_0 + \alpha_{i} + \beta_{j(i)}\\ \beta_{i{j}}&\sim&\mathcal{N}(0,\sigma_{B}^2)\\ \alpha_0, \alpha_i&\sim&\mathcal{N}(0,100000)\\ \sigma^2, \sigma_{B}^2&\sim&\mathcal{Cauchy}(0,25)\\ \end{array}$$ $$\begin{array}{rcl} y_{ijk}&\sim&\mathcal{N}(\mu_{ij},\sigma^2)\\ \mu_{ij} &=& \alpha\mathbf{X} + \beta_{j(i)}\\ \beta_{i{j}}&\sim&\mathcal{N}(0,\sigma_{B}^2)\\ \alpha&\sim&\mathcal{MVN}(0,100000)\\ \sigma^2, \sigma_{B}^2&\sim&\mathcal{Cauchy}(0,25)\\ \end{array}$$ $$\begin{array}{rcl} y_{ij}&\sim &\mathcal{N}(\beta_{j(i)}, \sigma^2)\\ \beta_{j(i)}&\sim&\mathcal{N}(\mu_{i},\sigma_{B}^2)\\ \mu_{i}&=&\alpha\mathbf{X}\\ \alpha_i&\sim&\mathcal{N}(0, 1000000)\\ \sigma^2, \sigma_{B}^2&\sim&\mathcal{Cauchy}(0,25)\\ \end{array}$$

The full parameterization, shows the effects parameterization in which there is an intercept ($\alpha_0$) and two treatment effects ($\alpha_i$, where $i$ is 1,2).

The matrix parameterization is a compressed notation, In this parameterization, 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 parameterization 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 3x3 matrix of zeros and 100 in the diagonals. $$\mu\sim\left[ \begin{array}{c} 0\\ 0\\ 0\\ \end{array} \right], \hspace{2em} \sigma^2\sim\left[ \begin{array}{ccc} 1000000&0&0\\ 0&1000000&0\\ 0&0&1000000\\ \end{array} \right]$$

In the Heirarchical parameterization, we are indicating two residual layers - one representing the variability in the observed data between individual observations (within sites) and the second representing the variability between site means (within the three treatments).

#### Full effects parameterization

modelString="
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- alpha0 + alpha[A[i]] + beta[site[i]]
}

#Priors
alpha0 ~ dnorm(0, 1.0E-6)
alpha[1] <- 0
for (i in 2:nA) {
alpha[i] ~ dnorm(0, 1.0E-6) #prior
}
for (i in 1:nSite) {
beta[i] ~ dnorm(0, tau.B) #prior
}
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.B <- pow(sigma.B,-2)
sigma.B <-z/sqrt(chSq.B)
z.B ~ dnorm(0, .0016)I(0,)
chSq.B ~ dgamma(0.5, 0.5)
}
"
writeLines(modelString,con="../downloads/BUGSscripts/tut9.2bS3.1.f.txt")

data.nest.list <- with(data.nest,
list(y=y,
site=as.numeric(Sites),
A=as.numeric(A),
n=nrow(data.nest),
nSite=length(levels(Sites)),
nA = length(levels(A))
)
)

params <- c("alpha0","alpha","sigma","sigma.B")
adaptSteps = 1000
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
rnorm(1)

[1] -0.5428819

jags.effects.f.time <- system.time(
data.nest.r2jags.f <- jags(data=data.nest.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/tut9.2bS3.1.f.txt",
n.chains=3,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 503

Initializing model

jags.effects.f.time

   user  system elapsed
7.812   0.008   7.874

print(data.nest.r2jags.f)

Inference for Bugs model at "../downloads/BUGSscripts/tut9.2bS3.1.f.txt", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
alpha[1]   0.000   0.000   0.000   0.000   0.000   0.000   0.000 1.000     1
alpha[2]  30.088   9.059  12.031  24.445  30.049  35.796  48.112 1.001  3000
alpha[3]  37.141   9.101  19.689  31.396  37.121  42.774  55.403 1.002  1700
alpha0    42.171   6.447  29.586  38.114  42.136  46.235  55.100 1.001  3000
sigma      4.428   0.273   3.940   4.241   4.415   4.598   4.996 1.001  3000
sigma.B   14.057   3.179   9.421  11.801  13.572  15.590  22.385 1.001  3000
deviance 869.642   6.183 859.552 865.156 868.846 873.323 883.745 1.002  1300

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 = 19.1 and DIC = 888.7
DIC is an estimate of expected predictive error (lower deviance is better).

data.nest.mcmc.list.f <- as.mcmc(data.nest.r2jags.f)

[1] "alpha0"   "alpha[1]" "alpha[2]" "alpha[3]" "deviance" "sigma"    "sigma.B"


#### Matrix parameterization

modelString="
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- inprod(alpha[],X[i,]) + beta[site[i]]
}

#Priors
alpha ~ dmnorm(a0,A0)
for (i in 1:nSite) {
beta[i] ~ dnorm(0, tau.B) #prior
}
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.B <- pow(sigma.B,-2)
sigma.B <-z/sqrt(chSq.B)
z.B ~ dnorm(0, .0016)I(0,)
chSq.B ~ dgamma(0.5, 0.5)

}
"
writeLines(modelString,con="../downloads/BUGSscripts/tut9.2bS3.1.m.txt")

A.Xmat <- model.matrix(~A,data.nest)
data.nest.list <- with(data.nest,
list(y=y,
site=as.numeric(Sites),
X=A.Xmat,
n=nrow(data.nest),
nSite=length(levels(Sites)),
nA = ncol(A.Xmat),
a0=rep(0,3), A0=diag(0,3)
)
)

params <- c("alpha","sigma","sigma.B")
adaptSteps = 1000
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
rnorm(1)

[1] 0.79071

jags.effects.m.time <- system.time(
data.nest.r2jags.m <- jags(data=data.nest.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/tut9.2bS3.1.m.txt",
n.chains=3,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 965

Initializing model

jags.effects.m.time

   user  system elapsed
9.577   0.008   9.670

print(data.nest.r2jags.m)

Inference for Bugs model at "../downloads/BUGSscripts/tut9.2bS3.1.m.txt", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
alpha[1]  42.362   6.448  29.963  38.129  42.351  46.617  54.771 1.003   900
alpha[2]  29.814   8.935  12.056  23.942  29.743  35.587  47.542 1.001  2300
alpha[3]  36.868   9.024  18.302  31.194  36.774  42.735  54.574 1.002  1400
sigma      4.431   0.270   3.945   4.242   4.417   4.603   5.009 1.001  2500
sigma.B   14.133   3.020   9.599  11.972  13.638  15.759  21.385 1.001  3000
deviance 869.516   5.925 860.027 865.256 868.874 873.104 882.656 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 = 17.6 and DIC = 887.1
DIC is an estimate of expected predictive error (lower deviance is better).

data.nest.mcmc.list.m <- as.mcmc(data.nest.r2jags.m)
Data.Nest.mcmc.list.m <- data.nest.mcmc.list.m

Or, more generally..
modelString="
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- inprod(alpha[],X[i,]) + inprod(beta[], Z[i,])
}

#Priors
alpha ~ dmnorm(a0,A0)
for (i in 1:nZ) {
beta[i] ~ dnorm(0, tau.B) #prior
}
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.B <- pow(sigma.B,-2)
sigma.B <-z/sqrt(chSq.B)
z.B ~ dnorm(0, .0016)I(0,)
chSq.B ~ dgamma(0.5, 0.5)

}
"

A.Xmat <- model.matrix(~A,data.nest)
Zmat <- model.matrix(~-1+Sites, data.nest)
data.nest.list <- with(data.nest,
list(y=y,
X=A.Xmat,
n=nrow(data.nest),
Z=Zmat, nZ=ncol(Zmat),
nA = ncol(A.Xmat),
a0=rep(0,3), A0=diag(0,3)
)
)

params <- c("alpha","sigma","sigma.B",'beta')
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
rnorm(1)

[1] 0.3485099

jags.effects.m2.time <- system.time(
data.nest.r2jags.m2 <- jags(data=data.nest.list,
inits=NULL,
parameters.to.save=params,
model.file=textConnection(modelString),
n.chains=3,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 3231

Initializing model

jags.effects.m2.time

   user  system elapsed
16.541   0.016  16.755

print(data.nest.r2jags.m2)

Inference for Bugs model at "5", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
alpha[1]  42.067   6.578  29.032  37.725  42.160  46.241  55.169 1.001  3000
alpha[2]  30.069   9.259  12.117  24.021  30.031  36.209  48.381 1.002  1800
alpha[3]  37.285   9.101  19.517  31.284  37.507  43.232  54.917 1.002  1500
beta[1]   -8.128   6.675 -21.114 -12.440  -8.222  -3.889   4.897 1.001  3000
beta[2]   -0.612   6.632 -13.657  -4.847  -0.702   3.688  12.406 1.001  3000
beta[3]  -11.405   6.683 -24.807 -15.689 -11.483  -7.171   1.555 1.001  3000
beta[4]   17.700   6.673   4.694  13.437  17.662  22.145  31.058 1.001  3000
beta[5]    3.490   6.698  -9.694  -0.894   3.387   8.082  16.528 1.001  3000
beta[6]  -11.674   6.603 -24.545 -15.902 -11.629  -7.462   1.011 1.002  1900
beta[7]    3.086   6.667  -9.762  -1.372   3.019   7.503  15.970 1.001  2300
beta[8]    9.543   6.607  -3.742   5.282   9.611  13.820  22.654 1.001  2500
beta[9]    2.287   6.676 -10.642  -2.088   2.300   6.640  15.307 1.001  2800
beta[10]  -3.286   6.602 -16.114  -7.658  -3.222   0.915   9.541 1.001  2400
beta[11]  18.438   6.367   5.868  14.225  18.382  22.586  31.068 1.002  2000
beta[12]   4.476   6.433  -8.237   0.294   4.428   8.738  17.024 1.002  1600
beta[13] -10.203   6.400 -22.966 -14.385 -10.263  -6.120   2.182 1.002  1700
beta[14] -27.562   6.431 -40.197 -31.830 -27.430 -23.309 -14.994 1.002  1800
beta[15]  14.658   6.402   2.338  10.408  14.640  18.816  27.198 1.001  2300
sigma      4.424   0.271   3.933   4.237   4.410   4.599   5.003 1.001  3000
sigma.B   14.088   2.992   9.553  11.970  13.675  15.744  21.031 1.001  2800
deviance 869.547   6.092 859.967 865.210 868.616 873.111 883.730 1.002  2000

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 = 18.5 and DIC = 888.1
DIC is an estimate of expected predictive error (lower deviance is better).

data.nest.mcmc.list.m2 <- as.mcmc(data.nest.r2jags.m2)
Data.Nest.mcmc.list.m2 <- data.nest.mcmc.list.m2


#### Hierarchical parameterization

modelString="
model {
#Likelihood (esimating site means (gamma.site)
for (i in 1:n) {
y[i]~dnorm(quad.means[i],tau)
quad.means[i] <- gamma.site[site[i]]
}
for (i in 1:s) {
gamma.site[i] ~ dnorm(site.means[i], tau.site)
site.means[i] <- inprod(beta[],A.Xmat[i,])
}
#Priors
for (i in 1:a) {
beta[i] ~ dnorm(0, 1.0E-6) #prior
}
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.B <- pow(sigma.B,-2)
sigma.B <-z/sqrt(chSq.B)
z.B ~ dnorm(0, .0016)I(0,)
chSq.B ~ dgamma(0.5, 0.5)

tau.site <- pow(sigma.site,-2)
sigma.site <-z/sqrt(chSq.site)
z.site ~ dnorm(0, .0016)I(0,)
chSq.site ~ dgamma(0.5, 0.5)
}
"
writeLines(modelString,con="../downloads/BUGSscripts/tut9.2bS3.1.txt")

A.Xmat <- model.matrix(~A,ddply(data.nest,~Sites,catcolwise(unique)))
data.nest.list <- with(data.nest,
list(y=y,
site=Sites,
A.Xmat= A.Xmat,
n=nrow(data.nest),
s=length(levels(Sites)),
a = ncol(A.Xmat)
)
)

params <- c("beta","sigma","sigma.site")
adaptSteps = 1000
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
rnorm(1)

[1] 0.3643516

jags.effects.simple.time <- system.time(
data.nest.r2jags <- jags(data=data.nest.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/tut9.2bS3.1.txt",
n.chains=3,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 406

Initializing model

print(data.nest.r2jags)

Inference for Bugs model at "../downloads/BUGSscripts/tut9.2bS3.1.txt", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]     42.291   6.619  29.102  38.006  42.265  46.554  55.150 1.001  3000
beta[2]     29.542   9.263  11.148  23.640  29.710  35.687  47.273 1.001  3000
beta[3]     36.933   9.273  18.407  31.142  36.908  42.774  54.705 1.001  3000
sigma        4.428   0.278   3.925   4.233   4.409   4.614   5.013 1.001  3000
sigma.site  14.137   3.079   9.566  11.986  13.642  15.762  21.504 1.001  3000
deviance   869.595   6.056 859.795 865.211 868.824 873.347 882.948 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 = 18.3 and DIC = 887.9
DIC is an estimate of expected predictive error (lower deviance is better).

data.nest.mcmc.list <- as.mcmc(data.nest.r2jags)

modelString="
model {
#Likelihood (esimating site means (gamma.site)
for (i in 1:n) {
y[i]~dnorm(quad.means[i],tau)
quad.means[i] <- gamma.site[site[i]]
y.err[i]<- quad.means[i]-y[i]
}
for (i in 1:s) {
gamma.site[i] ~ dnorm(site.means[i], tau.site)
site.means[i] <- inprod(beta[],A.Xmat[i,])
site.err[i] <- site.means[i] - gamma.site[i]
}
#Priors
for (i in 1:a) {
beta[i] ~ dnorm(0, 1.0E-6) #prior
}
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.site <- pow(sigma.site,-2)
sigma.site <-z/sqrt(chSq.site)
z.site ~ dnorm(0, .0016)I(0,)
chSq.site ~ dgamma(0.5, 0.5)

sd.y <- sd(y.err)
sd.site <- sd(site.err)
sd.A <- sd(beta)
}
"
writeLines(modelString,con="../downloads/BUGSscripts/tut9.2bS3.1SD.txt")

A.Xmat <- model.matrix(~A,ddply(data.nest,~Sites,catcolwise(unique)))
data.nest.list <- with(data.nest,
list(y=y,
site=Sites,
A.Xmat= A.Xmat,
n=nrow(data.nest),
s=length(levels(Sites)),
a = ncol(A.Xmat)
)
)

params <- c("beta","sigma","sd.y",'sd.site','sd.A','sigma','sigma.site')
adaptSteps = 1000
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
rnorm(1)

[1] 1.075778

jags.SD.time <- system.time(
data.nest.r2jagsSD <- jags(data=data.nest.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/tut9.2bS3.1SD.txt",
n.chains=3,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 571

Initializing model

print(data.nest.r2jagsSD)

Inference for Bugs model at "../downloads/BUGSscripts/tut9.2bS3.1SD.txt", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]     42.356   6.247  30.043  38.368  42.410  46.124  55.080 1.001  3000
beta[2]     29.793   8.891  11.548  24.253  29.953  35.495  47.365 1.001  3000
beta[3]     37.054   8.860  19.147  31.318  37.084  42.866  54.862 1.002  3000
sd.A         9.506   5.273   1.628   5.674   8.766  12.360  22.348 1.002  1400
sd.site     13.662   1.167  12.182  12.905  13.383  14.090  16.702 1.002  1500
sd.y         4.378   0.084   4.245   4.317   4.368   4.426   4.573 1.001  2600
sigma        4.429   0.276   3.949   4.230   4.414   4.605   5.013 1.001  3000
sigma.site  13.989   3.014   9.501  11.876  13.487  15.624  21.078 1.002  1700
deviance   869.639   6.104 859.832 865.322 868.942 873.053 883.310 1.001  2800

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 = 18.6 and DIC = 888.3
DIC is an estimate of expected predictive error (lower deviance is better).

data.nest.mcmc.listSD <- as.mcmc(data.nest.r2jagsSD)


### Rstan

#### Cell means parameterization

rstanString="
data{
int n;
int nA;
int nB;
vector [n] y;
int A[n];
int B[n];
}

parameters{
real alpha[nA];
real<lower=0> sigma;
vector [nB] beta;
real<lower=0> sigma_B;
}

model{
real mu[n];

// Priors
alpha ~ normal( 0 , 100 );
beta ~ normal( 0 , sigma_B );
sigma_B ~ cauchy( 0 , 25 );
sigma ~ cauchy( 0 , 25 );

for ( i in 1:n ) {
mu[i] <- alpha[A[i]] + beta[B[i]];
}
y ~ normal( mu , sigma );
}
"

data.nest.list <- with(data.nest, list(y=y, A=as.numeric(A), B=as.numeric(Sites),
n=nrow(data.nest), nB=length(levels(Sites)),nA=length(levels(A))))
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
library(rstan)
rstan.c.time <- system.time(
data.nest.rstan.c <- stan(data=data.nest.list,
model_code=rstanString,
pars=c('alpha','sigma','sigma_B'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 1.41 seconds (Warm-up)
#                5.4 seconds (Sampling)
#                6.81 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 1.34 seconds (Warm-up)
#                4.83 seconds (Sampling)
#                6.17 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 1.3 seconds (Warm-up)
#                5.41 seconds (Sampling)
#                6.71 seconds (Total)

print(data.nest.rstan.c)

Inference for Stan model: rstanString.
3 chains, each with iter=13000; warmup=3000; thin=10;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha[1]   42.08    0.13 6.94   28.10   37.64   42.08   46.73   55.54  2721    1
alpha[2]   71.70    0.13 7.04   57.34   67.33   71.90   76.06   85.41  2929    1
alpha[3]   78.84    0.13 6.95   65.13   74.50   78.94   83.19   92.65  2829    1
sigma       4.41    0.01 0.27    3.92    4.23    4.40    4.59    4.95  2793    1
sigma_B    15.28    0.07 3.63   10.07   12.76   14.68   17.05   24.30  2715    1
lp__     -340.78    0.06 3.44 -348.51 -342.84 -340.33 -338.29 -335.37  2877    1

Samples were drawn using NUTS(diag_e) at Tue Jan 13 07:14:31 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

data.nest.rstan.c.df <-as.data.frame(extract(data.nest.rstan.c))
head(data.nest.rstan.c.df)

  alpha.1 alpha.2 alpha.3 sigma sigma_B   lp__
1   49.61   76.03   68.16 4.713   13.83 -343.9
2   38.30   78.58   78.08 4.527   11.79 -340.6
3   49.83   73.75   79.77 4.872   12.61 -336.2
4   41.54   50.96   92.55 3.949   20.76 -342.8
5   33.08   69.08   84.15 4.101   14.88 -336.3
6   39.11   78.25   78.48 4.427   10.71 -342.1

data.nest.mcmc.c<-rstan:::as.mcmc.list.stanfit(data.nest.rstan.c)

plyr:::adply(as.matrix(data.nest.rstan.c.df),2,MCMCsum)

       X1  Median      X0.     X25.    X50.     X75.    X100.    lower    upper lower.1  upper.1
1 alpha.1   42.08    9.533   37.639   42.08   46.729   75.866   28.980   56.283   35.90   48.855
2 alpha.2   71.90   44.736   67.325   71.90   76.061  103.790   57.346   85.439   65.43   78.292
3 alpha.3   78.94   48.500   74.501   78.94   83.189  108.715   64.980   92.419   72.22   85.211
4   sigma    4.40    3.621    4.228    4.40    4.585    5.835    3.904    4.927    4.10    4.626
5 sigma_B   14.68    7.428   12.757   14.68   17.055   34.375    9.190   22.720   11.18   17.284
6    lp__ -340.33 -355.469 -342.838 -340.33 -338.291 -332.218 -347.562 -334.755 -342.90 -336.451


#### Full effects parameterization

rstanString="
data{
int n;
int nB;
vector [n] y;
int A2[n];
int A3[n];
int B[n];
}

parameters{
real alpha0;
real alpha2;
real alpha3;
real<lower=0> sigma;
vector [nB] beta;
real<lower=0> sigma_B;
}

model{
real mu[n];

// Priors
alpha0 ~ normal( 0 , 100 );
alpha2 ~ normal( 0 , 100 );
alpha3 ~ normal( 0 , 100 );
beta ~ normal( 0 , sigma_B );
sigma_B ~ cauchy( 0 , 25 );
sigma ~ cauchy( 0 , 25 );

for ( i in 1:n ) {
mu[i] <- alpha0 + alpha2*A2[i] +
alpha3*A3[i] + beta[B[i]];
}
y ~ normal( mu , sigma );
}
"

A2 <- ifelse(data.nest$A=='a2',1,0) A3 <- ifelse(data.nest$A=='a3',1,0)
data.nest.list <- with(data.nest, list(y=y, A2=A2, A3=A3, B=as.numeric(Sites),
n=nrow(data.nest), nB=length(levels(Sites))))
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps + ceiling((numSavedSteps * thinSteps)/nChains)
library(rstan)
rstan.f.time <- system.time(
data.nest.rstan.f <- stan(data=data.nest.list,
model_code=rstanString,
pars=c('alpha0','alpha2','alpha3','sigma','sigma_B'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 3.28 seconds (Warm-up)
#                13.47 seconds (Sampling)
#                16.75 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 3.48 seconds (Warm-up)
#                14.72 seconds (Sampling)
#                18.2 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 3.23 seconds (Warm-up)
#                15.54 seconds (Sampling)
#                18.77 seconds (Total)

print(data.nest.rstan.f)

Inference for Stan model: rstanString.
3 chains, each with iter=13000; warmup=3000; thin=10;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha0    42.45    0.13 6.93   29.24   37.94   42.50   46.83   56.15  2899    1
alpha2    29.58    0.18 9.82   10.24   23.39   29.53   35.72   49.69  2897    1
alpha3    36.83    0.18 9.50   17.37   30.62   37.07   42.96   55.63  2925    1
sigma      4.40    0.00 0.27    3.91    4.22    4.40    4.57    4.95  3000    1
sigma_B   14.95    0.07 3.50    9.95   12.52   14.35   16.77   23.11  2704    1
lp__    -340.59    0.07 3.34 -347.82 -342.64 -340.25 -338.20 -335.05  2584    1

Samples were drawn using NUTS(diag_e) at Mon May  4 09:23:09 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

data.nest.rstan.f.df <-as.data.frame(extract(data.nest.rstan.f))
head(data.nest.rstan.f.df)

    alpha0   alpha2   alpha3    sigma   sigma_B      lp__
1 32.47415 39.49903 40.69339 4.435729 13.640733 -339.0126
2 25.72473 49.36534 48.50811 4.599241 25.359650 -346.7289
3 55.57984 25.72382 44.14519 4.099557 19.877447 -343.5320
4 37.53338 32.71673 32.61893 4.619255 13.468995 -341.5376
5 43.42154 30.63651 33.91730 4.343031  9.067916 -340.4652
6 31.66594 50.17470 44.25908 3.999964 13.984926 -344.7312

data.nest.mcmc.f<-rstan:::as.mcmc.list.stanfit(data.nest.rstan.f)

plyr:::adply(as.matrix(data.nest.rstan.f.df),2,MCMCsum)

       X1      Median         X0.        X25.        X50.        X75.       X100.       lower      upper     lower.1     upper.1
1  alpha0   42.501846    6.911381   37.940547   42.501846   46.825055   73.056778   29.233689   56.14772   36.086736   49.269184
2  alpha2   29.531136   -1.666353   23.392484   29.531136   35.715508   69.217756   10.215123   49.68790   20.143402   38.562762
3  alpha3   37.071971    1.197981   30.623994   37.071971   42.958094   80.889476   18.257879   56.04997   28.298563   46.016004
4   sigma    4.396666    3.431330    4.218224    4.396666    4.574869    5.656815    3.875447    4.90624    4.119459    4.651383
5 sigma_B   14.349013    6.790137   12.516834   14.349013   16.766142   42.674970    9.345433   21.81630   11.270147   17.121305
6    lp__ -340.250078 -358.072758 -342.640643 -340.250078 -338.196360 -331.776385 -347.300454 -334.91509 -342.944121 -336.533948


#### Matrix parameterization

rstanString="
data{
int n;
int nB;
int nA;
vector [n] y;
matrix [n,nA] X;
int B[n];
vector [nA] a0;
matrix [nA,nA] A0;
}

parameters{
vector [nA] alpha;
real<lower=0> sigma;
vector [nB] beta;
real<lower=0> sigma_B;
}

model{
real mu[n];

// Priors
//alpha ~ normal( 0 , 100 );
alpha ~ multi_normal(a0,A0);
beta ~ normal( 0 , sigma_B );
sigma_B ~ cauchy( 0 , 25);
sigma ~ cauchy( 0 , 25 );

for ( i in 1:n ) {
mu[i] <- dot_product(X[i],alpha) + beta[B[i]];
}
y ~ normal( mu , sigma );
}
"

X <- model.matrix(~A, data.nest)
nA <- ncol(X)
data.nest.list <- with(data.nest, list(y=y, X=X, B=as.numeric(Sites),
n=nrow(data.nest), nB=length(levels(Sites)), nA=nA,
a0=rep(0,nA), A0=diag(100000,nA)))
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
library(rstan)
rstan.m.time <- system.time(
data.nest.rstan.m <- stan(data=data.nest.list,
model_code=rstanString,
pars=c('alpha','sigma','sigma_B'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 4.55 seconds (Warm-up)
#                16.06 seconds (Sampling)
#                20.61 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 4.63 seconds (Warm-up)
#                24.03 seconds (Sampling)
#                28.66 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 4.26 seconds (Warm-up)
#                22.02 seconds (Sampling)
#                26.28 seconds (Total)

print(data.nest.rstan.m)

Inference for Stan model: rstanString.
3 chains, each with iter=13000; warmup=3000; thin=10;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
alpha[1]   42.25    0.13 6.87   28.24   37.98   42.33   46.62   55.54  2790    1
alpha[2]   29.88    0.18 9.80   10.29   23.69   29.84   36.21   49.29  3000    1
alpha[3]   37.02    0.18 9.61   18.11   30.73   37.00   43.07   56.38  3000    1
sigma       4.41    0.00 0.27    3.92    4.23    4.39    4.59    4.95  3000    1
sigma_B    14.88    0.06 3.27   10.00   12.61   14.33   16.61   22.68  2740    1
lp__     -340.48    0.06 3.37 -347.89 -342.52 -340.18 -338.08 -334.84  2952    1

Samples were drawn using NUTS(diag_e) at Thu Feb 19 08:32:30 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

data.nest.rstan.m.df <-as.data.frame(extract(data.nest.rstan.m))
head(data.nest.rstan.m.df)

  alpha.1 alpha.2 alpha.3 sigma sigma_B   lp__
1   37.90   35.27   38.67 4.306   12.75 -337.2
2   39.83   28.98   44.57 4.704   16.34 -337.6
3   47.66   18.27   41.29 4.575   10.54 -344.4
4   39.40   32.05   39.63 4.334   14.11 -337.9
5   44.49   34.05   37.88 4.588   15.83 -340.6
6   42.56   23.89   32.94 4.658   14.37 -339.4

data.nest.mcmc.m<-rstan:::as.mcmc.list.stanfit(data.nest.rstan.m)

plyr:::adply(as.matrix(data.nest.rstan.m.df),2,MCMCsum)

       X1   Median      X0.     X25.     X50.     X75.    X100.    lower    upper  lower.1  upper.1
1 alpha.1   42.331   15.703   37.976   42.331   46.616   70.261   28.782   55.863   36.425   49.449
2 alpha.2   29.843  -16.683   23.688   29.843   36.212   68.901   11.076   49.889   21.084   39.769
3 alpha.3   36.996   -5.127   30.726   36.996   43.066   80.704   17.129   55.169   26.566   44.635
4   sigma    4.394    3.617    4.226    4.394    4.589    5.455    3.890    4.916    4.146    4.672
5 sigma_B   14.327    7.738   12.613   14.327   16.608   38.417    9.569   21.283   11.197   16.928
6    lp__ -340.184 -354.122 -342.521 -340.184 -338.079 -332.052 -347.122 -334.425 -342.609 -336.099


#### Heirarchical parameterization

rstanString="
data{
int n;
int nA;
int nSites;
vector [n] y;
matrix [nSites,nA] X;
matrix [n,nSites] Z;
}

parameters{
vector[nA] beta;
vector[nSites] gamma;
real<lower=0> sigma;
real<lower=0> sigma_S;

}

model{
vector [n] mu_site;
vector [nSites] mu;

// Priors
beta ~ normal( 0 , 1000 );
gamma ~ normal( 0 , 1000 );
sigma ~ cauchy( 0 , 25 );
sigma_S~ cauchy( 0 , 25 );

mu_site <- Z*gamma;
y ~ normal( mu_site , sigma );
mu <- X*beta;
gamma ~ normal(mu, sigma_S);
}
"

dt.A <- ddply(data.nest,~Sites,catcolwise(unique))
X<-model.matrix(~A, dt.A)
Z<-model.matrix(~Sites-1, data.nest)
data.nest.list <- list(y=data.nest$y, X=X, Z=Z, n=nrow(data.nest), nSites=nrow(X),nA=ncol(X)) burnInSteps = 2000 nChains = 3 numSavedSteps = 3000 thinSteps = 10 nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains) library(rstan) rstan.time <- system.time( data.nest.rstan <- stan(data=data.nest.list, model_code=rstanString, pars=c('beta','sigma','sigma_S'), chains=nChains, iter=nIter, warmup=burnInSteps, thin=thinSteps, save_dso=TRUE ) )  TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW. SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1). Iteration: 1 / 12000 [ 0%] (Warmup) Iteration: 1200 / 12000 [ 10%] (Warmup) Iteration: 2001 / 12000 [ 16%] (Sampling) Iteration: 3200 / 12000 [ 26%] (Sampling) Iteration: 4400 / 12000 [ 36%] (Sampling) Iteration: 5600 / 12000 [ 46%] (Sampling) Iteration: 6800 / 12000 [ 56%] (Sampling) Iteration: 8000 / 12000 [ 66%] (Sampling) Iteration: 9200 / 12000 [ 76%] (Sampling) Iteration: 10400 / 12000 [ 86%] (Sampling) Iteration: 11600 / 12000 [ 96%] (Sampling) Iteration: 12000 / 12000 [100%] (Sampling) # Elapsed Time: 0.95 seconds (Warm-up) # 2.64 seconds (Sampling) # 3.59 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2). Iteration: 1 / 12000 [ 0%] (Warmup) Iteration: 1200 / 12000 [ 10%] (Warmup) Iteration: 2001 / 12000 [ 16%] (Sampling) Iteration: 3200 / 12000 [ 26%] (Sampling) Iteration: 4400 / 12000 [ 36%] (Sampling) Iteration: 5600 / 12000 [ 46%] (Sampling) Iteration: 6800 / 12000 [ 56%] (Sampling) Iteration: 8000 / 12000 [ 66%] (Sampling) Iteration: 9200 / 12000 [ 76%] (Sampling) Iteration: 10400 / 12000 [ 86%] (Sampling) Iteration: 11600 / 12000 [ 96%] (Sampling) Iteration: 12000 / 12000 [100%] (Sampling) # Elapsed Time: 0.84 seconds (Warm-up) # 3.51 seconds (Sampling) # 4.35 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3). Iteration: 1 / 12000 [ 0%] (Warmup) Iteration: 1200 / 12000 [ 10%] (Warmup) Iteration: 2001 / 12000 [ 16%] (Sampling) Iteration: 3200 / 12000 [ 26%] (Sampling) Iteration: 4400 / 12000 [ 36%] (Sampling) Iteration: 5600 / 12000 [ 46%] (Sampling) Iteration: 6800 / 12000 [ 56%] (Sampling) Iteration: 8000 / 12000 [ 66%] (Sampling) Iteration: 9200 / 12000 [ 76%] (Sampling) Iteration: 10400 / 12000 [ 86%] (Sampling) Iteration: 11600 / 12000 [ 96%] (Sampling) Iteration: 12000 / 12000 [100%] (Sampling) # Elapsed Time: 0.9 seconds (Warm-up) # 2.06 seconds (Sampling) # 2.96 seconds (Total)  print(data.nest.rstan)  Inference for Stan model: rstanString. 3 chains, each with iter=12000; warmup=2000; thin=10; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta[1] 42.25 0.13 6.89 28.35 37.88 42.33 46.61 56.50 3000 1 beta[2] 29.71 0.18 9.71 9.95 23.67 29.72 35.90 48.92 3000 1 beta[3] 36.93 0.18 9.79 17.08 30.50 36.93 43.37 56.33 3000 1 sigma 4.42 0.01 0.28 3.92 4.22 4.40 4.60 5.00 3000 1 sigma_S 14.78 0.06 3.34 9.90 12.32 14.21 16.68 22.71 2720 1 lp__ -340.53 0.06 3.40 -348.18 -342.64 -340.13 -338.06 -334.94 2912 1 Samples were drawn using NUTS(diag_e) at Thu Feb 19 08:35:35 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).  data.nest.rstan.df <-as.data.frame(extract(data.nest.rstan)) head(data.nest.rstan.df)   beta.1 beta.2 beta.3 sigma sigma_S lp__ 1 52.90 5.691 20.97 4.830 17.77 -343.7 2 40.98 18.733 49.22 4.913 18.30 -341.5 3 41.67 26.881 44.36 4.762 10.83 -341.9 4 43.24 28.760 43.06 4.190 14.72 -338.7 5 43.28 26.365 33.50 4.821 13.97 -343.7 6 49.65 27.695 28.97 4.356 13.20 -337.1  data.nest.mcmc<-rstan:::as.mcmc.list.stanfit(data.nest.rstan) plyr:::adply(as.matrix(data.nest.rstan.df),2,MCMCsum)   X1 Median X0. X25. X50. X75. X100. lower upper lower.1 upper.1 1 beta.1 42.328 14.535 37.88 42.328 46.611 72.173 26.915 54.617 36.422 49.209 2 beta.2 29.722 -4.106 23.67 29.722 35.901 77.658 9.729 48.599 20.820 39.048 3 beta.3 36.929 -4.409 30.50 36.929 43.372 80.227 15.942 55.042 27.826 46.240 4 sigma 4.399 3.504 4.22 4.399 4.601 5.586 3.924 4.995 4.127 4.666 5 sigma_S 14.214 6.880 12.32 14.214 16.675 36.266 9.141 21.190 11.013 16.944 6 lp__ -340.126 -362.898 -342.64 -340.126 -338.062 -331.824 -347.103 -334.244 -343.293 -336.816  rstanString=" data{ int n; int nA; int nSites; vector [n] y; matrix [nSites,nA] X; matrix [n,nSites] Z; } parameters{ vector[nA] beta; vector[nSites] gamma; real<lower=0> sigma; real<lower=0> sigma_S; } model{ vector [n] mu_site; vector [nSites] mu; // Priors beta ~ normal( 0 , 1000 ); gamma ~ normal( 0 , 1000 ); sigma ~ cauchy( 0 , 25 ); sigma_S~ cauchy( 0 , 25 ); mu_site <- Z*gamma; y ~ normal( mu_site , sigma ); mu <- X*beta; gamma ~ normal(mu, sigma_S); } generated quantities { vector [n] mu_site; vector [nSites] mu; vector [n] y_err; real sd_y; vector [nSites] mu_site_err; real sd_site; real sd_A; mu_site <- Z*gamma; y_err <- mu_site - y; sd_y <- sd(y_err); mu <- X*beta; mu_site_err <- mu - gamma; sd_site <- sd(mu_site_err); sd_A <- sd(beta); } "  Unlike JAGS, STAN's var() and sd() functions generate sample variance and standard deviation rather than population versions. Therefore, it is arguably better to calculate the population standard deviations within R from the MCMC samples. dt.A <- ddply(data.nest,~Sites,catcolwise(unique)) X<-model.matrix(~A, dt.A) Z<-model.matrix(~Sites-1, data.nest) data.nest.list <- list(y=data.nest$y, X=X, Z=Z, n=nrow(data.nest),
nSites=nrow(X),nA=ncol(X))
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)
library(rstan)
rstan.SD.time <- system.time(
data.nest.rstanSD <- stan(data=data.nest.list,
model_code=rstanString,
pars=c('beta','sigma','sigma_S','sd_A','sd_site','sd_y'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 0.93 seconds (Warm-up)
#                2.78 seconds (Sampling)
#                3.71 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 1.1 seconds (Warm-up)
#                2.41 seconds (Sampling)
#                3.51 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 1.1 seconds (Warm-up)
#                2.78 seconds (Sampling)
#                3.88 seconds (Total)

print(data.nest.rstanSD)

Inference for Stan model: rstanString.
3 chains, each with iter=13000; warmup=3000; thin=10;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
beta[1]   42.45    0.12 6.79   29.16   38.09   42.59   46.72   55.80  3000    1
beta[2]   29.81    0.18 9.63   10.51   23.40   30.04   35.83   49.00  2983    1
beta[3]   36.85    0.18 9.92   17.03   30.62   36.70   43.05   57.03  2956    1
sigma      4.42    0.00 0.27    3.93    4.24    4.41    4.60    4.99  3000    1
sigma_S   14.86    0.06 3.36    9.95   12.51   14.35   16.56   22.89  3000    1
sd_A      10.08    0.10 5.71    1.76    5.89    9.22   13.28   23.94  3000    1
sd_site   13.84    0.03 1.40   12.28   12.95   13.47   14.32   17.49  3000    1
sd_y       4.38    0.00 0.09    4.25    4.32    4.37    4.43    4.58  3000    1
lp__    -340.57    0.06 3.48 -348.45 -342.70 -340.19 -338.04 -334.92  2999    1

Samples were drawn using NUTS(diag_e) at Thu Feb 19 08:38:26 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

data.nest.rstanSD.df <-as.data.frame(extract(data.nest.rstanSD))
head(data.nest.rstanSD.df)

  beta.1 beta.2 beta.3 sigma sigma_S   sd_A sd_site  sd_y   lp__
1  37.64  25.04  40.55 4.460   16.63  8.243   13.51 4.312 -339.3
2  44.09  27.15  46.35 4.285   12.53 10.494   13.41 4.381 -339.2
3  33.73  41.89  40.89 3.930   12.94  4.448   13.10 4.355 -339.4
4  38.87  35.00  53.61 4.660   14.70  9.819   14.91 4.507 -344.8
5  51.41  15.90  31.59 4.721   16.01 17.793   13.59 4.519 -344.3
6  43.91  28.03  46.14 4.502   11.50  9.876   13.56 4.620 -349.1

data.nest.mcmcSD<-rstan:::as.mcmc.list.stanfit(data.nest.rstanSD)

plyr:::adply(as.matrix(data.nest.rstanSD.df),2,MCMCsum)

       X1   Median        X0.     X25.     X50.     X75.    X100.     lower    upper  lower.1  upper.1
1  beta.1   42.591   12.20757   38.093   42.591   46.722   67.341   28.7291   55.362   36.310   49.262
2  beta.2   30.043   -4.91213   23.402   30.043   35.833   64.149   11.0105   49.425   21.078   39.443
3  beta.3   36.697  -11.13186   30.624   36.697   43.054   97.367   16.7544   56.529   27.205   45.229
4   sigma    4.413    3.64025    4.241    4.413    4.595    5.439    3.9266    4.984    4.148    4.682
5 sigma_S   14.354    7.91922   12.510   14.354   16.561   37.159    9.0999   21.175   11.008   16.666
6    sd_A    9.219    0.09577    5.891    9.219   13.279   41.432    0.4164   20.659    3.394   13.617
7 sd_site   13.469   11.50089   12.947   13.469   14.320   28.962   12.0590   16.688   12.408   14.202
8    sd_y    4.366    4.18990    4.316    4.366    4.429    4.775    4.2227    4.540    4.276    4.433
9    lp__ -340.193 -356.63993 -342.696 -340.193 -338.037 -331.919 -347.2215 -334.098 -343.133 -336.460


#### $R^2$ and finite population standard deviations

modelString="
model {
#Likelihood (esimating site means (gamma.site)
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- gamma[site[i]] + inprod(beta[], X[i,])
y.err[i]<- mu[i]-y[i]
}
#Priors
for (i in 1:nSite) {
gamma[i] ~ dnorm(0, tau.site)
}
beta ~ dmnorm(a0,A0)
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.site <- pow(sigma.site,-2)
sigma.site <-z/sqrt(chSq.site)
z.site ~ dnorm(0, .0016)I(0,)
chSq.site ~ dgamma(0.5, 0.5)

sd.y <- sd(y.err)
sd.site <- sd(gamma)
}
"

A.Xmat <- model.matrix(~A,data.nest)
nX=ncol(A.Xmat)
data.nest.list <- with(data.nest,
list(y=y,
site=Sites,
X= A.Xmat,
n=nrow(data.nest),
nSite=length(levels(Sites)),
nX = ncol(A.Xmat),
a0=rep(0,nX), A0=diag(1.0E-06,nX)
)
)

params <- c("beta","sigma","sd.y",'sd.site','sigma','sigma.site')
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
jags.SD.time <- system.time(
data.nest.r2jagsSD <- jags(data=data.nest.list,
inits=NULL,
parameters.to.save=params,
model.file=textConnection(modelString),
n.chains=3,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 1119

Initializing model

jags.SD.time

   user  system elapsed
9.312   0.000   9.410

print(data.nest.r2jagsSD)

Inference for Bugs model at "5", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta[1]     42.493   6.446  29.226  38.399  42.479  46.558  55.119 1.001  3000
beta[2]     29.601   9.213  11.098  23.494  29.794  35.416  48.199 1.001  3000
beta[3]     36.647   8.929  19.397  30.733  36.608  42.360  53.886 1.008  3000
sd.site     13.691   1.216  12.243  12.905  13.378  14.139  16.829 1.001  3000
sd.y         4.378   0.084   4.249   4.315   4.368   4.427   4.580 1.001  3000
sigma        4.425   0.278   3.922   4.230   4.408   4.603   4.993 1.001  2800
sigma.site  14.029   2.986   9.526  11.930  13.601  15.633  21.261 1.001  3000
deviance   869.638   6.133 859.900 865.166 868.939 873.287 883.840 1.001  2100

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 = 18.8 and DIC = 888.4
DIC is an estimate of expected predictive error (lower deviance is better).

data.nest.mcmc.listSD <- as.mcmc(data.nest.r2jagsSD)

Xmat <- model.matrix(~A, data.nest)
coefs <- data.nest.r2jagsSD$BUGSoutput$sims.list[['beta']]
fitted <- coefs %*% t(Xmat)
X.var <- aaply(fitted,1,function(x){var(x)})
Z.var <- data.nest.r2jagsSD$BUGSoutput$sims.list[['sd.site']]^2
R.var <- data.nest.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.site <- (Z.var)/(X.var+Z.var+R.var)
R2.site <- data.frame(Mean=mean(R2.site), Median=median(R2.site), HPDinterval(as.mcmc(R2.site)))
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.site=R2.site, R2.marginal=R2.marginal, R2.res=R2.res, R2.conditional=R2.conditional)

                     Mean     Median      lower      upper
R2.site        0.40513614 0.37777481 0.26346438 0.62477617
R2.marginal    0.55316849 0.58096881 0.31506217 0.70640681
R2.res         0.04169536 0.04165901 0.02314323 0.06005702
R2.conditional 0.95830464 0.95834099 0.93994298 0.97685677


#### Graphical summary

newdata <- with(data.nest, data.frame(A=levels(A)))
Xmat <- model.matrix(~A, newdata)
coefs <- data.nest.r2jags.m$BUGSoutput$sims.list[['alpha']]
fit <- coefs %*% t(Xmat)
newdata <- cbind(newdata,
adply(fit, 2, function(x) {
data.frame(Mean=mean(x), Median=median(x), HPDinterval(as.mcmc(x)),
HPDinterval(as.mcmc(x), p=0.68))
})
)

library(ggplot2)
library(gridExtra)
library(grid)
p1 <- ggplot(newdata, aes(y=Median, x=A)) +
geom_errorbar(aes(ymin=lower, ymax=upper), width=0.01, size=1) +
geom_errorbar(aes(ymin=lower.1, ymax=upper.1), width=0, size=2) +
geom_point(size=4, shape=21, fill='white')+
scale_y_continuous('Y')+
scale_x_discrete('X')+
theme_classic()+
theme(axis.title.y=element_text(vjust=2, size=rel(1.25)),
axis.title.x=element_text(vjust=-2, size=rel(1.25)),
plot.margin=unit(c(0.5,0.5,2,2), 'lines')
)

data.nest.mcmc.listSD <- as.mcmc(data.nest.r2jagsSD)

Xmat <- model.matrix(~A, data.nest)
coefs <- data.nest.r2jagsSD$BUGSoutput$sims.list[['beta']]
fitted <- coefs %*% t(Xmat)
X.var <- aaply(fitted,1,function(x){var(x)})
Z.var <- data.nest.r2jagsSD$BUGSoutput$sims.list[['sd.site']]^2
R.var <- data.nest.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.site <- (Z.var)/(X.var+Z.var+R.var)
R2.site <- data.frame(Mean=mean(R2.site), Median=median(R2.site), HPDinterval(as.mcmc(R2.site)))
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)))

curdies.R2<-rbind(R2.site=R2.site, R2.marginal=R2.marginal, R2.res=R2.res, R2.conditional=R2.conditional)

curdies.R2$name <- factor(rownames(curdies.R2), levels=c('R2.conditional','R2.res','R2.site','R2.marginal'), labels=c('Total','Residuals','Sites','A')) #curdies.R2$Perc <- curdies.R2$median/sum(curdies.R2$median)
head(curdies.R2)

                     Mean     Median      lower      upper      name
R2.site        0.40513614 0.37777481 0.26346438 0.62477617     Sites
R2.marginal    0.55316849 0.58096881 0.31506217 0.70640681         A
R2.res         0.04169536 0.04165901 0.02314323 0.06005702 Residuals
R2.conditional 0.95830464 0.95834099 0.93994298 0.97685677     Total

p2<-ggplot(curdies.R2,aes(y=name, x=Median))+
geom_vline(xintercept=0,linetype="dashed")+
geom_hline(xintercept=0)+
#scale_x_continuous("Finite population \nvariance components (sd)")+
#geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
geom_point(size=3, shape=21, fill='white')+
# geom_text(aes(label=sprintf("(%4.1f%%)",Perc),vjust=-1))+
theme_classic()+
theme(axis.title.y=element_blank(),
axis.text.y=element_text(size=rel(1.2),hjust=1),
axis.title.x=element_text(vjust=-2, size=rel(1.25)),
plot.margin=unit(c(0.5,0.5,2,2), 'lines'))

grid.arrange(p1,p2,nrow=1)


#### BRMS (STAN)

library(brms)
data.nest.brm <- brm(y~A+(1|Sites), data=data.nest, family='gaussian',
prior=c(set_prior('normal(0,100)', class='b'),
set_prior('cauchy(0,5)', class='sd')),
n.chains=3, n.iter=2000, warmup=500, n.thin=2
)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 1, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 1, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 1, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 1, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 1, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 1, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 1, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 1.72613 seconds (Warm-up)
#                2.79666 seconds (Sampling)
#                4.52279 seconds (Total)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 2, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 2, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 2, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 2, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 2, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 2, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 2, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 2, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 1.53597 seconds (Warm-up)
#                2.79055 seconds (Sampling)
#                4.32652 seconds (Total)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 3, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 3, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 3, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 3, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 3, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 3, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 3, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 3, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 1.22493 seconds (Warm-up)
#                2.90738 seconds (Sampling)
#                4.13231 seconds (Total)

summary(data.nest.brm)

 Family: gaussian (identity)
Formula: y ~ A + (1 | Sites)
Data: data.nest (Number of observations: 150)
Samples: 3 chains, each with n.iter = 2000; n.warmup = 500; n.thin = 2;
total post-warmup samples = 2250
WAIC: 885.4

Random Effects:
~Sites (Number of levels: 15)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)    13.95      2.77     9.58    20.23        840 1.01

Fixed Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    41.92      6.30    29.52    54.55        759    1
Aa2          30.28      8.85    12.40    47.54        780    1
Aa3          37.38      9.14    19.16    55.96        731    1

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma(y)      4.4      0.27      3.9     4.95       1178    1

Samples were drawn using NUTS(diag_e). For each parameter, Eff.Sample is a
crude measure of effective sample size, and Rhat is the potential scale
reduction factor on split chains (at convergence, Rhat = 1).

brms:::stancode(data.nest.brm)

functions {
}
data {
int<lower=1> N;  # number of observations
vector[N] Y;  # response variable
int<lower=1> K;  # number of fixed effects
matrix[N, K] X;  # FE design matrix
# data for random effects of Sites
int<lower=1> J_1[N];  # RE levels
int<lower=1> N_1;  # number of levels
int<lower=1> K_1;  # number of REs
real Z_1[N];  # RE design matrix
}
transformed data {
}
parameters {
real b_Intercept;  # fixed effects Intercept
vector[K] b;  # fixed effects
vector[N_1] pre_1;  # unscaled REs
real<lower=0> sd_1;  # RE standard deviation
real<lower=0> sigma;  # residual SD
}
transformed parameters {
vector[N] eta;  # linear predictor
vector[N_1] r_1;  # REs
# compute linear predictor
eta <- X * b + b_Intercept;
r_1 <- sd_1 * (pre_1);  # scale REs
# if available add REs to linear predictor
for (n in 1:N) {
eta[n] <- eta[n] + Z_1[n] * r_1[J_1[n]];
}
}
model {
# prior specifications
b_Intercept ~ normal(0,100);
b ~ normal(0,100);
sd_1 ~ cauchy(0,5);
pre_1 ~ normal(0, 1);
sigma ~ cauchy(0, 21);
# likelihood contribution
Y ~ normal(eta, sigma);
}
generated quantities {
}

#brms:::standata(data.nest.brm)


### Quick glance Bayesian regression model specification

JAGS Full JAGS Matrix JAGS heirarchy Stan Cellmeans Stan Full Stan Matrix Stan hierarcical
Iterations 13000 13000 13000 13000 13000 13000 12000
Thinning 10 10 10 10 10 10 10
Samples 3000 3000 3000 3000 3000 3000 3000
Maximum ACF 0.0189 0.0169 0.0202 0.0531 0.0524 0.0462 0.0301
Intercept 42.14 [29.91 , 55.39] 42.35 [30.22 , 54.94] 42.27 [29.01 , 54.92] 42.08 [28.98 , 56.28] 42.50 [29.23 , 56.15] 42.33 [28.78 , 55.86] 42.33 [26.92 , 54.62]
alpha2 30.05 [12.79 , 48.63] 29.74 [13.51 , 48.63] 29.71 [10.96 , 46.82] 71.90 [57.35 , 85.44] 29.53 [10.22 , 49.69] 29.84 [11.08 , 49.89] 29.72 [9.73 , 48.60]
alpha3 37.12 [19.56 , 55.19] 36.77 [18.30 , 54.57] 36.91 [19.19 , 55.24] 78.94 [64.98 , 92.42] 37.07 [18.26 , 56.05] 37.00 [17.13 , 55.17] 36.93 [15.94 , 55.04]
sigma 4.41 [3.93 , 4.99] 4.42 [3.91 , 4.97] 4.41 [3.90 , 4.98] 4.40 [3.90 , 4.93] 4.40 [3.88 , 4.91] 4.39 [3.89 , 4.92] 4.40 [3.92 , 5.00]
sigma.site 13.57 [8.75 , 20.17] 13.64 [8.87 , 19.75] 13.64 [8.92 , 20.20] 14.68 [9.19 , 22.72] 14.35 [9.35 , 21.82] 14.33 [9.57 , 21.28] 14.21 [9.14 , 21.19]
Time (user self) 7.81 9.58 6.92 19.92 53.48 75.70 11.11
Time (system self) 0.01 0.01 0.01 0.01 0.57 0.04 0.02
Time (elapsed) 7.87 9.67 6.94 53.43 86.65 105.92 38.84
Time (user child) 0.00 0.00 0.00 33.03 29.31 28.77 26.60
Time (system child) 0.00 0.00 0.00 0.82 0.99 0.66 0.76
ModelJAGSrstan
Full effects parameterization
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- alpha0 + alpha[A[i]] + beta.site[site[i]] + beta.quad[quad[i]]
}

#Priors
alpha0 ~ dnorm(0, 1.0E-6)
alpha[1] <- 0
for (i in 2:nA) {
alpha[i] ~ dnorm(0, 1.0E-6) #prior
}
for (i in 1:nSite) {
beta.site[i] ~ dnorm(0, tau.Bs) #prior
}
for (i in 1:nQuad) {
beta.quad[i] ~ dnorm(0, tau.Bq) #prior
}
tau <- pow(sigma,-2)
sigma <-z/sqrt(chSq)

data{
int n;
int nB;
vector [n] y;
int A2[n];
int A3[n];
int B[n];
}

parameters{
real alpha0;
real alpha2;
real alpha3;
real<lower=0> sigma;
vector [nB] beta;
real<lower=0> sigma_B;
}

model{
real mu[n];

// Priors
alpha0 ~ normal( 0 , 100 );
alpha2 ~ normal( 0 , 100 );
alpha3 ~ normal( 0 , 100 );
beta ~ normal( 0 , sigma_B );
sigma_B ~ cauchy( 0 , 25 );
sigma ~ cauchy( 0 , 25 );

for ( i in 1:n ) {
mu[i] <- alpha0 + alpha2*A2[i] +
alpha3*A3[i] + beta[B[i]];
}
y ~ normal( mu , sigma );

Matrix parameterization
data{
int n;
int nB;
int nA;
vector [n] y;
matrix [n,nA] X;
int B[n];
vector [nA] a0;
matrix [nA,nA] A0;
}

parameters{
vector [nA] alpha;
real<lower=0> sigma;
vector [nB] beta;
real<lower=0> sigma_B;
}

model{
real mu[n];

// Priors
//alpha ~ normal( 0 , 100 );
alpha ~ multi_normal(a0,A0);
beta ~ normal( 0 , sigma_B );
sigma_B ~ cauchy( 0 , 25);
sigma ~ cauchy( 0 , 25 );

for ( i in 1:n ) {
mu[i] <- dot_product(X[i],alpha) + beta[B[i]];
}

Hierarchical parameterization
data{
int n;
int nA;
int nSites;
vector [n] y;
matrix [nSites,nA] X;
matrix [n,nSites] Z;
}

parameters{
vector[nA] beta;
vector[nSites] gamma;
real<lower=0> sigma;
real<lower=0> sigma_S;

}

model{
vector [n] mu_site;
vector [nSites] mu;

// Priors
beta ~ normal( 0 , 1000 );
gamma ~ normal( 0 , 1000 );
sigma ~ cauchy( 0 , 25 );
sigma_S~ cauchy( 0 , 25 );

mu_site <- Z*gamma;
y ~ normal( mu_site , sigma );
mu <- X*beta;
gamma ~ normal(mu, sigma_S);
}

Including finite-population standard deviations
data{
int n;
int nA;
int nSites;
vector [n] y;
matrix [nSites,nA] X;
matrix [n,nSites] Z;
}

parameters{
vector[nA] beta;
vector[nSites] gamma;
real<lower=0> sigma;
real<lower=0> sigma_S;

}

model{
vector [n] mu_site;
vector [nSites] mu;

// Priors
beta ~ normal( 0 , 1000 );
gamma ~ normal( 0 , 1000 );
sigma ~ cauchy( 0 , 25 );
sigma_S~ cauchy( 0 , 25 );

mu_site <- Z*gamma;
y ~ normal( mu_site , sigma );
mu <- X*beta;
gamma ~ normal(mu, sigma_S);
}

generated quantities {
vector [n] mu_site;
vector [nSites] mu;
vector [n] y_err;
real sd_y;
vector [nSites] mu_site_err;
real sd_site;
real sd_A;

mu_site <- Z*gamma;
y_err <- mu_site - y;
sd_y <- sd(y_err);

mu <- X*beta;
mu_site_err <- mu - gamma;
sd_site <- sd(mu_site_err);

sd_A <- sd(beta);
}


### Scenario and Data

Now imagine a similar experiment in which we intend to measure a response ($y$) to one of treatments (three levels; 'a1', 'a2' and 'a3'). As with the previous design, we decided to establish a nested design in which there are sub-replicate (1m Quadrats) within each Site.

In the current design, we have decided to further sub-replicate. Within each of the 5 Quadrats, we are going to randomly place 2x10cm pit traps. Now we have Sites nested within Treatments, Quadrats nested within Sites AND, Pits nested within Sites. The latter of these (Pits nested within Sites) are the observations ($y$). 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. Consequently, I have folded (toggled) this section away.

Random data incorporating the following properties
• the number of treatments = 3
• the number of sites per treatment = 5
• the total number of sites = 15
• the number of quadrats per site = 5
• the number of pits per quadrat = 2
• the mean of the treatments = 40, 70 and 80 respectively
• the variability (standard deviation) between sites of the same treatment = 15
• the variability (standard deviation) between quadrats within sites = 10
• the variability (standard deviation) between pits within quadrats = 5
set.seed(3)
nTreat <- 3
nSites <- 15
nSitesPerTreat <- nSites/nTreat
nQuads <- 5
nPits <- 2
site.sigma <- 10 # sd within between sites within treatment
quad.sigma <- 10
sigma <- 7.5
n <- nSites * nQuads * nPits
sites <- gl(n=nSites,n/nSites,n, lab=paste("site",1:nSites))
A <- gl(nTreat, n/nTreat, n, labels=c('a1','a2','a3'))
a.means <- c(40,70,80)

#A<-gl(nTreat,nSites/nTreat,nSites,labels=c('a1','a2','a3'))
a<-gl(nTreat,1,nTreat,labels=c('a1','a2','a3'))
a.X <- model.matrix(~a, expand.grid(a))
a.eff <- as.vector(solve(a.X,a.means))
site.means <- rnorm(nSites,a.X %*% a.eff,site.sigma)

A <- gl(nTreat,nSites/nTreat,nSites,labels=c('a1','a2','a3'))
A.X <- model.matrix(~A, expand.grid(A))
#a.X <- model.matrix(~A, expand.grid(A=gl(nTreat,nSites/nTreat,nSites,labels=c('a1','a2','a3'))))
site.means <- rnorm(nSites,A.X %*% a.eff,site.sigma)

SITES <- gl(nSites,(nSites*nQuads)/nSites,nSites*nQuads,labels=paste('site',1:nSites))
sites.X <- model.matrix(~SITES-1)
quad.means <- rnorm(nSites*nQuads,sites.X %*% site.means,quad.sigma)

#SITES <- gl(nSites,1,nSites,labels=paste('site',1:nSites))
#sites.X <- model.matrix(~SITES-1)
#quad.means <- rnorm(nSites*nQuads,sites.X %*% site.means,quad.sigma)

QUADS <- gl(nSites*nQuads,n/(nSites*nQuads),n,labels=paste('quad',1:(nSites*nQuads)))
quads.X <- model.matrix(~QUADS-1)
#quads.eff <- as.vector(solve(quads.X,quad.means))
#pit.means <- rnorm(n,quads.eff %*% t(quads.X),sigma)
pit.means <- rnorm(n,quads.X %*% quad.means,sigma)

PITS <- gl(nPits*nSites*nQuads,1, n, labels=paste('pit',1:(nPits*nSites*nQuads)))
data.nest1<-data.frame(Pits=PITS,Quads=QUADS,Sites=rep(SITES,each=2), A=rep(A,each=nQuads*nPits),y=pit.means)
A3 <- ifelse(data.nest1$A=='a3',1,0) data.nest.list <- with(data.nest1, list(y=y, A2=A2, A3=A3, Site=as.numeric(Sites), n=nrow(data.nest1), nSite=length(levels(Sites)), nQuad=length(levels(Quads)), Quad=as.numeric(Quads))) burnInSteps = 3000 nChains = 3 numSavedSteps = 3000 thinSteps = 10 nIter = burnInSteps + ceiling((numSavedSteps * thinSteps)/nChains) library(rstan) rstan.f.time <- system.time( data.nest1.rstan.f <- stan(data=data.nest.list, model_code=rstanString, pars=c('alpha0','alpha2','alpha3','sigma','sigma_Site', 'sigma_Quad'), chains=nChains, iter=nIter, warmup=burnInSteps, thin=thinSteps, save_dso=TRUE ) )  TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW. SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 1.76 seconds (Warm-up) # 5.5 seconds (Sampling) # 7.26 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 1.69 seconds (Warm-up) # 3.8 seconds (Sampling) # 5.49 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 1.83 seconds (Warm-up) # 5.73 seconds (Sampling) # 7.56 seconds (Total)  print(data.nest1.rstan.f)  Inference for Stan model: rstanString. 3 chains, each with iter=13000; warmup=3000; thin=10; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha0 38.52 0.04 2.28 33.99 37.00 38.57 39.91 43.09 2795 1.00 alpha2 26.38 0.06 3.25 19.80 24.36 26.48 28.39 32.63 2783 1.00 alpha3 45.40 0.06 3.22 38.94 43.44 45.44 47.41 51.60 2829 1.00 sigma 7.84 0.02 0.66 6.70 7.37 7.82 8.25 9.24 1609 1.00 sigma_Site 2.39 0.08 1.61 0.20 1.12 2.14 3.32 6.13 447 1.00 sigma_Quad 7.64 0.02 1.12 5.45 6.85 7.64 8.38 9.85 2098 1.00 lp__ -583.54 1.13 13.65 -607.82 -593.03 -584.72 -575.55 -550.38 145 1.01 Samples were drawn using NUTS(diag_e) at Sun Mar 8 15:12:28 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).  data.nest1.rstan.f.df <-as.data.frame(extract(data.nest1.rstan.f)) head(data.nest1.rstan.f.df)   alpha0 alpha2 alpha3 sigma sigma_Site sigma_Quad lp__ 1 37.75 27.12 46.91 7.732 0.7185 7.604 -563.0 2 38.05 25.86 48.26 7.777 2.4772 7.498 -585.3 3 36.21 29.77 45.46 8.557 2.1948 7.077 -590.4 4 37.73 33.21 49.94 8.704 6.7076 5.562 -601.9 5 36.58 29.90 47.96 7.692 0.9029 10.274 -582.6 6 38.30 29.33 47.19 8.029 3.3813 6.247 -593.5  data.nest1.mcmc.f<-rstan:::as.mcmc.list.stanfit(data.nest1.rstan.f) plyr:::adply(as.matrix(data.nest1.rstan.f.df),2,MCMCsum)   X1 Median X0. X25. X50. X75. X100. lower upper lower.1 upper.1 1 alpha0 38.575 28.2309 36.998 38.575 39.907 48.41 34.2376 43.251 36.4326 40.775 2 alpha2 26.482 13.6165 24.359 26.482 28.386 41.29 20.3399 33.092 23.0974 29.233 3 alpha3 45.438 31.8434 43.437 45.438 47.410 57.45 39.3895 51.914 42.5987 48.755 4 sigma 7.819 6.0881 7.373 7.819 8.247 10.59 6.6687 9.198 7.1035 8.334 5 sigma_Site 2.143 0.1809 1.118 2.143 3.325 12.33 0.1809 5.371 0.1809 2.963 6 sigma_Quad 7.639 3.3836 6.846 7.639 8.382 12.26 5.4516 9.851 6.6402 8.768 7 lp__ -584.722 -623.1653 -593.031 -584.722 -575.551 -540.56 -604.7332 -549.884 -598.8459 -573.657  ##### Matrix parameterization rstanString=" data{ int n; int nSite; int nQuad; int nA; vector [n] y; matrix [n,nA] X; int Site[n]; int Quad[n]; vector [nA] a0; matrix [nA,nA] A0; } parameters{ vector [nA] alpha; real<lower=0> sigma; vector [nSite] beta_Site; real<lower=0> sigma_Site; vector [nQuad] beta_Quad; real<lower=0> sigma_Quad; } model{ real mu[n]; // Priors //alpha ~ normal( 0 , 100 ); alpha ~ multi_normal(a0,A0); beta_Site ~ normal( 0 , sigma_Site ); sigma_Site ~ cauchy( 0 , 25); beta_Quad ~ normal( 0 , sigma_Quad ); sigma_Quad ~ cauchy( 0 , 25); sigma ~ cauchy( 0 , 25 ); for ( i in 1:n ) { mu[i] <- dot_product(X[i],alpha) + beta_Site[Site[i]] + beta_Quad[Quad[i]]; } y ~ normal( mu , sigma ); } "  X <- model.matrix(~A, data.nest) nA <- ncol(X) data.nest.list <- with(data.nest1, list(y=y, X=X, Site=as.numeric(Sites), Quad=as.numeric(Quads), n=nrow(data.nest1), nSite=length(levels(Sites)), nQuad=length(levels(Quads)), nA=nA, a0=rep(0,nA), A0=diag(100000,nA))) burnInSteps = 3000 nChains = 3 numSavedSteps = 3000 thinSteps = 10 nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains) library(rstan) rstan.m.time <- system.time( data.nest1.rstan.m <- stan(data=data.nest.list, model_code=rstanString, pars=c('alpha','sigma','sigma_Site', 'sigma_Quad'), chains=nChains, iter=nIter, warmup=burnInSteps, thin=thinSteps, save_dso=TRUE ) )  TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW. SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 3.19 seconds (Warm-up) # 9.42 seconds (Sampling) # 12.61 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 3.59 seconds (Warm-up) # 10.41 seconds (Sampling) # 14 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 3.97 seconds (Warm-up) # 28.8 seconds (Sampling) # 32.77 seconds (Total)  print(data.nest1.rstan.m)  Inference for Stan model: rstanString. 3 chains, each with iter=13000; warmup=3000; thin=10; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha[1] 38.53 0.04 2.30 33.94 37.00 38.52 40.01 42.98 2643 1 alpha[2] 26.27 0.06 3.26 19.84 24.10 26.23 28.33 32.59 3000 1 alpha[3] 45.41 0.06 3.25 38.98 43.26 45.49 47.46 51.61 3000 1 sigma 7.87 0.01 0.67 6.70 7.41 7.82 8.27 9.35 2749 1 sigma_Site 2.48 0.05 1.58 0.30 1.25 2.26 3.42 6.15 1177 1 sigma_Quad 7.64 0.02 1.14 5.51 6.89 7.63 8.40 9.86 2340 1 lp__ -584.35 0.41 13.01 -607.32 -593.42 -585.09 -576.81 -555.47 1005 1 Samples were drawn using NUTS(diag_e) at Sun Mar 8 15:20:14 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).  data.nest1.rstan.m.df <-as.data.frame(extract(data.nest1.rstan.m)) head(data.nest1.rstan.m.df)   alpha.1 alpha.2 alpha.3 sigma sigma_Site sigma_Quad lp__ 1 32.75 31.91 50.51 7.634 3.0127 7.464 -594.2 2 39.40 23.76 42.90 6.620 2.3567 8.884 -572.7 3 38.35 24.41 43.42 7.147 0.7239 9.753 -570.9 4 39.64 26.10 44.11 7.187 2.1584 8.311 -572.4 5 38.01 27.74 42.38 6.840 2.1693 9.145 -583.1 6 37.63 26.06 46.30 7.771 0.6094 8.379 -571.7  data.nest1.mcmc.m<-rstan:::as.mcmc.list.stanfit(data.nest1.rstan.m) plyr:::adply(as.matrix(data.nest1.rstan.m.df),2,MCMCsum)   X1 Median X0. X25. X50. X75. X100. lower upper lower.1 upper.1 1 alpha.1 38.525 30.21146 37.001 38.525 40.005 48.21 34.14305 43.134 36.2077 40.704 2 alpha.2 26.226 14.60060 24.104 26.226 28.333 39.52 20.32018 33.020 23.0563 29.346 3 alpha.3 45.488 32.75320 43.264 45.488 47.456 61.27 39.20062 51.789 42.4782 48.693 4 sigma 7.824 5.87588 7.407 7.824 8.272 10.88 6.58580 9.169 7.0994 8.378 5 sigma_Site 2.258 0.07567 1.246 2.258 3.425 13.30 0.07567 5.307 0.2434 3.118 6 sigma_Quad 7.626 1.84726 6.888 7.626 8.404 12.30 5.48287 9.823 6.6016 8.799 7 lp__ -585.088 -624.80265 -593.424 -585.088 -576.809 -527.77 -609.25382 -557.749 -599.5603 -575.088  ####$R^2$and finite population standard deviations If we use the JAGS matrix parameterization model from above, the JAGS model is already complete (as we defined the sd components in that model already). data.nest1.mcmc.listSD <- as.mcmc(data.nest1.r2jags.m) Xmat <- model.matrix(~A, data.nest1) coefs <- data.nest.r2jags.m$BUGSoutput$sims.list[['alpha']] fitted <- coefs %*% t(Xmat) X.var <- aaply(fitted,1,function(x){var(x)}) Z.var <- data.nest.r2jagsSD$BUGSoutput$sims.list[['sd.site']]^2 R.var <- data.nest.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.site <- (Z.var)/(X.var+Z.var+R.var) R2.site <- data.frame(Mean=mean(R2.site), Median=median(R2.site), HPDinterval(as.mcmc(R2.site))) 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.site=R2.site, R2.marginal=R2.marginal, R2.res=R2.res, R2.conditional=R2.conditional)   Mean Median lower upper R2.site 0.40326358 0.39162062 0.21710229 0.62151297 R2.marginal 0.55526156 0.56689160 0.32820803 0.76703797 R2.res 0.04147485 0.04007836 0.02189078 0.06277313 R2.conditional 0.95852515 0.95992164 0.93722687 0.97810922  #### BRMS (STAN) library(brms) data.nest1$fQuads <- interaction(data.nest1$Sites, data.nest1$Quads)
data.nest1.brm <- brm(y~A+(1|Sites) + (1|Quads), data=data.nest1, family='gaussian',
prior=c(set_prior('normal(0,100)', class='b'),
set_prior('cauchy(0,5)', class='sd')),
n.chains=3, n.iter=2000, warmup=500, n.thin=2
)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 1, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 1, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 1, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 1, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 1, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 1, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 1, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 1.08899 seconds (Warm-up)
#                1.704 seconds (Sampling)
#                2.793 seconds (Total)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 2, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 2, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 2, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 2, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 2, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 2, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 2, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 2, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.634045 seconds (Warm-up)
#                1.19444 seconds (Sampling)
#                1.82849 seconds (Total)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 3, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 3, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 3, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 3, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 3, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 3, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 3, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 3, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.6651 seconds (Warm-up)
#                0.884811 seconds (Sampling)
#                1.54991 seconds (Total)

summary(data.nest1.brm)

 Family: gaussian (identity)
Formula: y ~ A + (1 | Sites) + (1 | Quads)
Data: data.nest1 (Number of observations: 150)
Samples: 3 chains, each with n.iter = 2000; n.warmup = 500; n.thin = 2;
total post-warmup samples = 2250
WAIC: 1096.55

Random Effects:
~Quads (Number of levels: 75)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)      7.5      1.09     5.31     9.69        973    1

~Sites (Number of levels: 15)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     2.02       1.4     0.09     5.21        675    1

Fixed Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    38.49      2.21    34.17    42.76       1143    1
Aa2          26.19      3.09    20.28    32.27       1207    1
Aa3          45.47      3.10    39.51    51.49       1395    1

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma(y)     7.87      0.68     6.67     9.32        778    1

Samples were drawn using NUTS(diag_e). For each parameter, Eff.Sample is a
crude measure of effective sample size, and Rhat is the potential scale
reduction factor on split chains (at convergence, Rhat = 1).

brms:::stancode(data.nest1.brm)

functions {
}
data {
int<lower=1> N;  # number of observations
vector[N] Y;  # response variable
int<lower=1> K;  # number of fixed effects
matrix[N, K] X;  # FE design matrix
# data for random effects of Quads
int<lower=1> J_1[N];  # RE levels
int<lower=1> N_1;  # number of levels
int<lower=1> K_1;  # number of REs
real Z_1[N];  # RE design matrix
# data for random effects of Sites
int<lower=1> J_2[N];  # RE levels
int<lower=1> N_2;  # number of levels
int<lower=1> K_2;  # number of REs
real Z_2[N];  # RE design matrix
}
transformed data {
}
parameters {
real b_Intercept;  # fixed effects Intercept
vector[K] b;  # fixed effects
vector[N_1] pre_1;  # unscaled REs
real<lower=0> sd_1;  # RE standard deviation
vector[N_2] pre_2;  # unscaled REs
real<lower=0> sd_2;  # RE standard deviation
real<lower=0> sigma;  # residual SD
}
transformed parameters {
vector[N] eta;  # linear predictor
vector[N_1] r_1;  # REs
vector[N_2] r_2;  # REs
# compute linear predictor
eta <- X * b + b_Intercept;
r_1 <- sd_1 * (pre_1);  # scale REs
r_2 <- sd_2 * (pre_2);  # scale REs
# if available add REs to linear predictor
for (n in 1:N) {
eta[n] <- eta[n] + Z_1[n] * r_1[J_1[n]] + Z_2[n] * r_2[J_2[n]];
}
}
model {
# prior specifications
b_Intercept ~ normal(0,100);
b ~ normal(0,100);
sd_1 ~ cauchy(0,5);
pre_1 ~ normal(0, 1);
sd_2 ~ cauchy(0,5);
pre_2 ~ normal(0, 1);
sigma ~ cauchy(0, 22);
# likelihood contribution
Y ~ normal(eta, sigma);
}
generated quantities {
}

#brms:::standata(data.nest1.brm)


## Worked Examples

##### Nested ANOVA (Mixed effects) references
• McCarthy (2007) - Chpt ?
• Kery (2010) - Chpt ?
• Gelman & Hill (2007) - Chpt ?
• Logan (2010) - Chpt 12
• Quinn & Keough (2002) - Chpt 9

### Nested ANOVA - one between factor

In an unusually detailed preparation for an Environmental Effects Statement for a proposed discharge of dairy wastes into the Curdies River, in western Victoria, a team of stream ecologists wanted to describe the basic patterns of variation in a stream invertebrate thought to be sensitive to nutrient enrichment. As an indicator species, they focused on a small flatworm, Dugesia, and started by sampling populations of this worm at a range of scales. They sampled in two seasons, representing different flow regimes of the river - Winter and Summer. Within each season, they sampled three randomly chosen (well, haphazardly, because sites are nearly always chosen to be close to road access) sites. A total of six sites in all were visited, 3 in each season. At each site, they sampled six stones, and counted the number of flatworms on each stone.

Download Curdies data set
Format of curdies.csv data files
SEASONSITEDUGESIAS4DUGESIA
WINTER10.6480.897
........
WINTER21.0161.004
........
WINTER30.6890.991
........
SUMMER400
........

Each row represents a different stone
 SEASON Season in which flatworms were counted - fixed factor SITE Site from where flatworms were counted - nested within SEASON (random factor) DUGESIA Number of flatworms counted on a particular stone S4DUGESIA 4th root transformation of DUGESIA variable

The hierarchical nature of this design can be appreciated when we examine an illustration of the spatial and temporal scale of the measurements.

• Within each of the two Seasons, there were three separate Sites (Sites not repeatidly measured across Seasons).
• Within each of the Sites, six logs were selected (haphazardly) from which the number of flatworms were counted.
So the Logs (smallest sampling units) are the replicates for the Sites (six reps per Site) and the Site means are the replicates for the two Seasons (three replicates within each Season).

We construct the model in much the same way. That is:

• The observed number of flatworm per log (for a given Site) is equal to (modelled against) the mean number for that Site plus a random amount drawn from a normal distribution with a mean of zero and a certain degree of variability (precision).
Number of flatwormlog i = βSite j,i + εi where ε ∼ N(0,σ2)
• The mean number of flatworm per Site (within a given Season) is equal to (modelled against) the mean number for that Season plus a random amount drawn from a normal distribution with a means of zero and a certain degree of variability (precision).
μSite j = γSeason j + εj where ε ∼ N(0,σ2)

The SITE variable is conveniently already coded as a numeric series (starting at 1).

Open
the curdies data file.
Show code
curdies <- read.table('../downloads/data/curdies.csv', header=T, sep=',', strip.white=T)
head(curdies)

  SEASON SITE   DUGESIA   S4DUGES
1 WINTER    1 0.6476829 0.8970995
2 WINTER    1 6.0961516 1.5713175
3 WINTER    1 1.3105639 1.0699526
4 WINTER    1 1.7252788 1.1460797
5 WINTER    1 1.4593867 1.0991136
6 WINTER    1 1.0575610 1.0140897

1. Lets start by performing some exploratory data analysis to help us identify any issues with the data as well as helping to try to identify a possible model to fit.
• Boxplots for the main effect (Season)
Show code
library(dplyr)
curdies.ag <- curdies %>% group_by(SEASON,SITE) %>% summarize(DUGESIA=mean(DUGESIA))
curdies.ag

Source: local data frame [6 x 3]
Groups: SEASON

SEASON SITE   DUGESIA
1 SUMMER    4 0.4190947
2 SUMMER    5 0.2290862
3 SUMMER    6 0.1942443
4 WINTER    1 2.0494375
5 WINTER    2 4.1819078
6 WINTER    3 0.6782063

#OR equivalently with the plyr package
library(plyr)
ddply(curdies, ~SEASON+SITE, numcolwise(mean))

  SEASON SITE   DUGESIA   S4DUGES
1 SUMMER    4 0.4190947 0.3508213
2 SUMMER    5 0.2290862 0.1804622
3 SUMMER    6 0.1942443 0.3811223
4 WINTER    1 2.0494375 1.1329421
5 WINTER    2 4.1819078 1.2718698
6 WINTER    3 0.6782063 0.8678707

detach(package:plyr) #since dplyr and plyr do not play nicely together

library(ggplot2)
ggplot(curdies.ag, aes(y=DUGESIA, x=SEASON)) + geom_boxplot()

Is there any evidence of non-normality and/or heterogeneity (Y or N)
If so, assess whether a transformation (such as a forth-root transformation) will address the violations and then make the appropriate corrections
• Consider a transformation (such as forth-root transform) and redo boxplots for the main effect (Season)
Show code
curdies.ag <- curdies %>% group_by(SEASON,SITE) %>% summarize(DUGESIA=mean(DUGESIA^0.25))
curdies.ag

Source: local data frame [6 x 3]
Groups: SEASON

SEASON SITE   DUGESIA
1 SUMMER    4 0.3508213
2 SUMMER    5 0.1804622
3 SUMMER    6 0.3811223
4 WINTER    1 1.1329421
5 WINTER    2 1.2718698
6 WINTER    3 0.8678707

ggplot(curdies.ag, aes(y=DUGESIA, x=SEASON)) + geom_boxplot()

Does a forth-root transformation improve the data with respect to normality and homogeneity of variance (Y or N)?
2. Exploratory data analysis has indicated that the response variable could be normalized via a forth-root transformation.

3. Fit the hierarchical model. Use with effects (or matrix) parameterization for the 'fixed effects' and means (or matrix) parameterization for the 'random effects'. Note the following points:
• Non-informative priors for residual precision (standard deviations) should be defined as vague uniform distributions rather than gamma distributions. This prevents the chains from walking to close to parameter estimates of 0. Alternatively, Gelman recommends the use of weakly informative half-cauchy distributions for variances in hierarchical designs.
• Finite-population standard deviation for Site (the within group factor) should reflect standard deviation between Sites within season (not standard deviation between all Sites). Unfortunately, the BUGS language does not offer a mechanism to determine which indexes match certain criteria (for example which Sites correspond to which Season), nor does it allow variables to be redefined. One way to allow us to generate multiple standard deviations is to create an offset vector that can be used to indicate the range of Site numbers that are within which Season.
• As this is a multi-level design in which the Season 'replicates' are determined by the Site means, it is important that firstly, the number of levels of the within Group factor (Season) reflects this, and secondly that the factor levels of the within group factor (Season) correspond to the factor levels of the group factor (Site). In this case, the Sites are labeled 1-6, but the first three Sites are in Winter. R will naturally label Winter as the second level (as alphabetically W comes after S for summer). Thus without intervention or careful awareness, it will appear that the Winter means are higher than the Summer means, when it is actually the other way around. One way to check this, is create the data list and check that all the factor levels are listed in increasing numerical order.
View preparation code
sites<-ddply(curdies, ~SITE, catcolwise(unique))
offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1) offset  [1] 4 1 7  curdies.list <- with(curdies, list(dugesia=DUGESIA, season=as.numeric(SEASON), site=as.numeric(SITE), n=nrow(curdies), nSeason=length(levels(as.factor(SEASON))), nSite=length(levels(as.factor(SITE))), SiteOffset=offset ) ) curdies.list  $dugesia
[1]  0.6477  6.0962  1.3106  1.7253  1.4594  1.0576  1.0163 16.1968  1.1681  1.0243  2.0113  3.6746
[13]  0.6891  1.2191  1.1131  0.6569  0.1361  0.2547  0.0000  0.0000  0.9411  0.0000  0.0000  1.5735
[25]  1.3745  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.1329  0.0000  0.6580  0.3745

$season [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1$site
[1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6

$n [1] 36$nSeason
[1] 2

$nSite [1] 6$SiteOffset
[1] 4 1 7

## notice that The Season variable has as many entries as there were plots
## rather than having as many as there were Sites.
## Also notice that Winter is the second Season according to the factor levels
## Finally notice that the offset values are not in order

#season <- with(curdies,SEASON[match(unique(SITE),SITE)])
#season <- factor(season, levels=unique(season))

curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER"))
sites<-ddply(curdies, ~SITE, catcolwise(unique))
offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1) offset  [1] 1 4 7  curdies.list <- with(curdies, list(dugesia=DUGESIA, season=as.numeric(SEASON), site=as.numeric(SITE), n=nrow(curdies), nSeason=length(levels(SEASON)), nSite=length(levels(as.factor(SITE))), SiteOffset=offset ) ) curdies.list  $dugesia
[1]  0.6477  6.0962  1.3106  1.7253  1.4594  1.0576  1.0163 16.1968  1.1681  1.0243  2.0113  3.6746
[13]  0.6891  1.2191  1.1131  0.6569  0.1361  0.2547  0.0000  0.0000  0.9411  0.0000  0.0000  1.5735
[25]  1.3745  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.1329  0.0000  0.6580  0.3745

$season [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2$site
[1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6

$n [1] 36$nSeason
[1] 2

$nSite [1] 6$SiteOffset
[1] 1 4 7

## That is better

Show full effects parameterization JAGS code
modelString="
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- alpha0 + alpha[SEASON[i]] + beta[SITE[i]]
y.err[i] <- y[i] - mu[i]
}

#Priors
alpha0 ~ dnorm(0, 1.0E-6)
alpha[1] <- 0
for (i in 2:nSEASON) {
alpha[i] ~ dnorm(0, 1.0E-6) #prior
}
for (i in 1:nSITE) {
beta[i] ~ dnorm(0, tau.B) #prior
}

#Uniform uninformative prior on variance
tau <- 1 / (sigma * sigma)
tau.B <- 1 / (sigma.B * sigma.B)
sigma ~ dunif(0, 100)
sigma.B ~ dunif(0, 100)

for (i in 1:nSEASON) {
SD.Site[i] <- sd(beta[SiteOffset[i]:SiteOffset[i+1]-1])
}
sd.Season <- sd(alpha)
sd.Site <- mean(SD.Site) #sd(beta.site)
sd.Resid <- sd(y.err)
}
"
writeLines(modelString,con="../downloads/BUGSscripts/ws9.2bS1.2a.txt")

# Generate a data list
curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER"))
sites<-ddply(curdies, ~SITE, catcolwise(unique))
offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1) offset  [1] 1 4 7  curdies$SITE <- factor(curdies$SITE) curdies.list <- with(curdies, list(y=S4DUGES, SITE=as.numeric(SITE), SEASON=as.numeric(SEASON), n=nrow(curdies), nSITE=length(levels(SITE)), nSEASON = length(levels(SEASON)), SiteOffset=offset ) ) curdies.list  $y
[1] 0.8971 1.5713 1.0700 1.1461 1.0991 1.0141 1.0040 2.0061 1.0396 1.0060 1.1909 1.3845 0.9111 1.0508 1.0272 0.9003 0.6074 0.7104 0.0000 0.0000 0.9849 0.0000 0.0000 1.1200 1.0828 0.0000 0.0000 0.0000
[29] 0.0000 0.0000 0.0000 0.0000 0.6038 0.0000 0.9007 0.7823

$SITE [1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6$SEASON
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

$n [1] 36$nSITE
[1] 6

$nSEASON [1] 2$SiteOffset
[1] 1 4 7

# define parameters
params <- c("alpha0","alpha","sigma","sigma.B","sd.Season","sd.Site","sd.Resid","SD.Site")
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
curdies.r2jags.a <- jags(data=curdies.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/ws9.2bS1.2a.txt",
n.chains=nChains,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 176

Initializing model

print(curdies.r2jags.a)

Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2a.txt", fit using jags,
3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10
n.sims = 3000 iterations saved
mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
SD.Site[1]   0.116   0.091  0.003  0.042  0.098  0.172  0.331 1.006  1300
SD.Site[2]   0.095   0.073  0.003  0.038  0.079  0.135  0.273 1.007   770
alpha[1]     0.000   0.000  0.000  0.000  0.000  0.000  0.000 1.000     1
alpha[2]    -0.789   0.239 -1.261 -0.911 -0.788 -0.662 -0.338 1.002  2100
alpha0       1.093   0.169  0.770  1.010  1.090  1.175  1.426 1.001  2300
sd.Resid     0.387   0.015  0.366  0.377  0.386  0.394  0.425 1.002  3000
sd.Season    0.560   0.160  0.247  0.468  0.557  0.645  0.892 1.005  1200
sd.Site      0.105   0.071  0.004  0.048  0.098  0.151  0.269 1.008   870
sigma        0.404   0.055  0.316  0.366  0.399  0.433  0.529 1.001  3000
sigma.B      0.173   0.164  0.005  0.066  0.128  0.229  0.614 1.008   560
deviance    35.022   3.511 29.717 32.605 34.482 36.802 43.431 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 = 6.2 and DIC = 41.2
DIC is an estimate of expected predictive error (lower deviance is better).

curdies.mcmc.a <- as.mcmc(curdies.r2jags.a)

Show matrix parameterization JAGS code
modelString="
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- inprod(alpha[],X[i,]) + beta[SITE[i]]
y.err[i] <- y[i] - mu[i]
}

#Priors
alpha ~ dmnorm(a0, A0)
for (i in 1:nSITE) {
beta[i] ~ dnorm(0, tau.B) #prior
}

#Uniform uninformative prior on variance
tau <- 1 / (sigma * sigma)
tau.B <- 1 / (sigma.B * sigma.B)
sigma ~ dunif(0, 100)
sigma.B ~ dunif(0, 100)

# Finite-population standard deviations
for (i in 1:nSEASON) {
SD.Site[i] <- sd(beta[SiteOffset[i]:SiteOffset[i+1]-1])
}
seasons[1] <- alpha[1]
for (i in 2:nSEASON) {
seasons[i] <- alpha[1]+alpha[i]
}
sd.Season <- sd(seasons)
sd.Site <- mean(SD.Site)
sd.Resid <- sd(y.err)

}
"
writeLines(modelString,con="../downloads/BUGSscripts/ws9.2bS1.2b.txt")

# Generate a data list
curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER"))
sites<-ddply(curdies, ~SITE, catcolwise(unique))
offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1) offset  [1] 1 4 7  curdies$SITE <- factor(curdies$SITE) Xmat <- model.matrix(~SEASON, data=curdies) curdies.list <- with(curdies, list(y=S4DUGES, SITE=as.numeric(SITE), X=Xmat, n=nrow(curdies), nSITE=length(levels(SITE)), nSEASON=length(levels(SEASON)), a0=rep(0, ncol(Xmat)), A0=diag(0, ncol(Xmat)), SiteOffset=offset ) ) curdies.list  $y
[1] 0.8971 1.5713 1.0700 1.1461 1.0991 1.0141 1.0040 2.0061 1.0396 1.0060 1.1909 1.3845 0.9111 1.0508 1.0272 0.9003 0.6074 0.7104 0.0000 0.0000 0.9849 0.0000 0.0000 1.1200 1.0828 0.0000 0.0000 0.0000
[29] 0.0000 0.0000 0.0000 0.0000 0.6038 0.0000 0.9007 0.7823

$SITE [1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6$X
(Intercept) SEASONSUMMER
1            1            0
2            1            0
3            1            0
4            1            0
5            1            0
6            1            0
7            1            0
8            1            0
9            1            0
10           1            0
11           1            0
12           1            0
13           1            0
14           1            0
15           1            0
16           1            0
17           1            0
18           1            0
19           1            1
20           1            1
21           1            1
22           1            1
23           1            1
24           1            1
25           1            1
26           1            1
27           1            1
28           1            1
29           1            1
30           1            1
31           1            1
32           1            1
33           1            1
34           1            1
35           1            1
36           1            1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$SEASON [1] "contr.treatment"$n
[1] 36

$nSITE [1] 6$nSEASON
[1] 2

$a0 [1] 0 0$A0
[,1] [,2]
[1,]    0    0
[2,]    0    0

$SiteOffset [1] 1 4 7  # define parameters params <- c("alpha","sigma","sigma.B","sd.Season","sd.Site","sd.Resid","SD.Site") burnInSteps = 3000 nChains = 3 numSavedSteps = 3000 thinSteps = 10 nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains) library(R2jags) curdies.r2jags.b <- jags(data=curdies.list, inits=NULL, parameters.to.save=params, model.file="../downloads/BUGSscripts/ws9.2bS1.2b.txt", n.chains=nChains, n.iter=nIter, n.burnin=burnInSteps, n.thin=thinSteps )  Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 258 Initializing model  print(curdies.r2jags.b)  Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2b.txt", fit using jags, 3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10 n.sims = 3000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff SD.Site[1] 0.115 0.090 0.004 0.041 0.096 0.171 0.324 1.015 490 SD.Site[2] 0.092 0.071 0.003 0.036 0.077 0.134 0.265 1.023 370 alpha[1] 1.088 0.170 0.756 1.001 1.085 1.176 1.425 1.005 2400 alpha[2] -0.786 0.230 -1.245 -0.907 -0.786 -0.663 -0.327 1.002 3000 sd.Resid 0.388 0.015 0.366 0.378 0.386 0.394 0.424 1.001 3000 sd.Season 0.557 0.158 0.237 0.469 0.556 0.642 0.880 1.001 3000 sd.Site 0.104 0.070 0.004 0.046 0.094 0.151 0.258 1.028 330 sigma 0.403 0.054 0.315 0.366 0.396 0.434 0.522 1.001 3000 sigma.B 0.169 0.157 0.006 0.064 0.132 0.226 0.552 1.022 440 deviance 35.038 3.445 29.819 32.678 34.494 36.737 43.569 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 = 5.9 and DIC = 41.0 DIC is an estimate of expected predictive error (lower deviance is better).  curdies.mcmc.b <- as.mcmc(curdies.r2jags.b)  Hierarchical parameterization JAGS code modelString=" model { #Likelihood #Likelihood (esimating site means (gamma.site) for (i in 1:n) { y[i]~dnorm(quad.means[i],tau) quad.means[i] <- gamma.site[SITE[i]] y.err[i] <- y[i] - quad.means[i] } for (i in 1:nSITE) { gamma.site[i] ~ dnorm(site.means[i], tau.site) site.means[i] <- inprod(beta[],X[i,]) } #Priors for (i in 1:nX) { beta[i] ~ dnorm(0, 1.0E-6) #prior } sigma ~ dunif(0, 100) tau <- 1 / (sigma * sigma) sigma.site ~ dunif(0, 100) tau.site <- 1 / (sigma.site * sigma.site) # Finite-population standard deviations for (i in 1:nSEASON) { SD.Site[i] <- sd(gamma.site[SiteOffset[i]:SiteOffset[i+1]-1]) } seasons[1] <- beta[1] for (i in 2:nSEASON) { seasons[i] <- beta[1]+beta[i] } sd.Season <- sd(seasons) sd.Site <- mean(SD.Site) sd.Resid <- sd(y.err) } " writeLines(modelString,con="../downloads/BUGSscripts/ws9.2bS1.2c.txt") # Generate a data list curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER")) sites<-ddply(curdies, ~SITE, catcolwise(unique)) offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1)
offset

[1] 1 4 7

curdies$SITE <- factor(curdies$SITE)
Xmat <- model.matrix(~SEASON, data=ddply(curdies, ~SITE, catcolwise(unique)))
curdies.list <- with(curdies,
list(y=S4DUGES,
SITE=as.numeric(SITE),
X=Xmat,
nX=ncol(Xmat),
n=nrow(curdies),
nSITE=length(levels(SITE)),
nSEASON=length(levels(SEASON)),
SiteOffset=offset
)
)
curdies.list

$y [1] 0.8971 1.5713 1.0700 1.1461 1.0991 1.0141 1.0040 2.0061 1.0396 1.0060 1.1909 1.3845 0.9111 1.0508 1.0272 0.9003 0.6074 0.7104 0.0000 0.0000 0.9849 0.0000 0.0000 1.1200 1.0828 0.0000 0.0000 0.0000 [29] 0.0000 0.0000 0.0000 0.0000 0.6038 0.0000 0.9007 0.7823$SITE
[1] 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 6 6 6 6 6 6

$X (Intercept) SEASONSUMMER 1 1 0 2 1 0 3 1 0 4 1 1 5 1 1 6 1 1 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$SEASON
[1] "contr.treatment"

$nX [1] 2$n
[1] 36

$nSITE [1] 6$nSEASON
[1] 2

$SiteOffset [1] 1 4 7  # define parameters params <- c("beta","sigma","sigma.site","sd.Season","sd.Site","sd.Resid") burnInSteps = 3000 nChains = 3 numSavedSteps = 3000 thinSteps = 10 nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains) library(R2jags) curdies.r2jags.c <- jags(data=curdies.list, inits=NULL, parameters.to.save=params, model.file="../downloads/BUGSscripts/ws9.2bS1.2c.txt", n.chains=nChains, n.iter=nIter, n.burnin=burnInSteps, n.thin=thinSteps )  Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 156 Initializing model  print(curdies.r2jags.c)  Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2c.txt", fit using jags, 3 chains, each with 13000 iterations (first 3000 discarded), n.thin = 10 n.sims = 3000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff beta[1] 1.089 0.172 0.757 1.001 1.086 1.178 1.414 1.001 3000 beta[2] -0.783 0.241 -1.255 -0.905 -0.780 -0.662 -0.350 1.001 3000 sd.Resid 0.388 0.015 0.366 0.377 0.386 0.394 0.423 1.002 1800 sd.Season 0.556 0.162 0.254 0.468 0.552 0.640 0.887 1.006 3000 sd.Site 0.104 0.071 0.003 0.047 0.094 0.150 0.261 1.003 970 sigma 0.402 0.053 0.314 0.365 0.396 0.433 0.515 1.002 1400 sigma.site 0.179 0.187 0.004 0.061 0.134 0.235 0.654 1.002 1100 deviance 34.970 3.365 29.820 32.577 34.405 36.866 42.802 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 = 5.7 and DIC = 40.6 DIC is an estimate of expected predictive error (lower deviance is better).  curdies.mcmc.c <- as.mcmc(curdies.r2jags.c)  Matrix parameterization STAN code Unfortunately, STAN does not yet support vector indexing or the ability to multiply a matrix by a vector. A work around involves multiplication of a dummy matrix with the gamma vector.. rstanString= # Generate a data list curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER")) sites<-ddply(curdies, ~SITE, catcolwise(unique)) offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1)
offset

[1] 1 4 7

curdies$SITE <- factor(curdies$SITE)
Xmat <- model.matrix(~SEASON, data=curdies)
Zmat <- model.matrix(~-1+SITE, data=curdies)
curdies.list <- with(curdies,
list(y=S4DUGES,
Z=Zmat,
X=Xmat,
nX=ncol(Xmat),
nZ=ncol(Zmat),
n=nrow(curdies),
a0=rep(0,ncol(Xmat)), A0=diag(100000,ncol(Xmat)),
SiteOffset=offset,
m=matrix(rep(1,6),2)
)
)
curdies.list

$y [1] 0.8971 1.5713 1.0700 1.1461 1.0991 1.0141 1.0040 2.0061 1.0396 1.0060 1.1909 1.3845 0.9111 1.0508 1.0272 0.9003 0.6074 0.7104 0.0000 0.0000 0.9849 0.0000 0.0000 1.1200 1.0828 0.0000 0.0000 0.0000 [29] 0.0000 0.0000 0.0000 0.0000 0.6038 0.0000 0.9007 0.7823$Z
SITE1 SITE2 SITE3 SITE4 SITE5 SITE6
1      1     0     0     0     0     0
2      1     0     0     0     0     0
3      1     0     0     0     0     0
4      1     0     0     0     0     0
5      1     0     0     0     0     0
6      1     0     0     0     0     0
7      0     1     0     0     0     0
8      0     1     0     0     0     0
9      0     1     0     0     0     0
10     0     1     0     0     0     0
11     0     1     0     0     0     0
12     0     1     0     0     0     0
13     0     0     1     0     0     0
14     0     0     1     0     0     0
15     0     0     1     0     0     0
16     0     0     1     0     0     0
17     0     0     1     0     0     0
18     0     0     1     0     0     0
19     0     0     0     1     0     0
20     0     0     0     1     0     0
21     0     0     0     1     0     0
22     0     0     0     1     0     0
23     0     0     0     1     0     0
24     0     0     0     1     0     0
25     0     0     0     0     1     0
26     0     0     0     0     1     0
27     0     0     0     0     1     0
28     0     0     0     0     1     0
29     0     0     0     0     1     0
30     0     0     0     0     1     0
31     0     0     0     0     0     1
32     0     0     0     0     0     1
33     0     0     0     0     0     1
34     0     0     0     0     0     1
35     0     0     0     0     0     1
36     0     0     0     0     0     1
attr(,"assign")
[1] 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$SITE [1] "contr.treatment"$X
(Intercept) SEASONSUMMER
1            1            0
2            1            0
3            1            0
4            1            0
5            1            0
6            1            0
7            1            0
8            1            0
9            1            0
10           1            0
11           1            0
12           1            0
13           1            0
14           1            0
15           1            0
16           1            0
17           1            0
18           1            0
19           1            1
20           1            1
21           1            1
22           1            1
23           1            1
24           1            1
25           1            1
26           1            1
27           1            1
28           1            1
29           1            1
30           1            1
31           1            1
32           1            1
33           1            1
34           1            1
35           1            1
36           1            1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$SEASON [1] "contr.treatment"$nX
[1] 2

$nZ [1] 6$n
[1] 36

$a0 [1] 0 0$A0
[,1]  [,2]
[1,] 1e+05 0e+00
[2,] 0e+00 1e+05

$SiteOffset [1] 1 4 7$m
[,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1

# define parameters
burnInSteps = 3000
nChains = 3
numSavedSteps = 3000
thinSteps = 10
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(rstan)
curdies.rstan.a <- stan(data=curdies.list,
model_code=rstanString,
pars=c('beta','gamma','sigma','sigma_Z','sd_Resid','sd_Season','sd_Site'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 0.42 seconds (Warm-up)
#                1.39 seconds (Sampling)
#                1.81 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 0.45 seconds (Warm-up)
#                1.56 seconds (Sampling)
#                2.01 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:     1 / 13000 [  0%]  (Warmup)
Iteration:  1300 / 13000 [ 10%]  (Warmup)
Iteration:  2600 / 13000 [ 20%]  (Warmup)
Iteration:  3001 / 13000 [ 23%]  (Sampling)
Iteration:  4300 / 13000 [ 33%]  (Sampling)
Iteration:  5600 / 13000 [ 43%]  (Sampling)
Iteration:  6900 / 13000 [ 53%]  (Sampling)
Iteration:  8200 / 13000 [ 63%]  (Sampling)
Iteration:  9500 / 13000 [ 73%]  (Sampling)
Iteration: 10800 / 13000 [ 83%]  (Sampling)
Iteration: 12100 / 13000 [ 93%]  (Sampling)
Iteration: 13000 / 13000 [100%]  (Sampling)
#  Elapsed Time: 0.41 seconds (Warm-up)
#                1.42 seconds (Sampling)
#                1.83 seconds (Total)

print(curdies.rstan.a)

Inference for Stan model: rstanString.
3 chains, each with iter=13000; warmup=3000; thin=10;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
beta[1]    1.09    0.00 0.18  0.72  0.99  1.09  1.18  1.42  3000    1
beta[2]   -0.78    0.01 0.24 -1.22 -0.91 -0.78 -0.65 -0.30  2248    1
gamma[1]   0.03    0.00 0.17 -0.30 -0.05  0.01  0.09  0.42  3000    1
gamma[2]   0.08    0.00 0.18 -0.22 -0.01  0.04  0.15  0.51  2431    1
gamma[3]  -0.09    0.00 0.19 -0.53 -0.17 -0.06  0.00  0.19  2813    1
gamma[4]   0.02    0.00 0.17 -0.33 -0.05  0.01  0.09  0.36  3000    1
gamma[5]  -0.06    0.00 0.17 -0.46 -0.12 -0.03  0.02  0.23  2882    1
gamma[6]   0.03    0.00 0.17 -0.29 -0.04  0.02  0.10  0.40  2857    1
sigma      0.40    0.00 0.05  0.32  0.37  0.40  0.44  0.53  1447    1
sigma_Z    0.19    0.00 0.18  0.02  0.07  0.14  0.24  0.62  1463    1
sd_Resid   0.39    0.00 0.02  0.37  0.38  0.39  0.39  0.42  2299    1
sd_Season  0.55    0.00 0.16  0.22  0.46  0.55  0.64  0.87  2129    1
sd_Site    0.11    0.00 0.07  0.01  0.06  0.10  0.16  0.26   874    1
lp__      22.43    0.20 4.50 13.77 19.47 22.24 25.20 31.53   518    1

Samples were drawn using NUTS(diag_e) at Sun Mar  8 15:28:34 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

curdies.mcmc.d <- rstan:::as.mcmc.list.stanfit(curdies.rstan.a)
curdies.mcmc.df.d <- as.data.frame(extract(curdies.rstan.a))

SD.site <- NULL
sd.site <- curdies.mcmc.df.d[,3:8]
for (i in 1:length(levels(curdies$SEASON))) { SD.site <- c(SD.site,apply(sd.site[,offset[i]:(offset[i+1]-1)],1,sd)) } data.frame(mean=mean(SD.site), median=median(SD.site), HPDinterval(as.mcmc(SD.site)), HPDinterval(as.mcmc(SD.site),p=0.68))   mean median lower upper lower.1 upper.1 var1 0.1115 0.09236 0.003046 0.2739 0.003046 0.1377  Hierarchical parameterization STAN code rstanString=" data{ int n; int nZ; int nX; vector [n] y; matrix [nZ,nX] X; matrix [n,nZ] Z; } parameters{ vector [nX] beta; real<lower=0> sigma; vector [nZ] gamma; real<lower=0> sigma_Z; } transformed parameters { vector [n] mu_Z; vector [nZ] mu; mu_Z <- Z*gamma; mu <- X*beta; } model{ // Priors beta ~ normal(0, 10000); gamma ~ normal( 0 , 1000); sigma_Z ~ uniform( 0 , 100 ); sigma ~ uniform( 0 , 100 ); y ~ normal(mu_Z, sigma); gamma ~ normal( mu , sigma_Z ); } generated quantities { vector [n] y_err; real<lower=0> sd_Resid; y_err <- y - mu_Z; sd_Resid <- sd(y_err); } " # Generate a data list curdies$SITE <- factor(curdies$SITE) Xmat <- model.matrix(~SEASON, data=ddply(curdies, ~SITE, catcolwise(unique))) Zmat <- model.matrix(~-1+SITE, data=curdies) curdies.list <- with(curdies, list(y=S4DUGES, Z=Zmat, X=Xmat, nX=ncol(Xmat), nZ=ncol(Zmat), n=nrow(curdies) ) ) curdies.list  $y
[1] 0.8971 1.5713 1.0700 1.1461 1.0991 1.0141 1.0040 2.0061 1.0396 1.0060 1.1909 1.3845 0.9111 1.0508 1.0272 0.9003 0.6074 0.7104 0.0000 0.0000 0.9849 0.0000 0.0000 1.1200 1.0828 0.0000 0.0000 0.0000
[29] 0.0000 0.0000 0.0000 0.0000 0.6038 0.0000 0.9007 0.7823

$Z SITE1 SITE2 SITE3 SITE4 SITE5 SITE6 1 1 0 0 0 0 0 2 1 0 0 0 0 0 3 1 0 0 0 0 0 4 1 0 0 0 0 0 5 1 0 0 0 0 0 6 1 0 0 0 0 0 7 0 1 0 0 0 0 8 0 1 0 0 0 0 9 0 1 0 0 0 0 10 0 1 0 0 0 0 11 0 1 0 0 0 0 12 0 1 0 0 0 0 13 0 0 1 0 0 0 14 0 0 1 0 0 0 15 0 0 1 0 0 0 16 0 0 1 0 0 0 17 0 0 1 0 0 0 18 0 0 1 0 0 0 19 0 0 0 1 0 0 20 0 0 0 1 0 0 21 0 0 0 1 0 0 22 0 0 0 1 0 0 23 0 0 0 1 0 0 24 0 0 0 1 0 0 25 0 0 0 0 1 0 26 0 0 0 0 1 0 27 0 0 0 0 1 0 28 0 0 0 0 1 0 29 0 0 0 0 1 0 30 0 0 0 0 1 0 31 0 0 0 0 0 1 32 0 0 0 0 0 1 33 0 0 0 0 0 1 34 0 0 0 0 0 1 35 0 0 0 0 0 1 36 0 0 0 0 0 1 attr(,"assign") [1] 1 1 1 1 1 1 attr(,"contrasts") attr(,"contrasts")$SITE
[1] "contr.treatment"

$X (Intercept) SEASONSUMMER 1 1 0 2 1 0 3 1 0 4 1 1 5 1 1 6 1 1 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$SEASON
[1] "contr.treatment"

$nX [1] 2$nZ
[1] 6

$n [1] 36  # define parameters burnInSteps = 3000 nChains = 3 numSavedSteps = 3000 thinSteps = 10 nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains) library(rstan) curdies.rstan.b <- stan(data=curdies.list, model_code=rstanString, pars=c('beta','gamma','sigma','sigma_Z','sd_Resid'), chains=nChains, iter=nIter, warmup=burnInSteps, thin=thinSteps, save_dso=TRUE )  TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW. SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 0.48 seconds (Warm-up) # 2.24 seconds (Sampling) # 2.72 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 0.39 seconds (Warm-up) # 0.94 seconds (Sampling) # 1.33 seconds (Total) SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3). Iteration: 1 / 13000 [ 0%] (Warmup) Iteration: 1300 / 13000 [ 10%] (Warmup) Iteration: 2600 / 13000 [ 20%] (Warmup) Iteration: 3001 / 13000 [ 23%] (Sampling) Iteration: 4300 / 13000 [ 33%] (Sampling) Iteration: 5600 / 13000 [ 43%] (Sampling) Iteration: 6900 / 13000 [ 53%] (Sampling) Iteration: 8200 / 13000 [ 63%] (Sampling) Iteration: 9500 / 13000 [ 73%] (Sampling) Iteration: 10800 / 13000 [ 83%] (Sampling) Iteration: 12100 / 13000 [ 93%] (Sampling) Iteration: 13000 / 13000 [100%] (Sampling) # Elapsed Time: 0.44 seconds (Warm-up) # 0.98 seconds (Sampling) # 1.42 seconds (Total)  print(curdies.rstan.b)  Inference for Stan model: rstanString. 3 chains, each with iter=13000; warmup=3000; thin=10; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta[1] 1.09 0.00 0.17 0.77 1.00 1.09 1.18 1.44 2750 1.00 beta[2] -0.78 0.01 0.25 -1.25 -0.91 -0.79 -0.65 -0.26 2189 1.00 gamma[1] 1.11 0.00 0.13 0.85 1.01 1.11 1.19 1.37 1590 1.00 gamma[2] 1.17 0.00 0.14 0.92 1.07 1.16 1.25 1.47 1851 1.00 gamma[3] 0.99 0.00 0.14 0.68 0.91 1.00 1.09 1.26 1739 1.00 gamma[4] 0.33 0.00 0.13 0.08 0.24 0.33 0.42 0.60 1059 1.00 gamma[5] 0.25 0.00 0.13 -0.03 0.17 0.26 0.34 0.50 1228 1.00 gamma[6] 0.34 0.00 0.13 0.09 0.25 0.34 0.43 0.60 1176 1.00 sigma 0.40 0.00 0.05 0.31 0.36 0.39 0.43 0.52 1665 1.00 sigma_Z 0.19 0.01 0.18 0.03 0.07 0.14 0.24 0.64 1065 1.00 sd_Resid 0.39 0.00 0.01 0.37 0.38 0.38 0.39 0.42 2307 1.00 lp__ 22.25 0.17 4.10 13.61 19.58 22.37 25.13 29.87 560 1.01 Samples were drawn using NUTS(diag_e) at Sun Mar 8 15:29:07 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).  curdies.mcmc.e <- rstan:::as.mcmc.list.stanfit(curdies.rstan.b) curdies.mcmc.df.e <- as.data.frame(extract(curdies.rstan.b))  #Finite-population standard deviations ## Seasons seasons <- NULL seasons <- cbind(curdies.mcmc.df.e[,'beta.1'], curdies.mcmc.df.e[,'beta.1']+curdies.mcmc.df.e[,'beta.2']) sd.season <- apply(seasons,1,sd) library(coda) data.frame(mean=mean(sd.season), median=median(sd.season), HPDinterval(as.mcmc(sd.season)), HPDinterval(as.mcmc(sd.season),p=0.68))   mean median lower upper lower.1 upper.1 var1 0.5536985 0.5585557 0.1999721 0.8859165 0.4131178 0.6909496  ## Sites SD.site <- NULL sd.site <- curdies.mcmc.df.e[,3:8] for (i in 1:length(levels(curdies$SEASON))) {
SD.site <-  c(SD.site,apply(sd.site[,offset[i]:(offset[i+1]-1)],1,sd))
}
data.frame(mean=mean(SD.site), median=median(SD.site), HPDinterval(as.mcmc(SD.site)), HPDinterval(as.mcmc(SD.site),p=0.68))

          mean     median       lower     upper    lower.1   upper.1
var1 0.1115323 0.09316411 0.001736477 0.2707392 0.01333439 0.1432832

BRM code
curdies.brm <- brm(S4DUGES ~ SEASON + (1|SITE), data=curdies,
family='gaussian',
prior=c(set_prior('normal(0,100)', class='b'),
set_prior('cauchy(0,5)', class='sd')),
n.chains=3, n.iter=2000, warmup=500, n.thin=2
)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 1, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 1, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 1, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 1, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 1, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 1, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 1, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.108542 seconds (Warm-up)
#                0.196446 seconds (Sampling)
#                0.304988 seconds (Total)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 2, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 2, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 2, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 2, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 2, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 2, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 2, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 2, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.104188 seconds (Warm-up)
#                0.295826 seconds (Sampling)
#                0.400014 seconds (Total)

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3, Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3, Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3, Iteration:  501 / 2000 [ 25%]  (Sampling)
Chain 3, Iteration:  700 / 2000 [ 35%]  (Sampling)
Chain 3, Iteration:  900 / 2000 [ 45%]  (Sampling)
Chain 3, Iteration: 1100 / 2000 [ 55%]  (Sampling)
Chain 3, Iteration: 1300 / 2000 [ 65%]  (Sampling)
Chain 3, Iteration: 1500 / 2000 [ 75%]  (Sampling)
Chain 3, Iteration: 1700 / 2000 [ 85%]  (Sampling)
Chain 3, Iteration: 1900 / 2000 [ 95%]  (Sampling)
Chain 3, Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.17552 seconds (Warm-up)
#                0.196082 seconds (Sampling)
#                0.371602 seconds (Total)

summary(curdies.brm)

 Family: gaussian (identity)
Formula: S4DUGES ~ SEASON + (1 | SITE)
Data: curdies (Number of observations: 36)
Samples: 3 chains, each with n.iter = 2000; n.warmup = 500; n.thin = 2;
total post-warmup samples = 2250
WAIC: 39.88

Random Effects:
~SITE (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.18      0.16     0.01     0.62        567    1

Fixed Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept        1.09      0.17     0.74     1.45        479 1.01
SEASONSUMMER    -0.79      0.22    -1.25    -0.34        713 1.00

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma(S4DUGES)      0.4      0.05     0.31     0.53       1726    1

Samples were drawn using NUTS(diag_e). For each parameter, Eff.Sample is a
crude measure of effective sample size, and Rhat is the potential scale
reduction factor on split chains (at convergence, Rhat = 1).

brms:::stancode(curdies.brm)

functions {
}
data {
int<lower=1> N;  # number of observations
vector[N] Y;  # response variable
int<lower=1> K;  # number of fixed effects
matrix[N, K] X;  # FE design matrix
# data for random effects of SITE
int<lower=1> J_1[N];  # RE levels
int<lower=1> N_1;  # number of levels
int<lower=1> K_1;  # number of REs
real Z_1[N];  # RE design matrix
}
transformed data {
}
parameters {
real b_Intercept;  # fixed effects Intercept
vector[K] b;  # fixed effects
vector[N_1] pre_1;  # unscaled REs
real<lower=0> sd_1;  # RE standard deviation
real<lower=0> sigma;  # residual SD
}
transformed parameters {
vector[N] eta;  # linear predictor
vector[N_1] r_1;  # REs
# compute linear predictor
eta <- X * b + b_Intercept;
r_1 <- sd_1 * (pre_1);  # scale REs
# if available add REs to linear predictor
for (n in 1:N) {
eta[n] <- eta[n] + Z_1[n] * r_1[J_1[n]];
}
}
model {
# prior specifications
b_Intercept ~ normal(0,100);
b ~ normal(0,100);
sd_1 ~ cauchy(0,5);
pre_1 ~ normal(0, 1);
sigma ~ cauchy(0, 5);
# likelihood contribution
Y ~ normal(eta, sigma);
}
generated quantities {
}

#brms:::standata(curdies.brm)

4. Before fully exploring the parameters, it is prudent to examine the convergence and mixing diagnostics. Chose either any of the parameterizations (they should yield much the same) - although it is sometimes useful to explore the different performances of effects vs matrix and JAGS vs STAN.
Full effects parameterization JAGS code
preds <- c("alpha0", "alpha[2]", "sigma", "sigma.B")
plot(curdies.mcmc.a)

autocorr.diag(curdies.mcmc.a[, preds])

             alpha0     alpha[2]       sigma      sigma.B
Lag 0    1.00000000  1.000000000  1.00000000 1.0000000000
Lag 10  -0.04819928 -0.047917862  0.03408172 0.4186761852
Lag 50  -0.02022011 -0.030419011  0.00106724 0.0190520320
Lag 100  0.02621796  0.022830915  0.01421471 0.0175050583
Lag 500  0.01966791  0.004624531 -0.01605507 0.0008104889

raftery.diag(curdies.mcmc.a)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

Matrix parameterization JAGS code
preds <- c("alpha[1]", "alpha[2]", "sigma", "sigma.B")
plot(curdies.mcmc.b)

autocorr.diag(curdies.mcmc.b[, preds])

           alpha[1]     alpha[2]        sigma      sigma.B
Lag 0    1.00000000  1.000000000  1.000000000  1.000000000
Lag 10   0.03371611  0.004028258  0.013521460  0.419591468
Lag 50  -0.01935447  0.017372840  0.006912629  0.010740829
Lag 100 -0.01690649 -0.011917936 -0.013496281 -0.012760873
Lag 500  0.02668688 -0.013412913  0.034986445  0.003179485

raftery.diag(curdies.mcmc.b)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

Hierarchical parameterization JAGS code
preds <- c("beta[1]", "beta[2]", "sigma", "sigma.site")
plot(curdies.mcmc.c)

autocorr.diag(curdies.mcmc.c[, preds])

             beta[1]      beta[2]        sigma   sigma.site
Lag 0    1.000000000  1.000000000  1.000000000  1.000000000
Lag 10   0.021287585  0.047841790 -0.001047102  0.559956291
Lag 50   0.018520042  0.010010499  0.018111468  0.147260934
Lag 100  0.005360734 -0.011870958  0.010250879  0.007888448
Lag 500 -0.011742350 -0.008742964  0.021524966 -0.016054557

raftery.diag(curdies.mcmc.c)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

Matrix parameterization STAN code
library(rstan)
grid.arrange(stan_trace(curdies.rstan.a, pars=c('beta','sigma'), ncol=1),
stan_dens(curdies.rstan.a, pars=c('beta','sigma'), separate_chains=TRUE, ncol=1),
ncol=2)

stan_ac(curdies.rstan.a, pars=c('beta','sigma'))

summary(do.call(rbind, args = get_sampler_params(curdies.rstan.a, inc_warmup = FALSE)),digits = 2)

 accept_stat__    stepsize__    treedepth__   n_leapfrog__ n_divergent__
Min.   :0.00   Min.   :0.22   Min.   :1.0   Min.   : 1    Min.   :0.000
1st Qu.:0.80   1st Qu.:0.22   1st Qu.:3.0   1st Qu.: 7    1st Qu.:0.000
Median :0.93   Median :0.23   Median :3.0   Median : 7    Median :0.000
Mean   :0.81   Mean   :0.22   Mean   :3.3   Mean   :11    Mean   :0.006
3rd Qu.:0.98   3rd Qu.:0.23   3rd Qu.:4.0   3rd Qu.:15    3rd Qu.:0.000
Max.   :1.00   Max.   :0.23   Max.   :7.0   Max.   :95    Max.   :1.000

stan_diag(curdies.rstan.a)

stan_diag(curdies.rstan.a, information='stepsize')

stan_diag(curdies.rstan.a, information='treedepth')

stan_diag(curdies.rstan.a, information='divergence')

library(gridExtra)
grid.arrange(stan_rhat(curdies.rstan.a)+theme_classic(8) ,
stan_ess(curdies.rstan.a)+theme_classic(8) ,
stan_mcse(curdies.rstan.a)+theme_classic(8) , ncol=2)

preds <- c("beta[1]", "beta[2]", "sigma", "sigma_Z")
plot(curdies.mcmc.d)

autocorr.diag(curdies.mcmc.d[, preds])

             beta[1]       beta[2]      sigma     sigma_Z
Lag 0    1.000000000  1.0000000000 1.00000000 1.000000000
Lag 10  -0.003829357 -0.0001842772 0.06968841 0.232121726
Lag 50  -0.005281874  0.0092503105 0.04862283 0.060227526
Lag 100  0.013319955  0.0488246158 0.01624396 0.014919291
Lag 500  0.003200964 -0.0245307715 0.03004036 0.001269244

raftery.diag(curdies.mcmc.d)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

Hierarchical parameterization STAN code
preds <- c("beta[1]", "beta[2]", "sigma", "sigma_Z")
plot(curdies.mcmc.e)

autocorr.diag(curdies.mcmc.e[, preds])

           beta[1]    beta[2]       sigma       sigma_Z
Lag 0   1.00000000 1.00000000  1.00000000  1.0000000000
Lag 10  0.04707077 0.07503569  0.14474725  0.1911481063
Lag 50  0.02265390 0.02383334  0.03022470  0.0512720335
Lag 100 0.02853164 0.01991633 -0.02099284  0.0055391709
Lag 500 0.03686456 0.02466837 -0.01936357 -0.0009199443

raftery.diag(curdies.mcmc.e)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

View code for BRM model
library(rstan)
grid.arrange(stan_trace(curdies.brm$fit, pars=c('b_Intercept','b_SEASONSUMMER', 'sigma_S4DUGES'), ncol=1), stan_dens(curdies.brm$fit, pars=c('b_Intercept','b_SEASONSUMMER',
'sigma_S4DUGES'), separate_chains=TRUE, ncol=1),
ncol=2)

stan_ac(stehman.brm$fit, pars=c('b_Intercept','b_SEASONSSUMMER', 'sigma_S4DUGES'))  Error in inherits(x, "stanreg"): object 'stehman.brm' not found  summary(do.call(rbind, args = get_sampler_params(curdies.brm$fit, inc_warmup = FALSE)),digits = 2)

 accept_stat__    stepsize__    treedepth__   n_leapfrog__ n_divergent__
Min.   :0.00   Min.   :0.15   Min.   :1.0   Min.   : 1    Min.   :0.0000
1st Qu.:0.90   1st Qu.:0.15   1st Qu.:4.0   1st Qu.:15    1st Qu.:0.0000
Median :0.97   Median :0.22   Median :4.0   Median :15    Median :0.0000
Mean   :0.90   Mean   :0.20   Mean   :4.1   Mean   :17    Mean   :0.0036
3rd Qu.:0.99   3rd Qu.:0.23   3rd Qu.:4.0   3rd Qu.:15    3rd Qu.:0.0000
Max.   :1.00   Max.   :0.23   Max.   :6.0   Max.   :63    Max.   :1.0000

stan_diag(curdies.brm$fit)  stan_diag(curdies.brm$fit, information='stepsize')

stan_diag(curdies.brm$fit, information='treedepth')  stan_diag(curdies.brm$fit, information='divergence')

library(gridExtra)
grid.arrange(stan_rhat(curdies.brm$fit)+theme_classic(8) , stan_ess(curdies.brm$fit)+theme_classic(8) ,
stan_mcse(curdies.brm$fit)+theme_classic(8) , ncol=2)  5. The variances in particular display substantial autocorrelation, perhaps it is worth exploring a weekly informative Half-Cauchy prior (scale = 25) rather than a non-informative flat prior (as suggested by Gelman). It might also be worth increasing the thinning rate from 10 to 100. Show full effects parameterization JAGS code modelString=" model { #Likelihood for (i in 1:n) { y[i]~dnorm(mu[i],tau) mu[i] <- alpha0 + alpha[SEASON[i]] + beta[SITE[i]] y.err[i] <- y[i]-mu[i] } #Priors alpha0 ~ dnorm(0, 1.0E-6) alpha[1] <- 0 for (i in 2:nSEASON) { alpha[i] ~ dnorm(0, 1.0E-6) #prior } for (i in 1:nSITE) { beta[i] ~ dnorm(0, tau.B) #prior } #Half-cauchy weekly informative prior (scale=25) tau <- pow(sigma,-2) sigma <- z/sqrt(chSq) z ~ dnorm(0, .0016)I(0,) chSq ~ dgamma(0.5, 0.5) tau.B <- pow(sigma.B,-2) sigma.B <- z.B/sqrt(chSq.B) z.B ~ dnorm(0, .0016)I(0,) chSq.B ~ dgamma(0.5, 0.5) #Finite-population standard deviations for (i in 1:nSEASON) { SD.Site[i] <- sd(beta[SiteOffset[i]:SiteOffset[i+1]-1]) } sd.Season <- sd(alpha) sd.Site <- mean(SD.Site) #sd(beta.site) sd.Resid <- sd(y.err) } " writeLines(modelString,con="../downloads/BUGSscripts/ws9.2bS1.2a.txt") # Generate a data list curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER")) sites<-ddply(curdies, ~SITE, catcolwise(unique)) offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1)
offset

[1] 1 4 7

curdies$SITE <- factor(curdies$SITE)
curdies.list <- with(curdies,
list(y=S4DUGES,
SITE=as.numeric(SITE),
SEASON=as.numeric(SEASON),
n=nrow(curdies),
nSITE=length(levels(SITE)),
nSEASON = length(levels(SEASON)),
SiteOffset=offset
)
)

# define parameters
params <- c("alpha0","alpha","sigma","sigma.B","sd.Season","sd.Site","sd.Resid","SD.Site")
burnInSteps = 6000
nChains = 3
numSavedSteps = 3000
thinSteps = 100
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
curdies.r2jags.a <- jags(data=curdies.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/ws9.2bS1.2a.txt",
n.chains=nChains,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 182

Initializing model

print(curdies.r2jags.a)

Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2a.txt", fit using jags,
3 chains, each with 106000 iterations (first 6000 discarded), n.thin = 100
n.sims = 3000 iterations saved
mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
SD.Site[1]   0.117   0.091  0.004  0.043  0.099  0.173  0.330 1.001  3000
SD.Site[2]   0.095   0.074  0.004  0.037  0.080  0.137  0.272 1.001  3000
alpha[1]     0.000   0.000  0.000  0.000  0.000  0.000  0.000 1.000     1
alpha[2]    -0.781   0.240 -1.198 -0.906 -0.786 -0.667 -0.306 1.001  3000
alpha0       1.091   0.179  0.753  1.008  1.091  1.181  1.405 1.004  3000
sd.Resid     0.387   0.015  0.366  0.377  0.385  0.394  0.423 1.001  3000
sd.Season    0.557   0.156  0.229  0.472  0.556  0.641  0.848 1.005  3000
sd.Site      0.106   0.071  0.005  0.050  0.096  0.153  0.260 1.001  3000
sigma        0.402   0.053  0.312  0.364  0.397  0.433  0.518 1.001  2900
sigma.B      0.179   0.196  0.006  0.067  0.130  0.227  0.641 1.001  3000
deviance    34.963   3.506 29.762 32.500 34.456 36.737 43.734 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 = 6.1 and DIC = 41.1
DIC is an estimate of expected predictive error (lower deviance is better).

curdies.mcmc.a <- as.mcmc(curdies.r2jags.a)
preds <- c("alpha0", "alpha[2]", "sigma", "sigma.B")
plot(curdies.mcmc.a)

autocorr.diag(curdies.mcmc.a[, preds])

               alpha0      alpha[2]        sigma      sigma.B
Lag 0     1.000000000  1.0000000000  1.000000000  1.000000000
Lag 100  -0.010755755 -0.0540961677 -0.021166646  0.086534511
Lag 500   0.002678608 -0.0009848874 -0.008121714 -0.022853723
Lag 1000  0.033051082  0.0481177854 -0.037233461 -0.011397839
Lag 5000 -0.015636563 -0.0303283480 -0.016482653 -0.004621383

raftery.diag(curdies.mcmc.a)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

Show matrix parameterization JAGS code
modelString="
model {
#Likelihood
for (i in 1:n) {
y[i]~dnorm(mu[i],tau)
mu[i] <- inprod(alpha[],X[i,]) + beta[SITE[i]]
y.err[i] <- y[i] - mu[i]
}

#Priors
alpha ~ dmnorm(a0, A0)
for (i in 1:nSITE) {
beta[i] ~ dnorm(0, tau.B) #prior
}

#Half-cauchy weekly informative prior (scale=25)
tau <- pow(sigma,-2)
sigma <- z/sqrt(chSq)
z ~ dnorm(0, .0016)I(0,)
chSq ~ dgamma(0.5, 0.5)

tau.B <- pow(sigma.B,-2)
sigma.B <- z.B/sqrt(chSq.B)
z.B ~ dnorm(0, .0016)I(0,)
chSq.B ~ dgamma(0.5, 0.5)

# Finite-population standard deviations
for (i in 1:nSEASON) {
SD.Site[i] <- sd(beta[SiteOffset[i]:SiteOffset[i+1]-1])
}
seasons[1] <- alpha[1]
for (i in 2:nSEASON) {
seasons[i] <- alpha[1]+alpha[i]
}
sd.Season <- sd(seasons)
sd.Site <- mean(SD.Site)
sd.Resid <- sd(y.err)

}
"
writeLines(modelString,con="../downloads/BUGSscripts/ws9.2bS1.2b.txt")

# Generate a data list
curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER"))
sites<-ddply(curdies, ~SITE, catcolwise(unique))
offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1) curdies$SITE <- factor(curdies$SITE) Xmat <- model.matrix(~SEASON, data=curdies) curdies.list <- with(curdies, list(y=S4DUGES, SITE=as.numeric(SITE), X=Xmat, n=nrow(curdies), nSEASON=length(levels(SEASON)), nSITE=length(levels(SITE)), a0=rep(0, ncol(Xmat)), A0=diag(0, ncol(Xmat)), SiteOffset=offset ) ) # define parameters params <- c("alpha","sigma","sigma.B","sd.Season","sd.Site","sd.Resid","SD.Site") burnInSteps = 6000 nChains = 3 numSavedSteps = 3000 thinSteps = 100 nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains) library(R2jags) curdies.r2jags.b <- jags(data=curdies.list, inits=NULL, parameters.to.save=params, model.file="../downloads/BUGSscripts/ws9.2bS1.2b.txt", n.chains=nChains, n.iter=nIter, n.burnin=burnInSteps, n.thin=thinSteps )  Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 264 Initializing model  print(curdies.r2jags.b)  Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2b.txt", fit using jags, 3 chains, each with 106000 iterations (first 6000 discarded), n.thin = 100 n.sims = 3000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff SD.Site[1] 0.117 0.091 0.003 0.043 0.099 0.174 0.329 1.001 3000 SD.Site[2] 0.093 0.074 0.003 0.034 0.077 0.137 0.275 1.001 3000 alpha[1] 1.089 0.173 0.729 1.004 1.090 1.181 1.422 1.002 3000 alpha[2] -0.782 0.256 -1.261 -0.912 -0.786 -0.660 -0.287 1.005 3000 sd.Resid 0.388 0.015 0.366 0.378 0.386 0.394 0.424 1.001 3000 sd.Season 0.558 0.165 0.219 0.467 0.556 0.645 0.899 1.001 3000 sd.Site 0.105 0.072 0.004 0.045 0.097 0.151 0.262 1.001 3000 sigma 0.401 0.052 0.314 0.364 0.397 0.434 0.520 1.001 3000 sigma.B 0.178 0.187 0.005 0.063 0.134 0.235 0.632 1.001 3000 deviance 35.048 3.356 30.013 32.736 34.495 36.882 43.176 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 = 5.6 and DIC = 40.7 DIC is an estimate of expected predictive error (lower deviance is better).  curdies.mcmc.b <- as.mcmc(curdies.r2jags.b) preds <- c("alpha[1]", "alpha[2]", "sigma", "sigma.B") plot(curdies.mcmc.b)  autocorr.diag(curdies.mcmc.b[, preds])   alpha[1] alpha[2] sigma sigma.B Lag 0 1.000000 1.000000 1.000000 1.000000 Lag 100 0.013155 -0.008749 -0.016552 0.060065 Lag 500 -0.001854 0.006032 -0.006479 -0.011143 Lag 1000 0.023152 -0.010042 0.012903 0.019126 Lag 5000 -0.003289 -0.009786 -0.032329 -0.004125  raftery.diag(curdies.mcmc.b)  [[1]] Quantile (q) = 0.025 Accuracy (r) = +/- 0.005 Probability (s) = 0.95 You need a sample size of at least 3746 with these values of q, r and s [[2]] Quantile (q) = 0.025 Accuracy (r) = +/- 0.005 Probability (s) = 0.95 You need a sample size of at least 3746 with these values of q, r and s [[3]] Quantile (q) = 0.025 Accuracy (r) = +/- 0.005 Probability (s) = 0.95 You need a sample size of at least 3746 with these values of q, r and s  Hierarchical parameterization JAGS code modelString= writeLines(modelString,con="../downloads/BUGSscripts/ws9.2bS1.2c.txt") # Generate a data list curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER")) sites<-ddply(curdies, ~SITE, catcolwise(unique)) offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1)
curdies$SITE <- factor(curdies$SITE)
Xmat <- model.matrix(~SEASON, data=ddply(curdies, ~SITE, catcolwise(unique)))
curdies.list <- with(curdies,
list(y=S4DUGES,
SITE=as.numeric(SITE),
X=Xmat,
nX=ncol(Xmat),
n=nrow(curdies),
nSITE=length(levels(SITE)),
nSEASON=length(levels(SEASON)),
SiteOffset=offset
)
)

# define parameters
params <- c("beta","sigma","sigma.site","sd.Season","sd.Site","sd.Resid")
burnInSteps = 6000
nChains = 3
numSavedSteps = 3000
thinSteps = 100
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(R2jags)
curdies.r2jags.c <- jags(data=curdies.list,
inits=NULL,
parameters.to.save=params,
model.file="../downloads/BUGSscripts/ws9.2bS1.2c.txt",
n.chains=nChains,
n.iter=nIter,
n.burnin=burnInSteps,
n.thin=thinSteps
)

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 162

Initializing model

print(curdies.r2jags.c)

Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2c.txt", fit using jags,
3 chains, each with 106000 iterations (first 6000 discarded), n.thin = 100
n.sims = 3000 iterations saved
mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
beta[1]      1.089   0.172  0.755  1.000  1.090  1.180  1.395 1.002  1400
beta[2]     -0.783   0.248 -1.219 -0.907 -0.790 -0.664 -0.324 1.006  1100
sd.Resid     0.387   0.015  0.366  0.378  0.386  0.393  0.423 1.001  3000
sd.Season    0.558   0.161  0.239  0.470  0.560  0.642  0.865 1.003  2200
sd.Site      0.107   0.071  0.005  0.047  0.098  0.154  0.262 1.003  1800
sigma        0.402   0.054  0.313  0.364  0.396  0.434  0.525 1.001  3000
sigma.site   0.176   0.183  0.008  0.063  0.131  0.232  0.603 1.001  3000
deviance    34.992   3.396 29.839 32.650 34.470 36.755 42.968 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 = 5.8 and DIC = 40.8
DIC is an estimate of expected predictive error (lower deviance is better).

curdies.mcmc.c <- as.mcmc(curdies.r2jags.c)
preds <- c("beta[1]", "beta[2]", "sigma", "sigma.site")
plot(curdies.mcmc.c)

autocorr.diag(curdies.mcmc.c[, preds])

           beta[1]   beta[2]     sigma sigma.site
Lag 0     1.000000  1.000000 1.0000000   1.000000
Lag 100   0.011328 -0.002593 0.0001509   0.024686
Lag 500  -0.012719 -0.028442 0.0072140  -0.003313
Lag 1000  0.007359  0.005217 0.0399172   0.021000
Lag 5000  0.003469 -0.012810 0.0291639   0.007214

raftery.diag(curdies.mcmc.c)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

Matrix parameterization STAN code
rstanString="
data{
int n;
int nZ;
int nX;
vector [n] y;
matrix [n,nX] X;
matrix [n,nZ] Z;
vector [nX] a0;
matrix [nX,nX] A0;
}

parameters{
vector [nX] beta;
real<lower=0> sigma;
vector [nZ] gamma;
real<lower=0> sigma_Z;
}
transformed parameters {
vector [n] mu;

mu <- X*beta + Z*gamma;
}
model{
// Priors
beta ~ multi_normal(a0,A0);
gamma ~ normal( 0 , sigma_Z );
sigma_Z ~ cauchy(0,25);
sigma ~ cauchy(0,25);

y ~ normal( mu , sigma );
}
generated quantities {
vector [n] y_err;
real<lower=0> sd_Resid;

y_err <- y - mu;
sd_Resid <- sd(y_err);

}
"

# Generate a data list
curdies$SITE <- factor(curdies$SITE)
Xmat <- model.matrix(~SEASON, data=curdies)
Zmat <- model.matrix(~-1+SITE, data=curdies)
curdies.list <- with(curdies,
list(y=S4DUGES,
Z=Zmat,
X=Xmat,
nX=ncol(Xmat),
nZ=ncol(Zmat),
n=nrow(curdies),
a0=rep(0,ncol(Xmat)), A0=diag(100000,ncol(Xmat))
)
)

# define parameters
burnInSteps = 6000
nChains = 3
numSavedSteps = 3000
thinSteps = 100
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(rstan)
curdies.rstan.a <- stan(data=curdies.list,
model_code=rstanString,
pars=c('beta','gamma','sigma','sigma_Z', 'sd_Resid'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:      1 / 106000 [  0%]  (Warmup)
Iteration:   6001 / 106000 [  5%]  (Sampling)
Iteration:  16600 / 106000 [ 15%]  (Sampling)
Iteration:  27200 / 106000 [ 25%]  (Sampling)
Iteration:  37800 / 106000 [ 35%]  (Sampling)
Iteration:  48400 / 106000 [ 45%]  (Sampling)
Iteration:  59000 / 106000 [ 55%]  (Sampling)
Iteration:  69600 / 106000 [ 65%]  (Sampling)
Iteration:  80200 / 106000 [ 75%]  (Sampling)
Iteration:  90800 / 106000 [ 85%]  (Sampling)
Iteration: 101400 / 106000 [ 95%]  (Sampling)
Iteration: 106000 / 106000 [100%]  (Sampling)
#  Elapsed Time: 0.81 seconds (Warm-up)
#                21.03 seconds (Sampling)
#                21.84 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:      1 / 106000 [  0%]  (Warmup)
Iteration:   6001 / 106000 [  5%]  (Sampling)
Iteration:  16600 / 106000 [ 15%]  (Sampling)
Iteration:  27200 / 106000 [ 25%]  (Sampling)
Iteration:  37800 / 106000 [ 35%]  (Sampling)
Iteration:  48400 / 106000 [ 45%]  (Sampling)
Iteration:  59000 / 106000 [ 55%]  (Sampling)
Iteration:  69600 / 106000 [ 65%]  (Sampling)
Iteration:  80200 / 106000 [ 75%]  (Sampling)
Iteration:  90800 / 106000 [ 85%]  (Sampling)
Iteration: 101400 / 106000 [ 95%]  (Sampling)
Iteration: 106000 / 106000 [100%]  (Sampling)
#  Elapsed Time: 0.85 seconds (Warm-up)
#                15.14 seconds (Sampling)
#                15.99 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:      1 / 106000 [  0%]  (Warmup)
Iteration:   6001 / 106000 [  5%]  (Sampling)
Iteration:  16600 / 106000 [ 15%]  (Sampling)
Iteration:  27200 / 106000 [ 25%]  (Sampling)
Iteration:  37800 / 106000 [ 35%]  (Sampling)
Iteration:  48400 / 106000 [ 45%]  (Sampling)
Iteration:  59000 / 106000 [ 55%]  (Sampling)
Iteration:  69600 / 106000 [ 65%]  (Sampling)
Iteration:  80200 / 106000 [ 75%]  (Sampling)
Iteration:  90800 / 106000 [ 85%]  (Sampling)
Iteration: 101400 / 106000 [ 95%]  (Sampling)
Iteration: 106000 / 106000 [100%]  (Sampling)
#  Elapsed Time: 0.83 seconds (Warm-up)
#                15.24 seconds (Sampling)
#                16.07 seconds (Total)

print(curdies.rstan.a)

Inference for Stan model: rstanString.
3 chains, each with iter=106000; warmup=6000; thin=100;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
beta[1]   1.09    0.00 0.18  0.77  1.00  1.09  1.18  1.43  2224    1
beta[2]  -0.79    0.01 0.25 -1.27 -0.92 -0.79 -0.66 -0.32  2202    1
gamma[1]  0.01    0.00 0.18 -0.31 -0.05  0.01  0.08  0.37  2644    1
gamma[2]  0.08    0.00 0.19 -0.24 -0.01  0.04  0.15  0.47  2786    1
gamma[3] -0.10    0.00 0.19 -0.53 -0.18 -0.06  0.00  0.18  2588    1
gamma[4]  0.02    0.00 0.18 -0.30 -0.04  0.01  0.09  0.37  3000    1
gamma[5] -0.05    0.00 0.18 -0.43 -0.12 -0.03  0.02  0.27  3000    1
gamma[6]  0.04    0.00 0.18 -0.27 -0.04  0.02  0.10  0.42  3000    1
sigma     0.40    0.00 0.05  0.31  0.37  0.40  0.44  0.52  2934    1
sigma_Z   0.19    0.00 0.21  0.02  0.07  0.14  0.24  0.68  2542    1
sd_Resid  0.39    0.00 0.02  0.37  0.38  0.39  0.39  0.43  1517    1
lp__     22.62    0.11 4.60 13.78 19.56 22.38 25.45 32.09  1605    1

Samples were drawn using NUTS(diag_e) at Sun Mar  8 15:31:25 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

curdies.mcmc.d <- rstan:::as.mcmc.list.stanfit(curdies.rstan.a)
curdies.mcmc.df.d <- as.data.frame(extract(curdies.rstan.a))

preds <- c("beta[1]", "beta[2]", "sigma", "sigma_Z")
plot(curdies.mcmc.d)

autocorr.diag(curdies.mcmc.d[, preds])

          beta[1]  beta[2]     sigma   sigma_Z
Lag 0    1.000000  1.00000  1.000000  1.000000
Lag 100  0.049602  0.03405  0.011961  0.040739
Lag 500  0.001999  0.01131 -0.003834 -0.009070
Lag 1000 0.001019 -0.01372  0.022460  0.027303
Lag 5000 0.023995  0.01447  0.001216 -0.007945

raftery.diag(curdies.mcmc.d)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

#Finite-population standard deviations
## Site
curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER"))
sites<-ddply(curdies, ~SITE, catcolwise(unique))

Error in eval(expr, envir, enclos): could not find function "ddply"

offset <- c(match(levels(sites$SEASON),sites$SEASON), length(sites$SEASON)+1) SD.site <- NULL sd.site <- curdies.mcmc.df.d[,3:8] for (i in 1:length(levels(curdies$SEASON))) {
SD.site <-  c(SD.site,apply(sd.site[,offset[i]:(offset[i+1]-1)],1,sd))
}
data.frame(mean=mean(SD.site), median=median(SD.site), HPDinterval(as.mcmc(SD.site)), HPDinterval(as.mcmc(SD.site),p=0.68))

          mean     median       lower     upper     lower.1   upper.1
var1 0.1091143 0.09009361 0.001057878 0.2720609 0.007046607 0.1374872

#Season
seasons <- NULL
seasons <- cbind(curdies.mcmc.df.d[,'beta.1'],
curdies.mcmc.df.d[,'beta.1']+curdies.mcmc.df.d[,'beta.2'])
sd.season <- apply(seasons,1,sd)
library(coda)
data.frame(mean=mean(sd.season), median=median(sd.season), HPDinterval(as.mcmc(sd.season)), HPDinterval(as.mcmc(sd.season),p=0.68))

          mean    median     lower     upper   lower.1  upper.1
var1 0.5604005 0.5600543 0.2238934 0.8789284 0.4244865 0.701144

Hierarchical parameterization STAN code
rstanString="
data{
int n;
int nZ;
int nX;
vector [n] y;
matrix [nZ,nX] X;
matrix [n,nZ] Z;
}

parameters{
vector [nX] beta;
real<lower=0> sigma;
vector [nZ] gamma;
real<lower=0> sigma_Z;
}
transformed parameters {
vector [n] mu_Z;
vector [nZ] mu;

mu_Z <- Z*gamma;
mu <- X*beta;
}
model{

// Priors
beta ~ normal(0, 10000);
gamma ~ normal( 0 , 1000);
sigma_Z ~ uniform( 0 , 100 );
sigma ~ uniform( 0 , 100 );

y ~ normal(mu_Z, sigma);
gamma ~ normal( mu , sigma_Z );
}
generated quantities {
vector [n] y_err;
real<lower=0> sd_Resid;

y_err <- y - mu_Z;
sd_Resid <- sd(y_err);

}
"

# Generate a data list
curdies$SITE <- factor(curdies$SITE)
Xmat <- model.matrix(~SEASON, data=ddply(curdies, ~SITE, catcolwise(unique)))
Zmat <- model.matrix(~-1+SITE, data=curdies)
curdies.list <- with(curdies,
list(y=S4DUGES,
Z=Zmat,
X=Xmat,
nX=ncol(Xmat),
nZ=ncol(Zmat),
n=nrow(curdies)
)
)
curdies.list

$y [1] 0.8971 1.5713 1.0700 1.1461 1.0991 1.0141 1.0040 2.0061 1.0396 1.0060 1.1909 1.3845 0.9111 1.0508 1.0272 0.9003 0.6074 0.7104 0.0000 0.0000 0.9849 0.0000 0.0000 1.1200 1.0828 0.0000 0.0000 0.0000 [29] 0.0000 0.0000 0.0000 0.0000 0.6038 0.0000 0.9007 0.7823$Z
SITE1 SITE2 SITE3 SITE4 SITE5 SITE6
1      1     0     0     0     0     0
2      1     0     0     0     0     0
3      1     0     0     0     0     0
4      1     0     0     0     0     0
5      1     0     0     0     0     0
6      1     0     0     0     0     0
7      0     1     0     0     0     0
8      0     1     0     0     0     0
9      0     1     0     0     0     0
10     0     1     0     0     0     0
11     0     1     0     0     0     0
12     0     1     0     0     0     0
13     0     0     1     0     0     0
14     0     0     1     0     0     0
15     0     0     1     0     0     0
16     0     0     1     0     0     0
17     0     0     1     0     0     0
18     0     0     1     0     0     0
19     0     0     0     1     0     0
20     0     0     0     1     0     0
21     0     0     0     1     0     0
22     0     0     0     1     0     0
23     0     0     0     1     0     0
24     0     0     0     1     0     0
25     0     0     0     0     1     0
26     0     0     0     0     1     0
27     0     0     0     0     1     0
28     0     0     0     0     1     0
29     0     0     0     0     1     0
30     0     0     0     0     1     0
31     0     0     0     0     0     1
32     0     0     0     0     0     1
33     0     0     0     0     0     1
34     0     0     0     0     0     1
35     0     0     0     0     0     1
36     0     0     0     0     0     1
attr(,"assign")
[1] 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$SITE [1] "contr.treatment"$X
(Intercept) SEASONSUMMER
1           1            0
2           1            0
3           1            0
4           1            1
5           1            1
6           1            1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$SEASON [1] "contr.treatment"$nX
[1] 2

$nZ [1] 6$n
[1] 36

# define parameters
burnInSteps = 6000
nChains = 3
numSavedSteps = 3000
thinSteps = 100
nIter = burnInSteps+ceiling((numSavedSteps * thinSteps)/nChains)

library(rstan)
curdies.rstan.b <- stan(data=curdies.list,
model_code=rstanString,
pars=c('beta','gamma','sigma','sigma_Z','sd_Resid'),
chains=nChains,
iter=nIter,
warmup=burnInSteps,
thin=thinSteps,
save_dso=TRUE
)

TRANSLATING MODEL 'rstanString' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'rstanString' NOW.

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 1).

Iteration:      1 / 106000 [  0%]  (Warmup)
Iteration:   6001 / 106000 [  5%]  (Sampling)
Iteration:  16600 / 106000 [ 15%]  (Sampling)
Iteration:  27200 / 106000 [ 25%]  (Sampling)
Iteration:  37800 / 106000 [ 35%]  (Sampling)
Iteration:  48400 / 106000 [ 45%]  (Sampling)
Iteration:  59000 / 106000 [ 55%]  (Sampling)
Iteration:  69600 / 106000 [ 65%]  (Sampling)
Iteration:  80200 / 106000 [ 75%]  (Sampling)
Iteration:  90800 / 106000 [ 85%]  (Sampling)
Iteration: 101400 / 106000 [ 95%]  (Sampling)
Iteration: 106000 / 106000 [100%]  (Sampling)
#  Elapsed Time: 0.7 seconds (Warm-up)
#                11.92 seconds (Sampling)
#                12.62 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 2).

Iteration:      1 / 106000 [  0%]  (Warmup)
Iteration:   6001 / 106000 [  5%]  (Sampling)
Iteration:  16600 / 106000 [ 15%]  (Sampling)
Iteration:  27200 / 106000 [ 25%]  (Sampling)
Iteration:  37800 / 106000 [ 35%]  (Sampling)
Iteration:  48400 / 106000 [ 45%]  (Sampling)
Iteration:  59000 / 106000 [ 55%]  (Sampling)
Iteration:  69600 / 106000 [ 65%]  (Sampling)
Iteration:  80200 / 106000 [ 75%]  (Sampling)
Iteration:  90800 / 106000 [ 85%]  (Sampling)
Iteration: 101400 / 106000 [ 95%]  (Sampling)
Iteration: 106000 / 106000 [100%]  (Sampling)
#  Elapsed Time: 0.73 seconds (Warm-up)
#                16.12 seconds (Sampling)
#                16.85 seconds (Total)

SAMPLING FOR MODEL 'rstanString' NOW (CHAIN 3).

Iteration:      1 / 106000 [  0%]  (Warmup)
Iteration:   6001 / 106000 [  5%]  (Sampling)
Iteration:  16600 / 106000 [ 15%]  (Sampling)
Iteration:  27200 / 106000 [ 25%]  (Sampling)
Iteration:  37800 / 106000 [ 35%]  (Sampling)
Iteration:  48400 / 106000 [ 45%]  (Sampling)
Iteration:  59000 / 106000 [ 55%]  (Sampling)
Iteration:  69600 / 106000 [ 65%]  (Sampling)
Iteration:  80200 / 106000 [ 75%]  (Sampling)
Iteration:  90800 / 106000 [ 85%]  (Sampling)
Iteration: 101400 / 106000 [ 95%]  (Sampling)
Iteration: 106000 / 106000 [100%]  (Sampling)
#  Elapsed Time: 0.74 seconds (Warm-up)
#                16.8 seconds (Sampling)
#                17.54 seconds (Total)

print(curdies.rstan.b)

Inference for Stan model: rstanString.
3 chains, each with iter=106000; warmup=6000; thin=100;
post-warmup draws per chain=1000, total post-warmup draws=3000.

mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
beta[1]   1.09    0.00 0.19  0.72  1.01  1.09  1.18  1.41  2882    1
beta[2]  -0.78    0.00 0.27 -1.22 -0.92 -0.79 -0.66 -0.33  2879    1
gamma[1]  1.11    0.00 0.13  0.85  1.03  1.11  1.19  1.37  2668    1
gamma[2]  1.17    0.00 0.14  0.91  1.08  1.17  1.26  1.47  2491    1
gamma[3]  0.99    0.00 0.15  0.69  0.90  1.00  1.09  1.25  2716    1
gamma[4]  0.32    0.00 0.13  0.06  0.24  0.32  0.41  0.60  2799    1
gamma[5]  0.25    0.00 0.14 -0.04  0.16  0.25  0.34  0.50  2904    1
gamma[6]  0.34    0.00 0.13  0.09  0.26  0.34  0.42  0.61  2844    1
sigma     0.40    0.00 0.05  0.31  0.37  0.40  0.43  0.52  2690    1
sigma_Z   0.19    0.00 0.20  0.03  0.08  0.14  0.25  0.64  2545    1
sd_Resid  0.39    0.00 0.01  0.37  0.38  0.38  0.39  0.42  3000    1
lp__     22.16    0.11 4.11 13.75 19.57 22.15 24.82 30.24  1372    1

Samples were drawn using NUTS(diag_e) at Sun Mar  8 15:32:40 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

curdies.mcmc.e <- rstan:::as.mcmc.list.stanfit(curdies.rstan.b)
curdies.mcmc.df.e <- as.data.frame(extract(curdies.rstan.b))

preds <- c("beta[1]", "beta[2]", "sigma", "sigma_Z")
plot(curdies.mcmc.e)

autocorr.diag(curdies.mcmc.e[, preds])

          beta[1]   beta[2]     sigma   sigma_Z
Lag 0    1.000000  1.000000  1.000000  1.000000
Lag 100  0.005039  0.020794  0.034445  0.039067
Lag 500  0.010435 -0.009136 -0.010605  0.001239
Lag 1000 0.008577  0.010613  0.015420 -0.015927
Lag 5000 0.015266  0.022877 -0.008721  0.003651

raftery.diag(curdies.mcmc.e)

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

You need a sample size of at least 3746 with these values of q, r and s

#Finite-population standard deviations
## Seasons
seasons <- NULL
seasons <- cbind(curdies.mcmc.df.e[,'beta.1'],
curdies.mcmc.df.e[,'beta.1']+curdies.mcmc.df.e[,'beta.2'])
sd.season <- apply(seasons,1,sd)
library(coda)
data.frame(mean=mean(sd.season), median=median(sd.season), HPDinterval(as.mcmc(sd.season)), HPDinterval(as.mcmc(sd.season),p=0.68))

          mean   median     lower     upper   lower.1   upper.1
var1 0.5614793 0.558442 0.2493118 0.8651797 0.4270324 0.6952875

## Sites
SD.site <- NULL
sd.site <- curdies.mcmc.df.e[,3:8]
for (i in 1:length(levels(curdies$SEASON))) { SD.site <- c(SD.site,apply(sd.site[,offset[i]:(offset[i+1]-1)],1,sd)) } data.frame(mean=mean(SD.site), median=median(SD.site), HPDinterval(as.mcmc(SD.site)), HPDinterval(as.mcmc(SD.site),p=0.68))   mean median lower upper lower.1 upper.1 var1 0.1134773 0.09522703 0.00400208 0.2749261 0.01109254 0.1420809  6. Explore the parameter estimates Full effects parameterization JAGS code library(R2jags) print(curdies.r2jags.a)  Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2a.txt", fit using jags, 3 chains, each with 106000 iterations (first 6000 discarded), n.thin = 100 n.sims = 3000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff SD.Site[1] 0.117 0.091 0.004 0.043 0.099 0.173 0.330 1.001 3000 SD.Site[2] 0.095 0.074 0.004 0.037 0.080 0.137 0.272 1.001 3000 alpha[1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 1 alpha[2] -0.781 0.240 -1.198 -0.906 -0.786 -0.667 -0.306 1.001 3000 alpha0 1.091 0.179 0.753 1.008 1.091 1.181 1.405 1.004 3000 sd.Resid 0.387 0.015 0.366 0.377 0.385 0.394 0.423 1.001 3000 sd.Season 0.557 0.156 0.229 0.472 0.556 0.641 0.848 1.005 3000 sd.Site 0.106 0.071 0.005 0.050 0.096 0.153 0.260 1.001 3000 sigma 0.402 0.053 0.312 0.364 0.397 0.433 0.518 1.001 2900 sigma.B 0.179 0.196 0.006 0.067 0.130 0.227 0.641 1.001 3000 deviance 34.963 3.506 29.762 32.500 34.456 36.737 43.734 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 = 6.1 and DIC = 41.1 DIC is an estimate of expected predictive error (lower deviance is better).  Full effects parameterization JAGS code print(curdies.r2jags.b)  Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2b.txt", fit using jags, 3 chains, each with 106000 iterations (first 6000 discarded), n.thin = 100 n.sims = 3000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff SD.Site[1] 0.117 0.091 0.003 0.043 0.099 0.174 0.329 1.001 3000 SD.Site[2] 0.093 0.074 0.003 0.034 0.077 0.137 0.275 1.001 3000 alpha[1] 1.089 0.173 0.729 1.004 1.090 1.181 1.422 1.002 3000 alpha[2] -0.782 0.256 -1.261 -0.912 -0.786 -0.660 -0.287 1.005 3000 sd.Resid 0.388 0.015 0.366 0.378 0.386 0.394 0.424 1.001 3000 sd.Season 0.558 0.165 0.219 0.467 0.556 0.645 0.899 1.001 3000 sd.Site 0.105 0.072 0.004 0.045 0.097 0.151 0.262 1.001 3000 sigma 0.401 0.052 0.314 0.364 0.397 0.434 0.520 1.001 3000 sigma.B 0.178 0.187 0.005 0.063 0.134 0.235 0.632 1.001 3000 deviance 35.048 3.356 30.013 32.736 34.495 36.882 43.176 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 = 5.6 and DIC = 40.7 DIC is an estimate of expected predictive error (lower deviance is better).  Hierarchical parameterization JAGS code print(curdies.r2jags.c)  Inference for Bugs model at "../downloads/BUGSscripts/ws9.2bS1.2c.txt", fit using jags, 3 chains, each with 106000 iterations (first 6000 discarded), n.thin = 100 n.sims = 3000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff beta[1] 1.089 0.172 0.755 1.000 1.090 1.180 1.395 1.002 1400 beta[2] -0.783 0.248 -1.219 -0.907 -0.790 -0.664 -0.324 1.006 1100 sd.Resid 0.387 0.015 0.366 0.378 0.386 0.393 0.423 1.001 3000 sd.Season 0.558 0.161 0.239 0.470 0.560 0.642 0.865 1.003 2200 sd.Site 0.107 0.071 0.005 0.047 0.098 0.154 0.262 1.003 1800 sigma 0.402 0.054 0.313 0.364 0.396 0.434 0.525 1.001 3000 sigma.site 0.176 0.183 0.008 0.063 0.131 0.232 0.603 1.001 3000 deviance 34.992 3.396 29.839 32.650 34.470 36.755 42.968 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 = 5.8 and DIC = 40.8 DIC is an estimate of expected predictive error (lower deviance is better).  Full matrix parameterization STAN code print(curdies.rstan.a)  Inference for Stan model: rstanString. 3 chains, each with iter=106000; warmup=6000; thin=100; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta[1] 1.09 0.00 0.18 0.77 1.00 1.09 1.18 1.43 2224 1 beta[2] -0.79 0.01 0.25 -1.27 -0.92 -0.79 -0.66 -0.32 2202 1 gamma[1] 0.01 0.00 0.18 -0.31 -0.05 0.01 0.08 0.37 2644 1 gamma[2] 0.08 0.00 0.19 -0.24 -0.01 0.04 0.15 0.47 2786 1 gamma[3] -0.10 0.00 0.19 -0.53 -0.18 -0.06 0.00 0.18 2588 1 gamma[4] 0.02 0.00 0.18 -0.30 -0.04 0.01 0.09 0.37 3000 1 gamma[5] -0.05 0.00 0.18 -0.43 -0.12 -0.03 0.02 0.27 3000 1 gamma[6] 0.04 0.00 0.18 -0.27 -0.04 0.02 0.10 0.42 3000 1 sigma 0.40 0.00 0.05 0.31 0.37 0.40 0.44 0.52 2934 1 sigma_Z 0.19 0.00 0.21 0.02 0.07 0.14 0.24 0.68 2542 1 sd_Resid 0.39 0.00 0.02 0.37 0.38 0.39 0.39 0.43 1517 1 lp__ 22.62 0.11 4.60 13.78 19.56 22.38 25.45 32.09 1605 1 Samples were drawn using NUTS(diag_e) at Sun Mar 8 15:31:25 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).  Hierarchical parameterization STAN code print(curdies.rstan.b)  Inference for Stan model: rstanString. 3 chains, each with iter=106000; warmup=6000; thin=100; post-warmup draws per chain=1000, total post-warmup draws=3000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta[1] 1.09 0.00 0.19 0.72 1.01 1.09 1.18 1.41 2882 1 beta[2] -0.78 0.00 0.27 -1.22 -0.92 -0.79 -0.66 -0.33 2879 1 gamma[1] 1.11 0.00 0.13 0.85 1.03 1.11 1.19 1.37 2668 1 gamma[2] 1.17 0.00 0.14 0.91 1.08 1.17 1.26 1.47 2491 1 gamma[3] 0.99 0.00 0.15 0.69 0.90 1.00 1.09 1.25 2716 1 gamma[4] 0.32 0.00 0.13 0.06 0.24 0.32 0.41 0.60 2799 1 gamma[5] 0.25 0.00 0.14 -0.04 0.16 0.25 0.34 0.50 2904 1 gamma[6] 0.34 0.00 0.13 0.09 0.26 0.34 0.42 0.61 2844 1 sigma 0.40 0.00 0.05 0.31 0.37 0.40 0.43 0.52 2690 1 sigma_Z 0.19 0.00 0.20 0.03 0.08 0.14 0.25 0.64 2545 1 sd_Resid 0.39 0.00 0.01 0.37 0.38 0.38 0.39 0.42 3000 1 lp__ 22.16 0.11 4.11 13.75 19.57 22.15 24.82 30.24 1372 1 Samples were drawn using NUTS(diag_e) at Sun Mar 8 15:32:40 2015. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).  BRM code summary(curdies.brm)   Family: gaussian (identity) Formula: S4DUGES ~ SEASON + (1 | SITE) Data: curdies (Number of observations: 36) Samples: 3 chains, each with n.iter = 2000; n.warmup = 500; n.thin = 2; total post-warmup samples = 2250 WAIC: 39.88 Random Effects: ~SITE (Number of levels: 6) Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sd(Intercept) 0.18 0.16 0.01 0.62 567 1 Fixed Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 1.09 0.17 0.74 1.45 479 1.01 SEASONSUMMER -0.79 0.22 -1.25 -0.34 713 1.00 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sigma(S4DUGES) 0.4 0.05 0.31 0.53 1726 1 Samples were drawn using NUTS(diag_e). For each parameter, Eff.Sample is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1).  7. How about that summary figure then? Matrix parameterization JAGS code preds <- c('beta[1]','beta[2]') coefs <- curdies.r2jags.c$BUGSoutput$sims.matrix[,preds] curdies$SEASON <- factor(curdies$SEASON, levels=c("WINTER","SUMMER")) newdata <- data.frame(SEASON=levels(curdies$SEASON))
Xmat <- model.matrix(~SEASON, newdata)
pred <- coefs %*% t(Xmat)
library(plyr)
newdata <- cbind(newdata, adply(pred, 2, function(x) {
data.frame(Median=median(x)^4, HPDinterval(as.mcmc(x))^4, HPDinterval(as.mcmc(x), p=0.68)^4)
}))
newdata

  SEASON X1   Median     lower  upper   lower.1 upper.1
1 WINTER  1 0.008805 1.235e-07 0.1335 0.0006134 0.03213
2 SUMMER  2 1.412530 3.190e-01 3.7407 0.8355808 2.23421

library(ggplot2)
p1 <- ggplot(newdata, aes(y=Median, x=SEASON)) +
geom_linerange(aes(ymin=lower, ymax=upper), show_guide=FALSE)+
geom_linerange(aes(ymin=lower.1, ymax=upper.1), size=2,show_guide=FALSE)+
geom_point(size=4, shape=21, fill="white")+
scale_y_sqrt('Number of Dugesia per stone')+
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'))

preds <- c('sd.Resid','sd.Site','sd.Season')
curdies.sd <- adply(curdies.r2jags.c$BUGSoutput$sims.matrix[,preds],2,function(x) {
data.frame(mean=mean(x), median=median(x), HPDinterval(as.mcmc(x)), HPDinterval(as.mcmc(x),p=0.68))
})
head(curdies.sd)

         X1   mean  median     lower  upper   lower.1 upper.1
1  sd.Resid 0.3874 0.38556 0.3641406 0.4154 0.3720475  0.3955
2   sd.Site 0.1065 0.09773 0.0001281 0.2356 0.0007675  0.1368
3 sd.Season 0.5578 0.55956 0.2393835 0.8659 0.4214354  0.6854

rownames(curdies.sd) <- c("Residuals", "Site", "Season")
curdies.sd$name <- factor(c("Residuals", "Site", "Season"), levels=c("Residuals", "Site", "Season")) curdies.sd$Perc <- curdies.sd$median/sum(curdies.sd$median)

p2<-ggplot(curdies.sd,aes(y=name, x=median))+
geom_vline(xintercept=0,linetype="dashed")+
geom_hline(xintercept=0)+
scale_x_continuous("Finite population \nvariance components (sd)")+
geom_errorbarh(aes(xmin=lower.1, xmax=upper.1), height=0, size=1)+
geom_errorbarh(aes(xmin=lower, xmax=upper), height=0, size=1.5)+
geom_point(size=3, shape=21, fill='white')+
geom_text(aes(label=sprintf("(%4.1f%%)",Perc),vjust=-1))+
theme_classic()+
theme(axis.title.y=element_blank(),
axis.text.y=element_text(size=rel(1.2),hjust=1))

library(gridExtra)
grid.arrange(p1,p2,ncol=2)


## References

[1] A. Gelman. "Prior distributions for variance parameters in hierarchical models". In: _Bayesian Analysis_ (2006), pp. 515-533.