Outcome-Specific Predictor Sets in a Joint Regression

Different Mean Models, One Joint Response

In a multivariate regression, each observational unit supplies a vector of outcomes. It is often reasonable for those outcomes to have different mean structures: one outcome may depend on treatment and age, while another depends only on age. Different predictor sets do not require treating the outcomes as unrelated. A joint model can specify a separate regression for each outcome while allowing their errors to covary within an observational unit.

The key is to keep two modelling choices distinct. The predictor set for outcome $j$ determines its conditional mean. The covariance matrix of the error vector describes how outcomes vary together after those means have been accounted for. A joint representation makes both choices explicit.

Definition: For $n$ observational units and $q$ outcomes, let $Y$ be the $n\times q$ response matrix. In an outcome-specific regression model, column $j$ has mean $X_j\beta_j$, where $X_j$ is an $n\times p_j$ design matrix and $\beta_j$ is a $p_j$-vector. Write $y_j$ for column $j$ of $Y$, and $e_j=y_j-X_j\beta_j$ for its error vector. The errors across outcomes for unit $i$ form a vector whose covariance matrix is $\Sigma$.

For the usual covariance model, observational units are independent, have a common error covariance matrix $\Sigma$, and satisfy $E[e_i]=0$, where $e_i$ is the $q$-vector of errors for unit $i$. Thus $\operatorname{Cov}(e_i)=\Sigma$, while $\operatorname{Cov}(e_i,e_k)=0$ for $i\ne k$. If a Gaussian model is assumed, the row error vectors are additionally multivariate normal. The predictor sets may differ even though the same $\Sigma$ describes the joint errors.

Two Equivalent Ways to Write the Mean Model

One representation gives each outcome its own design matrix. Use column-stacking vectorization, so $\operatorname{vec}(Y)$ stacks the first column of $Y$, then the second, and so on. Define

D=\operatorname{blockdiag}(X_1,\ldots,X_q),\qquad \beta= \begin{pmatrix} \beta_1\\ \vdots\\ \beta_q \end{pmatrix}.

Then $D$ is an $nq\times \sum_j p_j$ matrix, and the model mean is $D\beta$. If the error rows are independent with common covariance $\Sigma$, the stacked error covariance is $\Sigma\otimes I_n$, where $\otimes$ denotes the Kronecker product. Consequently, the joint model is

\operatorname{vec}(Y)=D\beta+\operatorname{vec}(E), \qquad \operatorname{Cov}\{\operatorname{vec}(E)\}=\Sigma\otimes I_n.

The off-diagonal elements of $\Sigma$ are not removed by using different $X_j$. For example, $\Sigma_{12}$ is the covariance between the errors for outcomes 1 and 2 on the same unit.

A second representation starts with a common list of candidate predictors. Let $X$ be an $n\times p$ design matrix and let $B$ be a $p\times q$ coefficient matrix, so the mean is $XB$. Predictors absent from an outcome's model are represented by zero entries in that outcome's column of $B$. More generally, restrictions can be encoded by writing $\operatorname{vec}(B)=R\gamma+b_0$, where $R$ and $b_0$ specify the allowed coefficient values. Substitution gives the stacked model

\operatorname{vec}(Y) =(I_q\otimes X)(R\gamma+b_0)+\operatorname{vec}(E).

Worked Example: Encoding an Outcome-Specific Predictor Set

Suppose outcome 1 depends on an intercept, age, and treatment, while outcome 2 depends on an intercept and age but not treatment. Use the common design $X=[\mathbf 1,\ a,\ t]$, where $a$ is the age column and $t$ is the treatment column. Write

B= \begin{pmatrix} \beta_{01}&\beta_{02}\\ \beta_{a1}&\beta_{a2}\\ \beta_{t1}&0 \end{pmatrix}.

The zero in row 3, column 2 imposes the required absence of a treatment effect for outcome 2. In stacked coefficient order, $\operatorname{vec}(B)=(\beta_{01},\beta_{a1},\beta_{t1},\beta_{02},\beta_{a2},0)^\top$. If $\gamma=(\beta_{01},\beta_{a1},\beta_{t1},\beta_{02},\beta_{a2})^\top$, this is $R\gamma$, where $R$ is the $6\times5$ matrix formed by taking the first five columns of $I_6$. Thus the restricted mean is $(I_2\otimes X)R\gamma$. The error covariance remains $\Sigma\otimes I_n$; the coefficient restriction changes the mean model, not the covariance assumption.

Least Squares With Separate Predictor Sets

For ordinary multivariate least squares, the criterion is the sum of squared residuals across all outcomes. If each $X_j$ has full column rank, minimizing this criterion gives the usual least-squares estimate for each outcome separately. The cross-outcome covariance does not enter this particular criterion, although it remains a feature of the joint error distribution.

Theorem: Suppose each $X_j$ has full column rank. The minimizer of $\sum_{j=1}^q\|y_j-X_j\beta_j\|^2$ is $\hat\beta_j=(X_j^\top X_j)^{-1}X_j^\top y_j$ for every $j$.

To prove the claim, note that the stacked least-squares criterion separates:

\|\operatorname{vec}(Y)-D\beta\|^2 =\sum_{j=1}^q\|y_j-X_j\beta_j\|^2.

Each term depends on only one coefficient vector. Differentiating its quadratic form with respect to $\beta_j$ gives the normal equation $X_j^\top X_j\hat\beta_j=X_j^\top y_j$. Full column rank makes $X_j^\top X_j$ invertible, yielding the stated estimator. Since every term is minimized separately, the collection of these estimates minimizes the full criterion. $\square$

This result concerns least squares, not every possible fitting criterion. In particular, under a Gaussian model with non-diagonal $\Sigma$, maximizing the joint likelihood over coefficients uses the quadratic form $\operatorname{tr}\{(Y-M)\Sigma^{-1}(Y-M)^\top\}$, where $M$ is the mean matrix. Its cross-outcome terms can couple the coefficient estimates when the predictor sets differ. Therefore, separate least-squares fitting should not automatically be called joint Gaussian maximum likelihood.

Worked Example: Fitting Two Different Mean Structures

Let $x=(-1,0,1)^\top$. Fit outcome 1 with an intercept and $x$, and outcome 2 with an intercept only. Take $y_1=(0,2,2)^\top$ and $y_2=(0,1,3)^\top$. The first design is $X_1=[\mathbf 1,x]$, and the second is $X_2=\mathbf 1$.

For outcome 1, $\mathbf 1^\top x=0$, $\mathbf 1^\top\mathbf 1=3$, and $x^\top x=2$. Also, $\mathbf 1^\top y_1=4$ and $x^\top y_1=2$. Therefore,

\hat\beta_1= \begin{pmatrix} 4/3\\ 1 \end{pmatrix}, \qquad \hat y_1= \begin{pmatrix} 1/3\\ 4/3\\ 7/3 \end{pmatrix}, \qquad e_1=y_1-\hat y_1= \begin{pmatrix} -1/3\\ 2/3\\ -1/3 \end{pmatrix}.

For outcome 2, the intercept-only estimate is its sample mean, $\hat\beta_2=(0+1+3)/3=4/3$. Thus $\hat y_2=(4/3,4/3,4/3)^\top$, and $e_2=(-4/3,-1/3,5/3)^\top$. The realized residual cross-product is

e_1^\top e_2 =(-1/3)(-4/3)+(2/3)(-1/3)+(-1/3)(5/3) =4/9-2/9-5/9=-1/3.

This calculation is consistent with the residual definition $e_j=y_j-\hat y_j$: the first residual for outcome 1 is $-1/3$, not $1/3$. The nonzero cross-product describes these fitted residual vectors; by itself, it is not the assumed population error covariance.

How Cross-Outcome Covariance Enters

Even when ordinary least squares estimates the coefficient vectors separately, their sampling variation is generally related. This follows because the outcomes' errors for the same observational unit may be correlated.

Theorem: Suppose $E[e_j]=0$, errors are independent across units, and for every unit $i$, $\operatorname{Cov}(e_{i,j},e_{i,k})=\sigma_{jk}$. If $X_j$ and $X_k$ have full column rank and the coefficient estimates are ordinary least squares, then $$ \operatorname{Cov}(\hat\beta_j,\hat\beta_k) =\sigma_{jk}(X_j^\top X_j)^{-1}X_j^\top X_k(X_k^\top X_k)^{-1}. $$

For the derivation, write $\hat\beta_j=\beta_j+A_j e_j$, where $A_j=(X_j^\top X_j)^{-1}X_j^\top$, and similarly define $A_k$. The errors have cross-covariance $\operatorname{Cov}(e_j,e_k)=\sigma_{jk}I_n$: entries from different units have covariance zero, and corresponding entries have covariance $\sigma_{jk}$. Applying the covariance rule for linear transformations gives

\begin{aligned} \operatorname{Cov}(\hat\beta_j,\hat\beta_k) &=A_j\operatorname{Cov}(e_j,e_k)A_k^\top\\ &=\sigma_{jk}(X_j^\top X_j)^{-1}X_j^\top X_k(X_k^\top X_k)^{-1}. \end{aligned}

This establishes the result. $\square$ The formula also shows that coefficient covariance depends on both the error covariance and the relationship between the two design matrices; it need not be nonzero for every pair of coefficient summaries.

Worked Example: Covariance Between Coefficient Estimates

Continue with $X_1=[\mathbf 1,x]$, $X_2=\mathbf 1$, and $x=(-1,0,1)^\top$. Suppose the error covariance between the two outcomes is $\sigma_{12}=0.5$. We have $X_1^\top X_1=\operatorname{diag}(3,2)$, $X_2^\top X_2=3$, and

X_1^\top X_2= \begin{pmatrix} 3\\ 0 \end{pmatrix}.

Substituting into the covariance formula gives

\begin{aligned} \operatorname{Cov}(\hat\beta_1,\hat\beta_2) &=0.5 \begin{pmatrix} 1/3&0\\ 0&1/2 \end{pmatrix} \begin{pmatrix} 3\\ 0 \end{pmatrix} \frac{1}{3}\\ &= \begin{pmatrix} 1/6\\ 0 \end{pmatrix}. \end{aligned}

The estimated intercepts have covariance $1/6$, while the estimated slope for outcome 1 has covariance zero with the intercept estimate for outcome 2 in this design. The zero follows from $x^\top\mathbf 1=0$, not from assuming independent outcome errors.

Interpreting the Joint Model Carefully

A convenient joint specification might take

\Sigma= \begin{pmatrix} 1&0.5\\ 0.5&4 \end{pmatrix}.

This is a valid covariance matrix: its diagonal entries are positive and its determinant is $1(4)-(0.5)^2=3.75>0$. It says that, after accounting for the respective mean models, the two outcomes have error variances 1 and 4 and error covariance 0.5 within a unit. It does not require the same predictors for both outcomes.

In the numerical least-squares example, the residual cross-product was $-1/3$, not $0.5$. There is no contradiction. A realized residual cross-product is calculated from one data set after fitting; $\sigma_{12}$ is a population covariance parameter in the error model. Fitting and uncertainty quantification must respect their distinction. In particular, with different designs, the residual vectors arise from different projections, so their cross-product should not be casually interpreted as an unbiased estimate of $\sigma_{12}$.

Common Pitfalls

Assuming all outcomes need identical predictors. A joint response means the outcomes are modelled together; it does not force their mean structures to match.
Confusing zero coefficients with zero error covariance. Omitting a predictor from an outcome sets a mean coefficient to zero. It says nothing by itself about the covariance between outcome errors.
Calling separate least squares the joint Gaussian maximum-likelihood fit. With correlated errors, the likelihood weights residuals jointly. Under differing predictor sets, its coefficient estimates can be coupled.
Reading a sample residual cross-product as the covariance parameter. A realized residual product is data-dependent and affected by fitting. It need not equal the population error covariance, even in sign.

Key takeaway: Outcome-specific designs describe which predictors enter each conditional mean. A stacked design or coefficient restrictions express those choices within one model, while $\Sigma\otimes I_n$ preserves the covariance among outcomes measured on the same unit.

Check Your Understanding

Work through these questions to check your understanding of outcome-specific predictor sets in a joint regression.

For three outcomes with design matrices $X_1,X_2,X_3$, what is the block structure of the stacked design under column-stacking vectorization?
How would you encode the absence of a treatment effect for outcome 2 when treatment is one column of a common design matrix?
Why does a nonzero off-diagonal element of $\Sigma$ not require the outcomes to share predictors?
State the assumptions needed for the displayed covariance formula for $\hat\beta_j$ and $\hat\beta_k$.
Why is a realized residual cross-product not automatically the population error covariance?