Sampling-based vs. Design-based Uncertainty in Regression Analysis

Sampling-based vs. Design-based Uncertainty in Regression Analysis

\baselineskip=20pt \setcounter{page}{1}

Introduction

In the dominant approach to inference in the social sciences uncertainty about population parameters is induced by random sampling from the population. Moreover, it is typically assumed that the sample comprises only a small fraction of the population of interest. This perspective is a natural and attractive one in many instances. For example, if one analyzes individual-level data from the U.S. Current Population Survey, the Panel Study of Income Dynamics, or the one percent public-use sample from the U.S. Census, it is natural to regard the sample as a small random subset of the population of interest. In many other settings, however, this sampling perspective is less appropriate. For example, manski2018right write “Random sampling assumptions, however, are not natural when considering states or counties as units of observation.” In this article, we provide an alternative framework for the interpretation of uncertainty in regression analysis regardless of whether a substantial fraction of the population or even the entire population is included in the sample. While our framework accommodates sampling-based uncertainty, it also takes into account design-based uncertainty, which arises when the parameter of interest is defined in terms of the unobserved outcomes that some units would attain under a certain intervention. Design-based uncertainty is often explicitly accounted for in the analysis of randomized experiments where it is the basis of randomization inference (neyman1923, neyman1923; rosenbaum_book, rosenbaum_book ; imbens2015causal, imbens2015causal), but it is rarely explicitly acknowledged in regression analyses or, more generally, in observational studies (exceptions include samii2012equivalencies, samii2012equivalencies; Freedman2008, Freedman2008; lu2016randomization, lu2016randomization; lin, lin).

To illustrate the differences between sampling-based inference and design-based inference, consider two simple examples. In the example of Table (ref), there is a finite population consisting of $n$ units with each unit characterized by a pair of variables, $Y_i$ and $Z_i$. Consider an estimand that is function of the full set of pairs $\{(Y_i,Z_i)\}^n_{i=1}$. Uncertainty about such an estimand arises when we observe the values $(Y_i,Z_i)$ only for a sample, that is, for a subset of the population. In Table (ref), inclusion of unit $i$ in a sample is coded by the binary variable $R_i\in\{0,1\}$. An estimator is a function of the observed data, $\{(R_i,R_iY_i,R_iZ_i)\}^n_{i=1}$. Sampling-based inference uses information about the process that determines the sampling indicators $R_1,\ldots, R_n$ to assess the variability of estimators across different samples. The second and third sets of columns in Table (ref) depicts such alternative samples.

Table (ref) depicts a different scenario in which we observe for each unit in the population the value of one of two potential outcome variables, either $ Y^*_i(1)$ or $Y^*_i(0)$, but not both. The binary variable $X_i\in\{0,1\}$ indicates which potential outcome we observe. Consider an estimand that is a function of the full set of triples $\{(Y^*_i(1),Y^*_i(0),X_i)\}^n_{i=1}$. As before, an estimator is a function of the observed data, the pairs $(X_i,Y_i)$, for $i=1,\ldots,n$, where $Y_i=Y_i^*(X_i)$ is the realized value. Design-based inference uses information about the process that determines the assignments $X_1,\ldots, X_n$ to assess the variability of estimators across different samples. The second and third sets of columns in Table (ref) depicts such alternative samples.

More generally, we can have missing data processes that combine features of these two examples, with some units not included in the sample at all, and with some of the variables not observed for the sampled units. Articulating both the exact nature of the estimand of interest and the source of uncertainty that makes an estimator stochastic is a crucial first step to valid inference. For this purpose, it will be useful to distinguish between descriptive estimands, where uncertainty stems solely from not observing all units in the population of interest, and causal estimands, where the uncertainty stems, at least partially, from unobservability of some of the potential outcomes.

The main formal contribution of this article is to generalize the results for the approximate variance for multiple linear regression estimators associated with the work by eicker, huber, and white1980robust,white1980using,white1982maximum, EHW from hereon, in two ways. First, our framework allows for sampling from a finite population, whereas the EHW results assume random sampling from an infinite population. Second, our framework explicitly takes into account design-based uncertainty. Incorporating these generalizations requires developing a new framework for regression analysis, nesting as special cases the neyman1923 analysis of randomized experiments with binary treatments, as well as the generalizations to randomized experiments with additional regressors in samii2012equivalencies, Freedman2008, and lin. We show that in large samples the widely used EHW robust standard errors are conservative, and only correct in special cases. Moreover, we show that the presence of attributes -- that is, characteristics of the units fixed in our repeated sampling thought experiments -- can be exploited to improve on the EHW variance estimator, and we propose variance estimators that do so. Another advantage of the formal separation into sampling-based and design-based uncertainty is that it allows us to clarify the distinction between the assumptions needed for internal and external validity (shadishcookcampbell, shadishcookcampbell; manski2013public, manski2013public; deaton2010, deaton2010) in terms of these two sources of uncertainty.

Our results are relevant in empirical settings where researchers have a random sample from a finite population and the ratio of the sample size to the population size is sufficiently large so that the proposed finite-population correction matters. Examples of such settings include large-scale experiments muralidharan2017experimentation, settings where the cost of data acquisition motivates the use of random samples keels2005fifteen, as well as analyses based on public-use census samples, like the 2010 Integrated Public Use Microdata Series (IPUMS) data (which is a 10 percent sample of the U.S. Census). More importantly in our view, our results are relevant in empirical settings where it is not natural to think of the data as a random sample from a well-defined population. Instead the researcher may have the entire population, e.g., states or counties as in the manski2018right quote, or the set of all visits to a website, or the researcher may have a convenience sample. In that case our design-based approach to uncertainty provides a coherent interpretation for sampling-based standard errors. It also provides methods that exploit the presence of attributes to calculate improved (i.e., less conservative) standard errors.

A Simple Example

In this section we set the stage for the problems discussed in the current article by discussing least squares estimation in a simple example with a single binary regressor. We make four points. First, we show how design-based uncertainty affects the variance of regression estimators. Second, we show that the standard Eicker-Huber-White (EHW) variance estimator remains conservative when we take into account design-based uncertainty. Third, we show that there is a simple finite-population correction to the EHW variance estimator for descriptive estimands but not for causal estimands. Fourth, we discuss the relation between the two sources of uncertainty and the notions of internal and external validity. Proofs of the results in this section are in a supplementary appendix.

We focus on a setting with a finite population of size $n$. We sample $N$ units from this population, with $R_i\in\{0,1\}$ indicating whether a unit was sampled ($R_i=1$) or not $(R_i=0$), so that $N=\sum_{i=1}^n R_i$. There is a single binary regressor, $X_i\in\{0,1\}$, and $n_x$ (resp.\ $N_x$) is the number of units in the population (resp.\ the sample) with $X_i=x$. Units could be U.S. states and the binary regressor $X_i$ could be an indicator for a state regulation, say the state having a right-to-carry law (RTC), as in manski2018right and donohue2017right. We view the regressor $X_i$ not as a fixed attribute or characteristic of each unit, but as a cause or policy variable whose value could have been different from the observed value. This generates missing data of the type shown in Table (ref), where only some of the states of the world are observed, implying that there is design-based uncertainty. Formally, using the Rubin causal model or potential outcome framework ( neyman1923, neyman1923; rubin1974estimating, rubin1974estimating; holland1986statistics, holland1986statistics; imbens2015causal, imbens2015causal), we postulate the existence of two potential outcomes for each unit, denoted by $Y_i^{*}(1) $ and $Y_i^{*}(0)$. For the RTC example, $Y_i^{*}(1) $ and $Y_i^{*}(0)$ could be state-level crime rates with and without RTC. The realized outcome is

which is the observed state-level crime rate in the RTC example.

In our setting, potential outcomes are viewed as non-stochastic attributes for unit $i$, irrespective of the realized value of $X_i$. They, as well as the additional observed attributes, $Z_i$ allow for in later sections, remain fixed in repeated sampling thought experiments, whereas $R_i$ and $X_i$ are stochastic and, as a result, so are the realized outcomes in the sample, $Y_i$. In the current section, we abstract from the presence of fixed observed attributes, which will play an important role in Section (ref). Let ${\boldsymbol{Y}}$, ${\boldsymbol{Y}}^{*}(1)$, ${\boldsymbol{Y}}^{*}(0)$, ${ \boldsymbol{R}}$, and ${\boldsymbol{X}}$ be the population $n$-vectors with $i$-th element equal to $Y_i$, $Y_i^{*}(1)$, $Y_i^{*}(0)$, $R_i$, and $X_i$ respectively. For sampled units (units with $R_i=1$) we observe $X_i$ and $Y_i$. For all units we observe $R_i$.

In general, estimands are functions of the full set of values $({ \boldsymbol{Y}}^{*}(1),{\boldsymbol{Y}}^{*}(0),{\boldsymbol{X}},{\boldsymbol{R}})$ for all units in the population, both those in the sample and those not in the sample. We consider two types of estimands, descriptive and causal. If an estimand can be written as a function of $({\boldsymbol{Y}},{\boldsymbol{X}})$, free of dependence on ${\boldsymbol{R}}$ and on the potential outcomes beyond the realized outcome, we label it a descriptive estimand. Intuitively a descriptive estimand is an estimand whose value would be known with certainty if we observe the realized values of all variables for all units in the population. If an estimand cannot be written as a function of $ ({\boldsymbol{Y}},{\boldsymbol{X}},{\boldsymbol{R}})$ because it depends on the potential outcomes ${\boldsymbol{Y}}^{*}(1)$ and ${\boldsymbol{Y}}^{*}(0)$, then we label it a causal estimand.\footnote{This does not partition the set of estimands. For example, there could be estimands that are functions of $({\boldsymbol{Y}},{\boldsymbol{X}},{\boldsymbol{R}})$, but not of $({\boldsymbol{Y}},{\boldsymbol{X}})$, although it is difficult to think of interesting ones that are.}

We now consider in our binary regressor example three closely related estimands, one descriptive and two causal:

and

We focus on the properties of a particular estimator:

This is also the least squares estimator of the coefficient on $X_i$ for the regression in the sample of $Y_i^{}$ on $X_i$ and a constant. There are two sources of randomness in this estimator: a sampling component arising from the randomness of $\boldsymbol{R}$ and a design component arising from the randomness of $\boldsymbol{X}$. We refer to the uncertainty generated by the randomness in the sampling component as sampling-based uncertainty, and the uncertainty generated by the design component as design-based uncertainty.

Next we consider the first two moments of $\widehat\theta$ under the following two assumptions:

We start by studying the first moment of the estimator, conditional on $ (N_1,N_0)$, and only for the cases where $N_1\geq 1$ and $N_0\geq 1$. We leave this latter conditioning implicit in the notation throughout this section. We also condition implicitly on the fixed potential outcomes ${\boldsymbol{Y}}^{*}(1)$ and ${\boldsymbol{Y}}^{*}(0)$. Taking the expectation only over the random sampling, or taking the expectation only over the random assignment, or over both, we find:

Next, we look at the variance of the estimator, maintaining both the random assignment and random sampling assumption. Define the population variances

and

We consider the variance of $ \widehat\theta$, as well as two conditional versions of this variance. We define the “sampling variance” conditional on ${\boldsymbol{X}}$, so that only the sampling uncertainty is taken into account. Analogously, we define the “design variance” conditional on ${\boldsymbol{R}}$, so that only the design uncertainty is taken into account. To make the different variances interpretable, we look at the expected value of the variances, taking the expectation both over the assignment and the sampling.

In this single binary regressor example the EHW variance estimator can be written as

and $\widehat S_0^2$ is defined analogously. Adjusting the degrees of freedom, using the modification proposed in mackinnon1985some specialized to this binary regressor example, we obtain $\widetilde V^\mathrm{ ehw}=\widehat S_1^2/N_1+\widehat S_0^2/N_0,$ which is identical to the variance estimator proposed by neyman1923, with the expectation of this modified EHW variance estimator $\widetilde V^\mathrm{ehw}$ (conditional on $N_0$ and $ N_1$) equal to the sampling variance in the infinite population case, $V^ \mathrm{sampling}(N_1,N_0,\infty,\infty)$. In the infinite population case the design-based uncertainty does not matter, so the EHW variance can be intepreted as implicitly taking into account design-uncertainty by focusing on the infinite population case.

We could also estimate the variance using resampling methods, which would give us variance estimates close to $\widehat{V}^{\mathrm{ehw}}$. To be precise, suppose we use the bootstrap where we draw $N_1$ bootstrap observations from the $N_1$ treated units and $N_0$ bootstrap units from the $N_0$ control units. In that case the bootstrap variance would in expectation (over the bootstrap replications) be equal to $\widehat{V}^{\mathrm{ehw}}$.

The appendix contains a Bayesian version of the analysis ofthe example from this section. Similar to the results in this section we show that, when the estimand of interest is defined for a finite population, the posterior variance depends not only on the sample sizes for treated and non-treated but also on the respective population sizes. Also similar to the analysis of this section, the posterior variance formula depends on whether the estimand is descriptive or causal.

The General Case

This section contains the main formal results in the article. We focus on a regression setting where we estimate a linear regression function for a scalar outcome and a number of regressors. The setting we consider here allows for the presence of two types of regressors. First, regressors that are causal, in the sense that they generate potential outcomes. Second, regressors that are attributes, in the sense that they are kept fixed for each unit in the thought experiment that provides the basis for inference. Which regressors are viewed as causal and which are viewed as attributes depends on the interpretation we wish to give to the regression estimates. If we wish to give a coefficient a causal interpretation, the corresponding regressor must be a cause. If a regressor is an attribute, the corresponding coefficient is simply estimating a population difference between subpopulations of units. For example, if we regress earnings on years of education, years of education may be the causal variable of interest. On the other hand, if we regress earnings on an indicator for participation in a job search program, age, and years of education, then the indicator for the program participation may be viewed as the causal variable of interest and age and years of education may be viewed as attributes. Because the repeated sampling thought experiment treats causes different from attributes, the variance of the regression estimator will depend on this designation.

Set Up

Consider a sequence of finite populations indexed by population size, $n$. Unit $i$ in population $n$ is characterized by a set of fixed attributes $ Z_{n,i}$ (including an intercept) and by a potential outcome function, $ Y_{n,i}^{*}(\cdot)$, which maps causes, $U_{n,i}$, into outcomes, with the realized outcome denoted by $ Y_{n,i}^{}=Y_{n,i}^{*}(U_{n,i})$. $Z_{n,i}$ and $U_{n,i}$ are real-valued column vectors, and $Y_{n,i}^{}$ is scalar. We do not place restrictions on the types of the variables: they can be continuous, discrete, or mixed. The attributes and potential outcome functions remain fixed in our repeated sampling thought experiments. The realized outcomes in the sample vary from sample to sample because the units in the sample and the values of the causal variables change.

There is a sequence of samples associated with the sequence of populations. We will use $R_{n,i}=1$ to indicate that unit $i$ of population $n$ is sampled, and $R_{n,i}=0$ to indicate that it is not sampled. For each unit in sample $ n$, we observe the triple, $(Y_{n,i},U_{n,i},Z_{n,i})$. Relative to Section (ref) we now allow for more complicated assignment mechanisms. In particular, we relax the assumption that the causes have identical distributions.

This assumption assumes independence of the treatment assignments. This is somewhat in contrast to the example in Section (ref), where we fixed the marginal distribution of the regressor, allowing us to obtain exact finite sample results. We do not need this here because we are focused on asymptotic results. We can allow for some dependence in the assignment mechanism, {\it e.g.,} clustering of the type analyzed in abadie2017should.

For what follows, it is convenient to work with a transformation $ X_{n,1},\ldots ,X_{n,n}$ of $U_{n,1},\ldots ,U_{n,n}$ that removes the correlation with the attributes. This can be accomplished in the following way. We assume that the population matrix $\sum_{i=1}^n Z_{n,i}Z_{n,i}^{\prime }$ is full-rank. Then, define

Later we formally make an assumption that will guarantee that this transformation is well-defined for large $n$. It is important to notice that, because $\Lambda_n Z_{n,i}$ is deterministic in our setting and $U_{n,1}, \ldots, U_{n,n}$ are i.n.i.d., the variables $X_{n,1}, \ldots, X_{n,n}$ are i.n.i.d.\ too.

For population $n$, let $\boldsymbol{Y}_n$, $\boldsymbol{X}_n$, $\boldsymbol{ Z}_n$, $\boldsymbol{R}_n$, and $\boldsymbol{Y}^*_n(\cdot)$ be matrices that collect outcomes, causes, attributes, sampling indicators, and potential outcome functions. We analyze the properties of the estimator $\widehat\theta_n$ obtained by minimizing least square errors in the sample:

The properties of the population regression residuals, $ e_{n,i}=Y_{n,i}-X_{n,i}^{\prime }\theta_n-Z_{n,i}^{\prime }\gamma_n$, depend on the exact nature of the estimands, $(\theta_n,\gamma_n)$. In what follows, we will consider alternative target parameters, which in turn will imply different properties for $e_{n,i}$. Notice also that, although the transformation in ((ref)) is typically unfeasible (because the values of $E[U_{n,i}]$ may not be known), $\widehat\theta_n$ is not affected by the transformation in the sense that the least squares estimators $(\widetilde\theta_n,\widetilde\gamma_n)$, defined as

satisfy $\widehat\theta_n=\widetilde\theta_n$ (although, in general, $ \widehat\gamma_n\neq \widetilde\gamma_n$). As a result, we can analyze the properties of $\widehat\theta_n$, focusing on the properties of the regression on $X_{n,1},\ldots, X_{n,n}$ instead of on $U_{n,1},\ldots, U_{n,n}$.

We assume random sampling, with some conditions on the sampling rate to ensure that the sample size increases with the population size.

The first part of Assumption (ref)$(ii)$ guarantees that as the population size increases, the (expected) sample size also increases. The second part of Assumption (ref)$(ii)$ allows for the possibility that, as $n$ increases, the sample size becomes a negligible fraction of the population size so that the EHW results, corresponding to $\rho=0$, are included as a special case of our results.

The next assumption is a regularity condition bounding moments.

Let

So $\Omega_n=E[W_n]$, where the expectation is taken over the distribution of $\boldsymbol{X}_n$. We also consider sample counterparts of $W_n$ and $\Omega_n$:

where $\widetilde\Omega_n =E[\widetilde W_n|{\boldsymbol{R}}_n]$. We will use superscripts to indicate submatrices. For example,

with analogous partitions for $\Omega_n$, $\widetilde W_n$, and $ \widetilde\Omega_n$. Notice that the transformation in ((ref)) implies that $\Omega_n^{XZ}$ and $\Omega_n^{ZX}$ are matrices will all zero entries.

We first obtain convergence results for the sample objects, $\widetilde W_n$ and $\widetilde \Omega_n$.

See appendix for proofs.

The next assumption imposes (deterministic) convergence of the expected value of the second moments in the population.

Descriptive and Causal Estimands

We now define the descriptive and causal estimands that generalize $\theta ^{ \mathrm{descr}}$, $\theta ^{\mathrm{causal},\mathrm{sample}}$, and $\theta^{ \mathrm{causal}}$ from Section (ref) to a regression context.

Causal estimands depend on the values of potential outcomes beyond the values that can be inferred from the realized outcomes. Given a sample, the only reason we may not be able to infer the exact value of a descriptive estimand is that we do not see all the units in the population. In contrast, even if we observe all units in a population, we are unable to infer the value of a causal estimand because its value depends on potential outcomes.

We define three estimands of interest, which under the conditions above exist with probability approaching one:

and

Alternatively, the estimands in ((ref)) to ((ref)) can be defined as the coefficients that correspond to the orthogonality conditions in terms of the residuals $ e_{n,i}=Y_{n,i}-X_{n,i}^{\prime }\theta_n-Z^{\prime }_{n,i}\gamma_n$,

for the descriptive, causal-sample, and causal estimands respectively. We will study the properties of the least squares estimator, $\widehat\theta_n$ , defined by

as an estimator of the parameters defined in equations ((ref)) to ((ref)).

Notice that $\theta^{\mathrm{causal}, \mathrm{sample}}_n$ and $\theta^{\mathrm{causal}}_n$ are causal estimands, while $\theta^{\mathrm{descr} }_n$ is not. However, the fact that an estimand is causal according to our definition does not imply it has an interpretation as an average causal effect. In Section (ref) we present conditions under which the regression estimand does have such an interpretation.

Causal Interpretations of the Estimands

By construction, the descriptive estimand can be interpreted as the set of coefficients of a population best linear predictor (least squares), e.g., goldberger1991course (goldberger1991course). A more challenging question concerns the interpretation of the two causal estimands, and in particular their relation to the potential outcome functions. In this section we investigate this question.

The next assumption generalizes random assignment to allow for some dependence of the assignment $U_{n,i}$ on the attributes $Z_{n,i}$.

Assumption (ref) looks very different from conventional exogeneity or unconfoundedness conditions, where the residuals are assumed to be (mean-) independent of the regressors, and so it merits some discussion. First, note that if the treatment $U_{n,i}$ is randomly assigned $E[U_{n,i}]$ is constant and Assumption (ref)$(i)$ and $(ii)$ are automatically satisfied as long as $Z_{n,i}$ includes an intercept.

Formally, Assumption (ref) relaxes the completely randomized assignment setting by allowing the distribution of $ U_{n,i}$ to depend on the attributes. However, this dependence is restricted in that the mean of $U_{n,i}$ is linear in $Z_{n,i}$. For example, Assumption (ref) holds automatically when $U_{n,1},\ldots, U_{n,n}$ are identically distributed and $Z_{n,i}$ contains a saturated set of indicators for all possible values of the attributes.

In the special case where the treatment is binary and $E[U_{n,i}]$ is the propensity score, the assumption amounts to combination of an unconfoundedness assumption that the treatment assignment does not depend on the potential outcomes and a linear model for the propensity score.

Later in this section, we will show that under a set of conditions that includes Assumption (ref), the two estimands $\theta_n^\mathrm{causal}$ and $ \theta^{\mathrm{causal},\mathrm{sample}}$ can be interpreted as weighted averages of unit-level causal effects. The connection between linearity in the propensity score (rosenbaum1983central), in our analysis represented by $E[U_{n,i}]=B_n Z_{n,i}$, and the interpretation of population regression coefficients as weighted averages of heterogeneous causal effects has been previously noticed in related contexts (see angrist1998, angrist1998; angristpischke , angristpischke; aronow2016, aronow2016; tymon2017, tymon2017).

In this formulation, any dependence of the potential outcomes $ Y_{n,i}^{*}(u) $ on observed or unobserved attributes is subsumed by $ \theta_{n,i}$ and $\xi_{n,i}$, which are non-stochastic. Each element of the vector $\theta_{n,i}$ represents the causal effect of increasing the corresponding value of $U_{n,i}$ by one unit.

The linearity in Assumption (ref) is a strong restriction in many settings. However, in some leading cases -- in particular, when the causal variable is binary or, more generally when the causal variable takes on only a finite number of values -- one can ensure that this assumption holds by including in $U_{n,i}$ indicator variables representing each but one of the possible values of the cause. With Assumption (ref) we are able to provide a more transparent interpretation of the regression estimator.

The Asymptotic Distribution of The Least Squares Estimator

In this section we present the main result of the article, describing the properties of the least squares estimator viewed as an estimator of the causal estimands and, separately, viewed as an estimator of the descriptive estimand. In contrast to Section (ref), we do not have exact results, relying instead on asymptotic results based on sequences of populations.

First, we define the population residuals, denoted by $\varepsilon _{n,i}$, relative to the population causal estimands,

Under the assumption that the $X_{n,i}$ are jointly independent (but not necessarily identically distributed), the $n$ products $X_{n,i} \varepsilon _{n,i}$ are also jointly independent but not identically distributed. Most importantly, in general the expectations $E[X_{n,i} \varepsilon _{n,i}] $ may vary across $i$, and need not all be zero. However, as shown in Section (ref), the averages of these expectations over the entire population are guaranteed to be zero by the definition of $(\theta^\mathrm{causal}_n,\gamma^\mathrm{causal}_n)$. Define the limits of the population variance,

and the expected outer product

The difference between $\Delta^\mathrm{ehw}$ and $\Delta^\mathrm{cond}$ is the limit of the average outer product of the means,

which is positive semidefinite. We assume existence of these limits.

The Variance Under Correct Specification

Consider a constant treatment effect assumption, which is required for a correct specification of a linear regression function as a function that describes potential outcomes.

This strengthens Assumption (ref) by requiring that the $ \theta_{n,i}$ do not vary by $i$.

Under Assumption (ref), Theorem (ref) implies that $\theta_n^\mathrm{causal}=\theta_n$ (although it need not be the case that $\theta^\mathrm{descr}=\theta_n$). Then, for

we obtain that equation ((ref)) holds for $\gamma_n^ \mathrm{causal}=\Lambda_n^{\prime }\theta_n+\lambda_n$ and $ \varepsilon_{n,i}=\xi_{n,i}-Z_{n,i}^{\prime }\lambda_n$. In this case, the residuals, $\varepsilon_{n,i}$, are non-stochastic. As a result, $ E[X_{n,i}\varepsilon_{n,i}]=E[X_{n,i}]\varepsilon_{n,i}=0$, which implies $ \Delta^\mu=\Delta^\mathrm{ehw}-\Delta^\mathrm{cond}=0$. This leads to the following result.

Notice that the result of the theorem applies also with $\theta_n^{\mathrm{ causal},\mathrm{sample}}$ replacing $\theta_n^{\mathrm{causal}}$ because the two parameter vectors are identical (with probability approaching one) under Assumption (ref).

Estimating the Variance

Now let us turn to the problem of estimating the variance for the descriptive and causal estimands. In what follows, we will use the shorthands $V^{{\rm causal}}=\Gamma ^{-1}\left(\rho\Delta^{\mathrm{cond} }+(1-\rho) \Delta^\mathrm{ehw} \right) \Gamma ^{-1}$, $V^{{\rm causal}, {\rm sample}}=\Gamma ^{-1}\Delta^{\mathrm{cond}} \Gamma ^{-1}$, $V^{\mathrm{descr}}=(1-\rho )\Gamma ^{-1}\Delta^\mathrm{ehw} \Gamma ^{-1}$, and $V^{{\rm ehw}}=\Gamma ^{-1}\Delta^\mathrm{ehw} \Gamma ^{-1}$. There are four components to the asymptotic variances, $\rho$, $\Gamma $, $\Delta^\mathrm{ehw}$ and $\Delta^\mathrm{ cond}$. The first three are straightforward to estimate. $\rho$ can be estimated as $\widehat\rho_n=N/n$, as long as the population size is known. To estimate $\Gamma $, first estimate $\Lambda_n$ as

Then one can estimate $\Gamma $ as the average of the matrix of outer products over the sample:

It is also straightforward to estimate $\Delta^\mathrm{ehw}$. First we estimate the residuals for the units in the sample, $\widehat{\varepsilon} _{n,i}=Y_{n,i}-(U_{n,i}-\widehat{\Lambda}_n Z_{n,i})^{\prime }\widehat\theta_n -Z_{n,i}^{\prime }\widehat\gamma_n,$ and then we estimate $ \Delta^\mathrm{ehw}$ as:

The EHW large sample variance, $V^\mathrm{ehw}$, is then estimated as

Let $\widehat V^{{\rm descr}}_n=(1-\widehat\rho_n) \widehat V^\mathrm{ehw}_n$. The result of Lemma (ref) immediately implies $\widehat V^{{\rm descr}}_n\overset{p}{\rightarrow} V^{{\rm descr}}$.

It is more challenging to estimate $V^{{\rm causal}}$ and $V^{{\rm causal},{\rm sample}}$ because they involve $\Delta^\mathrm{cond}$. Estimating $\Delta^{{\rm cond}}$ is complicated because of the same reason that complicates the estimation of the variance of the average treatment effect estimator in Section (ref). In that case there are three terms in the expression for the variance in equation ((ref)). The first two are straightforward to estimate, but the third one, $S^2_\theta/n$ cannot be estimated consistently because we do not observe both potential outcomes for the same units. Often, researchers use the conservative estimator based on ignoring $S^2_\theta/n$. If we proceed in the same fashion for the regression context of Section (ref), we obtain the conservative estimator $\widehat V^\mathrm{ehw}$, based on ignoring $\Delta^\mu$. We show, however, that in the presence of attributes we can improve the variance estimator. We build on abadie2008estimation, abadie2014inference, and fogarty2016regression who, in contexts different than the one studied in this article, have used the explanatory power of attributes to improve variance estimators. abadie2008estimation and abadie2014inference do so using nearest-neighbor techniques. Here we follow fogarty2016regression and apply linear regression techniques. The proposed estimator replaces the expectations $E[X_{n,i}\varepsilon_{n,i}]$, which cannot be consistently estimated, with predictors from a linear least squares projection of estimates of $X_{n,i}\varepsilon_{n,i}$ on the attributes, $Z_{n,i}$. Let $ \widehat X_{n,i}= U_{n,i}-\widehat\Lambda_n Z_{n,i}$, and

The matrix $\widehat G_n$ contains the coefficients of a least squares regression of $\widehat X_{n,i}\widehat\varepsilon_{n,i}$ on $Z_{n,i}$. The next assumption ensures convergence of $\widehat G_n$.

Consider now the following estimator,

which uses $\widehat G_nZ_{n,i}$ in lieu of a consistent estimator of $ E[X_{n,i}\varepsilon_{n,i}]$. Notice that we do not assume that $ E[X_{n,i}\varepsilon_{n,i}]$ is linear in $Z_{n,i}$. However, we will show that, as long as the attributes can linearly explain some of the variance in $\widehat X_{n,i}\widehat\varepsilon_{n,i}$, the estimator $\widehat \Delta^Z_n$ is smaller (in a matrix sense) than $\widehat\Delta^\mathrm{ehw} _n$. These results are provided in the following lemma.

Estimators of $V^{{\rm causal},{\rm sample}}$ and $V^{{\rm causal}}$ follow immediately from Lemma (ref) by replacing $\Delta^\mathrm{cond}$ with the estimate $\widehat\Delta^Z_n$ in the asymptotic variance formulas of Theorem (ref), leading to $\widehat{V}^{{\rm causal},{\rm sample}}_n=\widehat \Gamma ^{-1}_n\widehat\Delta^Z_n\widehat \Gamma ^{-1}_n$ for the estimation of $V^{{\rm causal}, {\rm sample}}$ and $\widehat V^{{\rm causal}}_n=\widehat\rho_n \widehat{V}^{{\rm causal},{\rm sample}}_n + (1-\widehat\rho_n)\widehat V^{{\rm ehw}}_n$ for the estimation of $V^{{\rm causal}}$. These estimators are not larger (and typically smaller) than $\widehat V^{{\rm ehw}}_n$ and they remain conservative in large samples.

Simulations

In this section, we use a simple data-generating process as well as simulations to illustrate the difference between the conventional EHW variance estimator and the variance estimators proposed in this article. We focus on the case of a single causal variable, $X_{n,i}$. In addition to the causal variable, the simulations employ an outcome variable, $Y_{n,i}$, and a vector of attributes, $Z_{n,i}$, which consists of a constant equal to one and $k$ values drawn independently from the standard normal distribution. The potential outcome function has the form in equation ((ref)). Population values of $\theta_{n,i}$ are generated as independent draws from a normal distribution with mean $Z_{n,i}'\psi$, where $\psi=(0,\psi_1,\ldots,\psi_k)'$ is a $k+1$ vector, and variance $\sigma^2_\theta$. Population values of $\xi_{n,i}$ and $U_{n,i}$ are generated as independent draws from a normal distribution with mean zero and variance one. Because in this data-generating process $E[U_{n,i}]=0$, it follows that $\Lambda_n$ is a row-vector of zeros and $X_{n,i}=U_{n,i}$. We use this data-generating process to produce a population of size $n$. For this data-generating process it can be shown that $\Gamma =1$, $\Delta^{{\rm ehw}}=1+3(\psi'\psi+\sigma^2_\theta)$ and $\Delta^{{\rm cond}}=1+2(\psi'\psi+\sigma^2_\theta)$, and $\Delta^Z=1+2\psi'\psi+3\sigma^2_{\theta}$ with probability one. In each simulation repetition, we sample units at random with probability $\rho$ from the population. As a result, the sample size $N$ is random with $E[N]=n\rho$. For each sample we estimate $\widehat\theta_n$ by least squares (as in equation ((ref))) and a number of variance estimators.

In Table (ref) we report the results of the simulations. We consider seven designs. The first column reports the basic design, with $\rho=0.01$ and $n=100000$, so the average sample size is 1000. In this design, there is one stochastic regressor, so $k=1$, and the distribution of the treatment effect, $\theta_{n,i}$, is given by parameter values $\psi=(0,2)'$ and $\sigma^2_\theta=1$. The remaining designs in the second to seventh columns are variations of the basic design in the first column. In the second design, we increase the dimensionality of $Z_{n,i}$ used for estimation from two to ten. Still, in this design $\psi$ has all entries equal to zero except for $\psi_1=2$, so only the first stochastic regressor matters for the distribution of $\theta_{n,i}$. In the next design, we change the population size to 10000, so that the average sample size is 100. In the fourth design, we change the population size to 1000 and the sampling rate, $\rho$, to one. In the fifth design, we impose $\psi'\psi=0$, which makes the treatment effect unrelated to the regressors, $Z_{n,i}$. In the sixth design, we set $\sigma^2_\theta=0$, which removes the stochastic part of the treatment effect. In the last design, $\psi'\psi=0$ and $\sigma^2_\theta=0$, so the treatment effect is constant. The first panel of Table (ref) provides the parameters of each of the seven simulation designs.

The second panel of Table (ref) reports the standard deviations of $(\widehat\theta_n-\theta_n^\mathrm{descr})$, $(\widehat\theta_n-\theta_n^{\mathrm{causal},\mathrm{sample}})$ and $(\widehat\theta_n-\theta_n^\mathrm{causal})$ across simulation iterations. The remaining panels report feasible standard errors based on the estimators of Section (ref) as well as bootstrap standard errors, along with coverage rates of the corresponding 95 percent confidence intervals. We employ 50000 iterations for the simulations and 1000 bootstrap samples. The coverage rates in each of the panels of the table are based on the adding and subtracting 1.96 times the standard errors in the first row of that panel.

For the basic design in the first column of Table (ref), $\rho=0.01$ is very small, and EHW and bootstrap standard errors provide accurate estimates of the standard deviations of $(\widehat\theta_n-\theta_n^{{\rm descr}})$ and $(\widehat\theta_n-\theta_n^{{\rm causal}})$. However, the standard deviation of $(\widehat\theta_n-\theta_n^{{\rm causal},{\rm sample}})$ is substantially smaller than that of $(\widehat\theta_n-\theta_n^{\rm descr})$ and $(\widehat\theta_n-\theta_n^{\rm causal})$, and the EHW and bootstrap variance estimators are very conservative for the sample average causal effect, $\theta^{{\rm causal},{\rm sample}}_n$. The variance estimator based on $\widehat V^{{\rm causal},{\rm sample}}_n$ is substantially smaller, and still has more than correct coverage for $\theta^{{\rm causal},{\rm sample}}$. Increasing the number of regressors in the second design leaves the result virtually unaffected. The same patterns of results appear in the third column, albeit with less precise variance estimators due to much smaller sample sizes. In the fourth design, $\rho=1$ and, as predicted by the results in section (ref), EHW standard errors greatly overestimate the variability of $(\widehat\theta_n-\theta_n^\mathrm{descr})$. The same is true for bootstrap standard errors. In the fifth, design we go back to the small sampling rate, $\rho=0.01$ and this time $\psi'\psi=0$, so regressors do not explain variation in treatment effects, and $\Delta^Z=\Delta^{{\rm ehw}}$. As suggested by the results in sections (ref) and (ref), all variance estimators produce similar results in this design. In the sixth design, where regressors explain all the variation in treatment effects, $\Delta^{Z}=\Delta^{{\rm cond}}$ and standard errors based on $\widehat V_n^{{\rm causal},{\rm sample}}$ and $\widehat V_n^{{\rm causal}}$ closely approximate the standard deviations of $(\widehat\theta_n-\theta_n^{{\rm causal},{\rm sample}})$ and $(\widehat\theta_n-\theta_n^{{\rm causal}})$, respectively. In the final design with a constant treatment effect, all the variances are similar.

Conclusion

In this article we study the interpretation of standard errors in regression analysis when the assumption that the sample is drawn randomly from a much larger population of interest is not appropriate. We base our results on a potential outcome framework, where the estimands of interest may be descriptive or causal, and we provide a coherent interpretation for standard errors that allows for uncertainty coming from both random sampling and from conditional random assignment. The standard errors estimators proposed in this article may be different from the conventional ones, and they may vary depending on {\em (i)} the specific nature of the estimand of interest (i.e., descriptive or causal), {\em (ii)} the fraction of the population represented in the sample, and {\em (iii)} the extent to which measured attributes explain variation in treatment effects.

In the current article we focus exclusively on linear regression models. The concerns we raise in this article arise in many other settings and for other kinds of hypotheses, and the implications would need to be worked out for those settings. Thus, we see this article as a first step in broader research program.

\parindent 0cm\setcounter{equation}{0} \parskip 5pt {A.\arabic{equation}} \setcounter{lemma}{0}{A.\arabic{lemma}} \setcounter{theorem}{0}{A.\arabic{theorem}} \centerline{\sc Appendix}

\vskip0.5cm

\centerline{\sc I. A Bayesian Approach}

\vskip0.5cm

Given that we are advocating for a different conceptual approach to modeling inference, it is useful to look at the problem from more than one perspective. In this section we consider a Bayesian perspective and re-analyze the example from Section (ref). Viewing the problem from a Bayesian perspective reinforces the point that formally modeling the population and the sampling process leads to the conclusion that inference is different for descriptive and causal questions. Note that in this discussion the notation will necessarily be slightly different from the rest of the article; notation and assumptions introduced in this subsection apply only within this subsection.

Define $\boldsymbol{Y}^*_n(1)$, $\boldsymbol{Y}^*_n(0)$ to be the $n$ vectors with typical elements $Y^*_i(1)$ and $Y^*_i(0)$, respectively. We view the $n$ -vectors $\boldsymbol{Y}^*_n(1)$, $\boldsymbol{Y}^*_n(0)$, $\boldsymbol{R}_n$ , and $\boldsymbol{X}_n$ as random variables, some observed and some unobserved. We assume the rows of the $n\times 4$ matrix $[\boldsymbol{Y}^* _n(1),\boldsymbol{Y}^*_n(0), \boldsymbol{R}_n, \boldsymbol{X}_n]$ are exchangeable. Then, by appealing to DeFinetti's theorem, we model this, with no essential loss of generality (for large $n)$ as the product of $n$ independent and identically distributed random quadruples $ (Y^*_i(1),Y^*_i(0),R_i,X_i)$ given some unknown parameter $\beta$:

Inference then proceeds by specifying a prior distribution for $\beta $, say $p(\beta )$. To make this specific, consider the following model. Let $X_i$ and $ R_i$ have Binomial distributions with parameters $q$ and $\rho$,

The pairs $(Y^*_i(1),Y^*_i(0))$ are assumed to be jointly normally distributed:

so that the full parameter vector is $\beta=(q,\rho,\mu_1,\mu_0 ,\sigma^2_1,\sigma^2_0,\kappa)$.

We change the observational scheme slightly from Section (ref) to allow for the analytic derivation of posterior distributions. We assume that for all units in the population we observe the pair $(R_i,X_i)$, and for units with $R_i=1$ we observe the outcome $Y_i=Y^*_i(X_i)$. Define $ \widetilde Y_i=R_i Y_i$, so for all units in the population we observe the triple $(R_i,X_i,\widetilde Y_i)$. Let $\boldsymbol{R}_n$, $\boldsymbol{X}_n$ , and $\widetilde{\boldsymbol{Y}}_n$ be the $n$ vectors of these variables. $ \bar Y_1 $ denotes the average of $Y_i$ in the subpopulation with $R_i=1$ and $X_i=1$, and $\bar Y_0$ denotes the average of $Y_i$ in the subpopulation with $R_i=1 $ and $X_i=0$.

The descriptive estimand is

The causal estimand is

It is interesting to compare these estimands to an additional estimand, the super-population average treatment effect,

In general these three estimands are distinct, with their own posterior distributions, but in some cases, notably when $n$ is large, the three posterior distributions are similar.

It is instructive to consider a very simple case where analytic solutions for the posterior distribution for $\theta^\mathrm{descr}_n$, $\theta^ \mathrm{causal}_n$, and $\theta^\mathrm{causal}$ are available. Suppose $ \sigma^2_1 $, $\sigma^2_0$, $\kappa$ and $q$ are known, so that the only unknown parameters are the two means $\mu_1$ and $\mu_0$. Finally, let us use independent, diffuse (improper), prior distributions for $\mu_1$ and $ \mu_0$.

Then, a standard result is that the posterior distribution for $(\mu_1 ,\mu_0)$ given $({{\boldsymbol{R}}}_n,{\boldsymbol{X}}_n, \widetilde{{\boldsymbol{Y}}} _n) $ is

where $N_1$ is the number of units with $R_i=1$ and $X_i=1$, and $N_0$ is the number of units with $R_i=1$ and $X_i=0$. This directly leads to the posterior distribution for $\theta^\mathrm{causal}$:

A longer calculation leads to the posterior distribution for the descriptive estimand:

The implied posterior interval for $\theta^\mathrm{descr}_n $ is very similar to the corresponding confidence interval based on the normal approximation to the sampling distribution for $\bar Y_1 -\bar Y_0$. If $n_1$ and $n_0$ are large, this posterior distribution is close to the posterior distribution of the causal estimand. If, on the other hand, $N_1=n_1$ and $ N_0=n_0$, then the posterior distribution of the descriptive estimand becomes degenerate and centered at $\bar Y_1 -\bar Y_0$.

A somewhat longer calculation for $\theta^\mathrm{causal}_n$ leads to

Consider the special case of constant treatment effects, where $ Y_i(1)-Y_i(0)=\mu_1-\mu_0$. Then, $\kappa=1$, and $\sigma_1=\sigma_0$, and the posterior distribution of $\theta^\mathrm{causal}_n$ is the same as the posterior distribution of $\theta^\mathrm{causal}$. The same posterior distribution arises in the limit if $n$ goes to infinity, regardless of the values of $\kappa$, $\sigma_1$, and $\sigma_0$.

To summarize, if the population is large, relative to the sample, the posterior distributions of $\theta^\mathrm{descr}_n$, $\theta^\mathrm{causal}_n$ and $ \theta^\mathrm{causal}$ agree. However, if the population is small, the three posterior distributions differ, and the researcher needs to be precise in defining the estimand. In such cases, simply focusing on the super-population estimand $\theta^\mathrm{causal} =\mu_1 -\mu_0$ is arguably not appropriate, and the posterior inferences for such estimands will differ from those for other estimands such as $\theta^\mathrm{causal}_n$ or $\theta^ \mathrm{descr}_n$. \vskip0.5cm

\centerline{\sc II. Proofs}

Proof of Lemma (ref): See supplementary appendix.$\square$

Proof of Theorem (ref): For $n$ large enough $ \sum_{i=1}^n Z_{n,i}Z_{n,i}^{\prime }$ is full rank and $\Lambda_n$ exists, so $\Omega^{ZX}_n=0$. This implies

Moreover, for $n$ large enough, $\Lambda_n=B_n$, which implies $E[X_{n,i}]=0$ , $\widetilde \Omega^{XZ}_n=0$, and

with probability approaching one. Now,

implies the results.$\square$

{\bf Proof of Theorem (ref)}: Let $\nabla Y^*_{n,i}()$ be the gradient of $Y^*_{n,i}()$. By the mean value theorem there exist sets $\mathcal T_{n,i}\subseteq [0,1]$ such that for any $t_{n,i}\in\mathcal T_{n,i}$, we have $Y^*_{n,i}(U_{n,i})=Y^*_{n,i}(B_nZ_{n,i})+X_{n,i}'\nabla Y^*_{n,i}(B_nZ_{n,i}+t_{n,i}X_{n,i})$. We define $\varphi_{n,i}=\nabla Y^*_{n,i}(v_{n,i})$, where $v_{n,i}=B_nZ_{n,i}+\bar t_{n,i}X_{n,i}$ and $\bar t_{n,i}=\sup\,\mathcal T_{n,i}$. Now, $E[X_{n,i}Y_{n,i}]=E[X_{n,i}]Y^*_{n,i}(B_nZ_{n,i})+E[X_{n,i}'\varphi_{n,i}]=E[X_{n,i}'\varphi_{n,i}]$. The rest of the proof is as for Theorem (ref). $\square$

The following lemma will be useful for establishing asymptotic normality.

Proof: Notice that

and

Now the continuous mapping theorem implies

As a result, it is enough to prove

Let

Consider $n$ large enough so $s_n^2>0$. Notice that, for $i=1,\ldots, n$,

and

Therefore,

Using $\rho_n\leq \rho_n^{1/(2+\delta)}$, $|\mu_{n,i}|^{2+\delta}\leq E[|V_{n,i}|^{2+\delta}]$, and Minkowski's inequality, we obtain:

Applying Liapunov's theorem (see, e.g., davidson1994stochastic, davidson1994stochastic), we obtain

Now, the result of the lemma follows from equation ((ref)) and from $s_n/\sqrt{\sigma^2+(1-\rho)\kappa^2}\rightarrow 1$.$\square$

\noindentProof of Lemma (ref): To prove $(i)$, consider $ V_{n,i}=a^{\prime }X_{n,i}\varepsilon_{n,i}$ for $a\in \mathbb{R}^k$. We will verify the conditions Lemma (ref). Notice that,

By Minkowski's inequality and Assumption (ref), the right-hand side of last equation is bounded. In addition,

Let $a\neq 0$. Then,

This implies

Using the Cramer-Wold device, this implies

The proofs of $(ii)$ and $(iii)$ are similar.$\square$

Proof of Theorem (ref) : To prove $(i)$, notice that

is invertible with probability approaching one. Then,

Therefore,

where

Because (i) $\Omega _{n}^{XZ}=0$, (ii) the first term of $r_{n}$ is $ o_{p}(1) $, and (iii) $(1/\sqrt{N})\sum_{i=1}^{n}R_{n,i}X_{n,i}\varepsilon _{n,i}$ is $O_{p}(1)$ (under the conditions stated above), it follows that

if we can show

We can write this standardized sum as

As shown in Lemma A.1, $\left( n\rho _{n}/N\right) ^{1/2}\overset{p}{ \rightarrow }1$. Therefore, it suffices to show

This expression has zero mean because $R_{n,i}$ is independent of $ \varepsilon _{n,i}$ and

We can study each element of the vector separately. By Chebyshev's inequality it suffices to show that the variances are bounded. Consider the $ j^{th}$ element. Then, by independence across $i$,

where the last equality holds because $E[R_{n,i}] =\rho _{n}$ and $Z_{n,i,j}$ is nonrandom. Both $Z_{n,i,j}^{2}$ and $E[\varepsilon _{n,i}^{2}]$ are bounded by Assumption (ref), and so this completes the proof. The proofs of $(ii)$ and $(iii)$ are analogous.\newline \strut $\square $

\noindentProof of Theorem (ref): The result follows directly $E[X_{n,i}\varepsilon_{n,i}]=0$. $\square$

\noindentProof of Lemma (ref): First, notice that (with probability approaching one) $\Lambda_n$ exists and it is equal to $ B_n $. This implies,

which converges to zero in probability by Lemma (ref) and Assumption (ref). Direct calculations yield

Now, Lemma (ref) and Assumption (ref) imply $\widehat \Gamma _n\overset{p}{\rightarrow} \Gamma $, where $\Gamma $ is full rank. Theorem (ref) diretclty implies $ \widehat\theta_n-\theta^\mathrm{causal}_n\overset{p}{\rightarrow} 0$. $ \widehat\gamma_n-\gamma^\mathrm{causal}_n\overset{p}{\rightarrow} 0$ follows from Lemma (ref). Let

Let $\alpha$ be a multi-index of dimension equal to the length of $ T_{n,i}=(Y_{n,i} : X_{n,i}^{\prime }: Z_{n,i}^{\prime })$. In addition, let

Using the same argument as in the proof of Lemma (ref) and given that Assumption (ref) holds with $\delta=4$, it follows that $\widetilde T_n^\alpha-\Psi_n^\alpha\overset{p}{\rightarrow} 0$ for $|\alpha|\leq 4$. This result directly implies $\widetilde\Delta^\mathrm{ehw}_n-\Delta^\mathrm{ ehw}_n\overset{p}{\rightarrow} 0$. By the same argument plus convergence of $ \widehat\theta_n$ and $\widehat\gamma_n$, it follows that $\widehat\Delta^ \mathrm{ehw}_n-\breve\Delta^\mathrm{ehw}_n\overset{p}{\rightarrow} 0$ and $ \breve\Delta^\mathrm{ehw}_n-\widetilde\Delta^\mathrm{ehw}_n\overset{p}{ \rightarrow} 0$. Now, the result follows from $\widehat\Delta^\mathrm{ehw} _n-\Delta^\mathrm{ehw}=(\widehat\Delta^\mathrm{ehw}_n-\breve\Delta^\mathrm{ ehw}_n)+(\breve\Delta^\mathrm{ehw}_n-\widetilde\Delta^\mathrm{ehw} _n)+(\widetilde\Delta^\mathrm{ehw}_n-\Delta^\mathrm{ehw}_n)+(\Delta^\mathrm{ ehw}_n-\Delta^\mathrm{ehw})\overset{p}{\rightarrow} 0$, where the last difference goes to zero by Assumption (ref).$\square$

\noindentProof of Lemma (ref): Notice that,

so that $\widehat \Delta^Z_n$ is no larger than $\widehat \Delta^\mathrm{ehw} _n$ in a matrix sense.

Let

be the expected value of $\widehat G_n$. Under the assumptions of Lemma (ref) and using the same argument as in the proof of that lemma, we obtain $\widehat G_n-G_n\overset{p}{\rightarrow} 0$. Therefore, $ \widehat \Delta^\mathrm{proj}_n -\Delta^\mathrm{proj}_n\overset{p}{ \rightarrow} 0$, where

Moreover, $\widehat \Delta^Z_n -\Delta^Z_n\overset{p}{\rightarrow} 0$, where $\Delta^Z_n = \Delta^\mathrm{ehw}_n-\Delta^\mathrm{proj}_n$ and

Let

Notice that

Let $\boldsymbol{A}_n$ and $\boldsymbol{D}_n$ be the matrices with $i$-th rows equal to $E[\varepsilon_{n,i}X_{n,i}^{\prime }]/\sqrt n$ and $ Z_{n,i}^{\prime }/\sqrt n$, respectively. Let $\boldsymbol{I}_n$ be the identity matrix of size $n$. Then,

which is positive semi-definite. Because $\Delta^\mathrm{cond}_n=\Delta^ \mathrm{ehw}_n-\Delta^\mu_n$, we obtain,

where the inequalities are to be understood in a matrix sense. Now, it follow from Assumption (ref) that $G_n$ and, therefore, $ \Delta^\mathrm{proj}_n$ and $\Delta^Z_n$ have limits. Then,

where $\Delta^\mathrm{cond}$, $\Delta^Z$, and $\Delta^\mathrm{ehw}$ are the limits of $\Delta^\mathrm{cond}_n$, $\Delta^Z_n$, and $\Delta^\mathrm{ehw}_n$ , respectively.$\square$