Mixed linear model's residuals diagnosing I am just learning mixed models. 
After fitting a mixed linear model with random intercepts:
$$Y_{ij} = \beta_0 + \beta_1 X_{ij} + b_i + \epsilon_{ij}$$
I study the fitted vs residual plot. It is noted that predictions were made on the training data without the inclusion of random subject-level intercepts into the result.
I observe:


*

*residuals are not stationary mean-wise 

*probably they violate homoskedasticity as well


Q: is this the result of a poor choice of fixed predictors / random predictor(s) / or both?
This is to show poputation and subject specific models with an intercept and one slope (powered by my expert task knowledge).

Population means correspnd pretty exactly with abline (which should be the case indeed), but the subject level models vary greatly. Below I provide data and approaches to the data I made.
Script #1 - ground zero:
## load libs

library(data.table)
library(ggplot2)
library(magrittr)
library(lme4)

# My real data
train_dat <- as.data.table(
     read.csv(
          url(
               'https://raw.githubusercontent.com/alexmosc/ds_lectures/master/train_dat.csv'
          )
     )
)

names(train_dat)

cat(
     'random levels: '
     , train_dat[, length(unique(teams))]
)


## Gaussian mixed model

desmat <- 
     train_dat[
          , .(
               measurements,
               score_diff_delta, 
               handicap_pred,
               score_diff
          )
          ] %>%
     as.matrix %>%
     scale %>%   ### scaling all numeric variables
     as.data.table %>%
     .[, teams := train_dat[, as.factor(teams)]]

rm(mixed.lmer)

mixed.lmer <- lmer(score_diff_delta ~ handicap_pred+score_diff + (1 | teams)
                   , data = desmat)

This model does not fit at all, resulting in error reading: 

Error in eval_f(x, ...) : Downdated VtV is not positive definite

After running certain checks I do not see obvious problems:
plot(density(desmat$score_diff_delta))
plot(density(desmat$handicap_pred))
plot(density(desmat$score_diff))
sapply(desmat, function(x) sum(is.na(x)))
plot(desmat$handicap_pred - desmat$score_diff, desmat$score_diff_delta)

Script #2: additional tries
rm(mixed.lmer)

mixed.lmer <- lmer(score_diff_delta ~ handicap_pred+score_diff + (handicap_pred+score_diff | teams)
                   , data = desmat)

results in 

boundary (singular) fit: see ?isSingular

Residuals for the model with zeroed random effects are not good:
y_hat <- predict(mixed.lmer, re.form = NA) # predict without random effects


Script #3: something that works better:
rm(mixed.lmer)

mixed.lmer <- lmer(score_diff_delta ~ handicap_pred+score_diff + (handicap_pred+score_diff-1 | teams)
                   , data = desmat)

residual vs y_hat

y vs y_hat

That works much better, but I don't understand why a random intercept causes the issues.
Adding more slopes can improve the model, but still without the random intecept.
Below are older rather chaotic attempts at this problem:

Experiment script:
## load libs

library(data.table)
library(ggplot2)
library(magrittr)
#library(lme4)
library(lqmm)

# My real data
train_dat <- as.data.table(
     read.csv(
          url(
               'https://raw.githubusercontent.com/alexmosc/ds_lectures/master/train_dat.csv'
               )
          )
     )

names(train_dat)

cat(
     'random levels: '
     , train_dat[, length(unique(teams))]
)


## Laplace mixed model


# design matrix

frml <- formula(
     score_diff_delta ~ # target
          measurements * # time variable
          (
               handicap_pred+ 
                    coef_diff+ 
                    score_diff+ 
                    start_coef1_win1+
                    start_coef2_win2+
                    diff_coef1+ 
                    diff_coef2
          ) # in-game dynamics
)

qs <- 0.5

qmem <- 
     lqmm(
          fixed = frml, 
          random = ~ 1, #handicap_pred + coef_diff + score_diff, ### experiments with random slopes
          group = teams, 
          covariance = "pdDiag", 
          tau = qs,
          nK = 7, 
          type = "normal", 
          rule = 1, 
          data = train_dat,
          fit = TRUE,
          control = lqmmControl(method = 'gs')
     )

# fixed coefs

qmem_coefs <- coefficients(qmem)

# rand coefs
if(F)
{
     rnf <- lqmm::ranef(qmem)
     rnf$teams <- rownames(rnf)
     setDT(rnf)
     colnames(rnf)[colnames(rnf) != 'teams'] <- paste0(colnames(rnf)[colnames(rnf) != 'teams'], '_rnf')
     setkey(rnf, teams)
}

## residual test

if(F)
{
     train_dat <- train_dat[rnf]

     desmat <-
          model.matrix(
               frml
               , data = train_dat
          )

     y_hat <- desmat %*% qmem_coefs %>% as.vector %>% `+`(train_dat$rnf_int)
} # matrix calculation of estimates with respect to random factors

y_hat <- predict(object = qmem, level = 1) ### default predicting with respect to random factors (!!!! level = 1 !!!)

y <- train_dat[, score_diff_delta]

dt <- data.table(
     y = y
     , y_hat = y_hat
) %>%
     .[, resids := y - y_hat]

ggplot(dt) +
     geom_point(
          aes(
               x = y_hat
               , y = resids
          )
          , size = 2
          , alpha = 0.2
     ) +
     geom_smooth(
          aes(
               x = y_hat
               , y = resids
          )
     ) +
     geom_hline(
          aes(
               yintercept = 0
          )
          , color = 'red'
     ) +
     theme_minimal()

ggplot(dt) +
     geom_point(
          aes(
               x = y_hat
               , y = y
          )
          , size = 2
          , alpha = 0.2
     ) +
     geom_smooth(
          aes(
               x = y_hat
               , y = y
          )
     ) +
     geom_abline(color = 'red') +
     theme_minimal()


It is not really the answer answering my questions, but I need to add some information for further research.
I did the same analysis of residuals, but this time I used a Gaussian mixed linear model thanks to package lme4.
I used almost the same formula (only one term was excluded because it was highly correlated with the other one and caused errors). The same change in the lqmm model resulted in the plot identical to what I got in the first time.
With lmer I am getting a drastically different outcome, where the residuals are just fine:

It makes me more confused this time, as I thought the median VS mean would result in approximately the same fits...
Even stranger to me is the fact that when predicting without addition of the random effects, y_hat <- predict(mixed.lmer, re.form = NA), the residuals become a real mess. This is frustrating since the model should work well on the data with completely new random levels:

## load libs

library(data.table)
library(ggplot2)
library(magrittr)
library(lme4)

# My real data
train_dat <- as.data.table(
     read.csv(
          url(
               'https://raw.githubusercontent.com/alexmosc/ds_lectures/master/train_dat.csv'
               )
          )
     )

names(train_dat)

cat(
     'random levels: '
     , train_dat[, length(unique(teams))]
)


## Gaussian mixed model

desmat <- 
     train_dat[
          , .(
               measurements,
               score_diff_delta, 
               handicap_pred,
               coef_diff,
               score_diff,
               start_coef1_win1,
               start_coef2_win2,
               diff_coef1,
               diff_coef2,
               handicap_diff,
               final_score_diff
               )
          ] %>%
     as.matrix %>%
     #scale %>%   ### scaling all numeric variables
     as.data.table %>%
     .[, teams := train_dat[, teams]]

mixed.lmer <- lmer(score_diff_delta ~ 
                        measurements * (
                             handicap_pred+ 
                                  coef_diff+ 
                                  #score_diff+ 
                                  start_coef1_win1+
                                  start_coef2_win2+
                                  diff_coef1+ 
                                  diff_coef2
                        ) + (1|teams)
                   , data = desmat)

summary(mixed.lmer)

plot(mixed.lmer)

y_hat <- predict(mixed.lmer, re.form = NULL) # predict using random effects

y <- desmat$score_diff_delta

resids <- y - y_hat

sd(resids)

#plot(y_hat, resids)
#plot(y_hat, y)

dt <- data.table(
     y = y
     , y_hat = y_hat
     , resids =resids
)

ggplot(dt) +
     geom_point(
          aes(
               x = y_hat
               , y = resids
          )
          , size = 2
          , alpha = 0.2
     ) +
     geom_smooth(
          aes(
               x = y_hat
               , y = resids
          )
     ) +
     geom_hline(
          aes(
               yintercept = 0
          )
          , color = 'red'
     ) +
     theme_minimal()

 A: There could multiple reasons why you observe this systematic behavior. A couple are:


*

*You have modeled appropriately the relationship between $Y$ and $X$. For example, it could be that $X$ is nonlinearly associated with $Y$. You could investigate that by relaxing the linearity assumption either using polynomials or splines.

*You have missing data in $Y$ which are of the missing at random type. In this case the model may perfectly correct but you see a systematic pattern in the residuals because of the missing data. I.e., because the observed data are selected sample from your target population. 

A: So, I figured where the problem is with high probability. The real issue is my data...
I created target variable $Y$ so that for each random subject the dependancy on covariates was 100%.
Call:
lm(formula = score_diff_delta ~ score_diff, data = ddt)

Residuals:
       Min         1Q     Median         3Q        Max 
-3.202e-15 -2.600e-18  7.460e-17  7.460e-17  6.077e-16 

Coefficients:
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  1.600e+01  8.434e-16  1.897e+16   <2e-16 ***
score_diff  -1.000e+00  6.253e-17 -1.599e+16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.625e-16 on 35 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 2.557e+32 on 1 and 35 DF,  p-value: < 2.2e-16

Warning message:
In summary.lm(.) : essentially perfect fit: summary may be unreliable

It caused the problems with convergence, I believe. That was a suprise to find because I thought (without good understanding of mixed models) that various subjects would contribute to a mix of covariates / targets (like with minibatching for the neural network) and the fit would not be nearly 100%. It appears the model builds estimates for each level of the random factor and them plunges into this singularity, and creates some strange patterns of coefficients. So this is totally my desing flaw. I will need to rethink my approach (Mixed effect predictive model where dependent variable is constant for each random instance) to $Y$ first.
