Solved – The ‘best’ model selected with AICc have lower $R^2$ -square than the full/global model

aicmixed modelr-squared

I have used the R lme function (nlme package) to construct linear mixed models, with a single random effect (as a random intercept) and a varIdent variance structure on a fixed effect (that is a factor).

The data is from a commercial fishery and I am looking to see which are the most important variables that need to be reported for this fishery. I have 3 biological variables (total length, weight, proportion of females) and two operational variables (effort and season) (unfortunately I can't post my data as it's confidential and I haven't managed to reproduce my problem with dummy data).

I have used the AICc to select the 'best' models from a list of a priori models and I am now calculating the R-squared values from each of the 'best' models (using the code from https://github.com/jslefche/rsquared.glmer/blob/master/rsquaredglmm.R from the blog post http://jonlefcheck.net/2013/03/13/r2-for-linear-mixed-effects-models/ and taken from Nakagawa and Schielzeth (2013) r2 for GLMMs). I have 6 models with a AICc < 10 and I've used that as my cut off as it gives me a balanced number of models to the calculate the relative importance of each variable of interest using the sum Akaike's weights (the purpose of the modelling).

However, the R-squared values are lower for my 'best' models (Akaike weights between 0.5-0.8; marginal $R^2$ = 0.48) than the model with all my fixed effects (Akaike weights 0.00; marginal $R2^$ = 0.89).

I've re-run the model selection using BIC and I get the same problem. My question is: Are most 'best' models the models that have a lower $R^2$ value (and lower AICc) or have higher R2 values but higher AICc values?

Best Answer

Is your goal model parsimony or the predictive power of the model? If parsimony, then use AIC, if predictive power then $R^2$. Usually the answer is similar, but if you are comparing models with very similar $R^2$ or a number of low quality predictors the answers can be different. This is why in regular regression we tend to look at adjusted $R^2$ rather than just $R^2$, that is, because the adjusted value penalizes $R^2$ to adjust for the variance one might expect to be explained by chance if a predictor was not really effective at all. As the author says in the blog post "although I should note that [$R^2$ is] a poor tool for model selection, since it almost always favors the most complex models".

P.S. If you are interested in the fixed effects relative to the overall variance (regardless of nesting factor), then you probably want to be looking at the marginal $R^2$.

Related Solutions

Solved – How to fix an overdispersion in a Poisson GLMM with glmer function in R

Pulling out the answer from the discussion in the comments: If you create an obs_effect variable (observation-level random effect) with a unique value for each observation (say, 1:nrow(dat)), then you can incorporate overdispersion in the model by fitting

m2 <- glmer(dependent.var ~  cat.var1 * cat.var2 * contin.var + (1|obs_effect) + (1|group),
           data = dat, family = "poisson")

You also state in the comments that your problem is that the residual plot is triangle-shaped, which I interpret as the variability of the residuals increases with the predicted value. Depending on what kind of residuals you are plotting, this might mean nothing (the variability of observed - fitted should increase with fitted), or might mean that you have a problem other than overdispersion, which does not show up in a residual plot.

Reference: Harrison, X.A., 2014. Using observation-level random effects to model overdispersion in count data in ecology and evolution. PeerJ 2, e616. doi:10.7717/peerj.616

Solved – Calculating $R^2$ in mixed models using Nakagawa & Schielzeth’s (2013) R2glmm method

I am answering by pasting Douglas Bates's reply in the R-Sig-ME mailing list, on 17 Dec 2014 on the question of how to calculate an $R^2$ statistic for generalized linear mixed models, which I believe is required reading for anyone interested in such a thing. Bates is the original author of the lme4package for R and co-author of nlme, as well as co-author of a well-known book on mixed models, and CV will benefit from having the text in an answer, rather than just a link to it.

I must admit to getting a little twitchy when people speak of the "R2 for GLMMs". R2 for a linear model is well-defined and has many desirable properties. For other models one can define different quantities that reflect some but not all of these properties. But this is not calculating an R2 in the sense of obtaining a number having all the properties that the R2 for linear models does. Usually there are several different ways that such a quantity could be defined. Especially for GLMs and GLMMs before you can define "proportion of response variance explained" you first need to define what you mean by "response variance". The whole point of GLMs and GLMMs is that a simple sum of squares of deviations does not meaningfully reflect the variability in the response because the variance of an individual response depends on its mean.

Confusion about what constitutes R2 or degrees of freedom of any of the other quantities associated with linear models as applied to other models comes from confusing the formula with the concept. Although formulas are derived from models the derivation often involves quite sophisticated mathematics. To avoid a potentially confusing derivation and just "cut to the chase" it is easier to present the formulas. But the formula is not the concept. Generalizing a formula is not equivalent to generalizing the concept. And those formulas are almost never used in practice, especially for generalized linear models, analysis of variance and random effects. I have a "meta-theorem" that the only quantity actually calculated according to the formulas given in introductory texts is the sample mean.

It may seem that I am being a grumpy old man about this, and perhaps I am, but the danger is that people expect an "R2-like" quantity to have all the properties of an R2 for linear models. It can't. There is no way to generalize all the properties to a much more complicated model like a GLMM.

I was once on the committee reviewing a thesis proposal for Ph.D. candidacy. The proposal was to examine I think 9 different formulas that could be considered ways of computing an R2 for a nonlinear regression model to decide which one was "best". Of course, this would be done through a simulation study with only a couple of different models and only a few different sets of parameter values for each. My suggestion that this was an entirely meaningless exercise was not greeted warmly.

Best Answer

Related Solutions

Solved – How to fix an overdispersion in a Poisson GLMM with glmer function in R

Solved – Calculating $R^2$ in mixed models using Nakagawa & Schielzeth’s (2013) R2glmm method

Related Question