Fixed Effects Individual Slopes: Accounting and Testing for Heterogeneous Effects in Panel Data or Other Multilevel Models

In this section, we will derive the bias of the conventional FE estimator analytically. For simplicity, we use only a single covariate x and one individual slope variable w. Assume the true DGP is

where x_i and w_i are $T\times 1$, and ${\text{α}}_{i1}$ and ${\text{α}}_{i2}$ are individual-specific scalars. In full matrix notation, we can write the DGP as

where x and w are $NT\times 1$ vectors, ${\mathbf{\alpha }}_1$ and ${\mathbf{\alpha }}_2$ are stacked $NT\times 1$ vectors of subvectors ${\mathit{i}}_T{\text{α}}_{i1}$ and ${\mathit{i}}_T{\text{α}}_{i2}$, and $\circ$ denotes the Hadamard (element-wise) product. Now suppose, we would erroneously assume the DGP was

where the scalar $\text{γ}$ (note the missing subscript) denotes a constant effect of w_i on y_i over all units i. In full matrix notation, we can write this as

This represents a DGP of a model with homogeneous effects of w_i on y_i across all units (no individual slopes), which equals the DGP assumed in case of a conventional FE model with additional control variables as given by equation (6).

Why would the FE estimator be biased? As mentioned above, FE requires parallel slopes. Otherwise strict exogeneity does not hold. Define ${\text{α}}_{i2}={\overline{\text{α}}}_2+{\ddot{\text{α}}}_{2i}$, where ${\overline{\text{α}}}_2$ constitutes the overall mean of ${\text{α}}_{2i}$, for example, a common trend in the population, and ${\ddot{\text{α}}}_{2i}$ are individual-specific deviations from the mean. The parallel slopes condition would be $E({\text{α}}_{2i}|{\mathit{x}}_i,{\mathit{w}}_i)={\overline{\text{α}}}_2$. In the case with w_i being calendar time, this means that individual outcome trajectories may depend on unobservables ${\text{α}}_{2i}$, as long as the average slope in the population does not differ by values of covariates, notably x_i. With the definition of ${\text{α}}_{2i}$ above, we can state the parallel slopes condition also as $E({\ddot{\text{α}}}_{2i}|{\mathit{x}}_i,{\mathit{w}}_i)=0$, which implies $E({\mathit{x}}_i{\ddot{\text{α}}}_{2i})=0$. That is, individual deviations from the common trend may not be correlated with the causal variable.

Thus, estimating a DGP like equation (9) by FE as in equation (6), we receive

with ${\ddot{\mathbf{\xi }}}_i={\ddot{\mathit{w}}}_i{\ddot{\text{α}}}_{2i}+{\ddot{\mathbf{\epsilon }}}_i$. Since individual deviations from the common trend are contained in the idiosyncratic error term of the model and we estimate the demeaned equation (11) by OLS, a necessary condition for strict exogeneity, that is, for $E({\text{ξ}}_{it}|{\mathit{x}}_i,{\mathit{w}}_i,{\text{α}}_{i1})=0$, is $E\left({\ddot{\mathit{x}}}_{it}{\ddot{\text{ξ}}}_{it}\right)=E\left[\right({\mathit{x}}_{it}-{\overline{\mathit{x}}}_i\left)\right({\text{ξ}}_{it}-{\overline{\text{ξ}}}_i\left)\right]=0$. Hence, even if ${\text{∊}}_{it}$ is strictly exogeneous in the true model (7), strict exogeneity may be violated by FE estimation of a model that specifies a homogeneous effect for w_i. It must be violated if the parallel slopes assumption does not hold.

To show the bias, we apply the conventional FE estimator controlling for w to the DGP specified in equation (8). We use the Frisch–Waugh theorem (Frisch and Waugh 1933; Lovell 1963) and define $\mathbf{M}=\mathbf{I}-\ddot{\mathit{w}}{({\ddot{\mathit{w}}}^{⊺}\ddot{\mathit{w}})}^{-1}{\ddot{\mathit{w}}}^{⊺}$ to receive an estimate of the parameter of interest:

As ${\ddot{\mathit{i}}}_{NT}={\mathit{i}}_{NT}-{\overline{\mathit{i}}}_{NT}=0$, we can eliminate the time-constant part of the unobserved heterogeneity. Moreover, by assuming strict exogeneity $\text{E \hspace{0.05em} \hspace{0.05em}}({\ddot{\mathit{x}}}^{⊺}\ddot{\mathbf{\epsilon }})=0$ conditional on $\ddot{\mathit{x}}$, $\ddot{\mathit{w}}$, ${\mathbf{\alpha }}_1$, and ${\mathbf{\alpha }}_2$, equation (12) further simplifies to:

As above, we split up ${\mathbf{\alpha }}_2={\overline{\alpha }}_2+{\ddot{\mathbf{\alpha }}}_2$, where ${\overline{\text{α}}}_2$ constitutes the overall mean of ${\mathbf{\alpha }}_2$, and ${\ddot{\mathbf{\alpha }}}_2$ the individual-specific deviation from the mean. Then, we can rewrite equation (13) as

Now suppose, we specify a relationship between x and w of the form:

where $\mathbf{\delta }$ is a stacked $NT\times 1$ vector of the subvector ${\mathit{i}}_T{\text{δ}}_i$, $\overline{\delta }$ is the overall mean, $\ddot{\mathbf{\delta }}$ is the deviation from the mean, and v is an independent random vector. Note that ${\text{δ}}_i$ specifies how strongly the slope variable w is associated with our explanatory variable x. As for ${\text{α}}_{2i}$, we may view ${\text{δ}}_i$ as unobserved effect.

If we substitute for x in equation (14), we receive:

Recall that $\mathit{M}(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)$ equals the residual vector from regressing $(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)$ on $\ddot{\mathit{w}}$, and the expected value from the regression $(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)=\ddot{\mathit{w}}\text{λ}+\mathbf{\epsilon }$ is given by $\text{E \hspace{0.05em} \hspace{0.05em}}(\widehat{\text{λ}})={\overline{\ddot{\mathbf{\alpha }}}}_2=0$, as ${\ddot{\mathbf{\alpha }}}_2$ is the vector of the overall-demeaned ${\mathbf{\alpha }}_2$ values. It follows that the residual vector $\mathit{M}(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)=(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)-\ddot{\mathit{w}}\text{E \hspace{0.05em} \hspace{0.05em}}(\widehat{\text{λ}})=(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)$. The same argument applies to the first part of the denominator of equation (16), and we can thus simplify equation (16) to:

Equation (17) demonstrates that the conventional FE model produces biased point estimates if $\text{Cov \hspace{0.05em} \hspace{0.05em}}(\ddot{\mathbf{\delta }},{\ddot{\mathbf{\alpha }}}_2)\ne 0$. Intuitively, this means that FE is biased if the same individuals who experience a stronger influence of w_i on the dependent variable y_i also experience a stronger influence of w_i on the independent variable x_i. For instance, imagine that those individuals who have steeper wage earnings over the life course are also those who invest more in occupational training over the life cycle. Using a conventional FE model, we would conclude that additional training increases earnings, though that is due to the fact that more motivated people participate in more training and experience steeper earning trajectories—even without any further training (age is the slope variable). As noted earlier, we may view both $\mathbf{\delta }$ and ${\mathbf{\alpha }}_2$ as unobserved variables. Obviously, if some unobservable (e.g., motivation) drives both w_i and x_i, the covariance of $\mathbf{\delta }$ and ${\mathbf{\alpha }}_2$ will not be zero. Similar conclusions apply with dichotomous treatment variables: V. Ludwig and Brüderl (2018) demonstrate that those individuals who have a higher wage growth tend to marry at higher rates (work experience is the slope variable), which in the conventional FE model shows up in a spurious “marriage wage premium.” However, the slope variables are not restricted to timing or life-course variables but can also be some other characteristic in general multilevel scenarios. For instance, some families might have lower abilities to support their children’s cognitive development, thereby having a higher correlation between children’s starting abilities and later test score outcomes. At the same time, those families might be more selective in choosing children with lower starting abilities to enroll in preschool programs, which means a higher negative correlation between children’s starting abilities and the probability of preschool enrolment (e.g., Behrman and Rosenzweig 2004; Griliches 1979). FE models would then underestimate the effect of preschool enrolment on test score outcomes, as the models suffer from selection on starting abilities (starting ability is the slope variable).

As the previous section shows, several quite realistic situations exist in which the conventional FE estimator might yield biased estimates of the true effects. To deal with the problem, several authors (Brüderl and Ludwig 2015; Frees 2001; Lemieux 1998; Polachek and Kim 1994; Wooldridge 2010) have proposed a generalized FE model allowing for unit-specific slopes: the FEIS estimator.

The DGP for the FEIS estimator is given by:

or, respectively:

where y_i is $T\times 1$, X_i is $T\times K$, and ${\mathbf{\epsilon }}_i$ is $T\times 1$. W_i is a $T\times J$ matrix of slope variables, and ${\mathbf{\alpha }}_i$ a $J\times 1$ vector of individual-specific slope parameters, for J slope parameters including a constant term. If W_i includes a constant term, ${\mathbf{\alpha }}_i$ contains the ${\alpha }_{i1}$ known from the conventional FE estimator. Thus, the conventional FE estimator constitutes a specific case of the FEIS estimator where W_i contains only a constant term, or ${\mathit{W}}_i={\mathbf{i}}_T$, as this would reduce equation (19) to (2). Note that estimating equation (19) with W_i containing a constant plus additional slope variables would involve the estimation of N individual regressions of y_i on W_i or a least square dummy variable model with $NJ$ fixed effects. However, as with the conventional FE, the same can be obtained by running OLS on transformed data. We specify the “residual maker” matrix ${\mathit{M}}_i={\mathit{I}}_T-{\mathit{W}}_i{({\mathit{W}}_i^{⊺}{\mathit{W}}_i)}^{-1}{\mathit{W}}_i^{⊺}$ and apply OLS to the transformed data:

where ${\tilde{\mathit{y}}}_i$, ${\tilde{\mathit{X}}}_i$, and ${\tilde{\mathbf{\epsilon }}}_i$ are the residuals of regressing y_i, each column-vector of X_i, and ${\mathbf{\epsilon }}_i$ on W_i. In full matrix notation, the model can be written as

where we define $\tilde{\mathit{M}}={\mathit{I}}_N\otimes {\mathit{M}}_i$, where $\otimes$ is the Kronecker (block-wise) product and ${\mathit{M}}_i={\mathit{I}}_T-{\mathit{W}}_i{\left({\mathit{W}}_i{ }^{⊺}{\mathit{W}}_i\right)}^{-1}{\mathit{W}}_i{ }^{⊺}$. $\tilde{\mathit{M}}$ constitutes a stacked block-diagonal matrix with block-diagonal elements equal to M_i and off-block-diagonal elements equal to zero. The resulting estimator is given by:

To see why premultiplying by $\tilde{\mathit{M}}$ in equation (23) eliminates the bias in equation (16), we will now only consider the last term of equation (16), which is decisive for the bias: $\mathit{M}(\ddot{\mathit{w}}\circ {\ddot{\mathbf{\alpha }}}_2)$. From estimating equation (23) instead of (12), we would receive $\tilde{\mathit{M}}({\mathbf{\alpha }}_1+\mathit{w}\circ {\mathbf{\alpha }}_2)$. Now recall that W always contains a constant intercept plus the slope variable(s). So for each group i, $\tilde{\mathit{M}}({\mathbf{\alpha }}_1+\mathit{w}\circ {\mathbf{\alpha }}_2)$ contains the residual vector of the following regression:

where we get $\text{E \hspace{0.05em} \hspace{0.05em}}({\widehat{\text{λ}}}_{i1})={\text{α}}_{i1}$ and $\text{E \hspace{0.05em} \hspace{0.05em}}({\widehat{\text{λ}}}_{i2})={\text{α}}_{i2}$, which is the individual-specific mean and the individual-specific slope parameter. Thus, the stacked $NT\times 1$ vector ${\widehat{\text{λ}}}_1$ of subvectors ${\mathit{i}}_T{\widehat{\text{λ}}}_{i1}$ equals the stacked vector ${\mathbf{\alpha }}_1$, and the stacked $NT\times 1$ vector ${\widehat{\mathbf{\lambda }}}_2$ equals the stacked vector ${\mathbf{\alpha }}_2$. Consequently, $\tilde{M}({\mathbf{\alpha }}_1+\mathit{w}\circ {\mathbf{\alpha }}_2)=({\mathbf{\alpha }}_1+\mathit{w}\circ {\ddot{\mathbf{\alpha }}}_2)-({\widehat{\mathbf{\lambda }}}_1+\mathit{w}\circ {\widehat{\mathbf{\lambda }}}_2)=0$, which eliminates the bias in equation (16) and provides a consistent estimator of $\text{β}$.

As shown previously, the conventional FE estimator requires the strict exogeneity assumption $\text{E \hspace{0.05em} \hspace{0.05em}}({\epsilon }_{it}|{\mathit{x}}_i,{\mathit{w}}_i,{\text{α}}_{i1})=0$, that is, conditional on w_i and the individual constant ${\alpha }_{i1}$, to hold for consistency. Note that ${\text{α}}_{i1}$ is a single scalar parameter for each group. In contrast, the FEIS estimator requires the strict exogeneity assumption of the form $\text{E \hspace{0.05em} \hspace{0.05em}}({\epsilon }_{it}|{\mathit{x}}_i,{\mathit{w}}_i,{\mathbf{\alpha }}_i)=0$, where ${\mathbf{\alpha }}_i={\text{α}}_{i1}\mathrm{,...},{\text{α}}_{iJ}$, that is, conditional on w_i and all J individual slopes of the vector ${\mathbf{\alpha }}_i$ (including the constant ${\text{α}}_{i1}$), to hold. Thus, consistency of FEIS relies on weaker assumptions than consistency of conventional FE models. Still, for consistent estimation of ${\widehat{\mathbf{\beta }}}_{\text{FEIS}}$, we need to impose the rank condition $\text{rank \hspace{0.05em} \hspace{0.05em} E }\left({\tilde{\mathit{X}}}_i^{⊺}{\tilde{\mathit{X}}}_i\right)=K$, where ${\tilde{\mathit{X}}}_i={\mathit{M}}_i{\mathit{X}}_i$ and $\text{rank \hspace{0.05em} \hspace{0.05em}}(\mathit{M})=T-J$ for K covariates, T time periods, and J slope parameters (including the constant). Thus, we need $T>J$, that is, the number of time periods to be larger than the number of specified slope parameters. While the conventional FE requires $T\ge 2$, the FEIS with one slope parameter (additional to the constant) requires $T\ge 3$, which in unbalanced data may lead to a problem of sample selection (for a discussion, see Limitations section). However, with large N, we should still be able to obtain a consistent estimate of ${\widehat{\mathbf{\beta }}}_{\text{FEIS}}$ if some observations have $\text{rank \hspace{0.05em} \hspace{0.05em}}\left({\tilde{\mathit{X}}}_i^{⊺}{\tilde{\mathit{X}}}_i\right)<K$ because of lacking within-variance on single covariates (Wooldridge 2010:379).

As has been shown by Mundlak (1978) and Chamberlain (1982), we can also derive the parameter estimates of the conventional FE model by estimating a correlated random effects (CRE) model. The CRE is given by:

where ${\overline{\mathit{X}}}_i$ and ${\overline{\mathit{W}}}_i$ contain the individual-specific means of the variables X_i and W_i. If we estimate the CRE specification by generalized least squares (GLS), the parameter vectors $\mathbf{\beta }$ and $\mathbf{\delta }$ represent estimates of the within effects—identical to the FE coefficients—and the vectors $\mathbf{\gamma }$ and $\mathbf{\theta }$ give us an estimate of the between effect minus the within effect (e.g., Allison 2009; Arellano 1993; Wooldridge 2010). Consistency of the RE estimator is based on the assumption that $\mathbf{\gamma }=\mathbf{\theta }=0$. Thus, we can perform a Wald test on the hypothesis ${\text{H}}_0:\mathbf{\gamma }=\mathbf{\theta }=0$ to obtain a regression-based version of the Hausman test statistic.

An advantage of the ART is that we can easily extend it to test for consistency of the FE estimator or more precisely for comparing FE against FEIS estimates. Therefore, we need to compute the individual-specific predicted values ${\widehat{\mathit{X}}}_i={\mathit{W}}_i{({\mathit{W}}_i^{⊺}{\mathit{W}}_i)}^{-1}{\mathit{W}}_i^{⊺}{\mathit{X}}_i$ based on the slope variables in W_i and extend the CRE in equation (27) to:

Analogue to the discussion above, the parameter vector $\mathbf{\beta }$ of the extended CRE is identical to the parameters estimated by FEIS. Consistency of the FE estimator is based on the assumption that $\mathbf{\rho }=0$. Thus, we can perform a Wald test of ${\text{H}}_0:\mathbf{\rho }=0$ (the FE estimator is consistent). If the Null is not rejected, $\mathbf{\rho }=0$ and we can omit the term ${\widehat{\mathit{X}}}_i\mathbf{\rho }$, making the conventional FE estimator more efficient than the FEIS estimator. Note that the terms ${\overline{\mathit{X}}}_i\mathbf{\gamma }$, ${\overline{\mathit{W}}}_i\mathbf{\theta }$, and ${\mathit{W}}_i\mathbf{\delta }$ are actually not necessary in equation (28) to obtain the FEIS coefficients in $\widehat{\mathbf{\beta }}$ ($\widehat{\mathbf{\beta }}={\widehat{\mathbf{\beta }}}_{\text{FEIS}}$ if we include ${\widehat{\mathit{X}}}_i$ as additional covariates). However, they are necessary if we want to compare FEIS against FE, as equation (28) only reduces to the FE estimator conditional on $\mathbf{\rho }=0$ if we include the individual-specific means as in equation (27).

The test discussed in the previous section comes with the advantage that we can specify standard errors which are robust to arbitrary forms of heteroscedasticity and serial correlation (e.g., Arellano 1987, 1993). This is not possible with the standard Hausman test, but it is necessary if the errors are correlated within clusters (units) and thus observations are not independent of each other. However, computation of cluster-robust standard errors is based on the assumption that the number of clusters approaches infinity (Stock and Watson 2008). Cameron, Gelbach, and Miller (2008) show by simulation that relying on cluster-robust standard errors can lead to downwardly biased standard errors and an overrejection of ${\text{H}}_0$ when the number of clusters is small (N cluster $<30$).

To overcome this problem, Cameron et al. (2008) and Cameron and Miller (2015) proposed a BSHT. The idea is to perform pairwise-clustered bootstrapping by randomly selecting $N^{\star }$ clusters with replacement from the original sample of N clusters and to repeat this over R replications. Estimating FEIS and FE models in each replication provides R estimates of ${\mathbf{\beta }}_1$ and ${\mathbf{\beta }}_0$, which can be used to estimate ${\widehat{\mathit{V}}}_{{\widehat{\mathbf{\beta }}}_1-{\widehat{\mathbf{\beta }}}_0}$. While the ART takes the variance–covariance of the two estimators into account by assuming an identical error variance for FEIS and FE, the BSHT relies on an estimate of the variance–covariance. The bootstrapped estimate of the variance–covariance matrix is given by

where ${\widehat{\mathbf{\beta }}}_{b\mathrm{,1}}^{\star }$ and ${\widehat{\mathbf{\beta }}}_{b\mathrm{,0}}^{\star }$ are the estimated coefficients for the consistent (${\widehat{\mathbf{\beta }}}_1$) and the efficient (${\widehat{\mathbf{\beta }}}_0$) model returned by single bootstrap replications, and ${\overline{\mathbf{\beta }}}_{b\mathrm{,1}}^{\star }$ and ${\overline{\mathbf{\beta }}}_{b\mathrm{,0}}^{\star }$ denote the mean values of the coefficients over all R replications. The resulting variance–covariance matrix can then be used to compute the Hausman test statistic of equation (26). This is equal to the method called “pairs cluster bootstrap-se” in the simulations by Cameron et al. (2008). Both methods presented in this section can be used to test for misspecification in conventional FE models.

In the simulations, we set up balanced panel data sets with $N=300$ units and $T=10$ time points, leading to a total of 3,000 observations. The data are generated by the following DGP:

where $\mathbf{\epsilon }$ and v are independent Gaussian random vectors $\mathbf{\epsilon },\mathbf{\nu }\sim N(0,1)$. The scalar parameter $\text{θ}$ specifies the correlation between x and the time-constant “unobserved” heterogeneity ${\mathbf{\alpha }}_1$, and the $NT\times 1$ parameter vector $\mathbf{\delta }$ the correlation between x and the slope variable w.

As described earlier, ${\mathbf{\alpha }}_1$, ${\mathbf{\alpha }}_2$, and $\mathbf{\delta }$ are stacked vectors of the individual-specific subvectors ${\mathit{i}}_T{\text{α}}_{\mathit{i}1}$, ${\mathit{i}}_T{\text{α}}_{\mathit{i}2}$, and ${\mathit{i}}_T{\delta }_i$. Further, we define ${\mathbf{\alpha }}_1$ as a normally distributed variable with mean ${\text{μ}}_{{\alpha }_{\text{1}}}$ and standard deviation ${\text{σ}}_{{\text{α}}_1}$. The values of random variable w are drawn from a normal distribution within each unit i, with mean ${\mathrm{\text{μ}}}_{\mathit{w}}$ and standard deviation ${\text{σ}}_w$. Finally, we draw ${\mathbf{\alpha }}_2$ and $\mathbf{\delta }$ from a bivariate normal distribution,

with N draws, which are then replicated T times for each unit i to get time-constant vectors ${\mathbf{\alpha }}_2$ and $\mathbf{\delta }$.¹

As shown earlier (see equation [17]), the theoretical bias of the standard FE estimator has four components: the covariance of $\ddot{\mathbf{\delta }}$ and $\ddot{\mathbf{\alpha }}$, and the variances of $\ddot{\mathbf{\delta }}$, $\ddot{\mathit{w}}$ and $\ddot{\mathbf{\nu }}$. In the simulations, we experiment with these components in that we vary parameters $\varphi \text{:= Cov \hspace{0.05em} \hspace{0.05em}}(\mathbf{\delta },{\mathbf{\alpha }}_2)$, ${\sigma }_{\text{δ}}^2:= \text{Var} (\mathbf{\text{δ}})$, ${\text{σ}}_w^2:= \text{Var} (\mathbf{w})$, and ${\text{σ}}_{\mathbf{\nu }}^2:= \text{Var} (\mathbf{\nu })$ systematically. For all our simulations, we set $\text{β}=1$. In the following, we present two sets of simulation results.

In setting (1), we vary $\varphi$ between −0.8 and 0.8 and set ${\sigma }_{\mathbf{\delta }}^2$, ${\sigma }_{\mathit{w}}^2$, and ${\sigma }_{\mathbf{\nu }}^2$ to fixed values. We are mainly interested in the performance of the FE and FEIS estimators, but we also include the conventional RE estimator because it is the main alternative to the FE model in applied research using panel data. The RE estimator is known to be biased if ${\mathbf{\alpha }}_1$ is related to x. In the initial setup, we therefore vary also $\text{θ}\in \{0,1\}$.

For each value of $\varphi$, we estimate FEIS, FE, and RE regression models of y on x and w with cluster-robust estimation of the covariance matrix. Similarly, the ART is performed using the robust error variance proposed by Wallace and Hussain (1969). We compute the mean bias of $\widehat{\beta }$ for each of the three estimators using $R=1,000$ replications. Since true $\beta$ remains fixed at 1 throughout, the mean of the absolute bias is equal to the mean relative bias. So, we will interpret the average bias as a percentage. We use $R_b=100$ replications to estimate the bootstrapped covariance matrix for the BSHT. To evaluate the size and power of the ART and BSHT in detecting biased estimates, we report rejection rates for the 95 percent confidence level, that is, the proportion of test statistics that yield a p value $<.05$.

In setting (2), we focus on the bias of the FE estimator and the performance of the ART and BSHT with respect to the bias. We set $\theta =0.5$ and expand the initial setting in that we vary all four parameters that go into the bias of FE. Specifically, we set

• $\varphi \in \left\{\mathrm{0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.65,0.8}\right\}$,

• ${\sigma }_{\delta }^2\in \{0.2,0.5,1,3,5\}$,

• ${\sigma }_w^2\in \{0.2,1,5\}$,

• ${\sigma }_v^2\in \{0.2,1,5\}$.

For each of the 360 resulting parameter combinations, we report the mean bias of the FE estimate and the rejection rates of the ART and BSHT comparing FE and FEIS coefficients over $R=1,000$ replications. This allows us to asses (1) the severity of the bias in FE estimates and (2) the sensitivity of our proposed test statistics across different situations.

In this example (V. Ludwig and Brüderl 2018), we use the data set of the original study, which consists of individual panel data from the National Longitudinal Survey of Youth 1979 (Bureau of Labor Statistics 2014).³ We investigate the “marital wage premium”: We analyze whether marriage leads to an increase in the hourly wage for men. By estimating FEIS models, we rule out the alternative that those men who eventually get married also show steeper wage growth over their career (even before marriage). Therefore, we specify actual work experience (in years) and squared work experience as slope variables.

Table 1 shows the regression results of the hourly wage rate on marital status (plus additional controls) for different model specifications: RE model, conventional FE model, a random intercept random slopes multilevel model (RIRS), and the FEIS model. Note that we include a random slopes model here, because the RIRS allows for individual-specific trends as well. As in the FEIS, random slopes are specified for work experience and work experience squared. In the RIRS model, however, the individual variation of the intercepts and slopes is assumed exogeneous. This assumption is not needed with the FEIS model.

Table 1. Example 1: Regression results for log hourly wage rates.

Table 1. Example 1: Regression results for log hourly wage rates.

View larger version

As can be seen in Table 1, the effect of marriage on the wage rate varies notably across the models. The RE and the FE return a significant effect of marriage on the log wage of 0.100 and 0.080, respectively. Similarly, the RIRS returns a significant effect of marriage, but at a slightly lower magnitude of 0.064. In comparison, the FEIS model reports an effect of only 0.005 which is far from significant. Note that the FE estimate exceeds the FEIS estimate by a factor of 16, and also the RIRS is nearly 14 times higher than the FEIS estimate. Consequently, according to the FE (as well as RE and RIRS), we would conclude that there really is a marital wage premium for men. However, the FEIS revises this conclusion by showing that the effect mostly stems from the fact that men with a steeper wage growth tend to marry earlier (those men with a stronger effect of experience on wage also exhibit a stronger effect of experience on marriage). In consequence, our replication underlines the conclusion of the original study by V. Ludwig and Brüderl (2018): The marital wage premium is mainly due to selection on wage growth rather than being a causal effect of marriage on men’s productivity.

This example illustrates how heterogeneous growth curves related to the covariates can drastically bias the conclusions drawn from conventional FE and RE models that are typically applied to panel data. Therefore, we propose to test for the presence of heterogeneous slopes in panel models by default.

Table 2 shows the results of the proposed ART and BSHT specification tests. Overall, both tests strongly reject the Null that the conventional FE and RE models provide consistent estimates. Regarding the first test comparing FEIS against FE, a highly significant ${\chi }^2$ value of 436.966 (581.808 in the bootstrapped version) indicates that estimates of the conventional FE are inconsistent because of heterogeneous slopes. Thus, we should use FEIS instead of FE. The second test of FE against RE is equivalent to the conventional Hausman test under the assumption of equal error variance for both models, but it does allow for clustered standard errors. The ART returns a highly significant ${\chi }^2$ and thus favors FE over an RE model.⁴ The third test statistic offers a direct comparison of FEIS against RE: Both the ART and the BSHT show a highly significant ${\chi }^2$, thereby indicating that we should reject the Null hypothesis of consistent estimates in the RE model. Taken together, when ignoring the possibility of heterogeneous slopes, we would lean toward using a conventional FE model, thereby relying on an effect of marriage which equals 16 times the effect obtained in a model accounting for individual-specific slopes.

Table 2. Example 1: Specification Tests.

Table 2. Example 1: Specification Tests.

View larger version

While the first example constitutes a textbook case for the application of FEIS with heterogenous growth curves, we will now turn to an example in which we apply the FEIS to a more general hierarchical data structure. The example concerns the impact of Head Start, a large-scale U.S. preschool program for disadvantaged families, on children’s cognitive performance. In general, there has been a large debate on whether Head Start participation induces cognitive benefits, and especially whether it generates long-term benefits on schooling and other outcomes (e.g., J. Ludwig and Miller 2007). We know from experimental evidence of other preschool programs that preschool interventions improve educational outcomes of children from disadvantaged families in the short and in the long run (Heckman 2006; Schweinhart et al. 2005). In addition, using experimental data of a randomly drawn subsample of Head Start participants, J. Ludwig and Phillips (2008) report short-term treatment effects of 0.15–0.35 standard deviations for a range of cognitive outcomes measured at age 3 or 4. However, findings with nonexperimental data for children participating in Head Start differ markedly in that effects seem to be much smaller, and effects fade-out over time (Duncan and Magnuson 2013). As we usually lack continuous information on children’s cognitive outcomes, most observational studies apply a within-family design (siblings estimators) to account for children’s individual preprogram conditions.

Below, we replicate and extend the study by Deming (2009). The author used data from the National Longitudinal Mother–Child Supplement, which surveyed the mothers of the NLSY79 (Bureau of Labor Statistics 2014) biannually from 1986 through 2004. When relying on simple OLS models, he does not find substantial short-term effects on cognitive outcomes. Deming (2009) suspects the estimates are biased downward. In the authors’ own words,

Essentially, using family fixed effects assumes that selection into Head Start among children from the same family is not related to unobserved characteristics of the children, which simultaneously correlate with school performance. This certainly is a strong assumption. Therefore, the author included an extensive set of 11 covariates to control for different pretreatment conditions of children within the same family. He also assembled an index of children’s starting conditions from pretreatment covariate values in which high values indicate overall favorable starting conditions. However, controlling for child-specific pretreatment conditions did not alter the results substantively. In the final (preferred) model, Deming (2009) finds that Head Start increases the test scores of children at ages 5–6 by 0.145 standard deviations, which meets the lower bound of results reported in J. Ludwig and Phillips (2008). Moreover, the effects further declines for test score results at older age.

However, the low (long term) returns may be the result of a problem of heterogenous slopes on the family level: If disadvantaged mothers are less able or willing to balance disparate starting conditions, their children’s starting conditions would be strongly (positively) linked to later outcomes. At the same time, children’s starting conditions may be more strongly (negatively) linked to participation in Head Start in these families (e.g., Behrman and Rosenzweig 2004; Griliches 1979). Perhaps disadvantaged mothers anticipate they will not be able to fully support their offspring given their restricted family resources. In this case, the treatment effect would be underestimated in a conventional family FE model which neglects family-specific heterogeneity in the effect of starting conditions.

Indeed, our reanalysis of Deming’s data set⁵ supports this argument. Model 1 in Table 3 perfectly replicates Deming’s main result. Even after controlling for family FE and an extensive set of covariates to capture children’s starting conditions, the effect of participation in Head Start is rather weak, ranging from 0.05 to 0.14 standard deviations of test performance. Model 2 estimates the same family FE model using a slightly smaller sample (families with fewer than three observations excluded) and a different specification of control variables, replacing the covariates for children’s starting conditions by Deming’s pretreatment index. The estimates of the treatment effect remain very similar in this model (albeit somewhat smaller compared to Model 1). Model 3 extends the specification by allowing for family-specific slopes on the pretreatment index. This more flexible FEIS model yields substantially larger treatment effects, in the range of 0.24–0.35 standard deviations. Unlike the FE results, this finding is in line with experimental evidence (J. Ludwig and Phillips 2008). Apparently, the effect of children’s pretreatment conditions on later test performance differs between families, and the heterogeneity in this effect is also linked to program enrolment.

Table 3. Example 2: Regression results for cognitive test scores.

Table 3. Example 2: Regression results for cognitive test scores.

View larger version

Despite the large difference in effects between the FE and FEIS estimates (models 2 and 3), a specification test indicates that the family-FE model is not biased. The ART returns a p value of $.219$, which seems to indicate that there is no need to specify family-specific slopes for the pretreatment index. However, the nonsignificant test statistic follow from construction of the estimation sample. As Deming (2009) notes, in his study, “the impact of Head Start is identified by comparing siblings in the same family who vary in their participation in preschool programs. Thus, if all three of a mother’s children were enrolled in Head Start, I cannot say anything about the effectiveness of the program for them” (p. 115).

Accordingly, we further restrict the sample to mothers with at least two children who differ in Head Start participation and reestimate our FE and FEIS models. Models 4 and 6 still provide similar, but somewhat stronger estimates for the Head Start treatment effect (compared to models 2 and 3). As in the first example, the estimates of the random slope multilevel model (Model 5) lie between FE and FEIS, but it is still far away from treatment effects reported in FEIS. Interestingly, the effect of visiting other Preschools also turns positive in the FEIS model. Based on model 6, the ART returns a p value of $.019$, which indicates that we should in fact specify heterogeneous slopes for the pretreatment index (Table 4).

Table 4. Example 2: Specification Tests.

Table 4. Example 2: Specification Tests.

View larger version

Why does the test work with the smaller sample? The reason is that, unless we apply the restriction to families with variation on the treatment variable, the CRE model estimated for the test suffers from collinearity. If children from the same family do not differ with respect to program participation, the family-specific predicted values for Head Start covariates are equal to the family-specific mean values. Since this holds for 61 percent of families contained in the full original sample, predicted and mean values are highly collinear in the CRE model, resulting in large standard errors and an insignificant test statistic for the ART.⁶ In terms of our simulation results, the power of the ART is insufficient when using the full sample because we run into a situation shown in Figure 3G with low variation of both the treatment variable (${\sigma }_w^2>>{\sigma }_v^2$) and the effect of the slope variable (${\sigma }_{\delta }^2$ is small).