What is this class?

• Part 3 of a practical introduction to fitting Bayesian multilevel models in R and Stan
• Uses the brms package (Bayesian Regression Models using Stan)

How to follow the course

• Slides and text version of lessons are online
• Fill in code in the worksheet (replace ... with code)
• You can always copy and paste code from text version of lesson if you fall behind

What we’ve learned in Lessons 1 and 2

• Relationship between prior, likelihood, and posterior
• How to sample from the posterior with Hamiltonian Monte Carlo
• How to use posterior distributions from our models to make predictions and test hypotheses
• GLMMs with non-normal response distributions and random intercepts and slopes

What will we learn in this lesson?

• How to fit a linear mixed model with residuals that are correlated in time (repeated measures in time)
• How to fit a linear mixed model to a split plot design, and add a spatial additive component to the model

Space and time

Space and time relationship: Example

• We are studying the effect of farm-level management practices on soil carbon
• The closer together two farms are, the more likely they are to share environmental conditions
• Adjacent farms are not statistically independent even if they make independent management decisions
• This is an issue with observational data and with designed experiments

Residuals are independent … right?

• The models we’ve seen so far assume the residuals of every data point are independent of all the others
• Even the mixed models assume this is true after the random effects are accounted for
• Patterns of variation in space and time can cause this assumption to be wrong

Two ways to model data correlated in space and time

• Within the framework of a GLMM, we can model a residual covariance structure
  • Data points are modeled as being more closely correlated, the closer they are to each other in space or time
• We can include an additive term, otherwise known as a “smoother” or “spline”
  • Models the response variable as depending on a smooth function of either space or time
  • Our model is now a GAMM (generalized additive mixed model)

• Example 1: a GLMM with residual covariance structure to model autocorrelation in time
• Example 2: a GAMM with an additive term to model autocorrelation in space

Both ways work for space or time, and you can even combine them in the same model!

Repeated measures within subjects

• Many datasets have the same subjects (plots, animals) measured at different time points
• We need to account for measurements from the same individual at different times not being independent

Load packages

• Same as before except ape for spatial autocorrelation statistics

library(brms)
library(agridat)
library(tidyverse)
library(emmeans)
library(tidybayes)
library(ggdist)
library(bayestestR)
library(gt)
library(ape)

theme_set(theme_bw())
colorblind_friendly <- scale_color_manual(values = unname(palette.colors(5)[-1]), aesthetics = c('colour', 'fill'))

options(mc.cores = 4, brms.backend = 'cmdstanr', brms.file_refit = 'on_change')

• Relevel factors so control groups are first in the factor level ordering

data(diggle.cow)

diggle.cow <- diggle.cow %>%
  mutate(iron = relevel(iron, ref = 'NoIron'),
         infect = relevel(infect, ref = 'NonInfected'))

• Plot each individual animal’s weight gain trend as a thin line, colored by treatment combination
• Plot median trend for each treatment combination as a thick line (stat_summary())

ggplot(diggle.cow, aes(x = day, y = weight, color = iron, linetype = infect, group = animal)) +
  geom_line() +
  stat_summary(aes(group = interaction(iron, infect)), fun = 'median', geom = 'line', linewidth = 2, show.legend = FALSE) +
  colorblind_friendly +
  labs(y = 'weight (kg)') +
  theme(legend.title = element_blank())

Scaling the variables

• Very large or small values cause numerical issues with Monte Carlo sampling
• Standardization, or z-transformation, done using scale()
• Also helpful for setting priors because the units are always the same (standard deviations)

diggle.cow$dayafterstart <- diggle.cow$day - min(diggle.cow$day)
diggle.cow$weight_scaled <- scale(diggle.cow$weight)

Fitting the model

• Fixed part of model formula: weight_scaled ~ dayafterstart * iron * infect
  • Continuous time effect, two categorical treatment effects, and interactions
  • How is the slope of weight versus time affected by iron and infection?
• But each of the 26 cows has her own unique time trend
• Our model needs to account for non-independence of measurements of the same cow on different days

Residual covariance matrix

• Models we’ve fit so far assumed conditionally independent residuals (independent after accounting for fixed and random effects)
• Those models still had a residual variance-covariance matrix \(\textbf{R} = \textbf{I}\sigma^2\)
  • \(n \times n\) matrix, \(n\) is number of time points
  • \(\sigma^2\) on every element of the diagonal: all residuals have same variance
  • 0 everywhere else: zero covariance means no relationship between any two residuals

Covariance matrix for correlated residuals

• Our big step is to put nonzero values into \(\textbf{R}\) to model residual autocorrelation
• There are many ways to model residual autocorrelation, we will cover three of the simplest ones today
  • Unstructured covariance
  • Compound symmetry
  • First-order autoregressive

Unstructured covariance

• Every pair of time points is modeled as having its own unique correlation
• Every single off-diagonal element of \(\textbf{R}\) (symmetrical) has its own covariance parameter \(\rho_{ij}\) that we have to estimate
• Most flexible because it has the most parameters

\[\textbf{R}=\begin{pmatrix}1 & \rho_{12} & \rho_{13} & \rho_{14} \\ \rho_{12} & 1 & \rho_{23} & \rho_{24} \\ \rho_{13} & \rho_{23} & 1 & \rho_{34} \\ \rho_{14} & \rho_{24} & \rho_{34} & 1 \end{pmatrix}\sigma^2\]

Compound symmetry

• Simplest model with the fewest parameters, but less flexible
• Only a single covariance parameter \(\rho\) is estimated for all pairs of time points
• Two residuals from the same individual are correlated the same, no matter how far apart in time they are

\[\textbf{R}=\begin{pmatrix}1 & \rho & \rho & \rho \\ \rho & 1 & \rho & \rho \\ \rho & \rho & 1 & \rho \\ \rho & \rho & \rho & 1 \end{pmatrix}\sigma^2\]

First-order autoregressive, or AR(1)

• Correlation between two time points within same individual is higher, the closer in time they are
• Like compound symmetry, only a single covariance parameter \(\rho\)
• Pairs of times 1 unit apart have correlation \(\rho\), 2 time units apart have \(\rho^2\), then \(\rho^3\) . . .

\[\textbf{R}=\begin{pmatrix}1 & \rho & \rho^2 & \rho^3 \\ \rho & 1 & \rho & \rho^2 \\ \rho^2 & \rho & 1 & \rho \\ \rho^3 & \rho^2 & \rho & 1 \end{pmatrix}\sigma^2\]

Choosing an appropriate covariance structure

• We will fit all three as proof of concept
• 23 time points, so there are 253 pairs
• Unstructured covariance model will probably be overparameterized and cause model fitting issues
• Autoregressive or compound symmetry?

Model syntax for residual covariance structures

• brms syntax: name(time = ..., gr = ...)
  • name: name of covariance function: unstr, cosy, ar, etc.
  • time: name of time variable column
  • gr: name of subject or grouping variable column

Model syntax for residual covariance structures

• Unstructured covariance: unstr(time = dayafterstart, gr = animal)
• Compound symmetry covariance: cosy(time = dayafterstart, gr = animal)
• Autoregressive AR(1) covariance: ar(time = dayafterstart, gr = animal, p = 1)
  • Increase p to get more complex autoregressive structures with “deeper” time lags

Comparing the models: posterior predictive check plots

pp_check(diggle_unstructured)
pp_check(diggle_compoundsymm)
pp_check(diggle_ar1)

Cross-validation two different ways

• Leave-one-out (LOO) cross-validation does not work well with unstructured covariance model and gives misleading results
• K-fold cross-validation works but takes a lot of computation time (each model needs to be refit K times)

kfold(diggle_unstructured, diggle_compoundsymm, diggle_ar1, folds = 'grouped', group = 'animal', K = 5, compare = TRUE)

Model comparisons:
                   elpd_diff se_diff
diggle_ar1             0.0       0.0 
diggle_unstructured  -41.3      19.0 
diggle_compoundsymm -295.1      24.5 

Model comparison results

• AR(1) model performs best
• Compound symmetry performs very poorly
• Weight gain of an animal over time is an autoregressive process (cow’s weight at time \(t\) depends on its weight at time \(t-1\))
• Compound symmetry might be better for other processes that are subject to random environmental fluctuations at each time (e.g., annual yield trends)

Model comparison results

• Unstructured covariance model does much better than it would in a frequentist context
  • Bayesian priors “regularize” models (shrink parameters toward zero) so even complex models are not too overfit

Bayesian p-values

• No strong evidence of effect of iron supplementation (but only \(n=7\) per group)
• We can also contrast the contrasts if we want

describe_posterior(emt_diggle$contrasts, ci = NULL, test = c('pd', 'p_map'))

Presentation ideas: How to report this analysis in a manuscript?

• Description of statistical model for methods section (including model comparison results)
• Verbal description of results
• Figure
• Table

Results section

There was little evidence for an effect of iron supplementation on the rate of weight gain. In uninfected animals the posterior estimate of the difference in mean weight gain rate between animals receiving supplemental iron and those that did not was -0.046 kg/day (95% CI [-0.134, 0.043]; p_MAP = 0.55), while in infected animals it was -0.038 kg/day (95% CI [-0.091, 0.013]; p_MAP = 0.37).

Manuscript-quality table

• Rescale slopes from standard deviations/day to kg/day units (multiply by sd(weight))
• 66% and 95% CIs: if the posterior distribution is roughly normal, they are roughly mean ± 1 and 2 SD

weight_slope_summary <- gather_emmeans_draws(emt_diggle$emtrends) %>%
  mutate(.value = .value * sd(diggle.cow$weight, na.rm = TRUE)) %>%
  median_qi(.width = c(0.66, 0.95)) %>%
  pivot_wider(names_from = .width, values_from = c(.lower, .upper)) %>%
  select(-c(.point, .interval))

gt(weight_slope_summary) %>%
  fmt_number(decimals = 3) %>%
  cols_merge(c(.lower_0.66, .upper_0.66), pattern = '[{1}, {2}]') %>% 
  cols_merge(c(.lower_0.95, .upper_0.95), pattern = '[{1}, {2}]') %>%
  cols_label(infect = 'tuberculosis infection', .value = 'weight gain rate (kg/day)', .lower_0.66 = '66% CI', .lower_0.95 = '95% CI') %>%
  tab_options(column_labels.font.weight = 'bold')

Explicit space and time terms

• Instead of modeling residual covariance we can model space or time directly
• It is not just a “nuisance” variable; we can make inference about it
• Simple way: just add linear covariate for space and/or time
• But the effect is rarely linear, especially with geospatial data

Generalized additive mixed models (GAMMs)

• Model the response variable as a function of some combination of unknown smooth functions of the covariates: spline or smoother term
• Many different functions can be used
• There is a built-in penalty term that prevents overfitting

GAM models fit with different k

GAMMs

• We will only scratch the surface of this complex topic today
• See “further reading” section to learn more

Meet the example data

Visualize layout of experiment

• Main plots are laid out adjacent to each other in a long line
• Total area of study: 10 rows by 56 beds

data(durban.splitplot)

ggplot(durban.splitplot, aes(x = bed, y = row, fill = block)) +
  geom_tile() +
  theme(aspect.ratio = 10/56, legend.position = 'bottom') +
  ggtitle('Blocks')

ggplot(durban.splitplot, aes(x = bed, y = row, fill = fung)) +
  geom_tile() +
  theme(aspect.ratio = 10/56, legend.position = 'bottom') +
  ggtitle('Main plots (fungicide treatment)')

ggplot(durban.splitplot, aes(x = bed, y = row, fill = gen)) +
  geom_tile() +
  theme(aspect.ratio = 10/56, legend.position = 'none') +
  ggtitle('Subplots (genotypes)')

Examining spatial trends

• Can’t we just account for the spatial trends using the block random effect?
• Heat map of spatial pattern of yield note single-color palette
• Blocks and main plots will only account for environmental factors that vary from east to west

ggplot(durban.splitplot, aes(x = bed, y = row, fill = yield)) +
  geom_tile() +
  scale_fill_distiller(palette = 'Oranges', direction = 1) +
  theme(aspect.ratio = 10/56, legend.position = 'bottom')

Examining spatial trends

• For demo purposes, randomly shuffle subplots within each main plot and make another heat map

set.seed(125)

durban.splitplot_shuffle <- durban.splitplot %>%
  group_by(block, fung) %>%
  mutate(yield_shuffled = sample(yield))

ggplot(durban.splitplot_shuffle, aes(x = bed, y = row, fill = yield_shuffled)) +
  geom_tile() +
  scale_fill_distiller(palette = 'Oranges', direction = 1) +
  theme(aspect.ratio = 10/56, legend.position = 'bottom')

Checking residuals for spatial autocorrelation

• If we see spatial pattern in the residuals, our fixed and random effects did not fully account for spatial variation in yield
• residuals() gives us mean and 95% quantiles of residuals across all posterior draws
• Extract first column (mean), append to the data, and create a heat map

resid_durban_lmm <- residuals(durban_lmm)

durban.splitplot <- mutate(durban.splitplot, resid_nospatial = resid_durban_lmm[,1])

ggplot(durban.splitplot, (aes(x = bed, y = row, fill = resid_nospatial))) +
  geom_tile() +
  scale_fill_distiller(palette = 'Greens', direction = 1) +
  theme(aspect.ratio = 10/56, legend.position = 'bottom')

Moran’s I

• Test statistic to measure strength of spatial autocorrelation
• Correlation between the difference between any two data points and the spatial distance between them
• Uses inverse distance: the higher the correlation, the more likely two close-together data points are to have a similar value
• Ranges from -1 to 1
  • 0 means no spatial pattern
  • 1 means perfect spatial correlation
• We don’t want high values of Moran’s I for residuals

Calculation of Moran’s I

• Calculate inverse distance matrix of all pairs of spatial points in the dataset
• Set diagonal of that matrix (distance between a point and itself) to 0
• Moran.I() from ape package gives us the statistic

durban_inv_dist <- 1 / as.matrix(dist(cbind(durban.splitplot$bed, durban.splitplot$row)))
diag(durban_inv_dist) <- 0

Moran.I(durban.splitplot$resid_nospatial, durban_inv_dist)

Comparing the models

• Some warnings about problematic observations, few enough that we can overlook them

durban_lmm <- add_criterion(durban_lmm, 'loo')
durban_gamm_k10 <- add_criterion(durban_gamm_k10, 'loo')
durban_gamm_k25 <- add_criterion(durban_gamm_k25, 'loo')
durban_gamm_k50 <- add_criterion(durban_gamm_k50, 'loo')
loo_compare(durban_lmm, durban_gamm_k10, durban_gamm_k25, durban_gamm_k50)

Checking residuals of spatial additive model

• Heat map of best-performing model (k = 50)

durban.splitplot <- mutate(durban.splitplot, resid_spatial = residuals(durban_gamm_k50)[,1])

ggplot(durban.splitplot, (aes(x = bed, y = row, fill = resid_spatial))) +
  geom_tile() +
  scale_fill_distiller(palette = 'Greens', direction = 1) +
  theme(aspect.ratio = 10/56, legend.position = 'bottom')

• Calculate Moran’s I statistic

Moran.I(durban.splitplot$resid_spatial, durban_inv_dist)

Impact of spatial additive term

• Check whether inference is different for non-spatial and spatial models
• Estimate marginal mean for each combination of genotype and fungicide for each model

emm_durban_nonspatial <- emmeans(durban_lmm, ~ gen + fung)
emm_durban_spatial <- emmeans(durban_gamm_k50, ~ gen + fung)

• Combine both estimates into one dataset for plotting
• Error bars are 95% highest-density credible intervals
• fct_reorder() sorts genotypes in increasing order of mean

means_durban_combined <- bind_rows(
  tibble(model = 'nonspatial', as_tibble(emm_durban_nonspatial)),
  tibble(model = 'spatial', as_tibble(emm_durban_spatial))
)

means_durban_combined %>%
  mutate(gen = fct_reorder(gen, emmean)) %>%
  ggplot(aes(x = gen, y = emmean, ymin = lower.HPD, ymax = upper.HPD, group = interaction(gen, model), color = model)) +
    facet_wrap(~ fung, ncol = 1) +
    geom_errorbar(position = position_dodge(width = .5)) +
    geom_point(position = position_dodge(width = .5), size = 2) +
    colorblind_friendly +
    theme(legend.position = 'bottom', axis.text.x = element_text(angle = 90))

Estimating marginal means

• emmeans() to get model predictions averaged over random effects of block and main plot, and over spatial smoother
• Results still on z-transformed scale

emms_bygenotype <- emmeans(durban_gamm_k50, pairwise ~ fung | gen)
emms_overall <- emmeans(durban_gamm_k50, pairwise ~ fung)

Presentation ideas: How to report this analysis in a manuscript?

• Description of model for methods section
• Table

We’ve learned . . .

• How to fit a model that assumes residuals are correlated in time, and compare different covariance structures
• How to diagnose spatial autocorrelation in residuals, and fit a GAMM with a spatial additive term to deal with the problem

You’re really starting to get the hang of this Bayes thing!

• See text version of lesson for further reading and useful resources
• Please provide feedback to help me improve this course!
  • Send feedback to quentin.read@usda.gov