Confidence interval for covariance in multivariate linear regression

Question

Let $X\in R^{n\times q}$ and $Y\in R^{n\times p}$ be our data matrices, and we assume that they are related by a linear model $Y = BX + \Xi$, where $B$ is fixed and unknown and each row $\Xi_i \sim \mathcal{N}(0, \Sigma)$ is iid drawn from a multivariate Gaussian with fixed but unknown covariance matrix $\Sigma$.

We want to estimate the covariance of the two dependent variables, controlled for the independent variable. If I am correct, this is represented by the off-diagonal elements of $\Sigma$.

We can estimate $\Sigma$ as $\hat \Sigma = \frac{1}{n-q-1}(Y-X\hat B)^T(Y-X\hat B)$, where $\hat B = (X^T X)^{-1}X^T Y$.

But I want to double check how confident I should be in this estimate, ie I want to compute a confidence interval around the off-diagonal parameter of $\hat \Sigma$.

How do I do that?

Answer 1

The bootstrap is by far the most common approach to this kind of question. It's common in practice but also well-supported by theory and well-accepted by statisticians. The percentile bootstrap (generate K bootstrap samples, compute the off-diagonal in each, and take the 2.5%ile and 97.%ile of bootstrap replicates as confidence bounds) is the simplest, most common, and a good choice. The "accelerated and bias-corrected" bootstrap has slightly better properties, as I understand, but you have to do some math. The relevant chapter here is a good starting point.

Beyond the bootstrap, I think Fisher did some work on a variance-stabilizing transformation to estimate a correlation coefficient using a normal approximation. Maybe it could be used in the linear model context here, I'm not sure. It's definitely not a standard question, to get a confidence interval for the error covariance in a multi-response linear model.

If the sample size really is 30, like in the comment, the bootstrap might still be OK, but if one really thinks the sample is too small for a bootstrap, then an analytical solution isn't going to help you, because it's going to rely on the assumption of normal errors which probably isn't justified. That's my intuition at least. So one could spend time thinking about an analytical solution but I don't feel it'd actually be a good idea. You could run simulations with reasonable values of $n,X,B,\Sigma$ given your problem and see how the bootstrap performs.

A better reason to avoid the bootstrap would be if the sample were extremely large and the computation had to be done often. Then the bootstrap is just computationally infeasible. But then you can rely heavily on asymptotics (LLN and CLT) and get some formula that way, maybe after doing some Taylor expansion.

Minor corrections: I think you want $Y = XB + \Xi$, and that we want to estimate the covariance of any pair of the $p$ dependent variables (or are considering the $p = 2$ case). The formula for $\hat{\Sigma}$ can have $1 / (n - q)$ instead of $1 / (n - q - 1)$ if we assume $X$ already includes any desired intercept (as the formulas indicate).