The Effects of Misspecifying the Random Part of Multilevel Models

Random Effects and Standard Errors

To understand how the misspecification of the random part of the model affects the standard errors for the fixed effects, it is helpful to consider the unrealistic case where the variance components are known, every Level 1 predictor has a random slope, and each group has sufficient size to calculate ordinary least squares (OLS) estimates. Following Raudenbush and Bryk (2002, Chapter 3), the Level-1 model for group j in an example with P Level 1 predictors can be expressed as

where Y_j is a vector of n_i outcome scores. X_j is a n_i by (P + 1) matrix where the first column is a series of ones to represent the intercept and the other columns are the Level 1 predictors. β_j is a vector of unknown parameters and r_j are the random errors with a mean of zero and a variance of σ². The ordinary least squares (OLS) estimator for the unknown parameters is

The error variance matrix for ${\widehat{\beta }}_j$ in estimating β_j is denoted as V_j and is given by

Finally, pre-multiplying both sides of Equation 2 by ${\left(X_j^T{X}_j\right)}^{-1}{X}_j^T$produces a simpler model for ${\widehat{\beta }}_j$:

At Level 2, the model for β_j is

In Equation 6, $\gamma$ is a F by 1 vector of fixed effects, the matrix Z_j is a (P + 1) by F matrix that is stacked in block diagonal fashion. For example, in a model with one Level 2 variable Z_j predicting one Level 1 intercept and one Level 1 slope, Equation 6 would be:

Both u_0j and u_1j are the random effects representing the intercept and slope, respectively. They are assumed to be normally distributed with a mean vector of zeros and a variance covariance matrix Τ.

Substituting Equation 6 into Equation 5 creates a single combined model:

The dispersion of ${\widehat{\beta }}_j$ conditional on Z_j is a function of the error variance matrix at Level 1 and the random effects matrix at Level 2:

Finally, the dispersion matrix for the fixed effects (γ’s) based on generalized least squares is given by:

The square roots of the diagonals of V_γ are the standard errors. To trace the effects of Δ_j on the standard errors, it is helpful to consider a very simple model with two Level 1 predictors and no Level 2 predictors. The underlying concepts are identical for models with more fixed and random effects but are more difficult to demonstrate. Thus, in our simple model, Z is simply a 3 × 3 identity matrix and Equation 10 becomes:

standard errors for γ₀₀, γ_10, and γ₂₀ are simply the square roots of the diagonals of V_j. Importantly, the slope variance for one Level 1 predictor does not impact the standard error for the other Level 1 predictor. It follows that fixing the slope of a given predictor when it should vary should only affect the standard errors associated with that predictor.

In practice, the variance components are not known and need to be estimated with the fixed effects using full or restricted maximum likelihood estimation. Based on the above discussion, we would expect the bias in standard errors would be limited to those fixed effects that are associated with the random effect that is misspecified and that the amount of this bias should decrease as the amount of slope variance decreases. In addition, the amount of bias in standard errors should increase as group size increases. This is because, all things being equal, the elements of the error matrix V_j will decrease as group size increases. Thus, the Level 2 variances will have a greater relative impact on the standard errors of the regression coefficients.

Data Generation

We created the Level 2 variable Z by randomly sampling values from a standard normal distribution. We also generated values for the Level 2 random variables by randomly sampling values from a normal distribution with a mean of 0 and a standard deviation equal to the square root of the desired variance. For u_0j, the standard deviation was set to the value to produce an ICC for Y of .20. The value changed based on the amount of explained variance at each level. For u_1j, the standard deviation was set to the value to produce the desired slope variance. At Level 1, we randomly sampled values from a standard normal distribution to produce the Level 1 predictor X. In addition, the Level 1 error term, r_ij, was created by randomly drawing values from a normal distribution with a mean of 0 and a standard deviation set to produce the desired ICC for Y.

To generate the Y, we first created intercepts for each group by multiplying Z by 0.40 (the cross-level effect) and added u_0j. Similarly, we created the group slopes (β_1j) by adding the desired mean slope (i.e., 0.00, 0.10, 0.30, or 0.50) and the product of Z with the cross-level interaction effect to u_1j. Finally, we computed values for Y using the group intercept and slope values and the randomly generated values for X and r_ij.

Analyses

For each sample, we estimated two multilevel models using restricted maximum-likelihood estimation: one where the slope was fixed and one where the slope varied. For each model, we recorded the estimates for γ₀₁, γ₁₀, and γ₁₁, as well as their standard errors. We calculated the bias in standard errors by subtracting the standard error from the data generating model from that of the misspecified model. Thus, a positive value for bias reflects the standard errors for the misspecified model was higher than the data generating model. Cohen’s d statistics were also calculated to aid in interpreting the size of the bias. We used the guidelines (Cohen, 1988) of .20, .50, and .80 as indicative of small, medium, and large effects, respectively. We conducted ANOVA’s to determine the impact of the various conditions on the bias in standard errors. Because of the large number of samples, we focused on interpreting effect sizes (partial η²). All analyses were conducted using the LME4 (Bates, Maechler, Bolker, & Walker, 2015) package in the statistical program system R (R Core Team, 2018).

Bias in Standard Errors

Table 1 presents the average bias statistics for each condition. These statistics were calculated by collapsing across the other conditions. Misspecification of the random part of the model did not appear to impact γ₀₁. The Cohen’s d statistics ranged from -0.02 to 0.02. In addition, the bias values were small: They were all near .000.

Misspecification of the random part of the model did influence the standard errors for γ₁₀ and γ₁₁. Cohen’s d’s ranged from -1.14 to 0.09 for γ₁₀ and -1.08 to 0.09 to for γ_11. The smallest values were found for the no slope variance condition, and the biggest values was found for the largest group size condition. When there was no population slope variance, estimating a model with a random slope produced a bias of .001 (Cohen’s d = 0.09) in the standard errors for both γ₁₀ and γ₁₁. In contrast, misspecifying the slope as fixed when it should have been random produced biases ranging from -.019 to -.005 for both γ₁₀ and γ₁₁. Cohen’s d ranged from -1.03 to -0.28 for γ₁₀ and -0.96 to -0.26 for γ₁₁. As expected, bigger biases occurred for greater levels of slope variance.

ANOVA Results

To get a better sense of the effects of the study conditions on bias, we conducted separate ANOVA’s investigating the impact of the study characteristics on the bias in standard errors for the γ₁₀ and γ₁₁. We included only the conditions where there was slope variance to produce a fully crossed design. Table 2 presents the results. Results were generally consistent with those for both ANOVA’s. Not surprisingly, the biggest effect size was for the amount of slope variance (partial η² = .65 for γ₁₀ and .62 for γ₁₁). The number of groups (partial η² = .38 for γ₁₀ and .36 for γ₁₁) also had a sizable effect size. The partial η² for group size was .11 for γ₁₀ and .10 for γ₁₁. Based on Table 1, it appears that both larger numbers of groups and group size produced less bias. Note that as the number of groups increased, the pooled variance of bias decreased. Thus, the Cohen’s d for bias as a function of the number of groups did not exhibit a similar pattern as the raw bias. The slope mean and the amount of intercept – slope correlation did not have large effect sizes (partial η²s = .00). Finally, only the interaction between the number of groups and amount of slope variance had a sizeable effect (partial η²s = .11). Figure 1 displays the form of the interaction for γ₁₀. The effect for slope variance decreased as the number of groups increased.

Design

We used a 3 (number of groups) × 3 (average group size) × 4 (slope correlation) × 5 (slope variance effect size for X₁) × 2 (slope variance effect size for X₂) × 4 (predictor correlation) × 2 (between group variance for X₁) × 2 (between group variance for X₂) design. It was not fully crossed because there can be no slope correlation when either of the slopes have zero variance. The values for the number of groups, average group size, and slope variance effect size for X₁ were consistent with Study 1. The cross-level interaction and slope mean of X₁ were set at 0.20 and 0.30 respectively. We eliminated the other conditions for these variables because Study 1 demonstrated negligible differences in standard error bias across these conditions. In addition, we added another Level 1 predictor (X₂). We set the correlation between the two Level 1 predictors (r_X1X2) to be 0.00, 0.50, and 0.90. The slope mean for X₂ was also set to 0.30, and the slope variance was set to either 0.00 or 0.10. We set the correlation between the two slopes to be .00, .50, or .90. When the slope variance for X₁ or X₂ was 0.00, the correlation between slopes was fixed to zero. Since Level 1 predictors can also contain between group variance, we set the ICC for each predictor to be .00 or .50. The latter value represents a substantial ICC value. For each possible combination, we simulated 1000 samples resulting in 3,168,000 samples.

Data Generation

The Level 2 variables Z and u_0j were generated as in Study 1. For u_1j, the standard deviation was set to the value to produce the desired slope variance. For u_2j, the standard deviation was set to the value of either 0.00 or the square root of 0.10 depending on the condition. We transformed these values to produce the desired covariance matrix for random effects using procedures adapted from nonparametric bootstrapping (Goldstein, 2011). For X₁ and X₂, we randomly sampled values from a standard normal distribution to produce the desired within group variance. Between group variance for X₁ and X₂ was created by adding the desired between group variances to the X₁ and X₂ predictors. In addition, the Level 1 error term, r_ij was created as in Study 1. Values for Y were generated by multiplying the coefficients by the relevant predictors and adding the error terms.

Analysis

As in Study 1, we estimated multiple models. In Model 1, the slope for X₁ was fixed and the slope for X₂ was random. In Model 2, the X₁ slope was random and the X₂ slope was fixed. Finally, Model 3 specified both slopes as random. We compared the Model 1 standard errors for the fixed effects against those from Model 3 to determine the effects of misspecifying the slope of X₁ on the standard errors of γ₀₁, γ₁₀, γ₁₁, and γ₂₀. We also compared Model 2 with Model 3 to determine if misspecifying the slope for X₂ affected the standard errors for γ₁₁.

Comparing Model 1 to Model 3

Table 3 displays the bias in standard errors as a function of comparing Model 1 with Model 3. Overall, the results for the standard errors for γ_01, γ_10, and γ₁₁ were similar to the one-predictor case. Misspecifying the slope of X₁ had little effect on the standard errors of γ₀₁ with bias estimates ranging from .000 to .002 and Cohen’s d’s ranging from 0.01 to 0.07. There was considerable bias in the standard errors for both γ₁₀ and γ₁₁ when the slope of X₁ was incorrectly fixed to zero. Bias estimates ranged from -.016 to -.005 with Cohen’s d’s ranging from -0.58 to -0.16 for the γ₁₀ standard errors. Similarly, bias estimates for the γ₁₁ standard errors ranged from -.004 to -.016 with Cohen’s d’s ranging from -0.70 to -0.21. As with Study 1, the effects of misspecifying the τ₁₁ as random when it should have been fixed to zero were not as pronounced. Bias estimates for standard errors were .002 for γ₁₀ and .001 for γ₁₁. Cohen’s d’s were 0.08 and 0.07, respectively.

It is interesting to note that varying the amount of slope variance produced little bias in the standard errors for γ₂₀. In general, the amount of bias was similar to that for the standard errors of γ₀₁. In addition, the degree of correlation between slopes did not appear to impact the bias in the standard errors for γ₁₀ and γ₁₁. For example, when the correlation was .00, the bias in the standard errors for γ₁₀ was -.009 with a Cohen’s d of -0.28. When the correlation equaled .90, the bias equaled -.010 with a Cohen’s d of -0.32. The amount of between group variance in the Level 1 predictors had little effect on the bias in the standard errors. Bias estimates for the standard errors of γ₁₀ were -.010 with a Cohen’s d of -0.30 when there was no between group variance in X₁, and -.009 with a Cohen’s d of -0.29 when between group variance was 0.5. We found similar results for γ₁₁. The between group variance conditions for X₂ produced biases and effect sizes that mimicked the findings for X₁.

We conducted two ANOVA’s investigating the impact of the study characteristics on the bias in standard errors for the average γ₁₀ and γ₁₁. As with Study 1, only conditions where there was variance in the slopes were included. We limited the ANOVA to the main effects and all two-way interactions because of the large number of conditions. Table 4 presents the results. Results were generally consistent with those for both ANOVA’s. The largest effects were associated with the amount of slope variance (partial η² = .43 for γ₁₀ and .49 for γ₁₁) and the number of groups (partial η² = .27 for both γ₁₀ and γ₁₁). None of the two-way interactions a partial η²s greater than .06.

Comparing Model 2 to Model 3

Varying the slope variance for X₂ (τ₂₂) when it should have been fixed did not affect the standard errors for γ₁₀ or γ₁₁. The bias estimates were -.001 and .000 with a Cohen’s d of -0.03 and 0.01 for the standard errors for γ₁₀ and γ_11, respectively. Similarly, fixing the slope variance for X₂ (τ₂₂) when it should have been random also did not have a large effect. The bias estimates were .004 and .003 with a Cohen’s d of 0.11 and 0.12 for γ₁₀ and γ_11, respectively. This provides further evidence that the effects of misspecifying the random part of multilevel models are isolated to the standard errors for the fixed effects associated with the misspecified slope variance.