Exercise 1

  1. The conditional effect plots appear to have roughly similar inference compared to the model with autoregressive error covariance. The relative difference between the slopes of the iron and no-iron treatments is similar. The smoothers are very close to linear, so based on the plots I would argue that little is gained by fitting the GAMM. Your interpretation may vary!
  2. We use add_criterion() and loo_compare():
diggle_GAMM <- add_criterion(diggle_GAMM, 'loo')
loo_compare(diggle_GAMM, diggle_ar1)
##             elpd_diff se_diff
## diggle_ar1     0.0       0.0 
## diggle_GAMM -203.7      24.9

The AR1 model performs much better, providing additional evidence that there is not much benefit in introducing the added complexity of the GAMM in this case.

Exercise 2

  1. Here are the models, with default priors and model fitting options (you may experiment with varying the defaults). I’ve loaded them from prefit model objects.
broadbalk_cosy <- brm(grain ~ year + cosy(time = year, gr = plot), data = broadbalk.wheat, file = 'fits/broadbalk_cosy')
broadbalk_ar1 <- brm(grain ~ year + ar(time = year, gr = plot, p = 1), data = broadbalk.wheat, file = 'fits/broadbalk_ar1')
  1. Again we use add_criterion() and loo_compare().
broadbalk_cosy <- add_criterion(broadbalk_cosy, 'loo')
broadbalk_ar1 <- add_criterion(broadbalk_ar1, 'loo')
loo_compare(broadbalk_cosy, broadbalk_ar1)
##                elpd_diff se_diff
## broadbalk_cosy    0.0       0.0 
## broadbalk_ar1  -149.2      20.0

The compound symmetry covariance is a better fit than the autoregressive. This makes sense because the time trends are capturing yearly environmental fluctuations. Each year is more or less unique and the weather in one year is not highly correlated with the previous. Thus, the autoregressive structure that models years closer in time as more closely correlated is a poor fit.

  1. Use summary(broadbalk_cosy) or describe_posterior(broadbalk_cosy) to look at the parameter estimate for the slope and its 95% CI, or use conditional_effects(broadbalk_cosy) to see a plot. We can see that there is evidence for a declining trend in grain yield over time.
describe_posterior(broadbalk_cosy, test = 'p_map')
## Summary of Posterior Distribution
## 
## Parameter   |    Median |         95% CI | p (MAP) |  Rhat |     ESS
## --------------------------------------------------------------------
## (Intercept) |     19.75 | [16.87, 22.67] |  < .001 | 1.001 | 4246.00
## year        | -9.51e-03 | [-0.01, -0.01] |  < .001 | 1.001 | 4277.00

Exercise 3

  1. The heat map (fill color palette is optional) shows a strong trend for increased yield as you get closer to the upper right corner.
ggplot(burgueno.rowcol, aes(x = col, y = row, fill = yield)) + geom_tile() + scale_fill_distiller(palette = 'Oranges', direction = 1)

heat map of example yield data

  1. Here are the models with default priors and model fitting options (again you may experiment with changing those). I’ve loaded them from prefit model objects.
burgueno_nonspatial <- brm(yield ~ gen, data = burgueno.rowcol, file = 'fits/burgueno_nonspatial')
burgueno_spatial <- brm(yield ~ gen + s(col, row, k = 10, bs = 'gp'), data = burgueno.rowcol, file = 'fits/burgueno_spatial')
  1. The residuals of the model that ignores the spatial trend still have a similar spatial trend to the observed data. The residuals of the spatial model look random.
burgueno.rowcol$residuals_nonspatial <- residuals(burgueno_nonspatial)[,1]
burgueno.rowcol$residuals_spatial <- residuals(burgueno_spatial)[,1]

ggplot(burgueno.rowcol, aes(x = col, y = row, fill = residuals_nonspatial)) + geom_tile() + scale_fill_distiller(palette = 'Greens', direction = 1)

heat maps of residuals from two models fit to example yield data

ggplot(burgueno.rowcol, aes(x = col, y = row, fill = residuals_spatial)) + geom_tile() + scale_fill_distiller(palette = 'Greens', direction = 1)

heat maps of residuals from two models fit to example yield data

As a bonus, here is code to calculate Moran’s I for both the models. We see there is fairly strong spatial autocorrelation in the first model’s residuals but not the second.

burgueno_inv_dist <- 1 / as.matrix(dist(cbind(burgueno.rowcol$col, burgueno.rowcol$row)))
diag(burgueno_inv_dist) <- 0

Moran.I(burgueno.rowcol$residuals_nonspatial, burgueno_inv_dist)
## $observed
## [1] 0.1279599
## 
## $expected
## [1] -0.007874016
## 
## $sd
## [1] 0.009038702
## 
## $p.value
## [1] 0
Moran.I(burgueno.rowcol$residuals_spatial, burgueno_inv_dist)
## $observed
## [1] 0.009799946
## 
## $expected
## [1] -0.007874016
## 
## $sd
## [1] 0.009038264
## 
## $p.value
## [1] 0.05052883