2.1 Notation

We first set notation. To aid the reader, Appendix A.1 describe (i) symbols used in our notation (Table 1) and (ii) abbreviations we use relating to models (Table 2).

Let $g=1, \dots , G$ index the group. We index vectors belonging to group g with the subscript g, for example, the outcome for all units in group g is given by the vector $Y_g$ . The ith unit in group g is then indexed by $g[i]$ , for example, unit i in group g has outcome $Y_{g[i]}$ . This notation reminds the reader that unit i is contained in group g. Each group g has size $n_g$ with $N = \sum _{g=1}^{_G} n_g$ . Let $X_{g[i]}$ be the p-dimensional vector of covariates, including an intercept term, with an associated coefficient vector $\beta $ . One element of $X_{g[i]}$ in particular will be regarded as a treatment in settings described here. $X_g$ is the matrix of $X_{g[i]}$ for group g, and X is the matrix of $X_{g[i]}$ for the entire sample.

$$ \begin{align*} X_{g[i]} &= \begin{bmatrix} 1 \\ X_{g[i]}^{(1)} \\ \vdots \\ X_{g[i]}^{(p-1)} \end{bmatrix} \in \mathbb{R}^{p} \ , \ X_g = \begin{bmatrix} X_{g[1]}^{\top} \\ \vdots \\ X_{g[n_g]}^{\top} \end{bmatrix} \in \mathbb{R}^{n_g \times p} \ , \ X = \begin{bmatrix} X_{1} \\ \vdots \\ X_{_G} \end{bmatrix} \in \mathbb{R}^{N \times p} \ , \ \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{p-1} \end{bmatrix} \in \mathbb{R}^{p} \end{align*} $$

Similarly, $Z_{g[i]}$ is a d-dimensional vector of covariates, often containing a subset of the covariates in $X_{g[i]}$ , and possibly an intercept which will later function as an indicator of membership to group g. For each group g, the $Z_{g[i]}$ have an associated coefficient vector $\gamma _g$ . $Z_g$ is the matrix of $Z_{g[i]}$ for group g, Z is a block diagonal matrix of the $Z_g$ , and $\gamma $ stacks the set of $\gamma _g$ into a matrix.

$$ \begin{align*} Z_{g[i]} &= \begin{bmatrix} Z_{g[i]}^{(0)} \\ \vdots \\ Z_{g[i]}^{(d-1)} \end{bmatrix} \in \mathbb{R}^{d}, \ Z_g = \begin{bmatrix} Z_{g[1]}^{\top} \\ \vdots \\ Z_{g[n_g]}^{\top} \end{bmatrix} \in \mathbb{R}^{n_{g} \times d}, \ Z = \begin{bmatrix} Z_{1} & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & Z_G \end{bmatrix} \in \mathbb{R}^{N \times Gd}, \ \end{align*} $$

$$ \begin{align*} \gamma_g = \begin{bmatrix} \gamma_{0 g} \\ \vdots \\ \gamma_{(d-1) g} \end{bmatrix} \in \mathbb{R}^{d}, \ \gamma = \begin{bmatrix} \gamma_1 \\ \vdots \\ \gamma_G \end{bmatrix} \in \mathbb{R}^{Gd} \end{align*} $$

As noted, $Y_{g[i]}$ is the outcome of interest, and $\epsilon _{g[i]}$ is its associated residual/error term. Y and $\epsilon $ are $N \times 1$ vectors containing $Y_{g[i]}$ and $\epsilon _{g[i]}$ for the entire sample.

$$ \begin{align*} & Y_{g[i]} \in \mathbb{R} \ , \ Y_g = \begin{bmatrix} Y_{g[1]} \\ \vdots \\ Y_{g[n_g]} \end{bmatrix} \in \mathbb{R}^{n_g} \ , \ Y = \begin{bmatrix} Y_{1} \\ \vdots \\ Y_{G} \end{bmatrix} \in \mathbb{R}^{N} \ , \ \epsilon_{g[i]} \in \mathbb{R} \ , \ \epsilon_g = \begin{bmatrix} \epsilon_{g[1]} \\ \vdots \\ \epsilon_{g[n_g]} \end{bmatrix} \in \mathbb{R}^{n_g} \ , \ \epsilon = \begin{bmatrix} \epsilon_{1} \\ \vdots \\ \epsilon_{G} \end{bmatrix} \in \mathbb{R}^{N} \end{align*} $$

2.2 Fixed effect and multilevel models

Varying intercepts: Group fixed effects and random intercept models

We focus mainly on uses of FE and MLM that allow each group in the data a different intercept but no other group-varying coefficients.Footnote ⁴ Suppose data are generated according to the simple model, written in matrix form,

(1) $$ \begin{align} Y = X \beta + Z \gamma + \epsilon, \end{align} $$

where Z contains only indicators for membership to each group g, that is, each $Z_g$ is a vector of ones, meaning that only intercepts may differ between groups. While useful to write in this form because it shows the role of Z, the model is often expressed as

(2) $$ \begin{align} Y_{g[i]} = X^{\top}_{g[i]} \beta + \gamma_{g} + \epsilon_{g[i]}, \end{align} $$

where $\gamma _g$ are group-specific deviations from the overall intercept in $\beta $ . Although a model can be fitted regardless of the correlation between the residuals and $X_{g[i]}$ , investigators are usually interested in understanding the effect of a treatment variable (in $X_{g[i]}$ ) on $Y_{g[i]}$ , which requires the “conditional independence assumption,” ${\mathbb{E}}(\epsilon _{g[i]} \ | \ X, Z)=0$ . This is an identification assumption relating to the form of confounding, although it is also subject to model specification. We describe the types of permitted confounding that are used to justify this assumption in Section 2.3.

The distinction between FE and MLM for the varying intercepts model in Equation (2) relates to the assumptions they employ during estimation. Both FE and MLM regard $\beta $ as “fixed,” that is, nonrandom with no distributional assumption. However, FE and MLM diverge in their handling of $\gamma _g$ . In FE, the $\gamma _g$ are regarded as fixed, like $\beta $ . FE thus estimates $\beta $ and $\gamma $ by including indicator variables for group membership (in Z) as additional regressors in an OLS of Y on X, or by equivalent demeaning/partialing-out procedures. We refer to this as “group fixed effects” (Group-FE) here. In practice, one group indicator and its corresponding $\gamma _g$ must be dropped, or the intercept term must be dropped from $X_{g[i]}$ .

MLM, by contrast, treats the $\gamma _g$ as “random effects,” meaning that each $\gamma _g$ is estimated under the assumption that it may have a distribution, that is has nonzero variance rather than being a fixed quantity. This has implications for estimating both $\gamma $ and $\beta $ , as well as for constructing variance estimates, detailed in Section 2.4.Footnote ⁵ Specifically, we define the “random intercept” (RI) model,

(RI) $$ \begin{align} Y_{g[i]} = X^{\top}_{g[i]} \beta + \gamma_g + \epsilon_{g[i]}, \quad \gamma_g \ | \ X, Z \overset{iid}{\sim} N(0, \omega^2) \end{align} $$

where we also assume that for all $g,g'$ , and i. Additionally, the conditional independence assumption is elaborated to require a multivariate-normal distribution for the residual vectors in each group, $\epsilon _g \ | \ X, Z \overset {ind}{\sim } \mathcal {N}(0, \Sigma _g)$ where $\Sigma _g \in \mathbb{R}^{n_g \times n_g}$ is group g’s error covariance matrix. All $\gamma _g$ can be kept and estimated by this approach, together with an intercept in $\beta $ . It is also commonly assumed that the $\epsilon _{g[i]}$ are spherical (i.e., $\Sigma _g = \sigma ^2 I_{n_g}$ ), although we avoid that restriction unless otherwise noted.

General fixed effect and multilevel models

Although we focus mainly on Group-FE and RI models here, the claims made below pertain to FE and MLM in their more general form,

(3) $$ \begin{align} Y_{g[i]} = X^{\top}_{g[i]} \beta + Z_{g[i]}^{\top} \gamma_g + \epsilon_{g[i]}, \end{align} $$

where the $\gamma _g$ represent group-specific coefficients for the variables in $Z_g$ , which may now include variables other than group indicators as above. As before, FE estimates both $\beta $ and $\gamma $ by OLS regression of Y on X and Z. For identifiability, the covariates that $X_{g[i]}$ and $Z_{g[i]}$ share (or that result in perfect colinearity) are either dropped from $Z_{g[i]}$ for one g, or dropped from $X_{g[i]}$ .

When instead fit with random effects in MLM, the coefficients $\gamma _g$ are estimated with additional distributional assumptions, specifically

(MLM) $$ \begin{align} Y_{g[i]} = X^{\top}_{g[i]} \beta + Z_{g[i]}^{\top} \gamma_g + \epsilon_{g[i]}, \quad \gamma_g \ | \ X, Z \overset{iid}{\sim} \mathcal{N}(0, \Omega), \end{align} $$

where $\Omega \in \mathbb{R}^{d \times d}$ and we assume for all $g,g'$ , and i.Footnote ⁶ In addition, we define $ {{\rm var}}(\epsilon \ | \ X,Z) = \Sigma $ to be the ordered block diagonal matrix of $\Sigma _g$ , and $\Omega _{\text {block}} = {{\rm var}}(\gamma \ | \ X,Z)$ to be the block diagonal matrix of $\Omega$ ,

$$ \begin{align*} \Sigma = \begin{bmatrix} \Sigma_1 & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & \Sigma_G \end{bmatrix} \in \mathbb{R}^{N \times N}, \ \Omega_{\text{block}} = \begin{bmatrix} \Omega & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & \Omega \end{bmatrix} \in \mathbb{R}^{Gd \times Gd} \end{align*} $$

2.3 Identification: From nonparametric conditions to specification requirements

Let us assume there is no within-group confounding. In longitudinal settings, this is to say there are no unobserved time-varying confounders. Such an assumption ensures non-parametric identifiability, that is, the causal effect of the treatment can be identified conditionally on group and then averaged together as desired across groups.Footnote ⁷

For such nonparametric identification to be sufficient for unbiased estimation, however, would require the ability to condition on group non-parametrically (i.e., estimate the relationship between treatment and outcome within each group separately). In most settings the investigator is unable to do this, and so turns to modeling assumptions that instead “account for” group in a specific model. Because of this model dependence, identification of treatment effects then requires additional assumptions related to the specification of those models.

The traditional motive for Group-FE is an assumption that the only source of confounding is group-level confounding that takes the linear form of Equation (2): a constant, additive shift by group (i.e., the $\gamma _g$ ). The conditional independence assumption needed is ${\mathbb{E}}(\epsilon _{g[i]} \ | \ X, Z)=0$ (e.g., Greene Reference Greene2003). The Group-FE approach is powerful precisely because including Z fully purges the group-level intercepts, $\gamma _g$ , from this residual, thus removing confounding. So long as other (within-group) forms of confounding do not exist, this unbiasedly estimates $\beta $ .

The central concern with RI is that even under the conditional independence assumption, random effect estimation fails to remove this confounding. As we explain in Section 3.1, this is because RI does not fully account for the $\gamma _g$ , biasing the estimate of $\beta $ as well. Avoiding this bias requires additionally that the “uncorrelated random effects” assumption be true: that the $\gamma _g$ are uncorrelated with $X_{g[i]}$ , so that failure to account for $\gamma _g$ does not bias estimates of $\beta $ .Footnote ⁸ ${{\hspace {-0.0001pt}}}^{,}$ Footnote ⁹ This is problematic, since $\gamma _g$ is most important to account for when it is correlated with the treatment variable of interest in $X_{g[i]}$ and thus acts as a confounder. While this assumption, also referred to simply as the “random effects assumption” (Bell and Jones Reference Bell and Jones2015; Kim and Steiner Reference Kim and Steiner2019), is well-known in principle, it remains widely neglected in practice (see Section 3.1) despite attractive solutions (see Section 3.2).

2.4 Parameter estimation in MLM

We begin by briefly reviewing how MLM parameters are estimated.Footnote ¹⁰ Estimation proceeds in three steps: (1) estimate $\Omega $ and $\Sigma $ , (2) estimate $\beta $ using the estimate of $(\Omega , \Sigma )$ , and (3) estimate $\gamma $ using the estimate of $(\beta , \Omega , \Sigma )$ . The first two steps require the distribution of Y given X and Z. Because, in MLM, the $Z_{g[i]}^{\top } \gamma _g$ are mean-zero given X and Z, the $(Z_{g[i]}^{\top } \gamma _g + \epsilon _{g[i]})$ can be treated as combined mean-zero error terms. This allows one to formulate MLM on the sample-level as $Y = X \beta + \epsilon ^{*}$ where $\epsilon ^{*} = Z \gamma + \epsilon $ . Then, because $\gamma $ and $\epsilon $ are both normally distributed given X and Z, we have

(4) $$ \begin{align} Y \ | \ X, Z \sim \mathcal{N}(X \beta, V) \ \ \ \text{where} \ \ \ V = {{\rm var}}(\epsilon^{*} \ | \ X, Z) = Z \Omega_{\text{block}} Z^{\top} + \Sigma. \end{align} $$

The likelihood is then

(5) $$ \begin{align} L(\beta, \Omega, \Sigma \ | \ Y, X, Z) = p(Y \ | \ X, Z, \beta, \Omega, \Sigma) \propto |V| ^{-1/2} \text{exp} \left( - \frac{1}{2} (Y - X \beta)^{\top} V^{-1} (Y - X \beta) \right). \end{align} $$

Given a choice of $(\Omega , \Sigma )$ , which determines V, maximizing the likelihood for $\beta $ would yield

(6) $$ \begin{align} \hat{\beta} (\Omega, \Sigma) = (X^{\top} V^{-1} X)^{-1} X^{\top} V^{-1} Y. \end{align} $$

However, estimates of $\Omega $ and $\Sigma $ do not enjoy similarly simple closed solutions. Instead, they are typically found iteratively through either unrestricted maximum likelihood estimation or restricted maximum likelihood estimation, both using Equation (5). $\beta $ can then be estimated by plugging estimates of $\Sigma $ and $\Omega $ into Equation (6). Finally, $\gamma $ is estimated by maximizing its posterior probability given Y, X, Z, and the estimated $(\beta , \Omega , \Sigma )$ , that is,

(7) $$ \begin{align} \hat{\gamma} (\beta, \Omega, \Sigma) = {\arg\,\max}_{\gamma} \ p(\gamma \ | \ Y, X, Z, \beta, \Omega, \Sigma) = \Omega_{\text{block}} Z^{\top} V^{-1} (Y - X \beta). \end{align} $$

In this article, we focus solely on estimation and inference for $\beta $ and $\gamma $ . Because MLM’s estimates of these parameters are functions of $\Omega $ and $\Sigma $ , our results hold regardless of the choice between unrestricted or restricted maximum likelihood in estimating $\Omega $ and $\Sigma $ . For this reason, we refer to arbitrary MLM estimates of the parameters as $\hat {\Omega }_{\text {MLM}}$ , $\hat {\Sigma }_{\text {MLM}}$ , $\hat {\beta }_{\text {MLM}}$ , and $\hat {\gamma }_{\text {MLM}}$ .

Finally, treating $\hat {V}_{\text {MLM}}$ as fixed, the conditional variance of $\hat {\beta }_{\text {MLM}}$ is simply

(8) $$ \begin{align} {{\rm var}}(\hat{\beta}_{\text{MLM}} \ | \ X, Z) = (X^{\top} \hat{V}_{\text{MLM}}^{-1} X)^{-1} X^{\top} \hat{V}_{\text{MLM}}^{-1} {{\rm var}}(Y \ | \ X, Z) \hat{V}_{\text{MLM}}^{-1} X (X^{\top} \hat{V}_{\text{MLM}}^{-1} X)^{-1}, \end{align} $$

which is commonly simplified to $(X^{\top } \hat {V}_{\text {MLM}}^{-1} X)^{-1}$ as the standard model-based MLM variance estimator, by assuming that $\hat {V}_{\text {MLM}}$ is a consistent estimator of ${{\rm var}}(Y \ | \ X, Z)$ .

3.1 Random effects as regularization and the “incomplete conditioning” problem

The first piece shows an equivalence between the random effects employed in MLM and regularization, that is, a penalized fitting approach. We begin by introducing the “regularized fixed effects” (regFE) class of models. Suppose we are interested in optimal out-of-sample generalization. Then, instead of estimating Equation (2) by minimizing (in-sample) squared error, we minimize the squared error plus a penalty proportional to the sum of squared $\gamma _g$ values, known as $L_2$ or Tikhonov regularization:

(9) $$ \begin{align} (\hat{\beta}_{\text{regFE}}, \hat{\gamma}_{\text{regFE}}) = \underset{\beta, \gamma}{{\arg\,\min}} \bigg(\sum_{g=1}^{G} \sum_{i=1}^{n_g} [Y_{g[i]} - X^{\top}_{g[i]} \beta - \gamma_{g}]^2 + \lambda \sum_{g=1}^{G} \gamma_{g}^2\bigg), \end{align} $$

where $\lambda>0$ determines the extent of regularization on $\gamma _g$ —larger values shrink each $\gamma _g$ toward to 0, while smaller values yield estimated intercepts closer to Group-FE’s estimates—and may be chosen by various means such as some form of cross-validation.Footnote ¹¹ The key result is that for a particular choice of $\lambda $ , such a regularized model is equivalent to the RI model,

We omit the proof here, as Theorem 3.2 below offers a generalization. That MLM provides “shrinkage” or “partial-pooling” estimates is very well-known (e.g., Steenbergen and Jones Reference Steenbergen and Jones2002; Fitzmaurice, Laird, and Ware Reference Fitzmaurice, Laird and Ware2004; Luke Reference Luke2004; Gelman Reference Gelman2006), and in recent texts MLM has even been referred to as employing a “regularizing prior” (e.g., McElreath Reference McElreath2018). The equivalence between regularization and holding a prior on the regularized coefficients is also well known.Footnote ¹² That said, because we did not find any explicit formalization of the equivalence of estimation in MLM to regularization with a specific choice of $\lambda $ in any textbooks or articles we reviewed, we fill this gap.

We demonstrate with a simple simulation example, drawing 1,000 datasets from the following data generating process (DGP) with $G=25$ and $n_g = 15$ each time:

(10) $$ \begin{align} Y_{g[i]} = \beta_0 + \beta_1 X_{g[i]} + \gamma_g + \epsilon_{g[i]} \ \ \ \text{where} \ \ \ X_{g[i]} \overset{iid}{\sim} N(0, 1) \ , \ \gamma_{g} \overset{iid}{\sim} N(0, 1), \end{align} $$

where $\lambda $ for regFE is chosen by cross-validation and $(\hat {\sigma }^2_{\text {RI}}, \hat {\omega }_{\text {RI}}^2)$ for RI by restricted maximum likelihood.Footnote ¹³ We see in Figure 1 that the estimates of $\beta $ and $\gamma $ are nearly identical. They would be numerically equal if, instead of setting $\lambda $ in regFE by cross-validation, we chose $\lambda = \hat {\sigma }^2_{\text {RI}} / \hat {\omega }^2_{\text {RI}}$ as per Theorem 3.1. However, despite the different selection procedures, we also see from Figure 1 that the $\lambda $ from RI and regFE are quite similar. This is because while $\lambda = \hat {\sigma }^2_{\text {RI}}/\hat {\omega }^2_{\text {RI}}$ in RI is not made specifically with predictive accuracy in mind, the choice is a sensible one from a prediction point of view. To see this, note first that the portion of $Y_{g[i]}$ that is not explained by $X_{g[i]}$ is the sum of $\gamma _g$ and $\epsilon _{g[i]}$ , which have variances of $\omega ^2$ and $\sigma ^2$ , respectively. When, for example, the $\gamma _g$ have high variance in relation to that of $\epsilon _{g[i]}$ , $\lambda = \sigma ^2/\omega ^2$ will be small, which is appropriate since $\gamma _g$ will be helpful in predicting $Y_{g[i]}$ . The converse is true when the $\gamma _g$ have low variance relative to $\epsilon _{g[i]}$ , leading to a higher effective $\lambda $ and desirable shrinkage on $\gamma _g$ to avoid over-fitting.

Figure 1 (a) Boxplots of $\hat {\beta }_1$ from RI and regFE over 1,000 iterations. The dashed line indicates the true value for $\beta _1$ (i.e., $\beta _1=1$ ). The estimates from each method are correlated at $0.999$ . (b) Example for the estimated intercepts from the RI model ( $\hat {\gamma }_{\mathrm {RI}}$ ) or from regFE ( $\hat {\gamma }_{\mathrm {regFE}}$ ) in a single iteration of the simulation. The points correspond to the coordinates $(\hat {\gamma }_{g, \mathrm {RI}}, \hat {\gamma }_{g, \mathrm {regFE}})$ across g, and the dashed line indicates equality. The RI and regFE estimates have a correlation of $0.999$ . (c) Plot of $\lambda $ from RI and regFE over 1,000 iterations. At each iteration, $\lambda $ is chosen by cross-validation for regFE, and $\lambda = \hat {\sigma }^2_{\mathrm {RI}} / \hat {\omega }^2_{\mathrm {RI}}$ for RI. The dashed line indicates equality. The $\lambda $ from RI and regFE have means of 4.278 and 4.277, respectively, and are correlated at 0.910.

More generally, consider estimating Equation (3) with the regularized regression:

(11) $$ \begin{align} \begin{aligned} &(\hat{\beta}_{\text{regFE}}, \hat{\gamma}_{\text{regFE}}) = \underset{\beta, \gamma}{{\arg\,\min}} \bigg(\sum_{g=1}^{G} \sum_{i=1}^{n_g} [Y_{g[i]} - X^{\top}_{g[i]} \beta - Z_{g[i]}^{\top} \gamma_g]^2 + \sum_{g=1}^G \gamma_g^{\top} \Lambda \gamma_g \bigg), \\ \text{where} \ \ \ & \Lambda \in \mathbb{R}^{d \times d} \text{ is symetric and positive semi-definite}. \end{aligned} \end{align} $$

Again, $\Lambda $ may be obtained by various means, such as cross-validation.Footnote ¹⁴ The equivalence with MLM is then given in Theorem 3.2.

Proof of Theorem 3.2 is given in Appendix A.3 and proves Theorem 3.1 as a special case with $Z_{g[i]} = [1]$ . This equivalence is useful in comprehending several of MLM’s most central advantages, and limitations. First, the equivalence to regularization immediately explains why MLM can produce more accurate (out-of-sample) outcome predictions than does FE: regularization prevents the over-fitting that FE may allow, particularly in the case of small groups.Footnote ¹⁵

Further, MLM is sometimes understood to “wisely” adapt the level of shrinkage based on group size, but the comparison to regularization shows that such adaptation is not as sophisticated as it may appear. When MLM or regFE encounters a group that appears to have a very high or low mean relative to others, choosing an extreme $\gamma _g$ would decrease the sum of squared errors for that group, $\sum _{i=1}^{n_g} (Y_{g[i]} - X^{\top }_{g[i]} \beta - Z_{g[i]}^{\top } \gamma _g)^2$ . For a relatively small group, the savings in squared loss would be smaller relative to the additional “cost” paid through the regularization penalty, $ \gamma _g^{\top } \Lambda \gamma _g$ , and so $\gamma _g$ will be left relatively close to zero. By comparison, when MLM or regFE encounters a larger group, the savings in terms of squared loss would justify the added regularization penalty, so the choice of $\gamma _g$ minimizing the penalized sum of squared errors will be a more extreme one. Although this behavior appears to intelligently balance prior knowledge with allowing the data to speak, it is reproduced simply through regularization.

Another feature of MLM that can be understood through the regularization view is its ability to estimate coefficients for group-level variables even while including group-specific intercepts and slopes that would have prevented model identification under OLS. This occurs for the same reason that one can include more coefficients than observations in a “ridge regression.” Consider attempting to find $(\hat {\beta }, \hat {\gamma })$ purely by OLS,

(12) $$ \begin{align} (\hat{\beta}, \hat{\gamma}) = \underset{\beta, \gamma}{{\arg\,\min}} \left(\sum_{g=1}^{G} \sum_{i=1}^{n_g} [Y_{g[i]} - X^{\top}_{g[i]} \beta - Z_{g[i]}^{\top} \gamma_g]^2 \right) = \underset{\beta, \gamma}{{\arg\,\min}} \ \bigg| \bigg| Y - \begin{bmatrix} X & Z \end{bmatrix} \begin{bmatrix} \beta \\ \gamma \end{bmatrix} \bigg| \bigg|_2^2. \end{align} $$

This has no unique solution if $X_{g[i]}$ contains a group-level variable or its interaction with a variable in $Z_{g[i]}$ , as $[X \ Z]^{\top } [X \ Z]$ is then singular. However, by introducing the regularization penalty into the minimization problem, regFE essentially adds to $[X \ Z]^{\top } [X \ Z]$ a positive semi-definite matrix that allows the sum of the two matrices to be invertible.

Finally, this equivalence illuminates the main concern with MLM: bias when the random effects are correlated with $X_{g[i]}$ . Because the group-specific intercepts in a RI model are regularized, they do not achieve the values that would “fully absorb” group-specific confounding, leaving components unexplained that can instead be captured by biasing the coefficients on $X_{g[i]}$ . In Equation (2), where we hope to condition on group to absorb group-level confounders, random effects thus offer only “incomplete conditioning,” not fully accounting for the unobserved group-level variables’ influence on the outcome. To illustrate, consider the following DGP:

(DGP1) $$ \begin{align} & Y_{g[i]} = \beta_0 + \beta_1 X_{g[i]} + (W^{(1)}_g + W^{(2)}_g) + \epsilon_{g[i]} \nonumber \\ \text{where} \ \ \ & [W^{(1)}_g \ W^{(2)}_g]^{\top} \overset{iid}{\sim} \mathcal{N} (0, 2 I_2) \ , \ X_{g[i]} = W^{(1)}_g + N(0, 1)_{g[i]} \ , \ \epsilon_{g[i]} \overset{iid}{\sim} N (0, \sigma^2), \end{align} $$

where $X_{g[i]} \in \mathbb{R}$ is an observed observation-level variable and the $W^{(\ell )}_g$ are unobserved group-level covariates, with $W^{(1)}_g$ being a confounder. While Group-FE will unbiasedly estimate $\beta _1$ , a simple OLS regression of Y on X would produce a biased estimate of $\beta _1$ , having failed to account for $W^{(1)}_g$ . RI may seem more appealing than OLS because investigators might hope that the $\gamma _g$ will capture the group-level confounding. However, the shrinkage of $\gamma _g$ due to treating them as random effects yields estimated intercepts closer to zero than those obtained by Group-FE and required to account for $W^{(1)}_g$ . This leaves some of $W^{(1)}_g$ unabsorbed, allowing it to continue biasing $\hat {\beta }_1$ as in the OLS model. We call this “incomplete conditioning” because the intended analytical strategy required to estimate $\beta $ unbiasedly would condition on group, but the use of RI fails to achieve this.

Figure 2 illustrates this. Define bias and root mean square error (RMSE) as,

$$ \begin{align*} \text{Bias} = \frac{1}{M}\sum_{m=1}^{M} (\hat{\beta}^{(m)} - \beta), \quad \text{RMSE} = \sqrt{ \frac{1}{M} \sum_{m=1}^{M}(\hat{\beta}^{(m)}- \beta)^2} \end{align*} $$

where m indexes the number of iterations from 1 to M and $\hat {\beta }^{(m)}$ is the estimate from the mth iteration. Bias is large, as expected, for OLS as it fails to account for group-level confounding at all. RI makes almost no improvement when groups are small ( $n_g=5$ ), and only a partial improvement when the groups are quite large ( $n_g=50$ ). Note that unbiasedness of RI would require the absence of correlated random effects, that is no correlation between $X_{g[i]}$ and $W_g^{(1)}$ , which is to say an absence of group-level confounding. Had this been the case, OLS would also suffice, and would differ from RI only in its efficiency and how variance is estimated. By contrast, in the presence of such correlation, Group-FE effectively eliminates confounding bias at both group sizes. In terms of efficiency, while RI has the expected slight decrease in variance, its RMSE remains about twice that of Group-FE at either group size due to these biases. That RI’s average bias falls as $n_g$ rises may seem to be a cause for hope when one has a large enough dataset, and Lockwood, McCaffrey, et al. (Reference Lockwood and McCaffrey2007) describes the conditions under which this type of bias tends to 0 as $n_g \rightarrow \infty $ generally. However, in practice, there is no knowing if $n_g$ is large enough to ensure negligible bias in a given case, as this depends on the correlation between the covariates and the random effects. We also note that with group-level covariates, increasing $n_g$ may not alleviate bias (see Appendix A.4).

Figure 2 Comparison of estimates of $\beta _1$ from OLS, Group-FE, and RI in DGP1. Note: Results across 1000 iterations, each drawn from DGP 1 with $\beta _0 = \beta _1 = 1$ . The dashed-line represents the true $\beta _{\ell }$ . Due to correlated random effects, RI estimates are almost as biased as OLS estimates when group size is small (5). The bias is less severe but still appreciable at a group size of 50, and RMSE remains twice that of Group-FE.

3.2 Bias-corrected MLM

The second analytical piece is that while MLM’s bias problem with correlated random effects appears to be dire, there is a simple solution dating back to Mundlak (Reference Mundlak1978). With this fix, MLM unbiasedly estimates coefficients for variables below the level of grouping or clustering, like FE does, while retaining the ability to include group-level variables in the model and the superior predictive accuracy that arises from regularization/random effects. For RI, the fix requires adding to the model the group-level means of $X_{g[i]}$ (including any interactions or other nonlinear terms).

As a stepping stone, consider a similar approach that could be applied to OLS as a substitute for FE. We consider again DGP 1 from Section 3.1, but suppose we estimate the following model by OLS,

(13) $$ \begin{align} Y_{g[i]} = \beta_0 + \beta_1 X_{g[i]} + \alpha_1 \bar{X}_g + \epsilon_{g[i]}, \end{align} $$

where $\bar {X}_g =\frac {1}{n_g} \sum _{i=1}^{n_g} X_{g[i]}$ . Informally, adding $\bar {X}_g$ “soaks up” any contribution that $W_g^{(1)}$ could have made to $Y_{g[i]}$ that could be correlated with $X_{g[i]}$ , protecting the estimate of $\beta _1$ in the same way that FE would, and leaving an unbiased estimate of $\beta _1$ under the conditional independence assumption. Proof is given in Appendix A.5. This bias-curing effect of including $\bar {X}_g$ in OLS can similarly be applied to the RI model, and alleviates RI’s bias problem, with $\Sigma = \sigma ^2 I_N$ , in estimating $\beta _1$ . We omit the proof of this result, as we prove a more general claim below. In fact, an OLS model including $\bar {X}_g$ , the Group-FE model, and the RI model including $\bar {X}_g$ all produce exactly the same point estimates. Despite this equivalence, the RI model retains MLM’s greater (out-of-sample) predictive accuracy for the outcome compared to Group-FE—a benefit of the regularization imposed on $\gamma _g$ , as illustrated in Figure 3.

Figure 3 Outcome prediction error for debiased RI versus Group-FE. Note: Comparison of testing error for the predicted outcome (average standardized test MSE, $(N \mathbb{E}(\epsilon _{g[i]}^2))^{-1} \sum _{g, i} (Y_{g[i]} - \hat {Y}_{g[i]})^2$ ). The RI model with $\Sigma = \sigma ^2 I_N$ has been debiased by including $\bar {X}_g$ as a covariate, and shows lower testing error, especially when groups are smaller. Results are averaged across 1000 iterations, each drawn from DGP 1. Testing data are of the same size as the training data.

One can easily generalize this alteration to other models with varying intercepts but with multiple covariates or treatments: simply add to the RI model the group-level means of all included variables, $\bar {X}_g = \frac {1}{n_g} \sum _{i=1}^{n_g} X_{g[i]}$ , including any interactions or nonlinear transformations. If spherical observation-level errors are assumed, the coefficient estimates from this model for lower-level variables are exactly equal to those obtained by Group-FE. Historically, the inclusion of $\bar {X}_g$ was proposed by Mundlak (Reference Mundlak1978) and extended by Chamberlain (Reference Chamberlain1979, Reference Chamberlain1982), and is known in some econometrics-informed traditions as the “correlated random effects” (CRE) approach (Wooldridge Reference Wooldridge2010; Schunck Reference Schunck2013). Although originally motivated by imposing the assumption ${\mathbb{E}} (\gamma _g \ | \ X_{g}, Z_{g}) = \bar {X}_g^{\top } \alpha $ , we show unbiasedness without this assumption by showing its equivalence to “partialing out” of the group-level effects, which is in turn equivalent to Group-FE. A closely related approach is the hybrid model of Allison (Reference Allison2009), which group-demeans $X_{g[i]}$ in addition to including $\bar {X}_g$ . This model is an isomorphic variation on models that only include $\bar {X}_g$ without demeaning, producing the same coefficient of interest on X but altering the way in which coefficients are combined to interpret between-group differences (see Schunck Reference Schunck2013).

We now turn to the more general debiasing approach for cases allowing multiple random coefficients and not just group intercepts. We refer to this as “bias-corrected MLM” (bcMLM), and use this label throughout the remainder of the paper to cover the special case of bias-corrected RI as well.Footnote ¹⁶ We first take the projections of the “fixed effect variables” ( $X_g$ ), excluding the intercept, onto the random effect variables ( $Z_{g}$ ) within each group, $\tilde {X}_{g[i]} = Z_{g[i]}^{\top } (Z_g^{\top } Z_g)^{-1} Z_g^{\top } X_{g}$ .Footnote ¹⁷ These projections are then added to the regression as “fixed effect variables,”

(bcMLM) $$ \begin{align} Y_{g[i]} &= X_{g[i]}^{\top} \beta + \tilde{X}_{g[i]} \alpha + Z_{g[i]}^{\top} \gamma_g + \epsilon_{g[i]}, \quad \gamma_g \ | \ X, Z \overset{iid}{\sim} \mathcal{N}(0, \Omega) \end{align} $$

where $\beta $ and $\alpha $ are assumed fixed, and we continue to assume that $\epsilon _g \ | \ X, Z \overset {ind}{\sim } \mathcal {N}(0, \Sigma _g)$ and for all $g,g'$ , and i. When $\Sigma =\sigma ^2 I_N$ , and provided it is not overidentified (described below), bcMLM produces estimates for $\beta $ that are unbiased under the conditional independence assumption (Appendix A.8) and identical to FE estimates (Appendix A.9).

We make two remarks on this result. First, it provides unbiasedness only when spherical observation-level errors are assumed.Footnote ¹⁸ While the choice of variance-covariance matrix for the errors has no impact on point estimates under FE, the point estimates from MLM are sensitive to this choice. Second, the RI model including $\bar {X}_g$ noted above is a special case of bcMLM: if $Z_g = \vec {\mathbb {1}}_{n_g}$ , then $\tilde {X}_{g[i]} = [1] ( \vec {\mathbb {1}}_{n_g}^{\ \top } \vec {\mathbb {1}}_{n_g})^{-1} \vec {\mathbb {1}}_{n_g}^{\ \top } X_{g} = \bar {X}_g^{\top } $ .

3.3 Variance estimation

The third and final analytical piece is that the “default” MLM standard errors are based on narrow assumptions that are often inappropriate, but this too can be remedied. A commonly cited motive for using MLM is the claim that it ensures appropriate standard errors for multilevel data (e.g., Luke Reference Luke2004; Morgan and Kelly Reference Morgan and Kelly2017; Beazer and Blake Reference Beazer and Blake2018; Campbell and Ronfeldt Reference Campbell and Ronfeldt2018; Rueda Reference Rueda2018; Clements et al. Reference Clements, Sarama, Baroody, Joswick and Wolfe2019; Pardos-Prado and Xena Reference Pardos-Prado and Xena2019). This is only true when the sole source of heteroskedasticity or dependency between residuals arises due to the random effects included by $Z_{g[i]}$ and their contributions to $Y_{g[i]}$ . That a default standard error exists for MLM, and is the only choice in some software, does not imply it is always an appropriate choice.

Recall from Section 2.4 that the random effects contribution, $Z_{g[i]}^{\top } \gamma _g$ , and the idiosyncratic error, $\epsilon _{g[i]}$ , can be thought of as a single mean-zero error term $\epsilon ^{*}_{g[i]} = Z_{g[i]}^{\top } \gamma _g + \epsilon _{g[i]}$ in the model $Y_{g[i]} = X^{\top }_{g[i]} \beta + \epsilon ^{*}_{g[i]}$ . The dependency between two observations’ outcomes within the same group, conditional on X and Z, is then

(14) $$ \begin{align} {\text{cov}}(Y_{g[i]}, Y_{g[i']} \ | \ X, Z) = {\text{cov}}(\epsilon^{*}_{g[i]}, \epsilon^{*}_{g[i']} \ | \ X, Z) = Z_{g[i]}^{\top} \Omega Z_{g[i']} + {\text{cov}}(\epsilon_{g[i]}, \epsilon_{g[i']} \ | \ X, Z). \end{align} $$

MLM can be understood as a framework for structuring this covariance, by specifying which variables enter the models as random effects (i.e., $Z_{g[i]}$ ) and parameterizing $\Omega $ and $\Sigma $ . In the RI model, with $Z_{g[i]} = [1]$ , the combined error, $\epsilon ^*_{g[i]}$ , is $\gamma _g + \epsilon _{g[i]}$ . Under the conditional independence assumption and spherical errors, $\epsilon _{g[i]} \ | \ X, Z \overset {iid}{\sim } N(0, \sigma ^2)$ , the dependence structure is

(15) $$ \begin{align} {{\rm var}}(Y_{g[i]} \ | \ X, Z) = \omega^2 + \sigma^2 \quad\text{and}\quad{\text{cov}}(Y_{g[i]}, Y_{g[i']} \ | \ X, Z) = \omega^2 \ \ \ \text{for} \ \ \ i \neq i' \end{align} $$

In other words, $Y_{g[i]}$ is modeled as linear in $X_{g[i]}$ with error, but instead of independent observations, there is constant covariance between observations in the same group, and this covariance does not differ by group. This yields a compound symmetric covariance matrix for each group’s error vector, $\epsilon _g^{*} = Y_g - X_g \beta $ . If $Z_{g[i]} = [1 \ X^{(1)}_{g[i]}]^{\top }$ , then $\epsilon ^*_{g[i]}=\gamma _{0g} + \gamma _{1g} X^{(1)}_{g[i]} + \epsilon _{g[i]}$ . Maintaining that $\epsilon _{g[i]} \ | \ X, Z \overset {iid}{\sim } N(0, \sigma ^2)$ and that (conditionally on X and Z) $\gamma _{0g}$ and $\gamma _{1g}$ are drawn from $N(0, \omega _0^2)$ and $N(0, \omega _1^2)$ with covariance $\omega _{01}$ yields intragroup covariances of

(16) $$ \begin{align} \begin{aligned} {{\rm var}}(Y_{g[i]} \ | \ X, Z) &= \omega_0^2 + 2 X^{(1)}_{g[i]} \omega_{01} + \left[X^{(1)}_{g[i]}\right]^2 \omega_1^2 + \sigma^2 \\ \text{and} \ \ \ {\text{cov}}(Y_{g[i]}, Y_{g[i']} \ | \ X, Z) &= \omega_0^2 + \left[X^{(1)}_{g[i]} + X^{(1)}_{g[i']}\right] \omega_{01} + X^{(1)}_{g[i]} X^{(1)}_{g[i']} \omega_1^2 \ \ \ \text{for} \ \ \ i \neq i' \end{aligned} \end{align} $$

With the addition of more random effect variables, the variance structure becomes more complex. This complexity should not be equated with generality—the variance is still assumed to be a highly prescribed function of the data. To illustrate the potential for variance estimates with poor coverage, consider the following longitudinal DGP, where g indexes the person and $t=1, \dots , T$ indexes the time-point of the observation:

(DGP2) $$ \begin{align} & Y_{g[t]} = \beta_0 + \beta_1 X_{g[t]} + \beta_2 U_g + W_g + \epsilon_{g[t]} \nonumber \\ \text{where} \ \ \ & W_g \overset{iid}{\sim} N (0, 4), \nonumber\\& X_{g[t]} \sim N (0, 1) \quad\text{and} \quad {\text{cor}}(X_{g[t]}, X_{g[t+k]}) = (0.75)^{k},\nonumber\\ & U_g \overset{iid}{\sim} N (0, 1),\; \epsilon_{g[i]} \sim N (0, U_g^2 \sigma^2) \quad\text{and} \quad {\text{cor}}(\epsilon_{g[t]}, \epsilon_{g[t+k]}) = (0.75)^{k} \end{align} $$

Here, there is an observation-level variable $X_{g[t]}$ that is auto-correlated; a group-level variable $U_g$ ; and a random intercept $W_g$ . The observation-level error terms are auto-correlated as well with (heteroskedastic) variances that depend on $U_g$ . The correct dependence structure would be

(17) $$ \begin{align} {\text{cov}}(Y_{g[t]}, Y_{g[t+k]} \ | \ X, Z) = U_g^2 \sigma^2 (0.75)^{k} + 4 \end{align} $$

Using the “default” variance with RI, assuming $\Sigma =\sigma ^2 I_N$ , Figure 4 shows coverage rates for nominal 95% confidence intervals of $\beta _1$ and $\beta _2$ across draws from DGP 2. Typically we would focus interest on $\beta _1$ , motivated by an assumption of no within-group confounding. However, we also show results for $\beta _2$ because group-level variables were often of interest in the empirical works we reviewed. The RI standard errors are consistently too small for both $\beta _1$ and $\beta _2$ across all sample sizes. The undercoverage for $\beta _1$ worsens as the total number of time-periods (T) increases. Coverage improves for $\beta _2$ as T increases, but remains unsatisfactory even at $T=50$ . Similar undercoverage of the “default” MLM standard errors for RI has been noted by Bell, Fairbrother, and Jones (Reference Bell, Fairbrother and Jones2019), Heisig, Schaeffer, and Giesecke (Reference Heisig, Schaeffer and Giesecke2017), and Jacqmin-Gadda et al. (Reference Jacqmin-Gadda, Sibillot, Proust, Molina and Thiébaut2007).

Relaxing assumptions for MLM variance estimation

If the user has strong reason to believe errors (conditionally on the random effect contributions) are spherical, then the default MLM standard errors just described would be appropriate. However, such justifications are rarely offered. Fortunately, more flexible approaches can be employed in the MLM framework, relaxing this assumption. For example, with longitudinal data, it is possible to assume an AR(1) error structure for the $\epsilon _{g[t]}$ , in which ${\text {cor}}(\epsilon _{g[t]}, \epsilon _{g[t + k]}) = \rho ^{k}$ for $\rho \in (-1, 1)$ . Or, one may allow ${{\rm var}}(\epsilon _{g[i]} \ | \ X, Z) = \sigma _g^2$ , where $\sigma ^2_g$ can differ by group, to accommodate heteroskedasticity by group. One may also allow a common unstructured $\Sigma _g$ across g, which makes no assumptions on the intragroup covariance and assigns a separate parameter to each ${\text{cov}}(\epsilon _{g[i]}, \epsilon _{g[i']})$ .Footnote ²¹

Using such specialized error structures when estimating standard errors may be advisable when researchers can defend the corresponding assumptions. On the other hand, users often cannot claim to know the correct variance structure based on theory alone. Fortunately, cluster-robust standard errors (CRSE), popular in conjunction with FE, provide a useful low-assumption alternative by asking the user only to assume independence across groups while allowing within group covariance to be fully estimated, albeit at the cost of requiring more data.Footnote ²²

Note that both MLM and CRSEs operate as though we assume zero covariance of the residuals from units in different groups, sharing the assumption

(18) $$ \begin{align} {\text{cov}}(Y_{g[i]}, Y_{g'[i']} \ | \ X, Z) = \begin{cases}{\mathbb{E}} (\epsilon^{*}_{g[i]} \epsilon^{*}_{g'[i']} \ | \ X, Z) &\quad \text{if}\; g=g' \\0 &\quad \text{if}\; g \neq g' \end{cases} \end{align} $$

What differs is only how ${\mathbb{E}}(\epsilon ^{*}_{g[i]} \epsilon ^{*}_{g[i']} \ | \ X, Z)$ is determined. In MLM, this covariance is parametrized as described above (e.g., ${\mathbb{E}} [\epsilon ^{*}_{g[i]} \epsilon ^{*}_{g[i']} \ | \ X, Z] = \omega ^2 + \mathbb {1}_{\{i=i'\} } \sigma ^2$ in a RI model), and estimated by plugging in the parameter estimates. By contrast CRSEs are remarkable for the lack of structure they impose on these within-group covariances. For example, after estimating $\hat {\beta }$ with an OLS of Y on X, CRSEs would simply construct empirical covariance estimates

(19) $$ \begin{align} \hat{{\mathbb{E}}}_{\text{CRSE}}(\epsilon^{*}_{g[i]} \epsilon^{*}_{g'[i']} \ | \ X, Z) = \begin{cases} c \times \hat{e}_{g[i]} \hat{e}_{g'[i']} &\quad \text{if}\; g=g' \\ 0 &\quad \text{if}\; g \neq g' \end{cases} \end{align} $$

where $\hat {e}_{g[i]} = Y_{g[i]} - X_{g[i]}^{\top } \hat {\beta }$ and c is a scalar finite sample correction. In other words, while MLMs make a strict assumption on the within-group covariances, CRSEs impose among the weakest possible assumptions by simply employing an empirical estimate based on fitted residuals. Appendix A.14 provides a more detailed discussion of the CRSE structure.

Fortunately, nothing prevents MLM users from employing CRSE assumptions during variance estimation (see also Cameron and Miller Reference Cameron and Miller2015), estimating variance according to

(20) $$ \begin{align} \widehat{{{\rm var}}}_{\text{CRSE}} (\hat{\beta}_{\text{MLM}}) = c \times (X^{\top} \hat{V}_{\text{MLM}}^{-1} X)^{-1} X^{\top} \hat{V}_{\text{MLM}}^{-1} {\scriptstyle \begin{bmatrix} \hat{e}_1 \hat{e}_1^{\top} & & 0 \\ & \ddots & \\ 0 & & \hat{e}_G \hat{e}_G^{\top} \end{bmatrix}} \ \hat{V}_{\text{MLM}}^{-1} X (X^{\top}\hat{V}_{\text{MLM}}^{-1} X)^{-1} \end{align} $$

where $\hat {e}_{g} = Y_g - X_g \hat {\beta }_{\text {MLM}}$ and V is defined in Equation (4). We discuss the choice of c in Appendix A.15.

This leads to a suprising and useful equivalance. Mirroring the equivalence between estimates of $\beta $ from bcMLM with $\Sigma = \sigma ^2 I_N$ and FE, the CRSEs for $\beta $ from both models are also equal if both use the same c (as we recommend in Appendix A.15). This fact, proven in Appendix A.16, may be surprising since the point estimates of $\hat {Y}_{g[i]}$ differ between models, and thus the overall error variance differs. This equivalence avoids debate over which method is appropriate to estimate the coefficient and standard error on the covariate of interest, since the answers will be the same.

Figure 5 illustrates the performance of CRSE with MLM in DGP 2, in which RI with $\Sigma =\sigma ^2 I_N$ misspecifies the dependence structure. Confidence intervals for $\beta _1$ using CRSEs fixes the undercoverage seen above in Figure 4 using the conventional standard errors. Coverage remains poor for $\beta _2$ with $G=15$ , with some undercoverage remaining at $G=50$ .

Figure 4 Coverage rates under RI, assuming $\Sigma =\sigma ^2 I_N$ in DGP 2. Note: Coverage rates for 95% nominal confidence intervals (vertical axis) for $\beta _1$ (left) and $\beta _2$ (right). Results across 1000 iterations, each drawn from DGP 2 with $\beta _0 = \beta _1 = \beta _2 = 1$ . The dashed-line represents the target coverage rate of 0.95.

Figure 5 Coverage rates of 95% confidence intervals from RI, OLS, and Group-FE, all with CRSE, in DGP 2. Note: Results across 1000 iterations, each drawn from DGP 2, with $\beta _0 = \beta _1 = \beta _2 = 1$ .