Why are residuals correlated with prediction when using linear mixed model with low number of observations per group?

Question

I have a database composed of n = 1000 observations of a variable to explain (y), an explanatory variable (x). The observations are grouped according to a variable "group", and I have only few observations per "group" (around 2), while the number of different "groups" is high.

I fitted a linear mixed model (with R lmer function, from package lme4) explaining y as a function of x, using "group" as a random effect. I then checked if the residuals were normal and independant of the fitted values and of the random effects.

**At this stage I got a problem:

• the residuals show an upward-trend in relation to the fitted values
• the residuals show an upward-trend in relation to the random effect prediction, which contradicts the assumption 1 in Pinheiro and Bates (2000, p174) : "the within group errors are [...] independent of the random effect"**

Here, the residuals are computed as "observed value - fixed effect - random effect" and the fitted are computed as "fixed effect + random effect".

I suspect these strange patterns are because the number of observations per "group" is low. Indeed, they nearly disappear when I simulate data with more observations per group.

I include this test code below to illustrate the problem:

# generate data 
library(nlme)
library(ggplot2)

nb.group <- 500
list.group <- as.character(seq(1,nb.group, 1)) # anmes of the different groups
nb.obs <- 1000

alpha <- 0.3 # fixed effect
gamma <- 1.5 # sd of random effect
sigma <- 1 # sd of individual erro

db.RE <- data.frame( # generate value for random effect
  group = list.group,
  A = rnorm(nb.group, 0, gamma)
)

db <- data.frame(
  x = rgamma(nb.obs, 10, 5), # explicative variable
  group = sample(list.group, size = nb.obs, replace = TRUE), # vector of group
  E = rnorm(nb.obs, 0, sigma) # individual error
)

db <- merge(db, db.RE, by = "group") # attribute the value of random effect corresponding to each group

db$y <- alpha*db$x + db$A + db$E # y : independant variable


# calibrate with RE

model <- lme(fixed = as.formula(y ~ x), random = ~ 1|group, data = db)

# diagnosis 

summary(model)

db_pred <- cbind(predict(model, level = c(0:1)), res = residuals(model, level = c(0:1))) # get fitted values and residuals
db_ranef <- data.frame(group = rownames(ranef(model)), ranef(model)) # get random effect prediction
names(db_ranef) <- c("group", "ranef")
db_pred <- merge(db_pred, db_ranef, by = "group")

# plot residuals as a function of fitted (at the population level)
graph.1 <- ggplot(data = db_pred, aes(x = cut(predict.fixed,
                                   breaks = seq(min(predict.fixed), max(predict.fixed), length.out = 10)), y = res.group)) +
  geom_boxplot() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggtitle("residuals as a function of fitted (at the population level)")

# plot residuals as a function of fitted (at the individual level)
graph.2 <- ggplot(data = db_pred, aes(x = cut(predict.group,
                                   breaks = seq(min(predict.group), max(predict.group), length.out = 10)), y = res.group)) +
  geom_boxplot() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggtitle("residuals as a function of fitted (at the individual level)")

# qq plot of residuals 
graph.3 <- ggplot(data = db_pred, aes(sample = res.group)) +
  stat_qq() + stat_qq_line() +
  theme_bw() +
  ggtitle("qq plot of residuals")

# qq plot of random effect 
graph.4 <- ggplot(data = db_pred, aes(sample = ranef)) +
  stat_qq() + stat_qq_line() +
  theme_bw() +
  ggtitle("qq plot of random effect")

# plot residuals as a function of random effect prediction
graph.5 <- ggplot(data = db_pred, aes(x = cut(ranef,
                                   breaks = seq(min(ranef), max(ranef), length.out = 10)), y = res.group)) +
  geom_boxplot() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggtitle("residuals as a function of random effect prediction")


cowplot::plot_grid(plotlist = list(graph.1, graph.2, graph.3, graph.4, graph.5))

The corresponding diagnosis graphs are here : diagnosis graphs We observed two problems : (i) correlation between residuals and fitted value (including random effect) and (ii) correlation between residuals and random effect.

Is this correlation between residuals and fitted values (or random effects) expected when we have few observations per group, and could you explain to me why it should be the case?