Individual Causal Inference
Using Panel Data With Multiple Outcomes

Suppose that we observe $K$ outcomes in domain $\mathcal{K}=\{1,2,\dots,K\}$ for $N$ individuals or units over $T\geq 2$ time periods, where a domain refers to a collection of related outcomes that depend on the same set of observed covariates and unobserved characteristics. For example, health-related outcomes may be affected by observed covariates such as age, education, occupation and income, as well as unobserved individual characteristics such as genetic inheritance, health habits and risk preferences. Assume that the $N_{1}$ individuals in the treated group $\mathcal{T}$ receive the treatment at period $T_{0}+1\leq T$ and remain treated afterwards, while the $N_{0}=N-N_{1}$ individuals in the control group $\mathcal{C}$ remain untreated throughout the $T$ periods. Denoting the binary treatment status for individual $i$ at time $t$ as $D_{it}$, we have $D_{it}=1$ for $i\in\mathcal{T}$ and $t>T_{0}$, and $D_{it}=0$ otherwise.

Following the “Rubin Causal Model” (Rubin,, 1974), the treatment effect on outcome $k\in\mathcal{K}$ for individual $i$ at time $t$ is given by the difference between the treated and untreated potential outcomes

where $Y_{it,k}^{1}$ is the treated potential outcome, the outcome that we would observe for individual $i$ at time $t$ if $D_{it}=1$, and $Y_{it,k}^{0}$ is the untreated potential outcome, the outcome that we would observe if $D_{it}=0$. Instead of assuming the unconfoundedness condition, we characterise the two potential outcomes for individual $i$ at time $t$ and $k\in\mathcal{K}$ using the interactive fixed effects models:

where $\boldsymbol{X}_{it}$ is the $r\times 1$ vector of observed covariates unaffected by the treatment, $\boldsymbol{\mu}_{i}$ is the $f\times 1$ vector of unobserved individual characteristics, $\boldsymbol{\beta}_{t,k}^{1}$ and $\boldsymbol{\lambda}_{t,k}^{1}$ are the $r\times 1$ and $f\times 1$ vectors of coefficients of $\boldsymbol{X}_{it}$ and $\boldsymbol{\mu}_{i}$ respectively for the treated potential outcome, $\boldsymbol{\beta}_{t,k}^{0}$ and $\boldsymbol{\lambda}_{t,k}^{0}$ are the coefficients for the untreated potential outcome, and $\varepsilon_{it,k}^{1}$ and $\varepsilon_{it,k}^{0}$ are the idiosyncratic shocks.

The regularity conditions on the observed covariates and the unobserved individual characteristics are stated in Assumption 1, and the assumptions on the idiosyncratic shocks are given in Assumption 2.

Given the models for $Y_{it,k}^{1}$ and $Y_{it,k}^{0}$ in (2) and (3), the individual treatment effect is identified by the observed covariates and the unobserved individual characteristics, i.e., two persons with the same values for these underlying predictors have identical individual treatment effect. Denote the set of observed covariates and unobserved individual characteristics as $\boldsymbol{H}_{it}=\left[\boldsymbol{X}_{it}^{\prime}\enspace\boldsymbol{\mu}_{i}^{\prime}\right]^{\prime}$, then the individual treatment effect for individual $i$ with $\boldsymbol{H}_{it}=\boldsymbol{h}_{it}$ is given by

which may appear similar to the conditional average treatment effect, but is different by conditioning not only on the observed covariates, but also on the unobserved individual characteristics.⁵⁵5As we assume the parametric models for the potential outcomes in (2) and (3) for all individuals, there is no need to impose additional assumptions on the propensity distribution for the individual treatment effect to be identified on its full support.

Our goal is to estimate the individual treatment effects $\bar{\tau}_{it,k}$, $i=1,\dots,N$. Once we have estimated the individual treatment effects, the estimates of the average treatment effects for heterogeneous subgroups defined by some observed covariates, also known as the conditional average treatment effects, and the estimate for the average treatment effect for the sample or the population are also readily available using the average of the estimated individual treatment effects in the corresponding groups.

As $\boldsymbol{\mu}_{i}$ is not observed, a direct application of least squares estimation to estimate the models in (2) and (3) would suffer from omitted variables bias. Since we have multiple outcomes that depend on the same set of underlying predictors, and we observe the untreated potential outcomes for all individuals prior to the treatment, we can use these pretreatment outcomes to replace $\boldsymbol{\mu}_{i}$ in the models.⁶⁶6This is analogous to the first step of the approach in Hsiao et al., (2012), who predict the posttreatment outcomes using pretreatment outcomes in lieu of the unobserved time factors in a small $N$, big $T$ environment. Stacking the $K$ outcomes observed in $t\leq T_{0}$, we have

where $\boldsymbol{Y}_{it}^{0}$ and $\boldsymbol{\varepsilon}_{it}^{0}$ are $K\times 1$, $\boldsymbol{\beta}_{t}^{0}$ is $K\times r$, and $\boldsymbol{\lambda}_{t}^{0}$ is $K\times f$. Let $\mathcal{P}\subseteq\left\{1,\cdots,T_{0}\right\}$ be a set of $P$ pretreatment periods. We can further stack the outcomes over these periods to get

where $\boldsymbol{\delta}_{i}^{\mathcal{P}}=\left[\cdots\enspace\left(\boldsymbol{\beta}_{s}^{0}\boldsymbol{X}_{is}\right)^{\prime}\enspace\cdots\right]^{\prime}$ with $s\in\mathcal{P}$ is $KP\times 1$, $\boldsymbol{\lambda}^{\mathcal{P}}$ is $KP\times f$, and $\boldsymbol{\varepsilon}_{i}^{\mathcal{P}}$ is $KP\times 1$.

To be able to recover $\boldsymbol{\mu}_{i}$ from the covariates and outcomes observed in $\mathcal{P}$, we need the following full rank condition, which ensures that there is enough variation in the effects of the unobserved individual characteristics over time or across different outcomes.

Under Assumption 3, we can pre-multiply both sides of equation (6) by $({\boldsymbol{\lambda}^{\mathcal{P}}}^{\prime}\boldsymbol{\lambda}^{\mathcal{P}})^{-1}{\boldsymbol{\lambda}^{\mathcal{P}}}^{\prime}$ to obtain

Substituting (7) into $Y_{it,k}^{0}=\boldsymbol{X}_{it}^{\prime}\boldsymbol{\beta}_{t,k}^{0}+\boldsymbol{\mu}_{i}^{\prime}\boldsymbol{\lambda}_{t,k}^{0}+\varepsilon_{it,k}^{0}$, $t>T_{0}$, and with a little abuse on the notation by omitting the superscript $\mathcal{P}$ on the new coefficients and error term, we have

where

Let $Z=r(P+1)+KP$. If we denote the $Z\times 1$ vector of observables $[\boldsymbol{X}_{it}^{\prime}\cdots\boldsymbol{X}_{is}^{\prime}\cdots{\boldsymbol{Y}_{i}^{\mathcal{P}}}^{\prime}]^{\prime}$ as $\boldsymbol{Z}_{it}$, and the $Z\times 1$ vector of coefficients $[{\boldsymbol{\beta}_{t,k}^{0}}^{\prime}\cdots{\boldsymbol{\alpha}_{st,k}^{0}}^{\prime}\cdots{\boldsymbol{\gamma}_{t,k}^{0}}^{\prime}]^{\prime}$ as $\boldsymbol{\theta}_{t,k}^{0}$, then equation (8) can be abbreviated as

Similarly, substituting (7) into $Y_{it,k}^{1}=\boldsymbol{X}_{it}^{\prime}\boldsymbol{\beta}_{t,k}^{1}+\boldsymbol{\mu}_{i}^{\prime}\boldsymbol{\lambda}_{t,k}^{1}+\varepsilon_{it,k}^{1}$, $t>T_{0}$, we have

where

Under Assumption 2, we have $\mathbb{E}(e_{it,k}^{1}\mid\boldsymbol{H}_{it})=0$ and $\mathbb{E}(e_{it,k}^{0}\mid\boldsymbol{H}_{it})=0$. This suggests using $\widehat{\tau}_{it,k}=\boldsymbol{Z}_{it}^{\prime}(\widehat{\boldsymbol{\theta}}_{t,k}^{1}-\widehat{\boldsymbol{\theta}}_{t,k}^{0})$, where $\widehat{\boldsymbol{\theta}}_{t,k}^{1}$ and $\widehat{\boldsymbol{\theta}}_{t,k}^{0}$ are some estimators of ${\boldsymbol{\theta}}_{t,k}^{1}$ and ${\boldsymbol{\theta}}_{t,k}^{0}$, to estimate $\bar{\tau}_{it,k}$. Note, however, that the error terms $e_{it,k}^{1}$ and $e_{it,k}^{0}$ are correlated with the regressors, since $\boldsymbol{Z}_{it}$ contains $\boldsymbol{Y}_{i}^{\mathcal{P}}$ which is correlated with $\boldsymbol{\varepsilon}_{i}^{\mathcal{P}}$. This renders the OLS estimators biased and inconsistent, which can be seen as a classical measurement errors in variables problem.⁷⁷7This is noted in Ferman and Pinto, (2019) as well, who also suggested using pre-treatment outcomes as instrumental variables to deal with the problem. Our method is also related to the quasi-differencing approach in Holtz-Eakin et al., (1988) and the GMM approach in Ahn et al., (2013). While these studies focus on estimating the coefficients on the observed covariates, our focus is on estimating the individual treatment effects. We thus use the remaining outcomes as instrumental variables for $\boldsymbol{Y}_{i}^{\mathcal{P}}$ to consistently estimate ${\boldsymbol{\theta}}_{t,k}^{1}$ and ${\boldsymbol{\theta}}_{t,k}^{0}$ in each period, which would then allow us to obtain asymptotically unbiased estimates for the individual treatment effects.⁸⁸8We may construct the vectors of regressors and instruments differently under alternative assumptions on the dependence structure of the idiosyncratic shocks. For example, if the idiosyncratic shocks are correlated across time but are independent across outcomes, then we can split different outcomes into regressors and instruments. This would be similar to using the characteristics of similar products (Berry et al.,, 1995) or trading countries (see the Trade-weighted World Income instrument in Acemoglu et al.,, 2008) as instrumental variables. Incorporating more complex structures of the idiosyncratic shocks in the model is left for future research. Since the outcomes depend on about the same set of observed and unobserved individual characteristics, the remaining outcomes are strongly correlated with the outcomes included in $\boldsymbol{Y}_{i}^{\mathcal{P}}$. Additionally, given that the idiosyncratic shocks are independent across time, the remaining outcomes are not correlated with $e_{it,k}^{1}$ or $e_{it,k}^{0}$. Thus, both the relevance and exogeneity conditions are satisfied, and the remaining outcomes can serve as valid instrumental variables.

Let $\boldsymbol{R}_{it}=[\boldsymbol{X}_{it}^{\prime}\cdots\boldsymbol{X}_{is}^{\prime}\cdots{\boldsymbol{Y}_{i}^{-\mathcal{P}}}^{\prime}]^{\prime}$ be the $R\times 1$ vector of instruments, where the $(KT-KP-1)\times 1$ vector $\boldsymbol{Y}_{i}^{-\mathcal{P}}$ comprises the remaining pretreatment outcomes as well as the posttreatment outcomes other than $Y_{it,k}$.⁹⁹9In the special case of $T_{1}=1$ and $T_{0}=1$, we can include $K-1$ pretreatment outcomes as regressors, and use the posttreatment outcomes other than $Y_{it,k}$ as instruments so that $R\geq Z$. Stacking $\boldsymbol{Z}_{it}$, $\boldsymbol{R}_{it}$ and $\boldsymbol{Y}_{it,k}^{0}$ respectively over the $N_{0}$ untreated individuals, we obtain the $N_{0}\times Z$ matrix of regressors $\boldsymbol{Z}_{t}^{0}$, the $N_{0}\times R$ matrix of instruments $\boldsymbol{R}_{t}^{0}$ and the $N_{0}\times 1$ matrix of outcomes $\boldsymbol{Y}_{t,k}^{0}$ for the untreated individuals. We can obtain $\boldsymbol{Z}_{t}^{1}$, $\boldsymbol{R}_{t}^{1}$ and $\boldsymbol{Y}_{t,k}^{1}$ similarly for the $N_{1}$ treated individuals. The GMM estimator for the individual treatment effect $\bar{\tau}_{it,k}$ can then be constructed as

where

with $\boldsymbol{W}^{1}$ and $\boldsymbol{W}^{0}$ being some $R\times R$ positive definite matrices.

The following result shows that the bias of the GMM estimator for the individual treatment effect in (11) goes away as both the number of treated individuals and the number of untreated individuals become larger.

Once we have the estimates for the individual treatment effects, the average treatment effect $\tau_{t,k}=\mathbb{E}\left(\tau_{it,k}\right)$ can be conveniently estimated using the average of the estimated individual treatment effects $\widehat{\tau}_{t,k}=\frac{1}{N}\sum_{i=1}^{N}\widehat{\tau}_{it,k}$, which can be shown to be consistent.

To satisfy Assumption 3, we need the number of pretreatment outcomes that we include as regressors in the model to be at least as large as $f$. Including more pretreatment outcomes may increase the variance of the estimator by increasing the variances of $\widehat{\boldsymbol{\theta}}_{t,k}^{1}$ and $\widehat{\boldsymbol{\theta}}_{t,k}^{0}$, but may also reduce the variance of the estimator when the sample is large and the variances of $\widehat{\boldsymbol{\theta}}_{t,k}^{1}$ and $\widehat{\boldsymbol{\theta}}_{t,k}^{0}$ are small, since

in the prediction error converges in probability to 0 as $KP$ grows.¹⁰¹⁰10Consistency of the individual treatment effect estimator may also be shown by allowing both $N$ and $KP$ to grow, with restrictions on the relative growth rate, e.g., $\frac{KP}{\min\left(\sqrt{N_{1}},\sqrt{N_{0}}\right)}\rightarrow 0$. We do not pursue this path in this study, as the number of pretreatment outcomes in empirical microeconomics that we focus on is usually not large.

To select the number of pretreatment outcomes to include in the model, we follow a model selection procedure similar to that in Hsiao et al., (2012), where for each usable number of pretreatment outcomes, we construct many different models by including a random subset of the pretreatment outcomes as regressors and the remaining outcomes as instruments. We then estimate the models using GMM and obtain the leave-one-out prediction errors for all or a subsample of the individuals. The best set of pretreatment outcomes is chosen as the one that minimises the mean squared leave-one-out prediction error.¹¹¹¹11An alternative way to select the best set of pretreatment outcomes is to use information criteria such as GMM-BIC and GMM-AIC (Andrews,, 1999). To avoid the potential problem of post-selection inference, we may also randomly split the sample into two parts, where we select the best model on one part, and conduct inference on the other.

In addition to the models using only a subset of the pretreatment outcomes, we also consider averaging different models that use the same number of pretreatment outcomes. Since the estimators constructed using only a subset of the pretreatment outcomes are asymptotically unbiased, as long as the number of pretreatment outcomes is larger than $f$, this property is passed on to the averaged estimator. The averaged estimator may also be more efficient as it uses more information in the sample and reduces uncertainty caused by a small number of sample splits.¹²¹²12We stick with simple averaging in this paper. More flexible averaging scheme, e.g., with larger weights on those with smaller out of sample prediction errors, would be an interesting direction for future research. The leave-one-out prediction errors are also averaged over the models, and the best number of pretreatment outcomes to be used for the averaged estimator is similarly determined by minimising the mean squared leave-one-out prediction error.

To measure the conditional variance of the individual treatment effect estimator, $\text{Var}\left(\widehat{\tau}_{it,k}\mid\boldsymbol{H},\boldsymbol{D}\right)$, where $\boldsymbol{H}$ is the matrix of observed covariates and unobserved individual characteristics and $\boldsymbol{D}$ is the matrix of the treatment status for all individuals and all time periods in the sample, we follow Xu, (2017) and employ a parametric bootstrap procedure.

First, we apply our method to all outcomes in all periods to obtain $\widehat{Y}_{it,k}^{1}$ and $\widehat{e}_{it,k}^{1}$ for the treated individuals in the posttreatment periods, and $\widehat{Y}_{it,k}^{0}$ and $\widehat{e}_{it,k}^{0}$ for the untreated individuals in the posttreatment periods and for all individuals in the pretreatment periods. Note that the residuals $\widehat{e}_{it,k}^{1}$ and $\widehat{e}_{it,k}^{0}$ are estimates for $\varepsilon_{it,k}^{1}-{\boldsymbol{\gamma}_{t,k}^{1}}^{\prime}\boldsymbol{\varepsilon}_{i}^{\mathcal{P}}$ and $\varepsilon_{it,k}^{0}-{\boldsymbol{\gamma}_{t,k}^{0}}^{\prime}\boldsymbol{\varepsilon}_{i}^{\mathcal{P}}$, respectively, rather than the idiosyncratic shocks in the original model, $\varepsilon_{it,k}^{1}$ and $\varepsilon_{it,k}^{0}$. Thus, the variance of the individual treatment effect estimator tends to be overestimated using the parametric bootstrap by resampling these residuals, especially when the number of pretreatment outcomes is small.¹⁵¹⁵15See the discussion on equation (14). Correcting for this bias would be a necessary step for future research.

These fitted values of the outcomes can be stacked into a $TK\times 1$ vector $\widehat{\boldsymbol{Y}}_{i}$ for each individual, where $\widehat{\boldsymbol{Y}}_{i}$ for $i\in\mathcal{T}$ contains $\widehat{Y}_{it,k}^{1}$ in the posttreatment periods and $\widehat{Y}_{it,k}^{0}$ in the pretreatment periods, and $\widehat{\boldsymbol{Y}}_{i}$ for $i\in\mathcal{C}$ contains $\widehat{Y}_{it,k}^{0}$ in all periods. The $TK\times 1$ vector of residuals $\widehat{\boldsymbol{e}}_{i}$ can be obtained similarly.

We then start bootstrapping for $B$ rounds:

1. 1.

   In round $b\in\{1,\dots,B\}$, generate a bootstrapped sample as

   $\displaystyle\boldsymbol{Y}_{i}^{(b)}$

   $\displaystyle=\widehat{\boldsymbol{Y}}_{i}+\widehat{\boldsymbol{e}}_{i}^{(b)},\enspace\text{for all}\enspace i,$

   where $\widehat{\boldsymbol{e}}_{i}^{(b)}$ is randomly drawn from $\{\widehat{\boldsymbol{e}}_{i}\}_{i\in\mathcal{T}}$ for $i\in\mathcal{T}$, and from $\{\widehat{\boldsymbol{e}}_{i}\}_{i\in\mathcal{C}}$ for $i\in\mathcal{C}$.¹⁶¹⁶16Since the entire series of residuals over the $T$ periods and $K$ outcomes are resampled, correlation and heteroskedasticity across time and outcomes are preserved (Xu,, 2017).

2. 2.

   Construct $\widehat{\tau}_{it,k}^{(b)}$ for each $i$ using the above bootstrapped sample.

The variance for the individual treatment effect estimator is computed using the bootstrap estimates as

and the $100(1-\alpha)\%$ confidence intervals for $\bar{\tau}_{it,k}$, $i=1,\dots,N$ can be constructed as

where the superscript denotes the index of the bootstrap estimates in ascending order. Alternatively, we can use a normal approximation and construct the confidence intervals as

where $\Phi\left(\cdot\right)$ is the cumulative distribution function for the standard normal distribution, and $\widehat{\sigma}_{it,k}=\sqrt{\text{Var}\left(\widehat{\tau}_{it,k}\mid\boldsymbol{H},\boldsymbol{D}\right)}$.

The variance for the average treatment effect estimator $\widehat{\tau}_{t,k}=\frac{1}{N}\sum_{i=1}^{N}\widehat{\tau}_{it,k}$ and the confidence interval for the average treatment effect $\tau_{t,k}$ can be obtained in similar manners using the bootstrap estimates $\widehat{\tau}_{t,k}^{(b)}=\frac{1}{N}\sum_{i=1}^{N}\widehat{\tau}_{it,k}^{(b)}$, $b=1,\dots,B$.