plot() does not show all diagnostic plots for lme/lmer

tl;dr ?plot.merMod explains in quite a bit of detail how the plotting methods work for fits produced by [g]lmer ...

You can get at least the first three plots corresponding to plot.lm fairly easily:

fitted vs residual with smooth line added

plot(lmer_model, type=c("p","smooth"), col.line=1)

(it's harder to get the smooth and the zero line drawn in different colours)

scale-location plot

plot(lmer_model,
     sqrt(abs(resid(.)))~fitted(.),
     type=c("p","smooth"), col.line=1)

Q-Q plot

lattice::qqmath(lmer_model)

residuals vs leverage

plot(fm1, rstudent(.) ~ hatvalues(.))

(the Cook's distances can be computed via cooks.distance() but superimposing the contours of CD={0.5,1} isn't so easy ...)

historical note

The design and implementation of lme4 diagnostic plot methods differ from plot.lm, which is the canonical example in base R. Why? I don't know for sure, but this approach is derived from the nlme package, which predates R; the earliest version I could find is this page from the Wayback Machine (1998), which links to a copy of the user's guide for version 1.2, dated February 1995; that's three months before the first source-code release of R (via ftp) in June 1995.

• it uses lattice (derived from Trellis™ graphics) rather than base-R graphics
• although it doesn't automatically construct e.g. scale-location plots, it is more flexible. You can use formulas to show fitted or residual values vs parameters, facet, etc., e.g. plot(fm1,residuals(.)~Days|Subject)
• there are separate commands for plotting residuals etc. (plot) and Q-Q plots (qqnorm in nlme, qqmath in lme4)

I know that this is a 2-year-old question, but I was having the same issue (September/2022) and then I found Panel of Diagnostic Residual Plots and redres

resid_panel(mod1, smoother = TRUE, qqbands = TRUE)

which shows:

In R, plot is a generic function. This means that when you call plot, R examines the class of the object you have passed to the first argument and chooses the plotting method according to this class.

Let's take an example. Suppose I use the lm function to create a model. The resulting model object will have class "lm":

lm_model <- lm(Sepal.Length ~ Sepal.Width, data = iris)
class(lm_model)
#> [1] "lm"

That means that when I call plot(lm_model), R will see that I am calling plot on an object of class lm. Instead of trying to construct a basic xy plot as it would if I did plot(1:10), R now knows to call a plotting method that has been specifically written to plot objects of type "lm". In this case, it will dispatch the method stats:::plot.lm, which is a long function that takes the "lm" object and creates the 4 diagnostic plots.

Now let's see what we get when we create a model with lmer:

library(lme4)
lmer_model <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
class(lmer_model)
#> [1] "lmerMod"
#> attr(,"package")
#> [1] "lme4"

Our model is an object of type "lmerMod". When we call plot on this object, R looks up the correct method to plot an object of this class. Since it has a completely different structure from an object of class "lm", it wouldn't make sense to plot it with plot.lm, so the authors who created the lme4 package had to decide what the best way to plot an object of class "lmerMod" was. They wrote the method lme4:::plot.merMod, which draws the single plot you see when you call plot on your model.

Why is this? That's one for the authors to answer, but it seems the main reason is that they wanted a plot method that would cover GLMM, LMM and REML models. The diagnostic plots for lm don't make sense for all of these model types.

So the short answer is that there is no problem to "solve" as such; this is just not how "lmerMod" objects are plotted. If you have specific concerns about some aspects of your fit that can be answered by these diagnostic plots, you should examine these individually.

As long as we're adding answers, the performance package is now available:

library(lme4)
library(performance)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
check_model(fm1)