How to construct a confidence interval for the coefficients of a multivariate regression with dependence between dependent variables?

Question

Suppose we have two linear regression models $y_1=a+bx+\epsilon_1$ where $\mathbb[\epsilon_1]=\sigma_1$ and $y_2=c+dx+\epsilon_2$ where $\mathbb[\epsilon_2]=\sigma_2$. In other words, I am using the variable $x$ in two regression models. We can estimate parameters $a,b,c,d$ separately through OLS, based on historical data find the covariance matrix for $(\hat{a},\hat{b})^T$. If we use $\beta=(a,b)^T$, then, we will have:

$$ \hat{\beta}-\beta=\left(X^TX \right)^{-1}X^T\epsilon $$ Assuming that noises are sub-gaussian (or anything else), we can find the confidence interval for $||\hat{\beta}-\beta||$ using Cram´er–Chernoff method. Similarly, we can use a similar procedure to find a confidence bound for $||\hat{\nu}-\nu||$ where $\nu=(c,d)^T$.

However, these confidence bounds are for parameters of each regression model separately. I gonna have confidence for the unified vector where all parameters are considered together. That is,

$$ \Pr\left(|\hat{\mu}-\mu|\ge e \right)\le \delta $$ where $\mu=(a,b,c,d)$

Additional Note: Here, I will not use exactly the whole values of $x_i$ for the section regression and just a subset of them will be used. For example, the first regression uses 100 observations but the second model uses 10 observations. To add more information, we just observe $y_2$ conditioning on the value of $y_1$. If we use this model as a classification model, then we observe $y_2$ when $y_1=1$.

I would be thankful to know your opinion how to handle it.

Answer 1

Although your question asks about a model where the error terms are independent, I'm going to generalise to drop this assumption in order to give a slightly broader answer. This will allow me to show you the effect of that assumption on your analysis more easily, and I hope it will also give a bit more insight into the interactions occurring in the multivariate regression model.

Also, at the outset it is worth noting that the analytical results for the multivariate regression model turn out to be mostly identical to —or analogous to— the case where we have a single response variable. This is true of the OLS estimator, the hat matrix, the error variance estimator, and many of the other statistical results. I'll show some of this below.

---

For your problem, let's start by writing your overall model equation in matrix form:

$$\begin{bmatrix} y_{1,1} & y_{1,2} \\ y_{2,1} & y_{2,2} \\ \vdots & \vdots \\ y_{n,1} & y_{n,2} \end{bmatrix} = \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{bmatrix} \begin{bmatrix} a & c \\ b & d \end{bmatrix} + \begin{bmatrix} \epsilon_{1,1} & \epsilon_{1,2} \\ \epsilon_{2,1} & \epsilon_{2,2} \\ \vdots & \vdots \\ \epsilon_{n,1} & \epsilon_{n,2} \end{bmatrix} \quad \quad \quad \quad \quad \begin{bmatrix} \epsilon_{i,1} \\ \epsilon_{i,2} \end{bmatrix} \sim \text{IID N}(\mathbf{0}, \boldsymbol{\Sigma}).$$

Here I am using a general variance matrix $\boldsymbol{\Sigma}$ for the error variance, without assuming that the errors of the two parts are uncorrelated (but we can look at the specific case of interest to you later). Defining the design matrix $\mathbf{x}$ and appropriate vectors for the responses, coefficients and errors of the two parts, we can write this model equation in more succinct form as:

$$\begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix} = \mathbf{x} \begin{bmatrix} \boldsymbol{\beta}_1 & \boldsymbol{\beta}_2 \end{bmatrix} + \begin{bmatrix} \boldsymbol{\epsilon}_1 & \boldsymbol{\epsilon}_2 \end{bmatrix}.$$

In multivariate regression, it turns out that the OLS estimator for the overall model is the same as imposing OLS estimation on the two parts. The OLS estimator is given by:$^\dagger$

$$\begin{align} \begin{bmatrix} \hat{\boldsymbol{\beta}}_1 & \hat{\boldsymbol{\beta}}_2 \end{bmatrix} &= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix} \\[6pt] &= \begin{bmatrix} (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbf{y}_1 & (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbf{y}_2 \end{bmatrix}. \\[6pt] \end{align}$$

Assuming the model is correctly specified, the OLS estimator is an unbiased estimator of the true coefficients. If we write the coefficient estimators as a single vector $\hat{\boldsymbol{\beta}} = [\hat{a},\hat{b},\hat{c},\hat{d}]^\text{T}$, this vector has variance matrix:

$$\begin{align} \mathbb{V}(\hat{\boldsymbol{\beta}}) &= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \otimes \boldsymbol{\Sigma} = \begin{bmatrix} \sigma_1^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \\ \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \sigma_2^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \end{bmatrix}. \end{align}$$

Moreover, under the Gaussian version of the model the OLS estimator is normally distributed, so we have:

$$\hat{\boldsymbol{\beta}} - \boldsymbol{\beta} \sim \text{N} \Bigg( \begin{bmatrix} \mathbf{0} \\[4pt] \mathbf{0} \end{bmatrix}, \begin{bmatrix} \sigma_1^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \\ \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \sigma_2^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \end{bmatrix} \Bigg).$$

From the above form, we can see that the estimators for the two coefficient vectors in the model parts have covariance $\mathbb{C}(\hat{\boldsymbol{\beta}}_1, \hat{\boldsymbol{\beta}}_2) = \rho \sigma_{1} \sigma_{2} (\mathbf{x}^\text{T} \mathbf{x})^{-1}$. Consequently, in the case where the underlying errors for the two model parts are uncorrelated (i.e., when $\rho = 0$) the coefficient estimators for the two model parts are also uncorrelated.

In order to form a confidence set for $\boldsymbol{\beta}$ you also need to estimate the error variance matrix $\boldsymbol{\Sigma}$. The standard estimator is the adjusted MLE (adjusted to be unbiased), which is:

$$\hat{\boldsymbol{\Sigma}} = \frac{\begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix}^\text{T} (\mathbf{I}_n - \mathbf{x} (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T}) \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix}}{n-2}.$$

Using these results you can form a confidence set for $\boldsymbol{\beta}$, which will be a hyper-ellipse in four-dimensional Euclidean space.

---

$^\dagger$ It's also worth noting that the OLS estimator is equivalent to the MLE in the Gaussian case, just as in standard regression.