Why Fit Two Outcomes Together?
In many studies, each subject has several measured outcomes and the same predictors are relevant to each. For example, a clinical study might relate treatment group and baseline severity to both systolic and diastolic blood pressure. A shared-predictor regression gives a compact way to estimate the two outcome-specific coefficient sets while retaining their residual covariance in a single matrix.
The shared design matrix does not force the outcomes to have the same coefficients. It means only that the same predictor columns are used for each response. The coefficient matrix has one column per outcome, and the residual sum-of-squares-and-cross-products matrix records both the residual variation of each outcome and their residual covariation.
For the model considered here, \(X\) includes a column of ones for the intercept and one or more columns for shared predictors. If \(X\) has full column rank, the least-squares estimate is \(\hat B=(X'X)^{-1}X'Y\). As in the earlier material on multivariate OLS, the regression is computed for each response; the matrix notation displays the two fits together.
How the Shared-Design Estimate Arises
Least squares chooses \(B\) to minimize the total squared residuals across both responses. Equivalently, it minimizes the trace of the residual cross-product matrix:
Let \(\hat B=(X'X)^{-1}X'Y\) and \(\Delta=B-\hat B\). Then \(Y-XB=(Y-X\hat B)-X\Delta\), and \(X'(Y-X\hat B)=0\) by the normal equations. Expanding the trace therefore gives
Because \(X'X\) is positive definite when \(X\) has full column rank, the added term is nonnegative and is zero only when \(\Delta=0\). Thus \(\hat B\) is the unique least-squares estimate. The same \(X\) appears in both response columns, but each column of \(\hat B\) is estimated from its corresponding response values.
Worked Example: Fitting Two Outcomes With One Predictor
Worked Example: Fitting Two Outcomes With One Predictor
Suppose five subjects have predictor values \(x=-2,-1,0,1,2\). The two outcomes are \(Y_1\) and \(Y_2\), recorded in the table. Fit a regression with an intercept and the shared predictor \(x\); then calculate its fitted values, residual SSCP, and covariance estimates.
| Subject | \(x\) | \(Y_1\) | \(Y_2\) |
|---|---|---|---|
| 1 | -2 | 7 | 8 |
| 2 | -1 | 6 | 2 |
| 3 | 0 | 12 | 11 |
| 4 | 1 | 10 | 0 |
| 5 | 2 | 15 | 4 |
1. Set up the matrices. Use an intercept column followed by the predictor column:
2. Calculate the coefficient matrix. The predictor values sum to zero and their squared values sum to \(10\), so
The outcome column sums are \(50\) and \(25\). The predictor-weighted sums are \(\sum_i x_iY_{i1}=20\) and \(\sum_i x_iY_{i2}=-10\). Thus
The first row gives the estimated intercepts, \(10\) for \(Y_1\) and \(5\) for \(Y_2\). The second row gives the estimated slopes, \(2\) and \(-1\). For a one-unit increase in \(x\), the fitted mean of \(Y_1\) increases by \(2\), while that of \(Y_2\) decreases by \(1\).
3. Calculate fitted values and residuals. The fitted outcome values are \(10+2x_i\) and \(5-x_i\). Subtracting these fitted values from the observed values gives
4. Calculate the residual SSCP matrix. The diagonal entries are the residual sums of squares for the two outcomes; the off-diagonal entries are their residual cross-products:
The matrix is symmetric, as a cross-product matrix must be. Its positive determinant, \(14(70)-30^2=80\), also confirms that this particular SSCP matrix is positive definite.
From Residual SSCP to Error Covariance
Two covariance estimates commonly accompany this fit. Under independent, mean-zero errors with common covariance \(\Sigma\), and with \(X\) fixed and full rank, the Gaussian maximum-likelihood estimate is \(S_R/n\). The unbiased estimate is \(S_R/(n-p)\). Here \(n-p\) is the residual degrees of freedom, because \(p\) coefficient rows, including the intercept, have been fitted for each outcome.
To verify the expectation, note that \(M\) is symmetric and idempotent, \(R=M\mathcal{E}\), and \(R'R=\mathcal{E}'M\mathcal{E}\). Since the rows \(\varepsilon_i\) of \(\mathcal{E}\) are row vectors, expansion gives
For \(i\ne j\), independence and zero means give \(E[\varepsilon_i'\varepsilon_j]=0\). For \(i=j\), \(E[\varepsilon_i'\varepsilon_i]=\Sigma\). Consequently,
The last equality follows because \(M\) is a projection of rank \(n-p\), so its trace is \(n-p\). Dividing the expectation by \(n-p\) proves unbiasedness. Under Gaussian errors, maximizing the likelihood over the mean coefficients first gives the residual matrix \(R\); maximizing over the common covariance then gives \(R'R/n\), provided a nonsingular maximizing covariance exists. \(\square\)
Worked Example: Calculating Both Covariance Estimates
Continue with the fitted regression above. There are \(n=5\) observations and \(p=2\) design columns, so the residual degrees of freedom are \(5-2=3\). The maximum-likelihood estimate is
The unbiased estimate is
The estimates differ because the maximum-likelihood divisor is \(n=5\), whereas the unbiased divisor is the residual degrees of freedom, \(3\). Both have positive determinant: the first has determinant \(2.8(14)-6^2=3.2\), and the second has determinant \((14/3)(70/3)-10^2=80/9\). In this fit, both are valid positive-definite covariance matrices.
Worked Example: Interpreting the Residual Covariance
Use the unbiased covariance estimate from the fitted model to describe residual variation after accounting for \(x\). Its diagonal entries estimate the residual variances: \(14/3\) for \(Y_1\), and \(70/3\) for \(Y_2\). Its off-diagonal entry, \(10\), estimates the covariance between the two outcomes’ errors.
The corresponding estimated residual correlation is obtained by dividing the covariance by the product of the estimated residual standard deviations:
Thus, in this small data set, observations with a residual above the fitted value for one outcome tend also to have a residual above the fitted value for the other. This is a statement about residual association after fitting the shared predictor; it does not say that one outcome causes the other. The estimate is based on only three residual degrees of freedom, so it should not be treated as a precise description of a broader population.
Common Mistake: Mixing Up the Matrix Dimensions
A frequent error is to treat the row error \(\varepsilon_i\) as a column vector when expanding \(\mathcal{E}'M\mathcal{E}\). Here \(\varepsilon_i\) is a \(1\times 2\) row. The product \(\varepsilon_i'\varepsilon_j\) is therefore a \(2\times 2\) matrix, as required in the SSCP expansion. In contrast, \(\varepsilon_i\varepsilon_j'\) is a scalar and cannot be added as a term in that \(2\times 2\) sum. Keeping the row-vector convention consistent also gives \(\Sigma=E[\varepsilon_i'\varepsilon_i]\).
Another easy mistake is to divide the SSCP by the wrong quantity without naming the target. Dividing by \(n\) gives the Gaussian maximum-likelihood covariance estimate; dividing by \(n-p\) gives the unbiased covariance estimate under the stated error assumptions. Neither divisor should be substituted for the other without saying which estimate is being reported.
Check Your Understanding
Use the definitions and calculations above to check your understanding.
- For \(n\) observations, two responses, and an intercept plus two predictors, what are the dimensions of \(X\), \(Y\), and \(B\)?
- In the worked fit, what does the second row of \(\hat B\) say about the change in each fitted outcome for a one-unit increase in \(x\)?
- Why is \(R'R\) a \(2\times 2\) matrix, and what does each type of entry represent?
- Under the stated assumptions, why does the unbiased covariance estimate use \(n-p\) rather than \(n\)?
- For row-vector errors \(\varepsilon_i\), which product, \(\varepsilon_i'\varepsilon_j\) or \(\varepsilon_i\varepsilon_j'\), has the dimensions needed in the SSCP expansion?