# Strong linear pattern in residuals of random intercept model (simulated unbalanced data)

I am working with an unbalanced data set where most subjects have only 1 observation, but many have >1. When I fit a linear mixed model to these data (random intercept), there is a strong pattern in the residuals vs. fitted. I was able to duplicate the problem with simulated data (below).

The linear pattern seems to mostly come from individuals with only 1 observation.

This is because the random effect estimate (i.e., ranef(model))for those individuals with only 1 observation is approximately equal to their residual (i.e., resid(model)).

So my question is, how do I assess the model fit for this mixed effects random intercept model via residual diagnostics given these type of unbalanced data?

### Example data:

I borrowed the code below to simulate data from here.

This code block generates a balanced dataset with varying intercepts and 1 explanatory covariate. In the next block, I sub-sample this to create an unbalanced dataset, fit the model, and plot residuals vs. fitted.

#This code allows for random inercepts too, but I just set these parameters close to 0.
library(lme4)
form<-as.formula(c("~env+(env|ID)"))
N.ind<-2000
N.obs<-5
simdat <-data.frame(ID=factor(rep(1:N.ind,each=N.obs)),env=runif(N.ind*N.obs,0,1))
beta<-c(100,0.8)
names(beta)<-c("(Intercept)","env")

V.ind0 <- 1
V.inde <- 0.01
V.err0 <- 1
COVind0.ind1e <- 0.01*(sqrt(V.ind0*V.inde))

vcov<-matrix(c(V.ind0,COVind0.ind1e,COVind0.ind1e,V.inde), 2, 2)

theta<-c((chol(vcov)/sqrt(V.err0))[1,1],
(chol(vcov)/sqrt(V.err0))[1,2],
(chol(vcov)/sqrt(V.err0))[2,2])
names(theta)<-c("ID.(Intercept)","ID.env.(Intercept)","ID.env")

#set.seed(25)
response<-simulate(form,newdata=simdat,family=gaussian,
newparams=list(theta = theta,beta = beta, sigma = sqrt(V.err0)))
simdat$resp<-as.vector(response[,1]) summary(simdat)  Sub-sample these data to create unbalanced dataset. require(dplyr) simdat2<- simdat %>% sample_frac(0.25) simdat2$ID<- droplevels(simdat2$ID) simdat2 %>% group_by(ID) %>% summarise(obs.per.ID=n()) %>% group_by(obs.per.ID) %>% summarise(n=n()) model2<-lmer(resp~env+(1|ID),data=simdat2) summary(model2) plot(model2)  ## 1 Answer I don't think there is any problem with the pattern you see in the residuals: it comes from the fact that the estimate of the random effects for the subjects with only 1 observation are confounded with the per-observation noise term. (See also these related questions: Q1, Q2). Moreover, you can see in the example you give that despite the large number of subjects with only 1 observation, the estimated value of the slope ($0.82\pm0.08$) is close (less than 1 standard error) from the "true" value used to generate the data (which was$0.8\$). So if the goal was to estimate the (fixed-effect) slope, then the mixed-model would have provided a good estimate. (Note also that you get the same value if you fit a simple linear model with no random effect, lm(resp~env,data=simdat2)).

To asses the model residuals, since in this example you only have one continuous predictor, you could just plot the residuals with respect to the predictor This plot suggests that the assumptions of the linear model are respected: for example the variance of the residuals remains constant with respect to the values of the predictor, and their distribution seems symmetrical with respect to 0.

• If there was >1 predictor, just check residuals against each one separately? – Dave M Jan 31 '18 at 21:36
• Yes, that would tell you whether there is a problem of correlation in the residuals, and whether their variance remains or not constant with respect to changes in the continuous predictors. – matteo Feb 1 '18 at 10:43
• On a related note, in the example data there is a very large number of subjects for which n.obs=1. In this case you should take into account that all these subjects do not provide information for the estimation of the variance of the random effects, which may not be very reliable. Instead, estimates of fixed-effects should be OK. Also note that the model assumes that subjects with n.obs=1 and n.obs>1 comes from the same population, and only differ in the number of times they are "measured" (which may not be true if there is a specific reason for which they are "measured" more than once). – matteo Feb 1 '18 at 10:51