4
$\begingroup$

I've fitted two models and while checking the residuals vs. fitted plot I saw a pattern that I haven't seen before in my analysis. The thing is, although the data seems to be well distributed (rectangular shape), it's all centered on the left, I first thought of this as a classical 'funnel-shaped' example and the following interpretation as detected heteroscedasticity, but a friend had a different point of view. (we were discussing the distribution on the left side of the x-axis and if how we'd be concerned about that)

  • Plot 1 (model 1):

plot1

  • model1:
# Can SCORE be explained by students(ID)'s performance on three different test types (types) ?
 
a <- lmer(SCORE ~ type + (1|ID), data = df_1, REML = F)
performance::check_model(a)
  • Plot 2 (model 2):

enter image description here

  • model2:
b <- lmer(SCORE2 ~ type + (1|ID), data = df_2, REML = F)
check_model(b)

Both plots look somehow different from this 'default' example:

plot3

Question for discussion: data doens't seen to be funnel-shaped, but it's all centered on the left of the x-axis. How should we interpret this situation? (specially concerning model1)

  • data1:
structure(list(ID = c("p1", "p1", "p1", "p2", "p2", "p2", "p3", 
"p3", "p3", "p4", "p4", "p4", "p5", "p5", "p5", "p6", "p6", "p6", 
"p7", "p7", "p7", "p8", "p8", "p8", "p9", "p9", "p9", "p10", 
"p10", "p10", "p11", "p11", "p11", "p12", "p12", "p12", "p13", 
"p13", "p13", "p14", "p14", "p14", "p15", "p15", "p15", "p16", 
"p16", "p16", "p17", "p17", "p17", "p18", "p18", "p18", "p19", 
"p19", "p19", "p20", "p20", "p20", "p21", "p21", "p21", "p22", 
"p22", "p22", "p23", "p23", "p23", "p24", "p24", "p24", "p25", 
"p25", "p25", "p26", "p26", "p26", "p27", "p27", "p27", "p28", 
"p28", "p28", "p29", "p29", "p29", "p30", "p30", "p30", "p31", 
"p31", "p31", "p32", "p32", "p32", "p33", "p33", "p33", "p34", 
"p34", "p34", "p35", "p35", "p35", "p36", "p36", "p36", "p37", 
"p37", "p37", "p38", "p38", "p38", "p39", "p39", "p39", "p40", 
"p40", "p40", "p41", "p41", "p41", "p42", "p42", "p42", "p43", 
"p43", "p43", "p44", "p44", "p44", "p45", "p45", "p45", "p46", 
"p46", "p46", "p47", "p47", "p47", "p48", "p48", "p48", "p49", 
"p49", "p49", "p50", "p50", "p50", "p51", "p51", "p51", "p52", 
"p52", "p52", "p53", "p53", "p53", "p54", "p54", "p54", "p55", 
"p55", "p55", "p56", "p56", "p56", "p57", "p57", "p57", "p58", 
"p58", "p58", "p59", "p59", "p59", "p60", "p60", "p60", "p61", 
"p61", "p61", "p62", "p62", "p62", "p63", "p63", "p63"), SCORE = c(69.33, 
52.42, 60.2, 50.09, 57.93, 56.79, 55.71, 55.77, 46.15, 59.74, 
61.02, 63.29, 60.6, 54.57, 65.28, 65.55, 66.31, 56.26, 57.63, 
59.74, 56.26, 60.33, 57.14, 58.31, 63.22, 63.38, 67.69, 52.82, 
59.26, 49.87, 56.29, 55.34, 53.33, 50.59, 50.4, 57.28, 70.68, 
75.26, 77.31, 55.71, 54.02, 55.34, 57.5, 64.3, 63.58, 76.89, 
68.43, 54.95, 64.08, 58.63, 63.42, 54.88, 54.09, 66.15, 57.2, 
60.81, 53.33, 56.73, 64.3, 61.93, 66.05, 72.23, 58.82, 61.25, 
61.33, 57.2, 60, 59.18, 65.18, 52.35, 61.93, 50.82, 60.65, 62.5, 
60.27, 66.96, 62.1, 68.05, 55.41, 61.34, 59.8, 58.14, 65.76, 
65.31, 63.82, 55.86, 60.74, 65.34, 56.27, 65.74, 61.6, 59.45, 
56.76, 58.93, 52.14, 59.79, 56.28, 53.6, 52.77, 64.11, 60.74, 
60.38, 59.09, 62.34, 64.34, 59.86, 60.55, 59.53, 67.27, 71.85, 
62.85, 62.52, 66.06, 60.9, 56.89, 63.01, 55.41, 56.73, 55.1, 
55.95, 63.83, 64.28, 68.05, 58.8, 61.76, 60.87, 67.27, 63.22, 
67.05, 68.57, 60.33, 62.51, 66.31, 60.2, 62.2, 65.58, 61.25, 
68.21, 65.3, 69.05, 72.63, 64.15, 58.46, 66.42, 60.14, 58.41, 
60.66, 66.37, 61.33, 58.53, 62.02, 60.25, 62.38, 66.6, 66.56, 
66.22, 54.64, 57.35, 62.18, 65.63, 65.13, 65.55, 63.67, 62, 58.14, 
75.93, 76.14, 72.5, 70.81, 74.38, 81.11, 56.52, 59, 55.77, 66.25, 
65.2, 59.74, 91.05, 78.49, 77.4, 63.47, 53.33, 63.58, 64.11, 
62.46, 67.14, 64.59, 63.03, 61.77), type = structure(c(1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("control", 
"testA", "testB"), class = "factor")), row.names = c(NA, -189L
), class = c("tbl_df", "tbl", "data.frame"))
  • data2:
> dput(df_2)
structure(list(ID = c("p1", "p1", "p1", "p2", "p2", "p2", "p3", 
"p3", "p3", "p4", "p4", "p4", "p5", "p5", "p5", "p6", "p6", "p6", 
"p7", "p7", "p7", "p8", "p8", "p8", "p9", "p9", "p9", "p10", 
"p10", "p10", "p11", "p11", "p11", "p12", "p12", "p12", "p13", 
"p13", "p13", "p14", "p14", "p14", "p15", "p15", "p15", "p16", 
"p16", "p16", "p17", "p17", "p17", "p18", "p18", "p18", "p19", 
"p19", "p19", "p20", "p20", "p20", "p21", "p21", "p21", "p22", 
"p22", "p22", "p23", "p23", "p23", "p24", "p24", "p24", "p25", 
"p25", "p25", "p26", "p26", "p26", "p27", "p27", "p27", "p28", 
"p28", "p28", "p29", "p29", "p29", "p30", "p30", "p30", "p31", 
"p31", "p31", "p32", "p32", "p32", "p33", "p33", "p33", "p34", 
"p34", "p34", "p35", "p35", "p35", "p36", "p36", "p36", "p37", 
"p37", "p37", "p38", "p38", "p38", "p39", "p39", "p39", "p40", 
"p40", "p40", "p41", "p41", "p41", "p42", "p42", "p42", "p43", 
"p43", "p43", "p44", "p44", "p44", "p45", "p45", "p45", "p46", 
"p46", "p46", "p47", "p47", "p47", "p48", "p48", "p48", "p49", 
"p49", "p49", "p50", "p50", "p50", "p51", "p51", "p51", "p52", 
"p52", "p52", "p53", "p53", "p53", "p54", "p54", "p54", "p55", 
"p55", "p55", "p56", "p56", "p56", "p57", "p57", "p57", "p58", 
"p58", "p58", "p59", "p59", "p59", "p60", "p60", "p60", "p61", 
"p61", "p61", "p62", "p62", "p62", "p63", "p63", "p63"), type = structure(c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("control", 
"testA", "testB"), class = "factor"), SCORE2 = c(0.7, 0.6, -1.6, 
0.4, -0.6, -0.199999999999999, -0.300000000000001, -1.1, 0.300000000000001, 
-2.8, -0.7, -1, 0.5, 1.4, 0.300000000000001, -1.3, -2.1, -0.3, 
-1.8, -2.4, -2.3, -1.1, 0.0999999999999996, -0.300000000000001, 
1.4, -1.2, 0.2, -0.199999999999999, -0.199999999999999, 0.199999999999999, 
-1, -0.300000000000001, 0.2, -2.8, 1.1, 2.4, -2, -1.2, -1.8, 
0.300000000000001, 0.5, -0.6, 0.7, -2.1, -2.4, 4.4, 4.4, 3.7, 
-1.6, 0.3, 0.3, -1.8, 1, 0.2, -0.6, -0.699999999999999, -0.300000000000001, 
-3, -2.4, -2.1, -0.2, 1.3, -0.0999999999999996, -1.6, -1.5, -0.6, 
0.3, 0.699999999999999, 2.3, 0, -0.800000000000001, 0.199999999999999, 
-2.5, -1.3, 0.6, -2.4, -1.8, -1, 0.4, -2.6, -0.3, 3, 1.2, 1.8, 
-2.8, -2.5, -0.0999999999999996, 0.9, 0, -0.4, 0.8, 0, 1.3, -2.6, 
-0.199999999999999, 0.8, -1.6, -1.6, -0.5, -0.8, -0.5, 0, 1, 
1.5, 0.300000000000001, -1.2, 0.300000000000001, -0.0999999999999996, 
-1.8, -1.5, -1.4, -1, -0.3, 1, -0.8, -0.199999999999999, -1.2, 
-1.8, 0.0999999999999996, -2.6, 0.5, -0.2, 1.3, -1.2, -0.6, -0.3, 
0.7, 2.8, -1.8, -0.6, -1.2, 1.7, -3, -1.6, -1.9, -1, -1.7, 2, 
0, 0.4, 0.8, 3.7, 4, 3, -2, -0.300000000000001, 0.5, 2.2, 3.8, 
3, -0.6, -1, 0, -0.6, -0.5, 1.3, 0.6, 0.699999999999999, 2, 0.4, 
2.8, 1.1, 0.2, -0.4, 0.4, -0.5, -1, -1, 0, 0.7, 1, -0.199999999999999, 
-0.4, 0.199999999999999, 1.4, 2.8, 1.6, -1.1, -1.8, 1, -0.2, 
1.2, 1.4, 0, -0.8, 0.2, 0.4, 0.4, 1.3)), row.names = c(NA, -189L
), class = c("tbl_df", "tbl", "data.frame"))
$\endgroup$
7
  • 1
    $\begingroup$ The "bunched" nature of the plots indicates some right skew or clustering somewhere, but not the kind that indicates any violated assumptions. I'd suggest looking at more graphs/descriptives to understand better. $\endgroup$ Nov 24, 2022 at 17:55
  • $\begingroup$ hi, BigBendRegion , any thoughts on how I should investigate the "right skew or clustering somewhere" based on the plots? I've studied a lot, but I've never seen such a pattern before $\endgroup$ Nov 25, 2022 at 15:46
  • $\begingroup$ Have a look at the q-q plots and histograms of all the variables, x and y. $\endgroup$ Nov 25, 2022 at 20:03
  • 1
    $\begingroup$ I would also argue that the "default" plot stems from really basic linear regression. The more complex the model becomes the less "funnely" these curves can become. $\endgroup$
    – Janosch
    Dec 1, 2022 at 9:29
  • 2
    $\begingroup$ Is the skewness you refer to the fact that most residuals lie to the left of 70 as predicted value? If so, that is just a fact of your data. Most of your data points have a predicted conditional mean of less than 70, that does not have to do with the distribution of your residuals. $\endgroup$
    – Kuku
    Dec 1, 2022 at 14:51

1 Answer 1

6
+50
$\begingroup$

As @Kuku points out, this pattern in the fitted value stems from a feature in the raw test scores: there is a quartet of high performing students with average test scores well above the rest of the class. A similar thing occurs at the bottom; however, the two students with the lowest scores don't perform much worse than the majority, so we don't see a big "gap" in the fitted values.

It's easiest to see these patterns by making our own residuals plots in order to highlight points of interest. The R code to reproduce the figures is attached at the end.

p13, p56, p57 & p60 are the students with the highest scores. They do quite a bit better than the rest of the class. We see this clearly if we plot the fitted scores against the raw scores.

Importantly, the model predicts the average score of over/under-achievers well but it under-predicts the within-student variation in test scores. Take for example student p60: the model predicts about the same score on all three tests (79, 78.5, 78.4) while in reality this student did particular well on the control test (91, 78.5, 77.4). This is a property of the model itself, so it might be relevant to understand.

What are the fitted values $\widehat{y}_{j,t}$ for student $j$ on test $t$? $$ \begin{aligned} \widehat{y}_{j,\text{control}} &= \widehat{\alpha}_j + \widehat{\beta}_0 \\ \widehat{y}_{j,\text{testA}} &= \widehat{\alpha}_j + \widehat{\beta}_0 + \widehat{\beta}_A \\ \widehat{y}_{j,\text{testB}} &= \widehat{\alpha}_j + \widehat{\beta}_0 + \widehat{\beta}_B \end{aligned} $$ where $\alpha_j$ is the student's random effect, $\beta_0$ is the expected score on the control test and $\beta_A$ and $\beta_B$ are the expected difference in score (relative to the control) on tests A and B. The hats indicate that these are estimates of the model parameters.

Here are the estimates for $\beta_0,\beta_A$ and $\beta_B$:

tidy(a, "fixed")
#>   effect term        estimate std.error statistic
#> 1 fixed  (Intercept)   62.2       0.797    78.1  
#> 2 fixed  typetestA     -0.672     0.708    -0.949
#> 3 fixed  typetestB     -0.516     0.708    -0.728

What we learn is that the predicted within-student variability is at most 0.672 (that's the difference between the control and A tests). And it is the same for every student! In practice there is more variability in the students' performance but the differences (between tests within a student) are not consistently in the same direction. That's why the improvement of test A and B compared to the control is only -0.7 and -0.5 points. (And this may be a good thing since the "improvement" is actually a decrease of about half a point on average.)

Created on 2022-12-04 with reprex v2.0.2

library("broom.mixed")
library("lme4")
library("tidyverse")

a <- lmer(SCORE ~ type + (1 | ID), data = df_1, REML = FALSE)

tidy(a, "fixed")

df_1a <- augment(a) %>% mutate(
  across(ID, as.character),
  # Standardize the residuals
  .std.resid = .resid / sigma(a)
)

df_1a %>%
  filter(
    .fitted >= 70
  ) %>%
  select(
    ID, SCORE, .fitted
  )

IDs_70plus <- filter(df_1a, .fitted >= 70) %>% pull(ID)
IDs_55minus <- filter(df_1a, .fitted <= 55) %>% pull(ID)

df_1a <- df_1a %>%
  mutate(
    label = if_else(ID %in% c(IDs_70plus, IDs_55minus), ID, "o"),
    color = case_when(
      ID %in% IDs_70plus ~ "a",
      ID %in% IDs_55minus ~ "b",
      TRUE ~ "c"
    )
  )

make_plot <- function(y, title.y, title = "X vs Y") {
  df_1a %>%
    ggplot(
      aes(.fitted, {{ y }})
    ) +
    geom_smooth(
      color = "#D9D9D9",
      se = FALSE
    ) +
    geom_text(
      aes(
        label = label,
        color = color
      )
    ) +
    labs(
      x = "Fitted values",
      y = title.y,
      title = title
    ) +
    guides(
      color = "none"
    )
}

make_plot(
  SCORE,
  "Test scores",
  "Fitted values against students' test scores"
)
make_plot(
  sqrt(abs(.std.resid)),
  "sqrt( |Std. residuals| )",
  "Fitted values against standardized residuals"
)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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