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

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.

• Are $\epsilon_1$ and $\epsilon_2$ independent? Are they mean zero?
– Ben
Jun 23 at 5:03
• Yes, they are independent and zero mean
– Amin
Jun 23 at 5:43
• Let $U_1$ be an unbiased estimating equation for the first outcome and likewise define $U_2$ for the second outcome. Then the stacked estimating equations $(U_1(Y_1, X_1), I(Y_1=1)/P(Y_1=1 \mid X_1) U_2(Y_2, X_1, X_2))$ will be unbiased and should be asymptotically normal.
– Ben
Jun 27 at 19:07
• Thanks Ben. I am wondering how it is possible to use this model? Because we dont know $P(Y_1=1|X_1)$ and it should be learned too.
– Amin
Jun 27 at 21:08
• You have to estimate it, so for example in a parametric setting it will depend on e.g. $\gamma \in \mathbb{R}^p$, where the first component of the estimating equation (involving $U_1$) also depends on $\gamma$.
– Ben
Jun 27 at 21:16

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.

• Thanks Ben. I think it is the general model for the multivariate regression model. But my problem is a bit different because the output of the second model will be observed when the output of the first model is 1. For example, if $y_1$ show that whether a clinic accepts or reject a person and $y_2$ show whether the person is classified as ill or not. So, we will observe $y_2$ just if $y_1=1$ (there might be some causality, but this an example to show how it works).
– Amin
Jun 26 at 3:38
• Okay, I'll have more of a think about it and see if I can amend my answer to say something more useful.
– Ben
Jun 26 at 3:39