How to derive the variance of the mean of predictions from a linear regression model?

Question

The context is linear regression analysis for estimating a sample mean. Assume the usual multivariate linear model: $$Y = X\beta + \epsilon$$ with $X$ a $n \times p$ covariate matrix with intercept, $y$ a $n \times 1$ vector, and $\epsilon$ and $n \times 1$ independent and identically distributed vector of error terms with mean zero and variance $\sigma^2$. We estimate $\beta$ by OLS as $\hat{\beta}=(X^TX)^{-1}X^Ty$ from data $\{(y_1,x^T_1),...,(y_n,x^T_n)\}$.

Now suppose we have a new sample of $X_0$ drawn from a population distribution $p(X_0)$. Assume the same model applies in this population, but we do not observe $Y$. We are interested in estimating $E(Y)$ and in particular the variance of its estimator $\hat{E}(Y)$.

Since $$E(Y)=E_X(E(Y|X))=E_X( X\beta )$$ it seems an unbiased Monte-Carlo integration type of estimator for $E(Y)$ is $$\hat{E}(Y)= n_0^{-1} \sum_{i=1}^n x^T_{0i} \hat{\beta} = n_0^{-1} e^T X_0 \hat{\beta},$$ with $e$ the unit vector. How can I derive the variance of $\hat{E}(Y)$? At first I thought it is simply $$V(\hat{E}(Y)) = n_0^{-2} e^T X_0 V(\hat{\beta}) X_0^T e.$$

with well-known $V(\hat{\beta})=\sigma^2 (X^TX)^{-1}$. However, by simulation I realized this is not true. I think the problem here is that both $\hat{\beta}$ and $X_0$ are random, but I am unsure how to derive an estimator taking this into account.

Answer 1

I assume that $S_n\equiv\{(x_i^{\top},\epsilon_i)\}_{i-1}^n$ and $S_{0,m}\equiv\{(x_{0,i}^{\top},\epsilon_{0,i})\}_{i=1}^m$are independent copies of $(x^{\top},\epsilon)$ with $\mathsf{E}[\epsilon\mid x]=0$ a.s. Then, indeed, a reasonable estimator of $\mathsf{E}y$ is $$ \hat{\mathsf{E}}y=\hat{\mathsf{E}}x^{\top}\times\hat{\beta}, $$ where $$ \hat{\mathsf{E}}x=\frac{1}{m}\sum_{i=1}^m x_{0,i} \quad\text{and}\quad \hat{\beta}=\left(\frac{1}{n}\sum_{i=1}^n x_ix_i^{\top}\right)^{-1}\frac{1}{n}\sum_{i=1}^n x_iy_i. $$

Since $S_n$ and $S_{0,m}$ are independent, $$ \mathsf{E}[\hat{\mathsf{E}}y]=\mathsf{E}[\hat{\mathsf{E}}x^{\top}]\times \mathsf{E}[\hat{\beta}]=\mathsf{E}x^{\top}\beta, $$ and $$ \mathsf{E}[\hat{\mathsf{E}}y]^2=\mathsf{E}\left[\hat{\mathsf{E}}x^{\top}\left(\mathsf{E}\hat{\beta}\hat{\beta}^{\top}\right)\hat{\mathsf{E}}x\right]=\beta^{\top}W_m\beta+\sigma^2\operatorname{tr}(W_mV_n), $$ where $V_n=\mathsf{E}\left[\sum_{i=1}^n x_ix_i^{\top}\right]^{-1}$ and $W_m\equiv\mathsf{E}[\hat{\mathsf{E}}x\hat{\mathsf{E}}x^{\top}]=\mathsf{E}x\mathsf{E}x^{\top}+m^{-1}\operatorname{Var}(x)$ (we assume that $V_n$ exists and is finite). Consequently, $$ \bbox[cornsilk,5px] { \begin{align} \operatorname{Var}(\hat{\mathsf{E}}y)&=\mathsf{E}[\hat{\mathsf{E}}y]^2-(\mathsf{E}[\hat{\mathsf{E}}y])^2\\ &=\sigma^2\mathsf{E}x^{\top}V_n\mathsf{E}x+\frac{1}{m}\left(\beta^{\top}\operatorname{Var}(x)\beta+\sigma^2\operatorname{tr}(\operatorname{Var}(x)V_n)\right). \end{align} } $$ It follows that $\operatorname{Var}(\hat{\mathsf{E}}y)\to\sigma^2\mathsf{E}x^{\top}V_n\mathsf{E}x$ as $m\to\infty$.

---

If you are interested in the conditional variance of $\hat{\mathsf{E}}y$ given $X_n\equiv\{x_i\}_{i=1}^n$, then using the same reasoning one obtains $$ \bbox[cornsilk,5px] { \operatorname{Var}(\hat{\mathsf{E}}y\mid X_n)=\sigma^2\mathsf{E}x^{\top}\tilde{V}_n\mathsf{E}x+\frac{1}{m}\left(\beta^{\top}\operatorname{Var}(x)\beta+\sigma^2\operatorname{tr}(\operatorname{Var}(x)\tilde{V}_n)\right), } $$ where $\tilde{V}_n=\left(\sum_{i=1}^n x_ix_i^{\top}\right)^{-1}$. Similarly to the previous case, $$ \operatorname{Var}(\hat{\mathsf{E}}y\mid X_n)\to \Sigma_n\equiv\sigma^2\mathsf{E}x^{\top}\tilde{V}_n\mathsf{E}x $$ as $m\to\infty$. Finally, by the SLLN (assuming that $\mathsf{E}\|x\|^2<\infty$ and $\mathsf{E}xx^{\top}$ is positive definite), $$ n\Sigma_n\to \sigma^2\mathsf{E}x^{\top}\left(\mathsf{E}xx^{\top}\right)^{-1}\mathsf{E}x\quad\text{a.s.} $$

---

You may be interested in deriving the asymptotic distribution of $\hat{\mathsf{E}}y$. For each $k\in\mathbb{N}$, we estimate $\hat{\mathsf{E}}y$ using two samples $\{(y_i,x_i^{\top})^{\top}\}_{i=1}^{n_k}$ and $\{x_{0,i}\}_{i=1}^{m_k}$ with $n_k,m_k\to\infty$ as $k\to\infty$. Let $\alpha_k$ be s.t. $$ \frac{\alpha_k}{\sqrt{m_k}}\to c_1<\infty \quad\text{and}\quad \frac{\alpha_k}{\sqrt{n_k}}\to c_2<\infty. $$ Then by the WLLN and CLT (assuming that $\mathsf{E}\|x\|^2<\infty$ and $\mathsf{E}xx^{\top}$ is positive definite), using the fact that two samples are independent, \begin{align} \alpha_k\left(\hat{\mathsf{E}}y-\mathsf{E}y\right)&=\frac{\alpha_k}{\sqrt{m_k}}\hat{\beta}^{\top}\sqrt{m_k}\left(\hat{\mathsf{E}}x-\mathsf{E}x\right)+\frac{\alpha_k}{\sqrt{n_k}}\hat{\mathsf{E}}x^{\top}\sqrt{n_k}\left(\hat{\beta}-\beta\right) \\ &=c_1\beta^{\top}N_1+c_2\mathsf{E}x^{\top}N_2+o_p(1), \end{align} where $N_1\sim \mathcal{N}(0,\operatorname{Var}(x))$, $N_2\sim \mathcal{N}(0,\sigma^2(\mathsf{E}xx^{\top})^{-1})$, and $N_1$ and $N_2$ are independent. Therefore, $$ \bbox[cornsilk,5px] { \alpha_k\left(\hat{\mathsf{E}}y-\mathsf{E}y\right)\xrightarrow{d}\mathcal{N}\left(0,c_1^2\beta^{\top}\operatorname{Var}(x)\beta+c_2^2\sigma^2\mathsf{E}x^{\top}\left(\mathsf{E}xx^{\top}\right)^{-1}\mathsf{E}x\right). } $$