# Unexpected residuals plot of mixed linear model using lmer (lme4 package) in R

I have conducted an experiment with multiple (categorical) conditions per subject, and multiple subject measurements.

My data-frame in short: A subject has one property, is_frisian which is either 0 or 1 depending on the subject. And it is tested for two conditions, person and condition. The measurement variable is error, which is either 0 or 1.

My mixed linear model in R is:

> model <- lmer(error~is_frisian*condition*person+(1|subject_id), data=output)


However, the residuals plot of this model gives an unexpected (?) result.

I was taught that this plot should show randomly scattered points, and they should be normal distributed. When plotting the density of the fitted and the residuals, it shows a reasonable normal distribution. The lines you can see in the graph, however, how is this to be explained? And is this okay?

The only thing I could come up with is that the graph has two lines due to the categorical variables. The output variable error is either 0 or 1. But I do not have that much knowledge of the underlying system to confirm this. And then again, the lines also seem to have a low negative slope, is this then perhaps a problem?

UPDATE:

> model <- glmer(error~is_frisian*condition*person + (1|subject_id), data=output, family='binomial')
> binnedplot(fitted(model),resid(model))


Gives the following result:

FINAL EDIT:

The density-plots have been omitted, they have nothing to do with satisfaction of assumptions in this case. For a list of assumptions on logistic regression (when using family=binomial), see here at statisticssolutions.com

• First of all: The new model is not a mixed model, it has no random intercept like the first model for subject_id; is this intended? The densityplots are useless to check the residuals in this case, I think. – COOLSerdash Jul 7 '13 at 22:24
• Consider a slightly simpler case with a single continuous predictor: If your response data are 0-1, then your data are two parallel horizontal lines (at y=0 and y=1, both of the form y=c). If you then fit a linear term (a + bx) and subtract it to obtain the residual, what you have is r = y - a - bx. What happens to the line y=c? It becomes r = (c-a) - bx ... a straight line, and for the two values of c, two parallel lines -- so you got what you should have expected. (Does it make much sense to fit such a model? I don't think so.) – Glen_b Jul 7 '13 at 23:18
• @COOLSerdash Sorry for the confusion, somehow it got left out, the plots shown include the subject_id. – Yeti Jul 7 '13 at 23:19
• @Glen_b I see, it is what COOLSerdash illustrated in his answer (he suggested to use logistic regression, by setting family=binomial), what would be your suggestion? – Yeti Jul 7 '13 at 23:26
• The assumptions of logistic regression? That would be a whole new question, not something I can answer in a 3 line comment. – Glen_b Jul 7 '13 at 23:31

Your residual structure is totally expected with this model specification and an indication of an ill-specified model. What you basically are trying to do is to fit a linear line through points that can only take values of 0 and 1 on the $y$-axis.

Let's look at a simple example with arbitrarily generated variables:

#-----------------------------------------------------------------------------
# Generate random data for logistic regression
#-----------------------------------------------------------------------------

set.seed(123)
x <- rnorm(1000)
z <- 1 + 2*x
pr <- 1/(1+exp(-z))
y <- rbinom(1000,1, pr)

#-----------------------------------------------------------------------------
# Plot the data
#-----------------------------------------------------------------------------

par(bg="white", cex=1.2)
plot(y~x, las=1, ylim=c(-0.1, 1.3))

#-----------------------------------------------------------------------------
# Fit a linear regression (nonsensical) and plot the fit
#-----------------------------------------------------------------------------

linear.mod <- lm(y~x)
segments(-2.32146, 0, 1.24196, 1, col="steelblue", lwd=2)
segments(1.24196, 1, 100, 28.71447, col="red", lwd=2)
segments(-100, -27.41153, -2.32146, 0, col="red", lwd=2)


As you can see, a linear line is fitted through the data. One problem of this is that the line predicts outcomes that are outside the interval $[0,1]$ (illustrated by the red lines outside that interval). Let's have a look at the residuals:

#-----------------------------------------------------------------------------
#-----------------------------------------------------------------------------

x.y0 <- sample(which(y==0), 50, replace=F)
x.y1 <- sample(which(y==1), 50, replace=F)

pre <- predict(linear.mod)

segments(x[x.y0], y[x.y0], x[x.y0], pre[x.y0], col="red", lwd=2)
points(x[x.y0], y[x.y0], pch=16, col="red", las=1)

segments(x[x.y1], y[x.y1], x[x.y1], pre[x.y1], col="blue", lwd=2)
points(x[x.y1], y[x.y1], pch=16, col="blue", las=1)


I randomly picked some values to show the pattern. The red and blue lines are depicting the residuals, which is the difference between the predicted value of the line and the actual observed value (red and blue dots). The blue lines correspond to the residuals where $y=1$ whereas the red residuals correspond to the situation where $y=0$. Because the outcome can only be either 0 or 1, the residuals are simply the distances between the regression line and either 0 or 1. The residuals take exactly the form that you see in your data:

#-----------------------------------------------------------------------------
# Plot the residuals
#-----------------------------------------------------------------------------

res.linear <- residuals(linear.mod, type="response")

par(bg="white", cex=1.2)
plot(predict(linear.mod)[y==0], res.linear[y==0], las=1,
xlab="Fitted values", ylab = "Residuals",
ylim = max(abs(res.linear))*c(-1,1), xlim=c(-0.4, 1.6), col="red")
points(predict(linear.mod)[y==1], res.linear[y==1], col="blue")
abline(h = 0, lty = 2)


The colors correspond to the residuals shown above: the blue dots are the residuals where $y=1$ and the red dots are the residuals where $y=0$. In normal linear regression, the residuals are assumed to be approximately normally distributed. But in this case, the residuals can hardly be normal. They are binomial.

We need a transformation that transformes the probability, which is bound within $[0,1]$ into a variable that ranges over $(-\infty, \infty)$. One such transformation is the logit (this is not the only possibility: we could also use probit or the complementary log-log function). Let's fit a logistic regression with a logit-link and again plot the binned residuals (explained on page 97 by Gelman and Hill (2007)). Plotting the raw residuals vs. fitted values are generally not useful after logistic regression:

#-----------------------------------------------------------------------------
# Fit a logistic regression
#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------
# Plot the binned residuals as recommended by Gelman and Hill (2007)
#-----------------------------------------------------------------------------

library(arm)
par(bg="white", cex=1.2, las=1)
binnedplot(predict(glm.fit), resid(glm.fit), cex.pts=1, col.int="black")


The residuals in logistic regression can be define -$~$as in linear regression$~$- as observed minus expected values: $$\text{residual}_{i}=y_{i}-\mathrm{E}(y_{i}|X_{i})=y_{i}-\text{logit}^{-1}(X_{i}\beta)$$ Because the data $y_{i}$ are discrete, so are the residuals. In the plot above, the residuals are binned by dividing the data into categories based on their fitted values, and are then plotted against the average residual versus the average fitted value for each category (bin). The lines indicate $\pm2$ standard-error bounds, within which one we would expect about 95% of the binned residuals to fall, under the assumption that the model is true.

So the remedy for your immediate problem is to fit a mixed effects logistic regression by typing:

model <- glmer(error~is_frisian*condition*person+(1|subject_id),
data=output, family="binomial")


For a good introduction to mixed effects logistic regression in R, see here. For a good overview of diagnostics in linear and generalized linear models, see here.

• Suddenly, all pieces fit together. Thank you, the links seem very nice as well! – Yeti Jul 7 '13 at 16:43
• I have read most relevant parts of your links. I may fail to see the obvious, but why exactly is such a 'distinct pattern' perfectly okay and expected? Also, I wouldn't call the density-plots really normally distributed in my case, in the contrary to your test-data as given in your answer. It seems rather skewed, and there seem to be multiple 'hills'. Is this a problem? – Yeti Jul 7 '13 at 19:30
• @Yeti Regarding the pattern of the residuals in logistic regression: I have updated my answer. The density plots of what are not normal? The independent variables? They doesn't need to be normal at all. – COOLSerdash Jul 7 '13 at 20:00
• Thanks for the update, the binned plot seems to make more sense. I made density plots of the fitted and residuals values of the model, I thought that made sense, although who am I to say that, I am far from being a statistician. Also, considering the binnedplot, there is like one outlier all the way at -19, the rest are in a curved line ranging from -2 to 0. It seems a bit weird, comparing with your included figure. – Yeti Jul 7 '13 at 20:29
• @Yeti It would be best if you could include those plots into your question by editing it. Just from your description, it is hard to judge what's going on. – COOLSerdash Jul 7 '13 at 20:35