Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
66,206 characters · 19 sections · 71 citation commands
Double Machine Learning for Static Panel Models with Fixed Effects
Recent advances in the econometric literature on Machine Learning (ML) use the power of ML algorithms, widely used in data science for solving prediction problems, to enhance existing estimation procedures for treatment and other kinds of causal effect. Notable developments include novel ML algorithms for causal analysis such as Causal Trees by atheyimbens2016, Causal Forests by atheywager2018 and Generalised Random Forests by atheytibshirani2019. However, the key development, as far as this paper is concerned, is Double/Debiased Machine Learning (DML) by chernozhukov2018 wherein ML is used to learn nuisance functions with ex ante unknown functional forms, and the predicted values of these functions used to construct (orthogonalized) scores for the interest parameters from which consistent and asymptotically normal estimators can be obtained. DML is a very general estimation framework but there are limited examples of its application to panel data, notable examples of which include chang2020, klosin2022, and semenova2023.
In this paper, we develop and assess novel DML procedures for estimating treatment (or causal) effects from panel data with fixed effects. The procedures we propose are extensions of the correlated random effects (CRE), within-group (WG) and first-difference (FD) estimators commonly used for linear models to scenarios where the underlying model is non-linear. Specifically, these are based on an extension of the partially linear regression (PLR) model proposed by robinson1988 to panel data through the inclusion of time-varying predictors and unobserved individual heterogeneity (i.e.\ individual fixed effects).
Our methodological contribution is twofold and complementary to the recent work on causal panel data estimation using DML by chang2020 on difference-in-differences, klosin2022 and semenova2023 on high-dimensional treatment heterogeneity in model with fixed effects.
A first contribution is that the procedures we develop are not based on {\em ex ante} taking the nuisance functions or fixed effects to be accurately approximated by high-dimensional sparse functions. Thus, we do not focus solely on the Least Absolute Shrinkage and Selection Operator (LASSO) in the first stage of DML. This is in contrast to klosin2022 and semenova2023 who explicitly exploit sparsity to propose alternative procedures. While we do not rule out using LASSO (and indeed do so in both our simulation study and empirical application), our approach requires only minimal assumptions (like Lipschitz continuity) about the nuisance functions. This is in line with chang2020 to encourage researchers to adopt estimation strategies which involve the use of different ML algorithms. Such strategies could involve the use of {\em ensemble} learning such as the stacking of breiman1996, the super-learning (a weighted average the best predictions over different choices of ML algorithm) of laan2007, or selecting the best-performing learner on an application-by-application basis to ensure inference is based on the most accurate predictions. This is important because, in practice, many ML algorithms are difficult to tune and have been shown to perform very differently for different datasets in different contexts.
Second, on a computational note, we use block-k-fold cross-fitting, where the entire time series of the sampled unit is allocated to one fold to allow for possible serial correlation within each unit as is common with panel data. In this, we are aligned with the concurrent work by klosin2022 in not needing to rely on the weak dependence assumption to deal with dependent data as do semenova2023. Conversely, chang2020 does not require any additional assumption or change in the splitting algorithms because they use only one observation per unit after treatment, ruling out the presence of any fixed effects.
Our method for panel data models with individual fixed effects is general and particularly relevant for applied researchers. We provide new estimation tools within the existing DML framework for use on panel data. In doing so, we broaden the reach of DML to a large family of empirical problems for which the time dimension must be properly accounted for. Our focus, in the subsequent development, on the homogeneous treatment effect case is also because such models are widely used by applied researchers, but we also show how our procedures can be extended to heterogeneous treatment effects provided that the analyst is prepared to specify a (finite-dimensional) parametric model for this heterogeneity. More widely, we encourage researchers to use the procedures we propose in place of existing ones, or to test the robustness of their results - based on, say, linear models - to non-linearity.
We carry out a simulation study to assess our DML procedures and find large accuracy gains from using DML with flexible learners when the data generating process is highly non-linear (specifically, one involving a non-linear discontinuous function of the regressors); and, while ordinary least squares (OLS) {\em can} outperform ML when the data generating process is linear (as expected) and for some non-linear processes (specifically, for a smooth process that excludes interactions), it is not robust because the analyst never knows the form of the data generating process ex ante. Finally, we apply our new procedures to observational panel data by re-analyzing part of the study by fazio2023 on the effect of the introduction National Minimum Wage (NMW) in the United Kingdom (UK) on voting for conservative parties. We find evidence that the first-difference procedures are the most robust and perform well, and confirm the conclusion from our simulation study that ensemble learning strategies are crucial to obtaining robust results.
The remainder of the paper is structured as follows. Section (ref) provides an overview of the literature and places our novel contribution within it. Section (ref) motivates the partially linear panel regression model and the causal assumptions which must hold to ensure the target parameter can be interpreted as a causal effect. Section (ref) introduces the two approaches we take to handle the fixed-effects problem. Section (ref) formally introduces the DML estimation procedures. Section (ref) briefly discusses the Monte Carlo simulation results. Finally, Section (ref) illustrates an empirical application of the procedure and we make concluding remarks in Section (ref).
There is a growing body of econometrics literature on using ML for causal inference. One strand focuses on building or modifying existing learners to consistently estimate and make inferences about causal effects atheyimbens2016,atheywager2018,atheytibshirani2019,lechner2022,difrancesco2022. Another strand focuses on incorporating ML into traditional statistical estimators (e.g., least squares, generalised method of moments, maximum likelihood) to estimate causal effects more accurately belloni2014,belloni2016,chernozhukov2018,chang2020,chernozhukov2022,huber2023. This paper falls into the second strand.
Much of the ML literature in econometrics, including that for causal estimation, is built around penalized regression methods like LASSO, where the onus is on the analyst to specify a sufficiently rich data dictionary for the problem at hand. For example, belloni2016 proposed two-step post-cluster LASSO procedures for panel data with additive individual-specific heterogeneity in which potential control variables are selected using LASSO at stage one followed by the estimation of the homogeneous treatment effect at stage two.
Two recent papers propose estimation procedures based on LASSO for non-linear panel models. Both focus on different causal targets to ours:\ in the context of continuous treatments, klosin2022 proposed an estimator for the average partial effect based on a PLR model suitable for static panels, where a first-difference data transformation is used to handle omitted fixed time-invariant confounding; and semenova2023 propose a procedure for CATE estimation based on a PLR model for dynamic panels. Both of these procedures rely crucially on the key nuisance functions being well approximated by high-dimensional sparse functions:\ without assuming this {\em ex-ante}, their procedures do not follow. We propose alternative procedures that make no such {\em ex-ante} assumptions, and set these up within the DML framework developed by chernozhukov2018 who considered the PLR model at length but did not provide any panel examples.
The generality of our approach allows us to move beyond LASSO to other types of ML algorithm. LASSO is powerful but it requires the analyst to pre-specify a sufficiently rich data dictionary, and the memory required to store these in panels with many waves could be prohibitively large such that alternative learners are preferable. However, by staying general, we rule out solutions to the fixed-effects problem like that proposed by semenova2023 who incorporate the fixed effects into the LASSO data dictionary under a weak sparsity assumption. While recent work by kolesar2023 argues that, in the absence of prior knowledge to the contrary, this assumption is unlikely to hold, semenova2023 avoid this issue by decomposing the overall fixed effect into a fixed part, involving the individual-specific means of the (fixed) predictors, and a fixed-effect residual which {\em can} be incorporated into the dictionary under a weak sparsity assumption.
By not formulating the learning problem in terms of LASSO, we must consider alternative approaches so, instead, propose three estimation procedures that obviate entirely the requirement to model the fixed effects. The first two are based on transforming the outcome to induce a reduced-form model that does not depend on the fixed effects. Because we are agnostic, our first-difference procedure differs from that proposed by klosin2022 in ways we explain further on. We also consider a procedure based on a correlated random effects model similar to that used by wooldridge2020 for limited dependent variable panel models to convert the fixed effects into random effects.
Finally, we note that our work also fits into the literature that leverages the power of ML for causal analysis and policy evaluation. The value added of doubly-robust procedures has been explored in applied works by, for example, knaus2022double, bach2023, strittmatter2023, baiardi2024plough,baiardi2024. Because panel data are widely used in applied analyses, our proposed procedures for panel data models have the potential to attract the interest of applied researchers from various fields broadening the applicability of DML.
\setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{lemma}{0} \setcounter{proposition}{0}
Suppose the panel study collected information on each of $N$ individuals at each of the $t$ time periods, or waves. Let $\{(Y_{it},D_{it},\vector X_{it}) : \ t=1,\ldots,T\}_{i=1}^{N}$ be $N$ independent and identically distributed (iid) random vectors for individual $i$ across all $T$ waves, where $Y_{it} \in \mathcal{Y}$ is the outcome variable, $D_{it}\in{\cal D}$ a continuous or binary treatment variable (or intervention), and $\vector X_{it} = (X_{it,1},\dots,X_{it,p})'\in \mathcal{X}$ a $p\times 1$ vector of control (pre-determined) variables, usually including a constant term, able to capture time-varying confounding induced by non-random treatment selection.\footnote{Throughout, we use ${\bf v}^\prime$ to indicate the matrix transpose of arbitrary vector {\bf v} and take vectors to be column vectors unless explicitly stated to the contrary.} We denote the realizations of these random variables by $\{(y_{it}, d_{it}, \vector x_{it})\}$, respectively. For continuous $D_{it} \in \mathcal{D}\subset\mathbb{R}$, if $d_{it}\geq 0$ then a dose-response relationship is presumed to hold with $d_{it}=0$ indicating null treatment; otherwise, $D_{it}$ is taken to be centered around its mean $\mu_D$ such that $D_{it} \equiv D_{it}-\mu_D$. For binary $D_{it}\in\{0,1\}$, $d_{it}=0$ is taken to indicate the absence and $d_{it}=1$ the presence of treatment.
The non-linear model we propose is a simple extension of the partially linear model of robinson1988 to panel data. We then use a potential outcomes causal framework (e.g.\ rubin1974) to set out the assumptions under which its parameters can be interpreted causally. As pointed out by lechner2015, this allows us to interpret the resulting partially linear panel model as reduced-form without relying on specific parametric model for the data generating process. To this end, further define the set $Y_{it}(.)=\{Y_{it}(d):d\in {\cal D}\}$ containing all potential outcomes for individual $i$ at wave $t$, where $Y_{it}(d)$ is the realization of the outcome for individual $i$ at wave $t$ were the treatment level set to $d$, with one potential outcome for every possible value the treatment could take. The realizations of the wave $t$ potential outcomes are taken to occur before treatment selection at wave $t$, and are linked to the observed outcome by the {\em consistency assumption} that $Y_{it}(d_{it})=Y_{it}$ with the others latent {\em counterfactuals}.\footnote{The stable unit treatment value (SUTVA) assumption, that $Y_{it}(d)$ does not depend on the treatment assignments of any other individual, is implicitly taken to hold.} In the interval preceding wave $t$, it is also presumed that the realization of time-varying predictor $\vector X_{it}$ precedes that of $(Y_{it},D_{it})$.
Finally, define the sets $\xi_i$ of time-invariant heterogeneity terms influencing $(Y_{it},D_{it})$, and $L_{t-1}(W_i) = \{W_{i1},\dots,W_{it-1}\}$ the lags of a generic random variable $W_{it}$ available at wave $t$ such that $L_0(W_i)\equiv\varnothing$.
We extend the partially linear model proposed by robinson1988 to panel data by introducing fixed effect $\alpha^*_i$ to give the partially linear panel regression (PLPR) model
where $g_1$ is a non-linear nuisance function of $\vector X_{it}$ and $E[U_{it}|D_{it},\vector X_{it},\alpha^*_i]=0$ but $E[\alpha_i^*\mid D_{it},\vector X_{it}]\neq 0$. The target parameter $\theta_0$ is the average partial effect of continuous $D_{it}$ such that $d\theta_0=E[Y_{it}(d)-Y_{it}(0)]$, or the average treatment effect (ATE) $E[Y_{it}(1)-Y_{it}(0)]$ for binary treatments.
Following chernozhukov2018, we focus on the partialled-out PLPR (PO-PLPR) model
where $l_1$ and $m_1$ are nuisance functions, $\alpha_i$ is a fixed effect, $E[U_{it}\midV_{it},\vector X_{it},\alpha_i]=0$, and $V_{it}$ the residual of a non-linear additive noise treatment model depending on fixed effect $\gamma_i$ and satisfying $E[V_{it}|\vector X_{it},\gamma_i]=0$. The focus will be on the PO-PLPR rather than the PLPR model because, as shown below, $l_1$ is a conditional expectation over $Y_{it}$ whereas $g_1$ is a conditional expectation over counterfactual $Y_{it}(0)$ that must generally be learnt iteratively. Moreover, the orthogonalized score on which DML estimation of $\theta_0$ is based involves $V_{it}$ whether it is derived under PLPR or PO-PLPR (see Section (ref)).
We now detail the underlying assumptions required for $\theta_0$ to be feasibly estimable with a causal interpretation. Note that some of these assumptions are stated in terms of stochastic conditional independence ${\perp\!\!\!\perp}$ to simplify the development. Recalling that $\xi_i$ represents the omitted time-invariant confounding variables, the first two assumptions can be stated as follows:
Assumptions (ref)-(ref) together ensure that the panel is static and that any observed lag dependence is due to non-causal autocorrelation. Specifically, these assumptions ensure that the joint distribution of outcomes and treatments given the time-varying predictors and omitted time-invariant influences satisfies
where $p(.)$ denotes a density function for the (conditional or joint) distribution indicated by its arguments. Moreover, there is no need to model the distribution of the time-varying predictors, and the {\em initial conditions problem}, which would arise were the panel study to have started after the joint process began, is avoided.
To give $\theta_0$ a causal interpretation, we require that the data generating process for the potential outcomes and treatment satisfies the following conditions:
Assumption (ref) states that treatment selection at wave $t$ is (strongly) ignorable given $\vector X_{it}$ and latent $\xi_i$. Assumption (ref) requires that $Y_{it}(d)-Y_{it}(0)$ is constant or varies between individuals independently of $\vector X_{it}$ and $\xi_i$. Alternatively, under the weakly ignorable assumption that only $Y_{it}(0)\indep{}D_{it} \vector X_{it},\xi_i$, target parameter $d\theta_0=E[Y_{it}(d)-Y_{it}(0)\mid D_{it}=d]$. However, we assume throughout that the strong version holds.
The final assumption is that the combined effect of the confounding variables is additively separable into that of the observed confounding variables and that of the omitted time-invariant confounding variables.
Assumptions (ref) and (ref) imply that $E[Y_{it}\mid D_{it},\vector X_{it}, \xi_i]=E[Y_{it}(0)\mid\vector X_{it},\xi_i]+D_{it}\theta_0$. Hence, Assumption (ref)(ref) leads instantly to PLPR model (ref). Moreover, Assumptions (ref) and (ref) also imply that $E[Y_{it}(0)\mid \vector X_{it}, \xi_i]= E[Y_{it}\mid\vector X_{it},\xi_i]- E[D_{it}\mid\vector X_{it},\xi]\theta_0$, from which it follows from both Assumptions (ref)(ref) and (ref)(ref) that PO-PLPR model ((ref))-((ref)) holds because $E[Y_{it}\mid\vector X_{it},\xi_i]=l_1(\vector X_{it})+\alpha_i$, where $l_1(\vector X_{it})=g_1(\vector X_{it})+m_1(\vector X_{it})\theta_0$ and $\alpha_i=\alpha^*_i+\gamma_i\theta_0$. \\
In the description of each approach, below, we focus on procedures for the homogenous treatment effects case. However, all of these procedures can be extended to heterogeneous treatment effects provided the analyst is prepared to specify a (finite-dimensional) parametric heterogeneity model. A summary of how this is done can be found in Section (ref).
\setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{lemma}{0} \setcounter{proposition}{0}
We propose two approaches for estimating $\theta_0$ from PO-PLPR model (ref)-(ref). The presence of the time-invariant fixed effects means that only $l_1(\vector X_{it})+E[\alpha_i|\vector X_{it}]$ and $m_1(\vector X_{it})+E[\gamma_i|\vector X_{it}]$ can be learnt from the observed data directly through $E[Y_{it}\mid\vector X_{it}]$ and $E[D_{it}\mid\vector X_{it}]$, respectively, so the first challenge is to remove the fixed effects just as it is in the linear case. As such, the approaches we propose are adaptations of existing techniques for {\em linear} panel data models. The first is based on correlated random effects (CRE) and the second on two varieties of data transformation:\ first-difference (FD) and within-group (WG).\footnote{Under the random effects assumption, where $E[\alpha_i\mid\vector X_{it}]=E[\alpha_i]=0$ and $E[\gamma_i\mid\vector X_{it}]=E[\gamma_i]=0$, $l_1$ and $m_1$ can be straightforwardly learnt from the observed data, but we assert that this is unlikely to hold in practice and, hence, causal inference would not be credible. However, were the analyst prepared to make the random effects assumption, we note that sela2012 developed an algorithm for using tree-based learners to estimate non-causal partially linear regression models.}
Crucially, to understand our contribution, it is important to note that the use of CRE or data transformations does not in itself lead to a consistent estimator of $\theta_0$. There are further challenges, created by the non-linearity of the nuisance functions $l_1$ and $m_1$, to be overcome. Specifically, we do not solve these challanges by {\em a priori} taking $l_1$ and $m_1$ to be accurately approximated by high-dimensional sparse functions of the form $l_1(\vector X_{it})\approx \vector B^{\prime}_{it}\pmb{\psi}$ and $m_1(\vector X_{it})\approx \vector B^{\prime}_{it}\pmb{\phi}$, where vector ${\bf B}_{it}={\bf B}(\vector X_{it})$ is an analyst-specified high-dimensional data dictionary of linear and non-linear terms, $\pmb{\psi}$ and $\pmb{\phi}$ are vectors of parameters with ${\rm dim}(\pmb{\psi})>>||\pmb{\psi}||_0$ and ${\rm dim}(\pmb{\phi})>>||\pmb{\phi}||_0$, and $||{\bf v}||_0$ is the number of non-zero elements in arbitrary vector ${\bf v}$. The FD and CRE procedures respectively developed by klosin2022 and semenova2023 rely explicitly on such representations. However, our procedures are suitable for DML based on ensemble learning strategies.
Correlated random effects (CRE) models are extensions of fixed effect models in which the fixed effect is replaced by a model for its dependence on the predictor variables. CRE models can depend on any function of $\{\vector X_{it}:t=1,\dots,T\}$ (e.g.\ wooldridge2020) but practice tends to focus on {\em exchangeable} functions and especially functions of $\overline{\vector X}_i=T^{-1}\sum_{t=1}^T\vector X_{it}$ following mundlak1978 who first demonstrated the equivalence of CRE and fixed effects estimators when the effects of $\vector X_{it}$ in the fixed-effects model and the effects of ${\overline {\vector X}}_i$ in the model for the fixed effect are all linear.
In general, the fixed-effect model must be correctly specified, but a major advantage of the Mundlak device for linear panel models, where the effects of $\vector X_{it}$ are linear, is that least-squares estimators of the coefficients remain consistent even if the linear model for $E[\alpha_i|{\overline {\vector X}}_i]$ is mis-specfied, provided that $\alpha_i$ and ${\overline {\vector X}}_i$ are correlated. However, this robustness is lost for PO-PLPR model ((ref))-((ref)) because $l_1$ and $m_1$ are non-linear. This is also the case for the linear panel model with random coefficients for heterogeneous treatment effects of wooldridge2019 and the non-linear probit model of wooldridge2020.
In the spirit of wooldridge2020, we derive a CRE model based on PO-PLPR model ((ref))-((ref)). Suppose that the data generating process satisfies Assumptions (ref)-(ref) and the fixed effects follow non-linear additive noise models $\alpha_i= \omega_\alpha({\overline{\vector X}_i})+a_i$ and $\gamma_i = \omega_\gamma({\overline{\vector X}_i}) + c_i$, where $a_i$ and $c_i$ are random effects satisfying $(a_i,c_i)\ {\perp\!\!\!\perp}\ L_T(\vector X_i)$ and $a_i\ {\perp\!\!\!\perp}\ c_i$. Then the CRE model with nuisance parameters $\widetilde{l}_1$ and $\widetilde{m}_1$ is
where $E[U_{it} \midV_{it},\vector X_{it},{\overline{\vector X}_i},a_i\big]=E[V_{it}|\vector X_{it},{\overline{\vector X}_i},c_i]=0$, $\widetilde{l}_1(\vector X_{it},{\overline{\vector X}_i})=l_1(\vector X_{it}) + \omega_{\alpha}({\overline{\vector X}_i})$ and $\widetilde{m}_1(\vector X_{it},{\overline{\vector X}_i}) = m_1(\vector X_{it})+\omega_\gamma({\overline{\vector X}_i})$. The random effects simply capture autocorrelation between observations on the same individual. Crucially, $\widetilde{l}_1$ and $\widetilde{m}_1$ can be learnt directly from $E[Y_{it}\mid\vector X_{it}, {\overline{\vector X}_i}]$ and $E[D_{it}\mid\vector X_{it}, {\overline{\vector X}_i}]$, respectively.\footnote{semenova2023 also use a CRE model in their formulation for dynamic models. In our static panel case with only {\em fixed} $\vector X_{it}$ predictors, this would be equivalent to specifying $l_1(X_{it})+\alpha_i=\psi_0 \overline{X}_i + \alpha_i$ with $(\alpha_1,\dots,\alpha_N)$ assumed to be a {\em weakly sparse} vector following a Mundlak-style model which can be estimated using LASSO. Alternative models for weakly sparse fixed effects have been developed by kock2016 and kock2019 which can be estimated using {\em de-sparsified} LASSO.}
The second approach we consider follows more conventional techniques for panel data by transforming the data to remove entirely the fixed effects from the analysis.
Let $W_{it}$ be a generic random variable and $Q$ a panel data transformation operator such that $Q(W_{it})=Q_t(W_{i1},\dots,W_{iT})$ is a function of the random variable at wave $t$ and the remaining realizations for individual $i$. We consider two such transformations:\ the WG (within-group) or time-demeaning transformation, i.e., $Q(W_{it})=W_{it}-\overline{W}_i$, where $\overline{W}_{i}= T^{-1} \sum_{t=1}^T W_{it}$; and the FD (first-difference) transformation, i.e., $Q(W_{it}) = W_{it}-W_{it-1}$ for $t=2,\dots,T$.
The reduced-form model for $Q(Y_{it})$ and $Q(V_{it})$ under PO-PLPR model (ref)-(ref) is
which does not depend on fixed effects $\alpha_i$ and $\gamma_i$ because $Q(\alpha_i)=Q(\gamma_i)=0$.
The challenge of learning the transformed nuisance functions $Q({l}_1(\vector X_{it}))$ and $Q({m}_1(\vector X_{it}))$ is complicated by the non-linearity of $l_1$ and $m_1$ because, as set out previously, we do not {\em a priori} assume that the nuisance functions admit high-dimensional sparse representations of $l_1$ and $m_1$ where $Q(l_1(\vector X_{it})\approx Q({\bf B}^{\prime}_{it})\pmb{\psi}$ and $Q(m_1(\vector X_{it}))\approx Q({\bf B}^{\prime}_{it})\pmb{\phi}$ could be learnt directly as per the FD procedure proposed by klosin2022.
If Assumption (ref) is implausible or the inferential target of the analysis is heterogeneity of the treatment effects itself then further modelling is required. Completely relaxing Assumption (ref) results in $\theta_0d=E[Y_{it}(d)-Y_{it}(0)\mid\vector X_{it},\xi_i]$, that is, conditional average partial effects or conditional average treatment effects (CATE) dependent on all of the confounding variables. Focussing on the binary treatment case for simplicity, extending Assumption (ref)(ref) so that $E[Y_{it}(1)\mid\vector X_{it},\xi_i]$ is also additively separable would lead to an additively separable CATE $\theta_0=h_0(\vector X_{it})+\delta_i$, where $\delta_i=\delta(\xi_i)$ is what wooldridge2019 refers to as a {\em random coefficient}.
CRE model ((ref))-((ref)) can be specified in terms of CATEs if Assumption (ref) fails. Temporarily discarding vector notation, if $E[Y_{it}(1)-Y_{it}(0)\mid X_{it},\xi_i]=h_0(X_{it})$ then our procedure allows parametric model $h_0(X_{it})=S_{it}\theta_0$, where $S_{it}=S(X_{it})$ captures treatment effect heterogeneity though interaction $S_{it}V_{it}$ replacing $V_{it}$ in ((ref)). If the heterogeneity also depends on $\xi_i$ then, following wooldridge2019, the analyst can proceed by additionally specifying a further linear model for $E[\delta_i\mid \overline{X}_i]$ so that $S_{it}$ is extended to $S_{it}=S( X_{it},\overline{X}_i)$ with no interactions between $X_{it}$ and $\overline{X}_i$ included. As noted by wooldridge2019, this model must be correctly specified to induce a consistent estimator, despite having the appearance of the Mundlak device.
Finally, for the transformation models, the failure of Assumption (ref) also requires the analyst to specify parametric model $h_0(X_{it})=S_{it}\theta_0$. Then (ref) becomes $Q(Y_{it}) = Q(V_{it} S_{it})\theta_0 + Q\big({l}_1(\vector X_{it})\big) + Q(U_{it})$. In contrast to the CRE, where the analyst must assume $\delta_i=0$ or specify the model $E[\delta_i\mid{\overline{\vector X}}_i]$, random coefficient $\delta_i$ cancels out.
Our approach contrasts with those of klosin2022 and semenova2023, both of which are based on learning the high-dimensional sparse representation of $h_1$ using a suitably specified data dictionary. The former is a FD procedure that removes $\delta_i$ from consideration provided the CATE is additively separable as above, whereas the latter implicitly assumes that $\delta_i=0$.
\setcounter{equation}{0}
We now describe our DML procedures for the fixed-effects estimators presented in Section (ref). The objective is to make inferences on the target parameter $\theta_0$ given a suitable predictions of nuisance functions $\bm{\eta}_0$ based on the observed data $W_i = \{W_{it}: t=1,\dots,T\}$, where $W_{it} = \{Y_{it},D_{it},\vector X_{it}, \overline{\vector X}_i\}$ for CRE, and $W_{it} = \{Q(Y_{it}),Q(D_{it}),Q(\vector X_{it})\}$ for the transformation approach (noting that $W_{i1}\equiv\varnothing$ if $Q$ is the FD transformation).
\setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{lemma}{0} \setcounter{proposition}{0} The first stage of DML involves using a suitably chosen learner (or combination of learners) to predict flexibly the unknown nuisance parameters. The nuisance functions $\bm{\eta}_{0i}=\bm{\eta}_0(\vector X_{i1},\ldots,\vector X_{iT})$ vary between individuals and so the task is to learn $\bm{\eta}_0$. As discussed previously, the presence of the fixed effects and non-linear functional forms pose challenges which need to be addressed. Below we describe different procedures for learning $\bm{\eta}_0$.
The estimation of $\theta_0$ based on model ((ref))-((ref)) requires learning $\widetilde{l}_1(\vector X_{it},\overline{\vector X}_i)$ from the data $\{Y_{it},\vector X_{it},\overline{\vector X}_i:t=1,\dots,T\}_{i=1}^N$, and obtaining predicted residual $\widehat{V}_{it}$ to plug into Equation (ref). Simply learning $\widetilde{m}_1(\vector X_{it},\overline{\vector X}_i)$ from the data $\{D_{it},\vector X_{it},\overline{\vector X}_i:t=1,\dots,T\}_{i=1}^N$ and using $\widehat{V}_{it}=D_{it} - {\widehat{m}_1}(\vector X_{it},\overline{\vector X}_i)$ would ignore individual random effect $c_i$ and so introduce bias. Hence, we propose the following approach: \setcounter{bean}{0}
where $m^*_1(\vector X_{it},\overline{\vector X}_i,\overline{D}_{i})=E[D_{it}\mid \vector X_{it},\overline{\vector X}_i]+c_i$. To capture $c_i$, it is essential to include the individual-specific treatment mean $\overline{D}_{i}$ to ensure $V_{it}$ is included in the prediction of $m_1^*$. For example, if the conditional distribution $D_{i1},\dots,D_{iT}\mid\vector X_{i1},...,\vector X_{iT}$ is multivariate normal with $E[D_{it}|\vector X_{it},\overline{\vector X}_i] = \widetilde{m}_1(\vector X_{it},\overline{\vector X}_i)$ and a homoskedastic variance-covariance matrix with random effects structure (induced by marginalizing over $c_i$) then $E(D_{it}\mid\vector{X}_{it},\overline{\vector X}_i,\overline{D}_i)=\widetilde{m}_1(\vector X_{it},\overline{\vector X}_i)+c_i=m_1^*$ as required, but it is generally necessary to follow the steps above. More details are given in Section (ref) of the Online Supplementary Information.
\
First, we consider learning based on the transformed data alone. Generally, $Q\big({l}_1(\vector X_{it})\big)\ne {l}_1\big(Q(\vector X_{it})\big)$ and $Q\big({m}_1(\vector X_{it})\big)\ne {m}_1\big(Q(\vector X_{it})\big)$ unless $l_1$ and $m_1$ are linear. This means a consistent estimator cannot straightforwardly be constructed from the transformed data. However, we propose the following procedures:
(a) Approximate procedure. Approximate model ((ref))-((ref)) by
where $l_1$ and $m_1$ can be learnt from the transformed data $\{Q(Y_{it}),Q(\vector X_{it}): t=1,\dots,T\}_{i=1}^N$ and $\{Q(D_{it}),Q(\vector X_{it}): t=1,\dots,T\}_{i=1}^N$, respectively. This approach can produce good approximations of some non-linear nuisance functions but it is not robust and we were able to choose functions where it performed poorly.
(b) Exact procedure. Consider first the FD transformation $Q(Y_{it})=Y_{it}-Y_{it-1}$. Then under Assumptions (ref)-(ref), $E[Y_{it}-Y_{it-1} \mid \vector X_{it-1} , \vector X_{it}]=\Delta l_1(\vector X_{it-1} ,\vector X_{it})$ and $E[D_{it}-D_{it-1}\mid\vector X_{it-1} ,\vector X_{it}]=\Delta m_1(\vector X_{it-1} ,\vector X_{it})$, so that $\Delta l_1(\vector X_{it-1} ,\vector X_{it})$ can be learnt using $\{Y_{it}-Y_{it-1},\vector X_{it-1} , \vector X_{it}: t=2,\dots,T\}_{i=1}^N$ and $\Delta m_1(\vector X_{it-1} ,\vector X_{it})$ can be learnt using $\{D_{it}-D_{it-1}, \vector X_{it-1} , \vector X_{it}:t=2,\ldots,T\}_{i=1}^N$. However, for WG transformation $E[Y_{it}-\overline{Y}_{i}\mid L_T(\vector X_i) ]=l_1(\vector X_{it}) - T^{-1}\sum_{s=1}^T l_1(\vector X_{is})$ and $E[D_{it}-\overline{D}_{i}| L_T(\vector X_i)]=m_1(\vector X_{it}) - T^{-1}\sum_{s=1}^T m_1(\vector X_{is})$, the dimension of the learning problem is $O(T^2)$ times greater than that of the FD transformation and so unfeasible for non-trivial $T$.
(c) Hybrid procedure. Suppose the conditions set out in Section (ref) for model (ref)-(ref) hold. Then $\widetilde{l}_1(\vector X_{it},\overline{\vector X}_i)=l_1 (\vector X_{it})+\omega_\alpha(\overline{\vector X}_i)$ and $\widetilde{m}_1(\vector X_{it},\overline{\vector X}_i)=m_1(\vector X_{it})+\omega_\gamma(\overline{\vector X}_i)$ learnt from the sample data satisfy $Q(\widetilde{l}_1(\vector X_{it},\overline{\vector X}_i))=Q(l_1(\vector X_{it}))$ and $Q(\widetilde{m}_1(\vector X_{it},\overline{\vector X}_i)\big)=Q(m_1(\vector X_{it})$. The hybrid approach can be used with both WG and FD transformations.
We elaborate on the justification for these procedures in Sections (ref) and (ref) of the Online Supplementary Information.
The second stage of DML requires specifying an orthogonal score function for the structural parameter of interest to be solved after the nuisance functions have been plugged in. Following chernozhukov2018, we construct a generic Neyman orthogonal score function for panel data that accounts for (a) the presence of the unobserved individual heterogeneity, and (b) non-linearity of the nuisance functions. Its properties allow us to obtain valid inferences using DML algorithms, provided the ML algorithms converge at rate $N^{1/4}$. We wish to make inference on the target parameter $\theta_0$ given a suitable estimate of the nuisance parameter $\bm{\eta}_0$ from the data. We repurpose the definition of $W_{it}$ to define the entire dataset as $W = \cup^{N}_{i=1}W_i$, where $W_i=\{W_{it}: t=1,\dots,T\}$, $W_{it} = w\{Y_{it},D_{it},\vector X_{it}\}$ and $w$ is a function or transformation of the data (possibly the identity) chosen by the analyst to implement one of the estimation approaches discussed in Section (ref).
For individual $i$, let ${\bf r}_i=(r_{i1},...,r_{iT})^\prime$ be the $T$ residuals of the PO-PLPR model (either ((ref))-((ref)) or ((ref))-((ref))), $X_i=X(\vector X_{i1},\dots,\vector X_{iT})$ the appropriately chosen set of predictor variables, and $\Sigma(X_i)=E[{\bf r}_i{\bf r}'_i\mid X_i]$ the (potentially) heteroskedastic residual variance-covariance matrix. Note that $r_{it}$ differs from $U_{it}$ in that it potentially incorporates the model random effect and other sources of autocorrelation.
Then the Neyman orthogonal score has the form
where row vector ${\bf V}^\perp_i=(V_{i1}^\perp,\ldots, V^\perp_{iT})$ contains the {\em orthogonalized regressors} chosen to ensure Neyman orthogonality.
A derivation of (ref) and its generalization to estimating conditional average treatment effects (CATEs) in the presence of heterogeneous treatment effects is provided in Section (ref) of the Online Supplement.
Suppose that suitable ML algorithms are available to learn the nuisance functions and obtain estimate ${\widehat\bm{\eta}}$ on a pathway converging to $\bm{\eta}_0$ at rate $N^{1/4}$. We now set out the DML2 chernozhukov2018 algorithm we use to make post-regularized inference about $\theta_0$ based on Neyman-orthogonal score ((ref)).
DML uses cross-fitting to control the overfitting or regularization bias introduced by plugging in ${\widehat\bm{\eta}}$ to ((ref)). A key feature of our procedure is that we have adapted DML to use block-k-fold cross-fitting for panel data. This involves treating the time-series for each sample unit, $W_i$, as {\em k-fold}-sampling units when partitioning the sample into $K$ equi-sized folds. Because $W_i$ satisfies the independent and identically distributed assumption, whereas $W_{it}$ does not, block cross-fitting accounts for the within-individual autocorrelation structure and so avoids having to make further post-fitting adjustments for autocorrelation to the estimmated standard errors as did semenova2023.
Letting ${\cal W}=\{1,\ldots, N\}$ denote the indices of the sample units, our DML estimation procedure with block-k-fold cross-fitting can be set out as follows:
Because score ((ref)) is linear in $\theta$, our proposition is that $\widehat{\theta}^{DML}$ converges to normality according to the vectorized equivalent of chernozhukov2018 under the regularity conditions set out by chernozhukov2018.
\setcounter{equation}{0} The primary focus of our Monte Carlo simulation study is to assess the performance of the DML procedures we propose in terms of bias and precision for different ML algorithms. We generated data under DGPs satisfying the assumptions required for the PO-PLPR model to hold for three different pairs of nuisance functions:\ linear, non-linear (smoothly continuous with no interactions), and non-linear (discontinuous with interactions). These are described along with a detailed discussion of the results in Online Supplement (ref). The DGP satisfies the usual sparsity constraints because only two of the thirty predictors included in the analysis have non-zero effects.
Table (ref) provides a summary of the main results under the most complex non-linear DGP for four of the procedures introduced above (Panels A-C). For each procedure, standard OLS estimates are compared with DML using four different learners for the non-linear nuisance functions. Across all three procedures, OLS can exhibit large bias and DML-LASSO typically outperforms the other learners in terms of both bias and precision, with the exception of the WG estimator (Panel C). In practice, LASSO requires the analyst to specify a sufficiently rich dictionary of non-linear terms; had only linear terms been included, for example, its performance would have been closer to OLS. Gradient boosting outperforms the other tree-based approaches (CART and RF) the performance of which in terms of bias and precision (SE/SD) is far inferior. A closer investigation, described in Online Supplement (ref), found the sampling distributions of tree-based ${\hat{\theta}}_{DML}$ to be highly non-normal (see Figure (ref)) so that inferences based on first-order asymptotic results are unreliable. Because trees are sensitive to hyperparameter choice (tree depth, etc.), we experimented with an alternative strategy for hyperparameter tuning (described in Online Supplement (ref)) that led to normal sampling distributions but larger biases. An extensive discussion on hyperparameter tuning can be found in Section (ref) of the Online Supplement.
The superior performance of LASSO is indicated by the Root Mean Square Errors (RMSEs) of the estimator, model, and nuisance parameters. These results support the use of an ensemble learning strategy, in which the best-performing learner is chosen because, for our example, this would have selected LASSO and ensured reliable inference.
There is little to choose between the CRE, FD (Exact) and WG (Approximate) procedures. This is not unexpected because the DGP in the simulation study satisfies the additional assumptions for modelling the unobserved heterogeneity with Mundlak-like CRE, but the results are slightly different because these procedures involve learning different nuisance functions. More generally, the FD (Exact) estimator is the most robust procedure because it does not rely on any additional assumptions. In comparison, the WG (Approximation) method performed poorly no matter which learner was chosen because the approximation error induced by using the linearly transformed data (the time-demanded transformation) is large for the non-linear DGP.
\setcounter{equation}{0}
We reanalyse fazio2023's study of the impact of the National Minimum Wage (NMW) on voting behaviour in the UK. The data used in the original investigation come from the British Household Panel Survey (BHPS) comprising $4,927$ working individuals (aged 18-64) between waves 1 and 16 (from 1991 until 2007).\footnote{BHPS is a longitudinal survey study for British households that run from 1991 until 2009. The online replication package provides the instructions for data access and the codes to run the analysis in this section. DML estimation is conducted in R using the XTDML package in the replication package (or accessible in its latest version at \url{https://github.com/POLSEAN/XTDML} at the time of writing) which is built on DoubleML by bach2024. \nocite{bhps}} The treatment is measured by the question `Were you paid the minimum wage' asked to those responding at Wave 9 in 1999. The authors use least squares to estimate the average treatment effect on the treated (ATT) between the NMW policy on various outcomes, including voting for conservative political parties which is our outcome of interest. fazio2023 compared the results of four regression specifications with different sets of covariates and fixed effects to capture all potential confounders, showing that the estimates did not vary across these, hence, concluding that workers who were paid the NMW rate were more likely to support conservative parties.
We revisit Specification (2) of Table 5 from the original paper that includes all control variables, regions and wave fixed effects, by using DML procedure to fit the following PO-PLPR model
where $Vote_{it}$ is an indicator corresponding to one if respondent $i$ voted for conservative parties in wave $t$, and zero otherwise; $NMW_{it}$ is the treatment variable which is equal to one if the respondent's hourly pay increased due to the introduction of the NMW in 1999 and the respondent is observed in wave 9 onwards; $\vector X_{it}$ are the confounding variables whose functional form is ex ante unknown; and $\theta$ is the target parameter. The original base control variables included information about the age, education, marital status, household size, and income of other household members for respondent $i$ in wave $t$.
We add thirty additional variables to the original specification, providing information on demographic characteristics, socio-economic status, employment and work-related variables, and ideology of the respondents. The complete list of confounders is reported in Table (ref) of the Online Supplement (ref). We also increase the sample size by including two additional waves for years 2008 and 2009 (waves 17 and 18, respectively). Our final sample includes $9,922$ working individuals. The functional form of the confounding variables is learnt using four different base learners:\ LASSO with a dictionary of non-linear terms (i.e., polynomials of order three and interaction terms of each raw variable), CART, RF, and boosted trees. The hyperparameter tuning of the base learners used in the implementation is discussed in the Online Supplement (ref).
The point estimates of NMW effect on voting conservative parties are reported in Table (ref). Column (1) shows OLS estimates (i.e.\ using standard panel estimation based on linear models) while the remaining columns contain the results using DML with different learners. Panel A displays the CRE estimates, Panel B the FD (Exact) estimates, and Panel C the WG (Approximation) estimates. The regression equations include the individual means of all included control variables when the CRE is used, and one-period lagged variables of the controls when FD exact is used. Cluster-robust standard errors are reported in parenthesis with clustering at the respondent level.
The results exhibit considerable differences in the estimated effects between learners and estimators but, we argue, are consistent with those from a DGP which is either close to linear or smoothly non-linear as in our simulation study (see Online Supplement (ref)). The first difference we consider is between the FD (Exact) and CRE estimators. The former should be the most robust because it does not rely on the Mundlak-type models for fixed effects and also allows for different sets of omitted variables fixed only between distinct wave pairs $t-1$ and $t$ (which could be represented by adding parameters like $\alpha_{i(t-1,t)}$ to the models). The magnitudes of the FD (Exact) point estimates are (a) smaller for OLS and LASSO and (b) more stable across all learners than those based on CRE. The RMSEs for the outcome ($l$) and treatment ($m$) models for FD (Exact) are smaller than those for CRE, which indicates the learners are more effective at learning nuisance functions of $\vector X_{it}$ and $\mathbf{x}_{it-1}$ than of $\vector X_{it}$ and $\overline{\mathbf{x}}_i$ and so more likely to have errors bounded by $N^{1/4}$. Moreover, the consistency across OLS and the learners would suggest a DGP that is linear or smoothly non-linear. Second, the WG (Approximate) estimates are also consistent across learners and larger than those obtained using FD (Exact). The differences between FD (Exact) and WG (Approximate) are also commensurate with the additional robustness of the former. The CRE estimates obtained using tree-based learners, compared with all other estimates, are both larger (in absolute value) and less precise, with the DML-Boosting estimate particularly imprecise. However, the RMSEs indicate the tree-based algorithms do a better job learning the nuisance functions in $\vector X_{it}$ and ${\bar {\bf x}}_i$ than LASSO, and an ensemble strategy would have selected the estimate based on the random forest (DML-RF) for comparison with OLS and so pointed to a slightly larger but less significant estimate than the best estimate obtained using WG (Approximate). We finally reiterate the preference of FD (Exact) in this application, and the importance of comparing different procedures and learners in reaching this conclusion.
DML is a powerful tool for leveraging the power of ML for robust estimation of treatment effects, or policy-intervention. We provide estimation tools for applying DML when panel data are available which practitioners can use in place of existing ones or in a complementary way to test the robustness of their results to non-linearity. Moreover, the nature of static panels means that our procedures extend naturally to unbalanced panels provided that the non-response at each wave is conditionally independent of $Y_{it}$ given $\mathbf{X}_{it}$ and $\xi_i$. Further work on relaxing this assumption for CRE could build on wooldridge2019.
In general, although the three approaches we presented can be considered complementary to each other, the suitability of each will depend on the specific empirical example and the assumptions that the analyst is prepared to make. However, following our empirical studies, our recommendation for practice is to employ the FD (exact) approach because the ML algorithms we consider proved effective at estimating the nuisance functions involving predictors from waves $t$ and $t-1$, the random coefficient treatment effect heterogeneity of wooldridge2019 is automatically accounted for, and it places fewer assumptions on the the distribution of the fixed effects.
From our simulation study, we also recommend using multiple (base) learners for DML stage one in the context of ensemble learning. In our simulations, the LASSO with an extended dictionary of non-linear terms performed best across all scenarios. However, despite Lipschitz continuity and weak sparsity conditions holding, we found tree-based learners to perform poorly in terms of bias (of the point and interval estimates) and in terms of standard deviation and normality, even after making extensive efforts to implement adaptive hyperparameter tuning. We do not claim this to be a general result (there are many successful applications of tree-based learners in the literature), but we strongly recommend that analysts use an ensemble strategy involving multiple learners (the tree-based learners could perform better than the others in some scenarios) as is standard practice in other disciplines where ML and other AI-based methods are widely used.
Finally, our method is limited when there is treatment effect heterogeneity but the analyst still wishes to target the ATE rather than the CATE. Further work on extending these procedures, based on adapting the {\em interactive model} of chernozhukov2018 to panel data, would therefore be of great value for practice.
The authors are grateful to Andreas Alfons, Anna Baiardi, Thomas Cornelissen, Riccardo Di Francesco, Maria Grith, Omar Hussein, Damian Machlanski, Spyros Samothrakis, David Zentler-Munro, Wendun Wang, and the participants to the Annual MiSoC Workshop, RSS International Conference 2023, RCEA-ICEEF 2024, IPDC 2024, the internal seminars in ISER and in the Econometric Institute at the Erasmus University Rotterdam for the helpful comments and suggestions. We also thank the editor Christoph Rothe and three anonymous referees, whose comments significantly improved the final version of this paper. This research was funded by the UK Economic and Social Research Council award ES/S012486/1 (MiSoC). The authors also acknowledge the use of the High Performance Computing Facility (Ceres) and its associated support services at the University of Essex in the completion of this work.
\setcounter{section}{0}