I am trying out residual covariance structures in brms. I have fit unstr(), cosy(), and ar(p=1) models to a bunch of example datasets in the agridat package. I uniformly find that the unstr() models perform best in LOO cross-validation. This goes against my intuition because to me it seems like those models are really overparameterized. Their effective number of parameters is much much higher than for the cosy() and ar() models.
I know that the regularizing effect of priors means that overfitting is not as much of a problem in the Bayesian context. But I still find it hard to believe that the unstructured covariance models, which estimate a unique parameter for every off-diagonal element of the residual covariance matrix, are performing that much better. I was wondering if I am somehow not doing a legitimate comparison with loo_compare(). Or is my intuition that badly off? Thanks for any insight anyone can provide.
Example
Here is an example comparing models fit with brm() using loo_compare(). Note that it was necessary to include a slightly regularizing prior on the covariance matrix for the unstructured covariance model, as well as set init=0
. The result of loo_compare()
below shows that the unstructured covariance model is superior despite having an effective number of parameters p_loo
of almost 900. In addition the posterior predictive check plot looks worse for the unstructured covariance model (see images at bottom).
Code
library(brms)
library(agridat)
data(broadbalk.wheat)
broadbalk.wheat$year <- broadbalk.wheat$year - 1852
options(mc.cores=4, brms.backend='cmdstanr')
broadbalk_unstructured <- brm(grain ~ year + unstr(time = year, gr = plot),
prior = c(
prior(normal(0, 5), class = b),
prior(lkj_corr_cholesky(3), class = Lcortime)
),
init = 0,
data = broadbalk.wheat, seed = 1041)
broadbalk_compoundsymm <- brm(grain ~ year + cosy(time = year, gr = plot),
prior = c(
prior(normal(0, 5), class = b)
),
data = broadbalk.wheat, seed = 1042)
broadbalk_ar1 <- brm(grain ~ year + ar(time = year, gr = plot, p = 1),
prior = c(
prior(normal(0, 5), class = b)
),
data = broadbalk.wheat, seed = 1043)
broadbalk_unstructured <- add_criterion(broadbalk_unstructured, 'loo')
broadbalk_compoundsymm <- add_criterion(broadbalk_compoundsymm, 'loo')
broadbalk_ar1 <- add_criterion(broadbalk_ar1, 'loo')
print(loo_compare(broadbalk_unstructured, broadbalk_compoundsymm, broadbalk_ar1), simplify = FALSE)
Result
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
broadbalk_unstructured 0.0 0.0 408.7 9.0 889.3 8.0 -817.5 18.0
broadbalk_compoundsymm -1506.2 24.7 -1097.5 26.3 2.2 0.1 2195.0 52.5
broadbalk_ar1 -1655.4 24.8 -1246.7 25.2 3.9 0.2 2493.4 50.4