3
$\begingroup$

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).

enter image description here

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

enter image description here

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

enter image description here

y vs y_hat

enter image description here

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:

enter image description here

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:

enter image description here

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:

enter image description here

## 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()
$\endgroup$
  • 2
    $\begingroup$ This could be a case of omitted variable bias, probably a fixed effect. Can you edit your question and include more detail about the study design, the model and the data ? $\endgroup$ – Robert Long Jan 24 at 13:52
  • $\begingroup$ @RobertLong, hello, thank you. I added thw script (which uses a sample of real data from my repo). I think I used slightly different sample though as I don't see the same pattern again... I used a quantile mixed linear model because I need quantiles for practical reasons. The dataset contains all possible input variables, including those used in the experiment design. $\endgroup$ – Alexey Burnakov Jan 24 at 14:29
  • 1
    $\begingroup$ So the residual plot is from a model that hasn't converged? Ok you need to address that first. Tbh I am not surprised that it doesn't converge. Aside from the problem of probably having too many variables as random slopes, if you don't also include them as fixed then you are forcing their average effects to all be zero. As a first measure, fit them all as fixed only. $\endgroup$ – Robert Long Jan 28 at 18:59
  • 1
    $\begingroup$ Let me update my experiment and present it more clearly. I got your point about random slopes, thank you. $\endgroup$ – Alexey Burnakov Jan 28 at 19:47
  • 1
    $\begingroup$ Also, are you aware that when you use predict with re.form = NA you are setting all the random effects to zero ? Is that what you want ? $\endgroup$ – Robert Long Jan 28 at 21:04
2
$\begingroup$

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.
$\endgroup$
  • $\begingroup$ thank you! I do not experience missing data in my sample. As for the nonlinearity, given that I only use random intercepts from the random effects, does it mean that I definitely screwed the fixed effects? Should I expect good residual behaviour on fixed effects only? $\endgroup$ – Alexey Burnakov Jan 24 at 14:31
  • $\begingroup$ I found a model desing manually after many trials. It appears that intercept in the random block was causing the patterns - probably - because when it was included the model failed to converge. Everything is quite complicated. lmer( score_diff_delta ~ (handicap_pred + score_diff + diff_coef1 + diff_coef2 + total_score_relat + coef_diff) * measurements + (handicap_pred:score_diff:measurements -1 | teams) , data = desmat ) $\endgroup$ – Alexey Burnakov Jan 28 at 15:58
  • $\begingroup$ I am not sure if your answer touches on this behaviour, as I had to carefully craft non-linear terms as well, as you advised. By the way, is this a typo? "You have modeled appropriately the relationship...." Is not mising? $\endgroup$ – Alexey Burnakov Jan 28 at 15:59
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.