Do Schools Discriminate Against Night Owls? Inference on Clustered Data (part 3)

Clustered Standard errors

One of the assumptions of the linear regression model is that the observations are independent and that the error term is constant across observations. These assumptions are reflected in what the variance covariance matrix \(\Sigma\) should look like. Independence of observations means that the observations are uncorrelated with each other so that the off-diagonal elements are zero. The constant error term (or homoscedasticity) implies that terms on the diagonal are the same.

\[\Sigma = \begin{bmatrix} \sigma & 0 & 0\\ 0 & \sigma & 0 \\ 0 & 0 & \sigma \end{bmatrix}\]

If these assumptions are not met, our model will be misspecified. A common way to proceed is to not fix the problems in the model, but to carry on with the misspecified model and patch it up by using robust standard errors. These are also known as sandwich estimators for the standard errors, because a key trick in estimating the robust standard error is to use the accompanying matrices to reduce the n by n matrix \(\Sigma\) to a k by k matrix, where k is the number of variables. That means far less parameters need to be estimated. We forgo a discussion here of how this is done.

The whole approach sounds dodgy. Why not fix the problems in the model instead of carrying on with a misspecified model? Moreover, if the assumptions of independence and homoscedasticity are not met, what else could be wrong with the model? Indeed, often a lot more is wrong with the model, which can lead to radically different estimates of the independent variables in the model. See here for an example. Patching up standard errors does nothing to change these estimates. Indeed, some argue for using robust standard errors as a diagnostic device. If they differ too much from classical standard errors, it signals problems in model specification, that must be fixed.

In the case of multilevel data structures, robust standard errors can be used to correct for the correlation between observations. These do not necessarily bias the estimation of the explanatory variables in the model. However, robust standard errors shouldn’t be used to fix other specification problems in the model or if the researcher cares about effects at the group level.

A bunch of simulations described by Bosker and Snijders in Chapter 12 of the 2nd edition of their Multilevel Analysis book confirm that the standard errors for fixed effects in otherwise well specified models can be estimated adequately with robust standard errors. This sounds dodgy, but from these simulations it appears that if there are around 40 clusters, which are not correlated among themselves and fairly balanced in size, it is okay to carry on with the misspecified classical linear model to estimate fixed effects. It must be mentioned though that using robust standard errors reduces the efficiency of the estimation.

Applying robust standard errors in our case to estimate the effect of chronotype is iffy. To do so, we should take the highest level of the multilevel into account. Arguably, that is the level of the students (havo or vwo). However, these are just two clusters, which are not even well balanced. Alternatively, we can use classes as the highest level of clustering, which brings us up to 20 well balanced clusters, and adjust for the level (havo/vwo) in the model specification.

m1coeffs_std <- data.frame(summary(fit_lm)$coefficients)
coi_indices <- which(!startsWith(row.names(m1coeffs_std), 'Class'))

m1coeffs_cl <- coeftest(fit_lm, vcov = vcovCL, cluster = ~Class)
m1coeffs_cl[coi_indices,]

##                       Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept)      12.7677699513 2.772980e+00  4.6043502 4.150006e-06
## MSFsc            -0.0564239655 2.428547e-02 -2.3233626 2.016456e-02
## late_13_14       -0.0625315624 2.027324e-02 -3.0844379 2.040732e-03
## removed_13_14    -0.0835816275 1.760564e-02 -4.7474346 2.067101e-06
## sick_dur         -0.0174821718 3.793230e-03 -4.6087825 4.062569e-06
## sexmale          -0.1242332652 4.511758e-02 -2.7535445 5.897981e-03
## age               0.0697627483 4.606825e-02  1.5143348 1.299487e-01
## Topicchemistry    0.1376305381 2.136292e-01  0.6442496 5.194172e-01
## TopicDutch       -0.3756184337 7.653118e-02 -4.9080445 9.234285e-07
## TopicEnglish      0.1233631008 1.235213e-01  0.9987196 3.179364e-01
## Topicgeography   -0.2328525962 9.440785e-02 -2.4664537 1.364991e-02
## Topichistory     -0.3327347013 7.102643e-02 -4.6846606 2.813259e-06
## Topicmathematics -0.1008544870 1.090878e-01 -0.9245258 3.552181e-01
## Topicphysics     -0.7068510556 1.051582e-01 -6.7217900 1.818598e-11
## Period           -0.0481227796 2.003387e-02 -2.4020713 1.630697e-02
## stud_ID          -0.0000488132 2.026035e-05 -2.4092970 1.598768e-02

Note that the standard errors have increased across the board by about an order of magnitude, as expected. We will compare these results with those of the multilevel model.

Multilevel models

A more principled way to model the data is to have levels in the model that match the levels in the data generating process. That is, we treat grades as the smallest unit in the analysis, students as a first level and classes as a second level, where students are nested inside classes.

One approach is to set up models for the group means: the mean score of each student could be modeled based on chronotype, age, sex and so on and the mean scores of each class could be modeled based on its level, for example. After running these two models we could set up models inside each group to estimate for example the grades inside a class. This way we are already doing a kind of multilevel analysis. We are taking the correlations inside classes and for individuals into account after all, by setting up a model for each one of them. But we can do better.

If we are willing to assume that the classes are drawn from some (super-)population, then we can set up a probability distribution for them. Say, a normal distribution. We can learn from the data what the mean and the standard deviation of this distribution is. What good is there in setting up this distribution? The answer is that it allows us to learn something about the class mean we expect by moving from class to class. If we get an exceptional result for one class sample, we can ‘regularize’ it based on the distribution for class means we learned from the other classes. This makes intuitive sense from a Bayesian point of view - where a so called prior plays this role - but for those unfamiliar with this way of thinking I refer to the formula I discussed here. The same goes for the groups of individual students (confusingly, because grades are the lowest level of analysis, individual students are also groups).

The starting point for multilevel modeling is sometimes called the null model. We use lmer to run a basic model with two levels on the imputed data sets.

null_fml <- "Grade ~ 1 + (1 | stud_ID)"
fit_null <- lmerModList(null_fml, data = a.out$imputations)

fastdisp(fit_null)

## lmer(formula = Grade ~ 1 + (1 | stud_ID), data = d)
##             estimate std.error
## (Intercept)     6.60      0.03
## 
## Error terms:
##  Groups   Name        Std.Dev.
##  stud_ID  (Intercept) 0.62    
##  Residual             1.43    
## ---
## number of obs: 40890, groups: stud_ID, 523
## AIC = 146788---

We observe that the standard deviation of grades between students is 0.62 and within students 1.43. So most of the spread in grades happens ‘inside’ a student (because of the amount of learning or the topic, we may assume).

We must make a reservation here. If we turn the sd’s into variances, it would be tempting to add them up and call the sum the total variance. This is what the law of total variance in probability theory says after all, which we loosely rephrase as:

• total variance = variance within groups + variance between groups

We are not dealing with population variances however, but with samples, so that there is uncertainty. Some of the variance between groups could be due to factors at the individual level that are distributed unevenly among the groups. Since this is a sampling problem, it goes away as the size of the sample goes to infinity. For the formula for the total variance in the case of a sample, see Snijders and Bosker 2012, p. 20.

What all of this means is that the numbers we just gave for the spread of grades inside students and between students are conditional on a bunch of individual level variables that could mask the true variances in the population.

We can do yet more justice to the data generating process by adding a level for classes as well.

fit_null_tl <- lmerModList(Grade ~ 1 + (1 | Class/stud_ID),
                 data = a.out$imputations)
fastdisp(fit_null_tl)

## lmer(formula = Grade ~ 1 + (1 | Class/stud_ID), data = d)
##             estimate std.error
## (Intercept)     6.61      0.07
## 
## Error terms:
##  Groups        Name        Std.Dev.
##  stud_ID:Class (Intercept) 0.29    
##  Class         (Intercept) 0.55    
##  Residual                  1.43    
## ---
## number of obs: 40890, groups: stud_ID:Class, 528; Class, 20
## AIC = 146719---

In order to arrive at the author’s model we add a forth level, consisting of the Level (havo/vwo) of the students. It is debatable whether these levels are exchangeable though. That is, we need to believe clusters are similar, so that it makes sense to set up a probability distribution for them (i.e. to treat them as ‘random’, in the frequentist lingo). Of course, this objection would apply to treating classes as exchangeable just as much.

We also include the explanatory variables that the authors settled upon. These variables were featured in a model that was the winner of a competition of nine models, that was decided by looking at their AIC scores. Fitting several models is better than one model. But why nine models and not ten or twenty? In response to these questions, a new trend in running regressions is to specify all possible models and look at the specification curve of tens of thousands of specifications. I hope to pursue this line of analysis in a future post.

# Author's model; Four levels as authors, with Level havo/vwo as highest level - but only two categories!
fit_null_fl <- lmerModList(Grade ~ 1 + sex + MSFsc + late_14_15 + removed_14_15 + sick_dur + Period + subject_area + (1 | Level/Class/stud_ID),
                    data = a.out$imputations)
fastdisp(fit_null_fl, digits=3)

## lmer(formula = Grade ~ 1 + sex + MSFsc + late_14_15 + removed_14_15 + 
##     sick_dur + Period + subject_area + (1 | Level/Class/stud_ID), 
##     data = d)
##                        estimate std.error
## (Intercept)               7.129     0.155
## late_14_15               -0.009     0.007
## MSFsc                    -0.053     0.014
## Period                   -0.044     0.006
## removed_14_15            -0.022     0.007
## sexmale                  -0.138     0.036
## sick_dur                 -0.018     0.004
## subject_areascientific    0.066     0.016
## 
## Error terms:
##  Groups                Name        Std.Dev.
##  stud_ID:(Class:Level) (Intercept) 0.233   
##  Class:Level           (Intercept) 0.176   
##  Level                 (Intercept) 0.52    
##  Residual                          1.43    
## ---
## number of obs: 40890, groups: stud_ID:(Class:Level), 528; Class:Level, 20; Level, 2
## AIC = 146655---

Across the board, we note that the coefficient estimates are quite similar to the original classical linear model that we fitted with lm. However, the standard errors are again quite a bit bigger, as in the case of the robust standard errors. Compared to the latter case, the standard errors of the four-stage multilevel model are bigger overall. So robust standard errors come close, but underestimate uncertainty a little bit, which is consistent with simulations discussed in chapter 12 of Snijders and Bosker (2012).

By adding the variables preferred by the authors, we can investigate the interaction between STEM subjects and chronotype in a multilevel model, where the influence of classes and students is taken into account.

# the difference between imputations is very small, so we take the first for convenience
fit_2 <- lmer(Grade ~ 1 + Level + sex + late_14_15 + removed_14_15 + sick_dur + Period + subject_area*MSFsc + (1 | Level/Class/stud_ID),
              data = a.out$imputations[[1]])
fastdisp(fit_2)

## lmer(formula = Grade ~ 1 + Level + sex + late_14_15 + removed_14_15 + 
##     sick_dur + Period + subject_area * MSFsc + (1 | Level/Class/stud_ID), 
##     data = a.out$imputations[[1]])
##                              coef.est coef.se
## (Intercept)                   6.82     0.35  
## Levelvwo                      0.27     0.49  
## sexmale                      -0.11     0.04  
## late_14_15                   -0.02     0.01  
## removed_14_15                -0.03     0.01  
## sick_dur                     -0.02     0.00  
## Period                       -0.04     0.01  
## subject_areascientific        0.59     0.06  
## MSFsc                        -0.01     0.01  
## subject_areascientific:MSFsc -0.12     0.01  
## 
## Error terms:
##  Groups                Name        Std.Dev.
##  stud_ID:(Class:Level) (Intercept) 0.52    
##  Class:Level           (Intercept) 0.23    
##  Level                 (Intercept) 0.34    
##  Residual                          1.43    
## ---
## number of obs: 40890, groups: stud_ID:(Class:Level), 528; Class:Level, 20; Level, 2
## AIC = 146592

We observe that the negative effect of chronotype on grade is indeed more pronounced for STEM subjects.

Analysis of Group Differences

Even though the results of the multilevel model are very similar to those of the robust standard error approach, there is an advantage to the former. By modeling the levels in the data generating process, we can investigate differences between groups. In this spirit, we fit a model with Class as the highest level, since it has a decent amount of clusters.

fit_4 <- lmer(Grade ~ 1 + MSFsc + (1 + MSFsc | Class/stud_ID),
              data = a.out$imputations[[1]])
fastdisp(fit_4)

## lmer(formula = Grade ~ 1 + MSFsc + (1 + MSFsc | Class/stud_ID), 
##     data = a.out$imputations[[1]])
##             coef.est coef.se
## (Intercept)  6.90     0.12  
## MSFsc       -0.07     0.02  
## 
## Error terms:
##  Groups        Name        Std.Dev. Corr  
##  stud_ID:Class (Intercept) 0.65           
##                MSFsc       0.06     -0.61 
##  Class         (Intercept) 0.44           
##                MSFsc       0.04     -1.00 
##  Residual                  1.43           
## ---
## number of obs: 40890, groups: stud_ID:Class, 528; Class, 20
## AIC = 146712

The usual assumptions of the linear model carry over to the multilevel version. If we wanted to make predictions from the model, we would want to check the homoscedacity of the residuals. We have no such intentions though, so we won’t. For the inference on the random parameters (the intercepts and slopes), we do need to check if it is reasonable to impose a normal distribution on them.

# for intercept random
qqnorm(ranef(fit_4)$Class[,1] )
qqline(ranef(fit_4)$Class[,1], col = "red")

# for slopes random
qqnorm(ranef(fit_4)$Class[,2] )
qqline(ranef(fit_4)$Class[,2], col = "red")

The fit is not perfect, but it seems acceptable.

We will look now at the plots of the estimated intercepts and slopes for Class and Student, according to our multilevel model with Class as the highest level.

reEx <- REsim(fit_4)
p1 <- as.ggplot(plotREsim(reEx))
p1

We shape of these graphs reflects that there are 20 classes and about 500 students. Grayed out confidence intervals intersect with 0. This is true for most intercepts of the classes and students. However, the slopes for chronotype are what really matters. In a way, the bottom part of the diagram is terrible, because we cannot see the differences between the slopes of the classes. We can do better if we wanted, like so.

p2 <- as.ggplot(plotREsim(reEx, facet= list(groupFctr= "Class", term= "MSFsc")))
p2

But the instructive point of the first, faceted diagram is that all the slopes are significantly different from zero (since they are black) even though they are hugging the horizontal line at 0. In other words, we have so many observations that the confidence interval is tiny, so as to be invisible.

This raises the question of how much shrinkage takes place in the multilevel model. Multilevel regression tries to find the sweet spot between pooled and unpooled regression. If there are very few observations per group, then the estimates inside the groups will be very similar to the estimates for the entire population. On the other hand, if there are a lot of observations in the groups, then the grand parameters will have a hard time pulling group parameters towards them. The latter may be the case here. Do does multilevel regression even make a difference here, given the large number of observations?

To that end we first estimate the regression lines separately for every class. We then pull the estimates of the multilevel model out of the lme4 object and add the relevant fixed and random parts together to get the multilevel estimates for the regression line per group. We won’t show uncertainty estimates here, because what we care about here is the amount of pull towards the unpooled estimates. So let’s see what happens.

# CREATE DF OF UNPOOLED REGRESSION COEFS

#function of the model
MSF_model <- function(df) {
  lm(Grade ~ 1 + MSFsc, data = df)
}
# loop over the groups
unpooled_coefs <- list()
for(i in seq_along(levels(a.out$imputations[[1]]$Class))){
  unpooled_coefs[[i]] <- a.out$imputations[[1]] %>% 
    filter(Class %in% c(levels(Class)[i])) %>% 
    MSF_model %>% 
    tidy %>% 
    filter(term %in% c("(Intercept)","MSFsc")) %>% 
    dplyr::select(term,estimate)
}
# pull the class names
classes <- scores %>% pull(Class) %>% levels()
# stuff everything into a df
unpooled_coefs_df <- map_dfr(unpooled_coefs,`[`, c("term", "estimate")) %>% 
  group_by(term) %>%
  mutate(row = row_number()) %>%
  tidyr::pivot_wider(names_from = term, values_from = estimate) %>%
  dplyr::select(-row) %>% 
  mutate(classes = classes) %>% 
  dplyr::select(classes, everything())

# CREATE A DF OF MULTILEVEL REGRESSION COEFS

#pull slopes and intercepts
class_re <- ranef(fit_4)$Class[[2]]
total_sl <- summary(fit_4)$coefficients[2]
multilevel_slopes <- total_sl + class_re 

class_re_int <- ranef(fit_4)$Class[[1]]
total_sl_int <- summary(fit_4)$coefficients[1]
multilevel_intercepts <- class_re_int + total_sl_int 
# add together in df
classes <- scores %>% pull(Class) %>% levels()
df_multilevel <- tibble(classes,multilevel_intercepts, multilevel_slopes)


# PLOT BOTH REGRESSION LINES PER GROUP

# some text to tell the lines apart
dat_text <- data.frame(
  label = "Multilevel",
  classes   = "1h3"
)
dat_text2 <- data.frame(
  label = "Unpooled",
  classes   = "1h3"
)
# the actual plot
  ggplot() +
  scale_y_continuous(limits=c(4,8)) +
  scale_x_continuous(limits=c(0,8)) +
  geom_abline(data=df_multilevel,aes(
    slope=multilevel_slopes,
    intercept=multilevel_intercepts),
    color="darkblue") +
  geom_abline(data=unpooled_coefs_df,aes(
      slope=MSFsc,
      intercept=`(Intercept)`), color="pink") +
    facet_wrap(vars(classes)) +
    labs(x="chronotype", y="Grade") +
    ggtitle("Multilevel and unpooled regression lines per class") +
      geom_text(
      data    = dat_text,
      color = "darkblue",
      mapping = aes(x = -Inf, y = -Inf, label = label),
      hjust   = -0.1,
      vjust   = -1
    ) +
    geom_text(
      data    = dat_text2,
      color = "pink",
      mapping = aes(x = 6, y = 7.5, label = label),
    )

We note that the blue multilevel lines are like a sceptical version of the unpooled ‘naive’ pink lines. Whenever the separate group regressions (or pooled regression) are sloping up or down too wildly, the multilevel model regularizes the estimates towards the mean of the distribution of the slopes.

If we believe that classes are exchangeable and that the normal is an appropriate distribution of the random intercepts and slopes, then our speculations about the lack of motivation of havo 3 students seem unwarranted. However, the point of our speculations was that havo 3 classes are different from other classes in terms of motivation, so that they may not be exchangeable. So a multilevel model on its own can’t really settle the type of speculation we engaged in. One would have to measure students motivation is some way, if one though the speculation was worth pursuing. However, it isn’t, unless there is a more convincing scientific story about the effect of the interaction of motivation and chronotype on grades.