# Finding covariance between internally studentized residuals

Consider a linear regression model $$\boldsymbol y=X\boldsymbol\beta+\boldsymbol\varepsilon$$, where $$\boldsymbol y$$ is an $$n\times 1$$ response vector, $$X$$ is an $$n\times p$$ matrix of covariates (fixed), and the error vector $$\boldsymbol\varepsilon$$ is multivariate normal $$N_n(\boldsymbol 0,\sigma^2I)$$.

The $$i$$th internally studentized residual is

$$r_i=\frac{e_i}{\hat\sigma\sqrt{1-h_{ii}}}\,,$$

where $$e_i=y_i-\hat y_i$$ is the $$i$$th residual, $$h_{ii}$$ is the $$i$$th diagonal entry of the hat matrix $$H=X(X^TX)^{-1}X^T$$, and $$\hat\sigma^2=\frac{\boldsymbol e^T\boldsymbol e}{n-p}$$ is the usual unbiased estimator of $$\sigma^2$$.

I am trying to find $$\operatorname{Cov}(r_i,r_j)$$ without explicitly finding the joint distribution of $$(r_i,r_j)$$.

I can see that $$\operatorname E(r_i)=0$$, since $$\frac{e_i}{\lVert \boldsymbol e\rVert}$$ is symmetric about $$0$$.

Now, for $$\operatorname{Cov}(r_i,r_j)=\operatorname E(r_ir_j)$$, I need $$\operatorname{Cov}\left(\frac{e_i}{\lVert \boldsymbol e\rVert},\frac{e_j}{\lVert \boldsymbol e\rVert}\right)=\operatorname E\left[\frac{e_ie_j}{\boldsymbol e^T\boldsymbol e}\right]$$

Any suggestion in finding this quantity is welcome.

I know that $$\boldsymbol e\sim N_n(\boldsymbol 0,\sigma^2(I-H))$$ and $$\frac{r_i^2}{n-p}\sim \text{Beta}\left(\frac12,\frac{n-p-1}{2}\right)$$.

The answer is supposed to be the same as the covariance if $$\sigma$$ was not estimated:

$$\operatorname{Cov}(r_i,r_j)=\operatorname{Cov}\left(\frac{e_i}{\sigma\sqrt{1-h_{ii}}},\frac{e_j}{\sigma\sqrt{1-h_{jj}}}\right)=\frac{-h_{ij}}{\sqrt{1-h_{ii}}\sqrt{1-h_{jj}}}\,,$$

where $$h_{ij}$$ is the $$(i,j)$$th entry of $$H$$.

#### Preparations

Lemma 1. If $$X = (X_1, X_2, \ldots, X_d) \sim N_d(0, I_{(d)})$$, then \begin{align} & E\left[\frac{X_j}{\|X\|}\right] = 0, 1 \leq j \leq d, \\ & E\left[\frac{X_iX_j}{\|X\|^2}\right] = \begin{cases} 0, & 1 \leq i \neq j \leq d, \\ d^{-1}, & 1 \leq i = j \leq d. \end{cases} \end{align}

Lemma 2. Suppose an order $$n$$ matrix $$A$$ is idempotent and symmetric with $$\operatorname{rank}(A) = r$$, then there exists an order $$n$$ orthogonal matrix $$P$$ such that \begin{align} A = P\operatorname{diag}(I_{(r)}, 0)P'. \end{align}

#### Proof

Without loss of generality, assume $$\sigma^2 \equiv 1$$. I am also assuming you are only interested in the case $$i \neq j$$ (the case $$i = j$$ is similar).

Denote $$I_{(n)} - H$$ by $$\bar{H}$$, then $$\bar{H}$$ is symmetric, idempotent, and $$\operatorname{rank}(\bar{H}) = n - p =: r$$. By Lemma 2, there exists an order $$n$$ orthogonal matrix $$P$$ such that $$\bar{H} = P\operatorname{diag}(I_{(r)}, 0)P'$$. Denote $$P'\varepsilon$$ by $$X := (\xi, \eta)$$, where $$\xi \in \mathbb{R}^r$$ and $$\eta \in \mathbb{R}^{n - r}$$. Since $$P$$ is orthogonal and $$\varepsilon\sim N_n(0, I_{(n)})$$, it follows that $$X \sim N_n(0, I_{(n)})$$ and $$\xi \sim N_r(0, I_{(r)})$$. Denote the matrix formed by the first $$r$$ columns of $$P$$ by $$Q$$ and the matrix formed by the last $$n - r$$ columns of $$P$$ by $$R$$, then \begin{align} & e = \bar{H}\varepsilon = P\operatorname{diag}(I_{(r)}, 0)P'\varepsilon = \begin{bmatrix} Q & R \end{bmatrix}\begin{bmatrix} \xi \\ 0 \end{bmatrix} = Q\xi. \tag{1} \\ & e'e = \|e\|^2 = \xi'Q'Q\xi = \|\xi\|^2. \tag{2} \end{align}

Denote the $$i$$-th row and $$j$$-th row of $$Q$$ by $$\begin{bmatrix}q_{i1} & \cdots & q_{ir} \end{bmatrix}$$ and $$\begin{bmatrix}q_{j1} & \cdots & q_{jr} \end{bmatrix}$$ respectively, it then follows by $$(1)$$ and $$(2)$$ that \begin{align} E\left[\frac{e_ie_j}{\|e\|^2}\right] = E\left[\frac{\sum_{k = 1}^r q_{ik}\xi_k \cdot \sum_{l = 1}^r q_{jl}\xi_l}{\|\xi\|^2}\right], \tag{3} \end{align} which by Lemma 1 and $$\xi \sim N_r(0, I_{(r)})$$ becomes to \begin{align} \sum_{k = 1}^r q_{ik}q_{jk}E\left[\frac{\xi_k^2}{\|\xi\|^2} \right] = r^{-1}\sum_{k = 1}^r q_{ik}q_{jk}. \tag{4} \end{align}

Finally, $$\bar{H} = I_{(n)} - H = QQ'$$ implies that \begin{align} \sum_{k = 1}^r q_{ik}q_{jk} = \bar{H}(i, j) = -h_{ij}. \tag{5} \end{align} Combining $$(3), (4), (5)$$ yields: \begin{align} E\left[\frac{e_ie_j}{\|e\|^2}\right] = -r^{-1}h_{ij}, \end{align} which immediately gives \begin{align} \operatorname{Cov}(r_i, r_j) = \frac{-h_{ij}}{\sqrt{1 - h_{ii}}\sqrt{1 - h_{jj}}}. \end{align}

#### Proof of Lemmas

Lemma 2 is a standard linear algebra result. So I will only discuss how to prove Lemma 1, which is a quite interesting result. Lemma 1 is a corollary of the following well-known theorem of spherical distributions (cf. Theorem 1.5.6 in Aspects of Multivariate Statistical Theory by R. Muirhead):

Theorem. If $$X$$ has an $$m$$-variate spherical distribution with $$P(X = 0) = 0$$ and $$r = \|X\| = (X'X)^{1/2}, T(X) = \|X\|^{-1}X$$, then $$T(X)$$ is uniformly distributed on $$S_m$$ and $$T(X)$$ and $$r$$ are independent.

The proof to this theorem is not hard (but might be tedious) by laying out the spherical coordinates of $$X$$.

To apply this theorem to the proof of Lemma 1, write \begin{align} & X_j = \frac{X_j}{\|X\|} \times r, \tag{6} \\ & X_iX_j = \frac{X_iX_j}{\|X\|^2} \times r^2. \tag{7} \end{align}

Then taking expectations on both sides of $$(6), (7)$$ and applying the independence of multipliers (by Theorem) finishes the proof.