# Why do model selection (AIC and LOO) outcomes differ between ML and bayesian approaches

I am interested in understanding whether my continuous data (dput code at bottom for reproducibility) are fit better by a linear model (Gaussian distribution) or a gamma distributed model.

I typically use the lme4 package in R (maximum likelihood), but have been toying with the idea of using rstanarm (bayesian) a bit more.

First the ML models:

library(lme4)
library(rstanarm)
library(bayesplot)
library(DHARMa)

lmm <- lmer(Area ~ dB.s + Temp.s + (1 | SITE), data = SPt)
glmm <- glmer(Area ~ dB.s + Temp.s + (1 | SITE),
family = Gamma(link = "log"), data = SPt)
plot(simulateResiduals(lmm))


plot(simulateResiduals(glmm))


Just looking at the residual plots, the lmm looks like a much better fit to the data, and the AIC output suggests the same (lower AIC = better):

AIC(lmm,glmm)

df      AIC
lmm   5 1286.038
glmm  5 1294.297


Now let's try the rstanarm package.

S_lmm <- stan_lmer(Area ~ dB.s + Temp.s + (1 | SITE), data = SPt)
S_glmm <- stan_glmer(Area ~ dB.s + Temp.s + (1 | SITE),
family = Gamma(link = "log"), data = SPt)
## All Rhat values of both models are 1.0, indicating good mixing of the chains.


following the loo vignette: https://cran.r-project.org/web/packages/loo/vignettes/loo2-example.html

looL <- loo(S_lmm, save_psis = TRUE)
looL

Computed from 4000 by 97 log-likelihood matrix

Estimate   SE
elpd_loo   -655.5  6.3
p_loo         6.8  0.9
looic      1311.0 12.5
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).

looGL <- loo(S_glmm, save_psis = TRUE)
looGL

Computed from 4000 by 97 log-likelihood matrix

Estimate   SE
elpd_loo   -646.6  6.5
p_loo         7.0  0.9
looic      1293.3 12.9
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).


All Pareto k estimates are good, so I think I am okay to go ahead and compare these two:

loo_compare(looL,looGL)
elpd_diff se_diff
S_glmm  0.0       0.0
S_lmm  -8.9       3.5


Perhaps I am misinterpreting this, but this looks like the stan version of the glmm (Gamma) model (instead of the lmm above) has the best fit to the data (higher number = better in this case). Still following the vignette (link above), comparing LOO-PIT values to generated samples:

yrep <- posterior_predict(S_lmm)
ppc_loo_pit_overlay(SPt$$Area, yrep, lw = weights(looL$$psis_object))


The model may be a bit underdispersed here, having some excessive zeros, but without a lot of experience with these plots, I am not sure how poor this is.

yrep <- posterior_predict(S_glmm)
ppc_loo_pit_overlay(SPt$$Area, yrep, lw = weights(looGL$$psis_object))


The lower end of the glmm looks a bit better (not sure about the hump in the middle).

Am I interpreting these outputs correctly, that the lmm is a better fit with the ML methods (lme4) and the glmm is a better fit with the bayesian methods (rstanarm)? If so, why is this the case?

Data for reproducible example:

SPt<-structure(list(SITE = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 6L, 6L, 6L,
6L, 7L, 4L, 7L, 7L, 7L, 7L, 12L, 12L, 12L, 12L, 12L, 13L, 13L,
13L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 10L, 10L, 15L, 14L, 14L, 14L,
14L, 9L, 9L, 9L, 3L, 2L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 12L, 8L,
8L, 8L), .Label = c("CU0", "CU1", "CU2", "CU3", "CU4", "CW1",
"CW2", "FI10", "FI2", "FI4", "FI5", "FI6", "FI7", "FI8", "FI9"
), class = "factor"), dB.s = c(-0.756084718341438, -0.912553673339041,
-0.756084718341438, -0.912553673339041, -0.549545697744603, 0.00122502384695639,
0.551995745438517, 0.551995745438517, 0.539478229038709, -0.975141255338082,
-0.649685828943069, 0.138917704244847, 0.138917704244847, 0.138917704244847,
0.138917704244847, 0.138917704244847, 0.658394634836886, 0.658394634836886,
-1.16916275953511, -0.57458073054422, -0.57458073054422, -0.680979619942589,
-1.22549158333425, -0.893777398739329, -0.0676213163519883, -0.0676213163519883,
-0.19279648035007, -0.19279648035007, -0.0676213163519883, -0.19279648035007,
-0.130208898351029, -0.167761447550454, 1.17787156542893, 1.17787156542893,
1.17787156542893, 1.04643764323094, 1.04643764323094, 0.545736987238613,
0.545736987238613, 0.332939208441874, 0.332939208441874, 0.345456724841681,
0.345456724841681, 0.345456724841681, 0.345456724841681, 0.345456724841681,
0.345456724841681, 0.345456724841681, 0.345456724841681, -0.0551037999521799,
0.282869142842641, 0.282869142842641, -0.0989151073515087, 0.5269607126389,
0.0512950894461898, 0.35797424124149, 0.35797424124149, 0.501925679839284,
0.501925679839284, 0.126400187845038, 0.126400187845038, 0.126400187845038,
-0.00503373435294734, -0.249125304149207, -0.305454127948344,
-0.305454127948344, -0.180278963950262, -0.267901578748919, -0.480699357545658,
-0.480699357545658, -0.480699357545658, -0.480699357545658, -0.230349029549495,
-0.230349029549495, 0.401785548640819, 0.589548294637941, 0.0262600566465731,
0.470631888839763, 0.401785548640819, 0.126400187845038, 0.126400187845038,
0.00748378204686102, 0.639618360237174, 0.639618360237174, -1.63231086632801,
-0.0363275253524678, -0.0363275253524678, -0.174020205750358,
-0.148985172950741, -1.11283393573597, -0.236607787749399, 0.320421692042066,
-0.205313996749878, -0.286677853348631, 0.239057835443313, 0.239057835443313,
0.239057835443313), Temp.s = c(-0.198220313123015, -0.557916086050254,
-0.198220313123015, -0.557916086050254, 0.0874204477309688, 0.20379260956037,
-0.145323875927833, -0.145323875927833, -0.430964636781816, -0.56849537348929,
-0.251116750318197, -0.642550385562545, -0.642550385562545, -0.642550385562545,
-0.642550385562545, -0.642550385562545, -0.642550385562545, -0.642550385562545,
-0.684867535318691, -0.684867535318691, -0.684867535318691, -0.6002332358064,
-0.557916086050254, -0.557916086050254, -0.621391810684472, -0.621391810684472,
-0.462702499098926, -0.462702499098926, -0.621391810684472, -0.462702499098926,
0.457695508097244, -0.621391810684472, -0.631971098123509, -0.631971098123509,
-0.631971098123509, -0.737763972513873, -0.737763972513873, -0.864715421782311,
-0.864715421782311, 0.6481226819999, 0.6481226819999, 0.6481226819999,
0.6481226819999, 0.6481226819999, 0.6481226819999, 0.6481226819999,
0.6481226819999, 0.6481226819999, 0.6481226819999, 0.552909095048572,
0.806811993585447, 0.806811993585447, 0.859708430780629, 0.859708430780629,
0.32016477138977, 0.679860544317009, 0.679860544317009, 1.02897702980521,
1.02897702980521, 0.425957645780135, -0.219378888001088, 0.425957645780135,
-0.219378888001088, 0.595226244804718, 0.383640496023989, 0.383640496023989,
0.394219783463026, 0.119158310048078, 0.0451032979748232, 0.0451032979748232,
0.0451032979748232, 0.0451032979748232, -0.388647487025671, -0.388647487025671,
0.468274795536281, -0.0924274387326505, -0.304013187513379, 0.0556825854138594,
-0.293433900074343, -0.0289517140984317, -0.0289517140984317,
0.50001265785339, 0.0133654356577138, 0.0133654356577138, -0.293433900074343,
1.00781845492714, 1.00781845492714, 0.690439831756046, -0.251116750318197,
-0.103006726171687, 0.256689046755552, -0.325171762391452, 0.595226244804718,
-0.177061738244942, -0.430964636781816, -0.430964636781816, -0.430964636781816
), Area = c(205.849913383358, 1062.818937407, 337.910550123144,
535.843159145604, 317.637329654352, 773.064126911591, 301.068546437375,
168.859119942524, 254.242114406942, 378.11692072456, 364.829384950198,
298.957412667524, 228.018365593874, 292.623569178805, 447.169350303542,
926.007918050677, 481.372105142623, 131.743944565349, 831.000265719896,
891.805163211596, 501.640404306523, 248.286664757774, 839.867646604103,
860.135945768002, 339.494010995323, 204.709821555389, 648.585573244797,
507.974247795241, 392.381604126125, 784.129823903378, 618.893752757189,
374.963534532148, 175.732839060179, 141.371669411541, 270.176968208722,
61.0647072041516, 235.619449019234, 323.584043319749, 291.382718620453,
241.706284785565, 312.588469032184, 400.356713791849, 828.595062384308,
203.418124319939, 588.852273007237, 235.423099478385, 417.046424764045,
292.168116783851, 499.513231920777, 530.929158456675, 71.4712328691678,
385.630498228147, 223.053078404875, 477.522083345649, 758.498276301086,
554.294753817749, 172.787595947439, 501.869426410969, 670.730031541421,
199.491133502952, 100.074727121756, 286.670329640069, 228.018365593874,
240.331837999619, 386.415896391545, 829.380460547705, 354.214571692249,
360.497756999429, 298.254952550181, 599.062449131404, 461.61777053685,
647.953484802895, 697.989552456803, 304.024487458499, 527.787565803085,
637.743308678728, 435.768432023848, 238.761041672824, 431.968125930617,
686.43799480937, 223.053078404875, 619.679150920587, 312.588469032184,
133.517687777566, 1122.1376259541, 202.436376615692, 584.92528219025,
582.765437240907, 334.426936204349, 461.8141200777, 497.20671386442,
144.411631542787, 387.631221509586, 390.164758905073, 454.769962490004,
333.160167506605, 576.379757473404)), row.names = c(2L, 7L, 8L,
11L, 27L, 31L, 45L, 46L, 50L, 56L, 57L, 63L, 66L, 68L, 70L, 72L,
74L, 77L, 78L, 86L, 88L, 91L, 92L, 96L, 98L, 99L, 100L, 102L,
104L, 105L, 108L, 110L, 112L, 113L, 116L, 117L, 119L, 123L, 127L,
128L, 129L, 131L, 132L, 133L, 134L, 135L, 136L, 138L, 139L, 141L,
143L, 144L, 145L, 146L, 147L, 150L, 151L, 153L, 154L, 157L, 163L,
166L, 168L, 173L, 177L, 178L, 181L, 185L, 187L, 188L, 189L, 190L,
192L, 193L, 198L, 199L, 201L, 208L, 212L, 213L, 214L, 216L, 222L,
223L, 226L, 229L, 230L, 231L, 232L, 235L, 241L, 250L, 260L, 277L,
278L, 279L, 280L), class = "data.frame")

• I'm not sure where the simulateResiduals come from, but it looks as if you are evaluating the Gamma mixed model under a Gaussian QQplot. If that is true, of course the fit is going to look bad. Where is this function found? May 29 '20 at 21:53
• I apologize for missing that. I have edited library(DHARMa) into the code above, where simulateResiduals comes from. As far as I know, it is designed to be compatible with different distribution families, but I cannot really say how it works: cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html May 30 '20 at 0:24
• Thanks. I'm a little dubious of comparing models with different likelihoods via AIC. Have you seen that done before? I've seen AIC for model selection, but it has always been within the same family, just different parameters for the model. May 30 '20 at 0:56
• I have seen it done many times in ecology literature. Here is a link asking if it is okay stats.stackexchange.com/questions/31768/… some of the answers within say that it is a myth that they need to be the same family, but I am not sure. May 30 '20 at 1:02

Something strikes me as particularly odd when comparing different likelihoods via AIC.

Suppose I observed $$x=2$$. The log likelihood for a gaussian, gamma, and poisson each with mean and variance 1 is -0.91, -1, and -1. Should I assume this observation came from a gaussian simply because of the likelihood, ignoring details about the data generating process? I don't buy that.

In my own opinion, the choice of family comes (partially) prior to modelling. Considering that you are modelling a necessarily non-negative quantity, the choice of Gaussian is suspicious. The areas are large, perhaps large enough to warrant making the gaussian approximation (as is sometimes done with height. The probability of negative height under this model is negligibly small), but the residual variance of the model is nearly 200.

That means that when dB.s=1 (whatever that means, but it happens), 0 is almost 1 standard deviation away and so unphysical areas are not so improbable. In fact, calling simulate on lmm results in negative areas. That means that drawing samples from the distribution learned by your model results in sampling negative areas, which is clearly not physical. From this alone I would opt for the gamma were it my only other choice of family since it is supported on the non-negative reals (much like area).

This doesn't answer your question per se, but I think it does address something important. The choice of the family, in my own opinion and by the arguments I present here, is not something that is chosen in a data driven fashion, and it probably isn't something you select based on comparing the same models in two different modelling frameworks. Have a think about what your modelling and what the assumptions you are making. That, in part, should help with family selection and it won't rely on measures of goodness of fit.

• Thank you for this thoughtful answer. I generally do what you say when dealing with counts, or binomial data, but with positive continuous data, I often struggle getting gamma models to fit (or they have wild residual plots). RE "I would opt for the gamma were it my only other choice of family " - Is there another family that comes to mind for positive, continuous data, other than gamma? May 30 '20 at 2:17
• What about log-normal? Try taking the log of your areas and using lmer May 30 '20 at 2:24
• That would help make the distribution tails look more Gaussian, but it doesn't help with your statement: "Considering that you are modelling a necessarily non-negative quantity, the choice of Gaussian is suspicious." since all values would still be positive. I guess it would mean, however, that I couldn't get negative predictions though, only approaching 0, which is certainly an improvement. Is it appropriate in this case to center the response around the mean to force it to be Gaussian? May 30 '20 at 2:30
• Actually, it would. The appraoch I describe says that the log of the outcome is normal, which means that the outcome would be log normal. This implies that the outcome (under new realizations) would be strictly positive thus avoiding the criticisms I made. Do not alter the outcome in anyway, especially by centering the outcome. That destroys information about the intercept. May 30 '20 at 2:36
• Taking the log will not affect information in the same way that centering the outcome will. I suggest that before you embark on this model, you familiarize yourself with log-linear models and their interpretation. You may also be interested to read on the differences between log-linear models and gamma models. The latter assumes the coefficient of variation is constant, where the former does not. May 30 '20 at 2:53