Solved – Adding an observation level random term messes up residuals vs fitted plot. Why

glmmoverdispersion

I run a mixed effects generalized model for proportional data (response variable). I used binomial family and logit link function. I suffered from overdispersion so I added an observation level random effect to absorb part of that unexplained variability.

The problem is: now my residuals plot (resid vs. fitted) has a clear pattern. If I delete that new "observation level random term" the plot looks good again.

Why? Is that a common issue or it depends on my particular data? Is there a way I can solve it?

EDITED: I attache the "before-after" plots.

Best Answer

Thanks for updating your post, Charly. I played with some over-dispersed Poisson data to see the impact of adding an observation level effect in the glmer model on the plot of residual versus fitted values. Here is the R code:

# generate data like here: https://rpubs.com/INBOstats/OLRE

set.seed(324)
n.i <- 10
n.j <- 10
n.k <- 10
beta.0 <- 1
beta.1 <- 0.3
sigma.b <- 0.5
theta <- 5
dataset <- expand.grid(
X = seq_len(n.i),
b = seq_len(n.j),
Replicate = seq_len(n.k)
)
rf.b <- rnorm(n.j, mean = 0, sd = sigma.b)
dataset$eta <- beta.0 + beta.1 * dataset$X + rf.b[dataset$b]
dataset$mu <- exp(dataset$eta)
dataset$Y <- rnbinom(nrow(dataset), mu = dataset$mu, size = theta)
dataset$OLRE <- seq_len(nrow(dataset))


require(lme4)

m.1 <- glmer(Y ~ X + (1 | b), family=poisson(link="log"), data=dataset)

m.2 <- glmer(Y ~ X + (1 | b) + (1 | OLRE), family=poisson(link="log"),  
             data=dataset)

Note that model m.2 includes an observation level random effect to account for over-dispersion.

To diagnose the presence of over-dispersion in model m.1, we can use the command:

# check for over-dispersion:  
# values greater than 1.4 indicate over-dispersion

require(blmeco)

dispersion_glmer(m.1)

The value returned by dispersion_glmer is 2.204209, which is larger than the cut-off of 1.4 where we would start to suspect the presence of over-dispersion.

When applying dispersion_glmer to model m.2, we get a value of 1.023656:

dispersion_glmer(m.2)

Here is the R code for the plot of residuals (Pearson or deviance) versus fitted values:

par(mfrow=c(1,2))
plot(residuals(m.1, type="pearson") ~ fitted(m.1), col="darkgrey")
     abline(h=0, col="red")
plot(residuals(m.2, type="pearson") ~ fitted(m.2), col="darkgrey")
     abline(h=0, col="red")

par(mfrow=c(1,2))
plot(residuals(m.1, type="deviance") ~ fitted(m.1), col="darkgrey")
abline(h=0, col="red")
plot(residuals(m.2, type="deviance") ~ fitted(m.2), col="darkgrey")
abline(h=0, col="red")

As you can see, the Pearson residuals plot for the model m.2 (which includes an observation level random effect) looks horrendous compared to the plot for model m.1.

I am not showing the deviance residuals plot for m.2 as it looks about the same (that is, horrendous).

Here is the plot of fitted values versus observed response values for models m.1 and m.2:

par(mfrow=c(1,2))
plot(fitted(m.1) ~ dataset$Y, col="darkgrey", 
     xlim=c(0, 250), ylim=c(0, 250), 
     xlab="Y (response)", ylab="Fitted Values")
abline(a=0, b=1, col="red")
plot(fitted(m.2) ~ dataset$Y, col="darkgrey", 
     xlim=c(0, 250), ylim=c(0, 250), 
     xlab="Y (response)", ylab="Fitted Values")
abline(a=0, b=1, col="red")

The plot of fitted values versus actual response values seems to look better for model m.2.

We should check the summary corresponding to the two models:

summary(m.1)

summary(m.2)

As argued in https://rpubs.com/INBOstats/OLRE, large discrepancies between the fixed effect coefficients and especially the random effects variance for b would suggest that something may be off. (The extent of overdispersion present in the initial model would drive the extent of these discrepancies.)

Let's look at some diagnostic plots for the two models obtained with the Dharma package:

require(DHARMa)

fittedModel <- m.1
simulationOutput <- simulateResiduals(fittedModel = fittedModel)
plot(simulationOutput)

#----

fittedModel <- m.2
simulationOutput <- simulateResiduals(fittedModel = fittedModel)
plot(simulationOutput)

The diagnostic plots for model m.1 (especially the left panel) clearly shows overdispersion is an issue.

The diagnostic plot for model m.2 shows overdispersion is no longer an issue.

See https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html for more details on these types of plots.

Finally, let's do a posterior predictive check for the two models (i.e., plotting the fitted values obtained across simulated data sets constructed from each model over a histogram of the real response values Y), as explained at http://www.flutterbys.com.au/stats/ws/ws12.html:

range(dataset$Y)  # Actual response values Y range from 0 to 247

set.seed(1234567)
glmer.sim1 <- simulate(m.1, nsim = 1000)
glmer.sim2 <- simulate(m.2, nsim = 1000)

out <- matrix(NA, ncol = 2, nrow = 251)
cnt <- 0:250
for (i in 1:length(cnt)) {
 for (j in 1:2) {
     eval(parse(text = paste("out[i,", j, "] <- 
          mean(sapply(glmer.sim", j,",\nFUN = function(x) {\nsum(x == cnt  
      [i]) }))", sep = "")))
 }
}


plot(table(dataset$Y), ylab = "Frequency", xlab = "Y", lwd = 2, 
     col="darkgrey")
lines(x = 0:250, y = out[, 1], lwd = 2, lty = 2, col = "red")    
lines(x = 0:250, y = out[, 2], lwd = 2, lty = 2, col = "blue")

The resulting plot shows that both models are doing a good job at approximating the distribution of Y.

Of course, there are other predictive checks one could look at, including the centipede plot, which would show where the model with observation level random effect would fail (e.g., the model would tend to under-predict low values of Y): http://newprairiepress.org/cgi/viewcontent.cgi?article=1005&context=agstatconference.

This particular example shows that it is possible for the addition of an observation level random effect to worsen the appearance of the plot of residuals versus fitted values, while producing other diagnostic plots which look fine. I wonder if other people on this site may be able to add further insights into how one should proceed in this situation, other than to report what happens with each diagnostic plot when the correction for over-dispersion is used.

Related Solutions

Solved – How to assess the fit of a binomial GLMM fitted with lme4 (> 1.0)

Short answer since I don't have time for better: this is a challenging problem; binary data almost always requires some kind of binning or smoothing to assess goodness of fit. It was somewhat helpful to use fortify.lmerMod (from lme4, experimental) in conjunction with ggplot2 and particularly geom_smooth() to draw essentially the same residual-vs-fitted plot you have above, but with confidence intervals (I also narrowed the y limits a bit to zoom in on the (-5,5) region). That suggested some systematic variation that could be improved by tweaking the link function. (I also tried plotting residuals against the other predictors, but it wasn't too useful.)

I tried fitting the model with all 3-way interactions, but it wasn't much of an improvement either in deviance or in the shape of the smoothed residual curve.

Then I used this bit of brute force to try inverse-link functions of the form $(\mbox{logistic}(x))^\lambda$, for $\lambda$ ranging from 0.5 to 2.0:

## uses (fragile) internal C calls for speed; could use plogis(),
##  qlogis() for readability and stability instead
logitpower <- function(lambda) {
    L <- list(linkfun=function(mu)
              .Call(stats:::C_logit_link,mu^(1/lambda),PACKAGE="stats"),
              linkinv=function(eta)
              .Call(stats:::C_logit_linkinv,eta,PACKAGE="stats")^lambda,
              mu.eta=function(eta) {
                  mu <-  .Call(stats:::C_logit_linkinv,eta,PACKAGE="stats")
                  mu.eta <-  .Call(stats:::C_logit_mu_eta,eta,PACKAGE="stats")
                  lambda*mu^(lambda-1)*mu.eta
              },
              valideta = function(eta) TRUE ,
              name=paste0("logit-power(",lambda,")"))
    class(L) <- "link-glm"
    L
}

I found that a $\lambda$ of 0.75 was slightly better than the original model, although not significantly so -- I may have been overinterpreting the data.

Solved – Binomial GLMM: Model validation & ceiling effect

This looks fairly reasonable to me; I don't think the inference based on this model is likely to be far off. However, to take a more positive attitude, any deviation in your residuals also implies a chance to improve the model (i.e., there is further information that could be modeled).

Does the full model show the same deviations? That is, even though the variables you've discarded were non-significant, they might help address the (slight) pattern in the residuals.
You might be able to improve the model fit by modifying the link function (or equivalently transforming the predictor variables/linear predictor). In How to assess the fit of a binomial GLMM fitted with lme4 (> 1.0)? , where a similar pattern of residuals is discussed, I show how to construct a power-logit family of link functions that can be used for testing goodness-of-link and/or improving the model. (Existing goodness-of-link tests such as Pregibon's test use linearization and score tests to evaluate goodness of fit in an efficient way by comparing the existing fit to a family of link functions; the procedure at the linked question does the same thing in a much more brute-force way.) You might also find similar families of alternate link functions provided in the glmx package.

Best Answer

Related Solutions

Solved – How to assess the fit of a binomial GLMM fitted with lme4 (> 1.0)

Solved – Binomial GLMM: Model validation & ceiling effect

Related Question