2
$\begingroup$

I have been learning how to draw the model prediction on a scatter plot, and noticed a bit counter-intuitive result. I would greatly help if you could kindly explain how I am mistaken here.

Let me explain my confusion using "grouseticks" data set embedded in lme4 package. I fitted GLMM and GLM to the same data set. The only difference between the two models is an inclusion/exclusion of the random effect, BROOD.

# GLMM fitting
# Poisson distribution with log-link, BROOD as the random effect
library(lme4)
fitm <- glmer(TICKS ~ cHEIGHT + (1|BROOD), family=poisson(link=log), data=grouseticks)

# GLM fitting
# Poisson distribution with log-link
fit <- glm(TICKS ~ cHEIGHT, family=poisson(link=log), data=grouseticks)

Then I calculated Nakagawa & Schielzeth's (2013) marginal and conditional R-squares (R2m and R2c), which concerns variance explained by fixed factors and variance explained by both fixed and random factors, respectively.

# R-squares for GLMM
library(MuMIn)
r.squaredGLMM(fitm)

The output is:

             R2m    R2c
delta     0.3055 0.9363
lognormal 0.3070 0.9409
trigamma  0.3038 0.9310

which indicates that the addition of the random effect "BROOD" substantially improves the predictive power of the model. Am I right?

OK, here's the thing. I tried to confirm the above result visually. What I have done is:

# scatter plot
plot(TICKS ~ cHEIGHT, data= grouseticks)

# prepare newdata for predict() function
nd <- data.frame(cHEIGHT=c(-60:60))
pr <- predict(fit, newdata=nd, re.form=NA, type="response")
par(new=T)
lines(nd$cHEIGHT, pr, lwd=2, col="black")

# obtain GLMM predictions and draw the prediction curve
prm <- predict(fitm, newdata=nd, re.form=NA, type="response")
par(new=T)
lines(nd$cHEIGHT, prm, lwd=2, col="red")

...Then I get this. plot

I don't see such a big improvement in fitting from GLM (black) to GLMM (red). Rather, it seems to me that GLM predict the no. of TICKS better than GLMM when cHEIGHT is low.

I suspected that somehow I might be using predict() function in a wrong way, and tried to use the estimated parameters from GLM and GLMM to predict the TICKS, i.e.

pr <- exp(1.5446 - 0.0231*nd$cHEIGHT)
prm <- exp(0.5684 - 0.0252*nd$cHEIGHT)

The results were completely the same with those obtained from predict(). Am I correctly performing the prediction? If not, please tell me how I am wrong, thanks.

I also calculated simple correlation coefficients between the observed TICKS and the model predictions. The predictions from GLM correlated slightly better with the observed data (r=0.3454) than GLMM (r=0.3442), although the difference was trivial. I guess this is not the standard way to compare the goodness-of-fit, but it is still counter-intuitive to me.

Thus, I don't see how the addition of random effect "BROOD" improved the model fit, indicated by Nakagawa & Schielzeth's (2013) R2m and R2c.

Thank you very much for your kind help!

$\endgroup$
3
$\begingroup$

I think, there is a misunderstanding regarding the meaning of the R^2 values. From the help of the package, you get the following information: Marginal R_GLMM² represents the variance explained by the fixed effects, and is defined as: R_GLMM(m)² = (σ_f²) / (σ_f² + σ_α² + σ_ε²)

Conditional R_GLMM² is interpreted as a variance explained by the entire model, including both fixed and random effects, and is calculated according to the equation:

R_GLMM(c)² = (σ_f² + σ_α²) / (σ_f² + σ_α² + σ_ε²)

That is, the left column simply reflects how much of the total variation is explained by the fixed effects. Conceptually, this is what you would refer to as R^2 in a linear regression context as well.

The right column gives you the estimate when you count the intercept variance as explained variance (which, strictly speaking, it is not because you could continue by adding level-2 predictors that explain differences in the intercept).

Putting this aside, you would not see the improvement in fit in the plot you created because you are doing an unfair comparison, comparing the predicted values of the GLMM to the raw data. Rather, you would need to "control" for the random effects, that is, you would need to subtract the individual random intercept estimates from the raw values and then look at the fit of these "corrected" values.

Some illustrations of this have already been provided elsewhere in the context of the LMM, and they apply here as well: What would be an illustrative picture for linear mixed models?

Does this answer your question?

$\endgroup$
  • $\begingroup$ Thank you very much, MagunSimurgh! As for my misunderstanding about R^2: I wrote "R2m << R2c indicates that the addition of the random effect substantially improves the predictive power of the model." This sounds the same with what happens when an additional fixed effect improves the model fit to the raw data (I admit that's how I understood). But actually R2m << R2c means that the fit of the GLMM model to TICKS data adjusted for the random effect is better than the fit of the GLM model to the raw TICKS data. Is that what you suggested? $\endgroup$ – bbKZO Jul 30 '19 at 6:59
  • $\begingroup$ The link you suggested looks exactly what I was looking for. I will carefully read it later. $\endgroup$ – bbKZO Jul 30 '19 at 7:00
  • 1
    $\begingroup$ What you suggest is a fixed-effect approximation. That is, you dummy-code your level-2-categories and add the dummies as predictors. This is conceptually equivalent to the random intercept model. So far so good. However, the improvement in fit does not result from a better approximation of the fitting shape. Rather, the random effect reduces the noise by partialling out parts of it into the intercept variance. $\endgroup$ – StoryTeller0815 Jul 30 '19 at 7:06
3
$\begingroup$

I redrew the corrected figure, following the suggested link. The corrected R code is as follows.

RndEff <- unlist(ranef(fitm))
FixEff <- fixef(fitm)

for(i in 1:length(RndEff)) {
    prm <- exp(FixEff[1] + FixEff[2]*nd$cHEIGHT + RndEff[i])
    par(new=T)
    lines(nd$cHEIGHT, prm, lwd=.1, col="red")
}

Then I got this (the right graph is the corrected version).

enter image description here

Red curves represent the predictions for individual BROOD's.

$\endgroup$
  • 2
    $\begingroup$ Looks amazing! Do you now understand why this fits the data well? $\endgroup$ – StoryTeller0815 Jul 30 '19 at 13:18
  • 1
    $\begingroup$ Thanks, now I know how it works. The actual GLMM uses individual random intercept estimates as in the corrected version, while the previous version equals to the model with the average random intercept (=0). $\endgroup$ – bbKZO Jul 30 '19 at 13:53
  • 1
    $\begingroup$ Correct. Happy to hear it helped! $\endgroup$ – StoryTeller0815 Jul 30 '19 at 13:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.