General code to fit ANOVA models with JAGS and ‘rjags’

R
JAGS
Bayesian_stat
R-bloggers
Author

Andrea Onofri

Published

September 16, 2026

There is a quote I really love from Marc Kéry’s book, “Introduction to WinBUGS for Ecologists” (2010, p. 11):

“WinBUGS helps free the modeler in you.”

And it’s so true! Once you deeply understand a statistical model in all its moving parts, translating it into BUGS code feels remarkably logical and natural.

The downside? More often than not, the resulting code turns out to be very problem-specific. If you want to reuse it for a slightly different experimental setup, you usually end up doing a lot of tedious copy-pasting and editing.

Take ANOVA models with all their different “flavors”: one-way, two-way with interaction, nested, randomized blocks, and so on. These models are the bread and butter of agricultural research and genotype trials. While coding them individually in BUGS isn’t overly difficult, transitioning from one model structure to another requires manual code tweaks. And as anyone who has worked with BUGS/JAGS knows, spotting a subtle typo in BUGS code can easily eat up your entire afternoon…

For example, consider a field experiment where seven wheat genotypes were compared in a field trial laid out as a randomized complete block design with three replicates. The data comes from statforbiology, the companion R package for this blog. Let’s load the package and data:

library(statforbiology)
dataset <- getAgroData("WinterWheat2002")
head(dataset)
##   Plot Block Genotype Yield
## 1   57     A COLOSSEO  4.31
## 2   61     B COLOSSEO  4.73
## 3   11     C COLOSSEO  5.64
## 4   60     A    CRESO  3.99
## 5   10     B    CRESO  4.82
## 6   42     C    CRESO  4.17

This is a classic scenario where we want to fit an ANOVA model with Yield as the response variable, and Block and Genotype as explanatory fixed factors. In R, that model would be straightforward to fit, using the lm() function, while the genotype means with confidence intervals can be obtained by the emmeans() function.

# Fitting an ANOVA model with R
library(emmeans)
mod <- lm(Yield ~ Genotype + Block, data = dataset)
mm <- emmeans(mod, ~Genotype)
mm
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   4.89 0.186 14     4.49     5.29
##  CRESO      4.33 0.186 14     3.93     4.73
##  DUILIO     4.24 0.186 14     3.84     4.64
##  GRAZIA     4.34 0.186 14     3.94     4.74
##  IRIDE      4.96 0.186 14     4.56     5.36
##  SANCARLO   4.50 0.186 14     4.10     4.90
##  SIMETO     3.34 0.186 14     2.94     3.74
##  SOLEX      4.79 0.186 14     4.39     5.19
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95

The above coding is rather ‘reusable’ in the sense that, every other linear model can be fitted by simply changing the listing of effects and their relationships, namely the equations Yield ~ Genotype + Block and ~ Genotype in the calls to lm() and emmeans(), respectively.

If we want to fit the same model in the Bayesian platform with BUGS or one of its dialects, we could also use a similar coding, based on listing the predictors and their relationships (see for example Kéry, 2010). However, such coding would be specific to the experiment at hand and would not be easy to reuse for other linear models with different predictors. The question is: “Can we write a general BUGS/JAGS code that works for all ANOVA models without needing major edits every single time?”

After some tinkering, I came up with a neat solution! Since it took me a little while to get everything working smoothly, I thought I’d share it here for anyone interested in fitting ANOVA models within a Bayesian framework. We will be using JAGS (Just Another Gibbs Sampler) along with the rjags package in R (Plummer, 2019), but the coding can be easily adapted to other BUGS dialects.

Linear models in matrix notation

Apart from the ‘conceptual’ formula above, every linear model can also be specified in matrix notation, as:

\[ Y = X \, \beta + \varepsilon\] where \(Y\) is the vector of observed responses, \(\beta\) is the vector of estimated parameters and \(\epsilon\) is the vector of residuals, that are assumed to be independent, gaussian and homoscedastic. Consequently, we could as well say that the observed responses in themselves are sampled from a gaussian distribution, with expectation equal to \(X \beta\) and standard deviation equal to \(\sigma\).

\[ Y \sim N(X \beta, \sigma^2 I)\] Does this sound unfamiliar? Keep in mind that matrix \(X\) represents the ‘translation’ of our conceptual formula (‘~ Genotype + Block’) into a form that is understood by the fitting algorithm. If we change the dataset, we need to supply a different \(X\) matrix, but the specification of the model does not change, as it is completely general.

Likewise, the calculation of marginal means can be obtained via matrix multiplication, as:

\[M = K \beta \] where \(K\) is the so-called ‘linear combination matrix’. Also in this case, the expression is general and we can change the specification of what means we want to calculate by supplying a different \(K\) matrix.

For any fitted model, both the matrices \(X\) and \(K\) can be easily obtained with R: look at the code below to grasp the equivalence between the two model specifications (‘conceptual’ formula and in matrix form). Please, note the use of the %*% operator for matrix multiplication.

X <- model.matrix(mod)
# X <- model.matrix( ~ Block + Genotype, data = dataset) # equivalent
K <- mm@linfct
beta <- coef(mod)

# Fitted data
fitted(mod)
##        1        2        3        4        5        6        7        8 
## 4.675000 4.883750 5.121250 4.108333 4.317083 4.554583 4.021667 4.230417 
##        9       10       11       12       13       14       15       16 
## 4.467917 4.121667 4.330417 4.567917 4.745000 4.953750 5.191250 4.285000 
##       17       18       19       20       21       22       23       24 
## 4.493750 4.731250 3.121667 3.330417 3.567917 4.571667 4.780417 5.017917
as.numeric( X %*% beta )
##  [1] 4.675000 4.883750 5.121250 4.108333 4.317083 4.554583 4.021667 4.230417
##  [9] 4.467917 4.121667 4.330417 4.567917 4.745000 4.953750 5.191250 4.285000
## [17] 4.493750 4.731250 3.121667 3.330417 3.567917 4.571667 4.780417 5.017917

# Marginal means
mm
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   4.89 0.186 14     4.49     5.29
##  CRESO      4.33 0.186 14     3.93     4.73
##  DUILIO     4.24 0.186 14     3.84     4.64
##  GRAZIA     4.34 0.186 14     3.94     4.74
##  IRIDE      4.96 0.186 14     4.56     5.36
##  SANCARLO   4.50 0.186 14     4.10     4.90
##  SIMETO     3.34 0.186 14     2.94     3.74
##  SOLEX      4.79 0.186 14     4.39     5.19
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95
as.numeric (K %*% beta)
## [1] 4.893333 4.326667 4.240000 4.340000 4.963333 4.503333 3.340000 4.790000

Specification of a JAGS model

In JAGS (and other BUGS dialects), we can code any linear model using matrix notation, as long as we feed it the appropriate design matrix \(X\). Likewise, we can calculate means and, if necessary, pairwise differences by feeding it the appropriate linear combination matrix \(K\).

The main problem, for me, was to figure out a way to multiply matrices in JAGS. Lately, I discovered that it is necessary to use a for() loop and the inprod() function, which, for each ith observation, sums the products of all element in the ith row of \(X\) by the corresponding elements in the vector of estimated parameters \(\beta\), to obtain the expected value. Let’s see the full JAGS coding in the box below.

# Coding a JAGS model
modelSpec <- "
data {
n <- length(Y)
np <- dim(X)
nk <- dim(K)
}

model {
# Model
for (i in 1:n) {
   expected[i] <- inprod(X[i,], beta)
   Y[i] ~ dnorm(expected[i], tau)
  }

# Priors
beta[1] ~ dunif(0, 1000000)
for (i in 2:np[2]){
  beta[i] ~ dnorm(0, 0.000001)
  }
sigma ~ dunif(0, 100)

# Derived quantities (model specific)
tau <- 1 / ( sigma * sigma)

# Contrasts of interest
for (i in 1:nk[1]) {
   mu[i] <- inprod(K[i,], beta)
  }
}"
writeLines(modelSpec, con="ModelAOV.txt")

Let’s break down what’s going on here. The code consists of two main blocks enclosed in curly braces:

  1. A data block
  2. A model block

In the data block, we dynamically compute three dimensions: the number of observations (\(n\)), the number of parameters (np), and the number of contrasts (nk). These are derived automatically by counting the elements in \(Y\), the columns of the design matrix \(X\), and the rows of the linear combination matrix \(K\).

In the model block, we have three key components:

  1. Model Specification: We combine deterministic and stochastic statements. For each observation \(i\), the expected value is calculated by taking the inner product of row \(i\) of \(X\) and the vector of parameters \(\beta\), using the statement inprod(X[i,], beta). Then, we specify that Y[i] follows a normal distribution (dnorm) centered at expected[i] with precision tau (remember that BUGS/JAGS uses precision \(\tau = 1/\sigma^2\), rather than standard deviation \(\sigma\)).
  2. Priors: These express our expectations about model parameters before observing the data. Here, we use very vague priors: a wide uniform prior for the intercept (beta[1] ~ dunif(0, 1000000)), normal priors with mean 0 and near-zero precision for all other parameters in \(\beta\), and a uniform prior from 0 to 100 for the residual standard deviation (\(\sigma\)).
  3. Derived Quantities & Linear combinations of model parameters: We convert \(\sigma\) into precision \(\tau\), and we compute treatment means by multiplying the linear combination matrix \(K\) with the vector of estimated parameters \(\beta\), using inprod(K[i,], beta).

We save the string containing the JAGS coding into an external text file (‘ModelAOV.txt’) using writeLines().

Notice that this JAGS script is completely generic! It expects only three inputs from R: the response vector \(Y\), the design matrix \(X\), and the linear combination matrix \(K\).

Fitting the JAGS model in R

Now we can fit our model using the package ‘rjags’ (Plummer, 2019). Here are the steps:

  1. Extract the response vector Y;
  2. Build the design matrix X using model.matrix();
  3. Retrieve the linear combination matrix from the ‘emmeans’ object
  4. Pass everything to JAGS and run the MCMC sampler.
# Same as above
Y <- dataset$Yield
X <- model.matrix(mod)
K <- mm@linfct

Now, let’s launch the sampler.

library(rjags)

# Create lists
dataList <- list(Y = Y, X = X, K = K)
initList <- list(beta = beta, sigma = summary(mod)$sigma)

# Start sampler
mcmc <- jags.model("modelAOV.txt",
                   data = dataList, inits = initList,
                   n.chains = 4, n.adapt = 100)
## Compiling data graph
##    Resolving undeclared variables
##    Allocating nodes
##    Initializing
##    Reading data back into data table
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 24
##    Unobserved stochastic nodes: 11
##    Total graph size: 432
## 
## Initializing model

# Get samples
res <- coda.samples(mcmc, variable.names = c("beta", "sigma", "mu"),
                    n.iter = 1000)

out <- summary(window(res, start = 110))
res <- cbind(out$statistics[,1:2], out$quantiles[,c(1,5)])
res
##                 Mean         SD        2.5%       97.5%
## beta[1]   4.65398002 0.23929862  4.15701444  5.12817601
## beta[2]  -0.54962271 0.29848707 -1.13750071  0.05184453
## beta[3]  -0.63964792 0.30157657 -1.22861887 -0.02035331
## beta[4]  -0.54389958 0.30166191 -1.15118418  0.04889763
## beta[5]   0.08642982 0.29752857 -0.49133032  0.69177365
## beta[6]  -0.37144196 0.29525048 -0.95965461  0.22869691
## beta[7]  -1.53600068 0.30401159 -2.11513501 -0.93889532
## beta[8]  -0.08611270 0.30175883 -0.68332678  0.50944385
## beta[9]   0.21724170 0.18506336 -0.14352118  0.57796623
## beta[10]  0.45596275 0.18176129  0.08564017  0.80975117
## mu[1]     4.87838150 0.21174486  4.45695902  5.29274595
## mu[2]     4.32875879 0.20868876  3.92488597  4.74520264
## mu[3]     4.23873358 0.21225231  3.80199854  4.66074092
## mu[4]     4.33448192 0.21428484  3.90241528  4.76105391
## mu[5]     4.96481133 0.21409207  4.52864771  5.38338787
## mu[6]     4.50693954 0.20929064  4.09016766  4.93291796
## mu[7]     3.34238083 0.21750265  2.91838315  3.76680942
## mu[8]     4.79226880 0.21188674  4.39012645  5.22396670
## sigma     0.35795765 0.08036836  0.24191568  0.55593441

From the posterior distributions, we obtain posterior means, standard deviations, and 95% credible intervals. Because we used vague priors, the numerical results are almost identical to those from a standard frequentist ANOVA.

Reusing the EXACT same code for a multi-environment experiment

Here is where the magic happens! What if our genotype evaluation trial was conducted across multiple years?

In that case, we would want an ANOVA model with Year, Block (nested within Year), Genotype, and the Year by Genotype interaction.

Do we need to write a new JAGS script? Not at all! The ‘ModelAOV.txt’ file remains 100% untouched. All we need to do is update \(Y\), \(X\), and \(K\) in R:

# Loading the data
library(dplyr)
dataset <- getAgroData("WinterWheat")
dataset <- dataset |>
  mutate(across(c(Block, Year, Genotype), .fns = factor))
head(dataset)
##   Plot Block Genotype Yield Year
## 1    2     1 COLOSSEO  6.73 1996
## 2  110     2 COLOSSEO  6.96 1996
## 3  181     3 COLOSSEO  5.35 1996
## 4    2     1 COLOSSEO  6.26 1997
## 5  110     2 COLOSSEO  7.01 1997
## 6  181     3 COLOSSEO  6.11 1997

# Fit the model in R and calculate means
mod <- lm(Yield ~ Genotype*Year +  Block:Year, data = dataset)
mm <- emmeans(mod, ~Genotype|Year)
mm
## Year = 1996:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   6.35 0.223 98     5.90     6.79
##  CRESO      5.60 0.223 98     5.16     6.04
##  DUILIO     5.64 0.223 98     5.20     6.08
##  GRAZIA     6.27 0.223 98     5.82     6.71
##  IRIDE      6.04 0.223 98     5.60     6.48
##  SANCARLO   5.70 0.223 98     5.26     6.14
##  SIMETO     5.08 0.223 98     4.63     5.52
##  SOLEX      6.14 0.223 98     5.70     6.59
## 
## Year = 1997:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   6.46 0.223 98     6.02     6.90
##  CRESO      6.09 0.223 98     5.65     6.53
##  DUILIO     8.06 0.223 98     7.61     8.50
##  GRAZIA     6.73 0.223 98     6.29     7.18
##  IRIDE      7.72 0.223 98     7.28     8.17
##  SANCARLO   6.77 0.223 98     6.33     7.22
##  SIMETO     7.19 0.223 98     6.75     7.64
##  SOLEX      6.40 0.223 98     5.95     6.84
## 
## Year = 1998:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   6.70 0.223 98     6.25     7.14
##  CRESO      6.13 0.223 98     5.69     6.57
##  DUILIO     7.15 0.223 98     6.71     7.60
##  GRAZIA     6.35 0.223 98     5.91     6.80
##  IRIDE      6.39 0.223 98     5.95     6.84
##  SANCARLO   6.81 0.223 98     6.37     7.26
##  SIMETO     6.43 0.223 98     5.99     6.88
##  SOLEX      6.44 0.223 98     6.00     6.88
## 
## Year = 1999:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   6.98 0.223 98     6.54     7.43
##  CRESO      7.13 0.223 98     6.69     7.57
##  DUILIO     7.99 0.223 98     7.55     8.44
##  GRAZIA     6.84 0.223 98     6.40     7.28
##  IRIDE      7.99 0.223 98     7.55     8.43
##  SANCARLO   7.41 0.223 98     6.97     7.86
##  SIMETO     7.07 0.223 98     6.63     7.51
##  SOLEX      6.87 0.223 98     6.43     7.32
## 
## Year = 2000:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   6.44 0.223 98     6.00     6.88
##  CRESO      6.09 0.223 98     5.64     6.53
##  DUILIO     5.18 0.223 98     4.74     5.62
##  GRAZIA     4.75 0.223 98     4.31     5.19
##  IRIDE      6.05 0.223 98     5.61     6.49
##  SANCARLO   5.67 0.223 98     5.22     6.11
##  SIMETO     4.82 0.223 98     4.38     5.26
##  SOLEX      5.45 0.223 98     5.01     5.89
## 
## Year = 2001:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   7.07 0.223 98     6.63     7.52
##  CRESO      6.45 0.223 98     6.01     6.90
##  DUILIO     7.88 0.223 98     7.44     8.32
##  GRAZIA     7.29 0.223 98     6.85     7.74
##  IRIDE      7.71 0.223 98     7.26     8.15
##  SANCARLO   6.67 0.223 98     6.23     7.12
##  SIMETO     7.55 0.223 98     7.11     8.00
##  SOLEX      7.52 0.223 98     7.08     7.97
## 
## Year = 2002:
##  Genotype emmean    SE df lower.CL upper.CL
##  COLOSSEO   4.89 0.223 98     4.45     5.34
##  CRESO      4.33 0.223 98     3.88     4.77
##  DUILIO     4.24 0.223 98     3.80     4.68
##  GRAZIA     4.34 0.223 98     3.90     4.78
##  IRIDE      4.96 0.223 98     4.52     5.41
##  SANCARLO   4.50 0.223 98     4.06     4.95
##  SIMETO     3.34 0.223 98     2.90     3.78
##  SOLEX      4.79 0.223 98     4.35     5.23
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95

# Create input matrices
Y <- dataset$Yield
X <- model.matrix(mod)
K <- mm@linfct

# Create lists
dataList <- list(Y = Y, X = X, K = K)
initList <- list(beta = coef(mod), sigma = summary(mod)$sigma)

# Start sampler
mcmc <- jags.model("modelAOV.txt",
                   data = dataList, inits = initList,
                   n.chains = 4, n.adapt = 100)
## Compiling data graph
##    Resolving undeclared variables
##    Allocating nodes
##    Initializing
##    Reading data back into data table
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 168
##    Unobserved stochastic nodes: 71
##    Total graph size: 16380
## 
## Initializing model

# Get samples
res <- coda.samples(mcmc, variable.names = c("beta", "sigma", "mu"),
                    n.iter = 1000)

out <- summary(window(res, start = 110))
res <- cbind(out$statistics[,1:2], '50%'=out$quantiles[,3],
             out$quantiles[,c(1, 5)])
head(res)
##                Mean        SD          50%       2.5%       97.5%
## beta[1]  6.38802956 0.2237061  6.382177443  5.9553187  6.86054080
## beta[2] -0.66857825 0.2907187 -0.669801108 -1.2391331 -0.10177751
## beta[3] -0.65462205 0.2865742 -0.654855764 -1.2254214 -0.10066396
## beta[4] -0.01499206 0.2962239 -0.005706012 -0.5991771  0.59661314
## beta[5] -0.23425348 0.2992963 -0.238321565 -0.8121839  0.35368330
## beta[6] -0.56729738 0.2982795 -0.552967278 -1.1877757 -0.01023127
#....
tail(res)
##            Mean         SD       50%      2.5%     97.5%
## mu[52] 4.340052 0.22826887 4.3376727 3.9017644 4.7996293
## mu[53] 4.966680 0.22564790 4.9584444 4.5317184 5.4221167
## mu[54] 4.504849 0.22199908 4.5048840 4.0645888 4.9352030
## mu[55] 3.338719 0.22477808 3.3395523 2.8946365 3.7820630
## mu[56] 4.790833 0.22964059 4.7935006 4.3262642 5.2347674
## sigma  0.391495 0.02801552 0.3902889 0.3414186 0.4498655

Discovering how to leverage the inprod() function alongside R’s model.matrix() was a major lightbulb moment for me. This approach makes Bayesian model specification extremely flexible and opens the door to fitting mixed models effortlessly—something I’ll cover in future posts!

Thanks for reading, and happy coding!

(P.S. If you enjoyed this post, check out my book “Field Research Methods in Agriculture: An Introduction with R” via the link below!)

Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: andrea.onofri@unipg.it

Book cover

This post was originally published on 2020-12-23

References

  1. Kery, M., 2010. Introduction to WinBUGS for ecologists. A Bayesian approach to regression, ANOVA, mixed models and related analyses. Academic Press, Burlington, MA (USA).
  2. Plummer M. (2019). rjags: Bayesian Graphical Models using MCMC. R package version 4-10. https://CRAN.R-project.org/package=rjags