Tutorials › Multivariate Statistics (Graduate) › Worked Tests of Multivariate Regression Effects

Home  /  Tutorials  /  Hypothesis Testing in Linear Models

Hypothesis Testing in Linear Models · Tutorial 328 of 975

Worked Tests of Multivariate Regression Effects

Learn how to obtain and interpret likelihood-ratio, score, and Wald statistics from restricted and unrestricted multivariate regression fits.

Intermediate 12 min read

What You'll Learn

  • Set up a block test of predictor coefficients in a multivariate regression model.
  • Construct the hypothesis and error matrices from fitted models.
  • Derive the likelihood-ratio statistic from the profiled Gaussian likelihood.
  • Derive score and Wald statistics using the restricted and unrestricted fits.
  • Calculate all three statistics for a full block, a partial test, and a predictor contrast.
  • Distinguish the covariance estimate required by each statistic.

From a Regression Question to Three Test Statistics

A multivariate regression can ask whether a predictor affects any of several responses, or whether a group of predictors has an effect across the response vector. The same null hypothesis can be assessed from three different places in the fitting process: the likelihood comparison between restricted and unrestricted fits, the score at the restricted fit, and the distance from the unrestricted estimate to the null restriction. These yield the likelihood-ratio, score, and Wald statistics.

The calculations are related, but they are not interchangeable. The likelihood-ratio statistic compares residual determinants; the score statistic uses the restricted covariance estimate; and the Wald statistic uses the unrestricted covariance estimate. A small model makes both the similarities and differences visible.

Model and notation: Let \(Y\) be an \(n\times q\) response matrix and \(X\) an \(n\times p\) full-column-rank design matrix. Assume the rows of \(Y\) are independent with \(Y_i\sim N_q(x_i'B,\Sigma)\), where \(B\) is \(p\times q\) and \(\Sigma\) is positive definite. To test \(CB=0\), take \(C\) to be an \(r\times p\) full-row-rank matrix. Write \(E\) for the unrestricted residual cross-product matrix and \(E_0\) for the restricted one, and define \(H=E_0-E\). Assume both residual cross-product matrices are positive definite so that the Gaussian covariance maximum-likelihood estimates exist.

For full-rank \(X\), the unrestricted estimate is \(\hat B=(X'X)^{-1}X'Y\), with \(E=(Y-X\hat B)'(Y-X\hat B)\). For the restriction \(CB=0\), the additional residual sum of squares has the equivalent expression \(H=(C\hat B)'[C(X'X)^{-1}C']^{-1}(C\hat B)\), and \(E_0=E+H\). These identities let us calculate the test statistics directly from the fitted models. They also show why the test depends on the estimated restricted effect across all \(q\) responses, rather than on a separate univariate calculation for each response.

How the Three Statistics Arise

For a residual cross-product matrix \(S\), the Gaussian likelihood is maximized over the covariance matrix at \(\hat\Sigma=S/n\). Substituting that estimate into the likelihood gives a profile likelihood proportional to \(|S|^{-n/2}\). Thus the likelihood-ratio statistic is naturally expressed using the determinant ratio \(|E|/|E_0|\).

Likelihood-ratio statistic: Define \(\Lambda=|E|/|E_0|\), the ratio of the unrestricted to restricted residual determinants. The maximized restricted-to-unrestricted likelihood ratio is \(\Lambda^{n/2}\), so the conventional deviance statistic is \(-2\log(\Lambda^{n/2})=n\log(|E_0|/|E|)\).

The score statistic asks how strongly the likelihood points away from the restriction, evaluated at the restricted fit. After accounting for the unrestricted nuisance coefficients, let \(Z\) denote the tested design directions remaining after projection off the nuisance design. Then \(H=R_0'P_ZR_0\), where \(R_0\) is the restricted residual matrix and \(P_Z\) is the projection onto those tested directions. Under the null fit, the score for the tested coefficient matrix is \(Z'R_0\hat\Sigma_0^{-1}\), where \(\hat\Sigma_0=E_0/n\). Its information matrix is the Kronecker product of \(Z'Z\) and \(\hat\Sigma_0^{-1}\). The score quadratic form is therefore

$$ S_{\mathrm{score}}=\operatorname{tr}(\hat\Sigma_0^{-1}H) =n\,\operatorname{tr}(E_0^{-1}H). $$

This is the score test with nuisance parameters estimated under the restriction. In particular, the covariance estimate in this expression comes from the restricted fit, not the unrestricted one.

For the Wald statistic, use the unrestricted estimate of the tested coefficient matrix. In the residualized coordinates, its covariance is determined by the tested design cross-product and the unrestricted covariance estimate \(\hat\Sigma=E/n\). The Wald quadratic form is the fitted effect measured in that estimated covariance metric. Since the corresponding hypothesis cross-product is \(H\), it becomes

$$ S_{\mathrm{Wald}}=\operatorname{tr}(\hat\Sigma^{-1}H) =n\,\operatorname{tr}(E^{-1}H). $$

The factor \(n\) in both trace formulas is a consequence of using maximum-likelihood covariance estimates \(E_0/n\) and \(E/n\). Replacing these with residual covariance estimates divided by residual degrees of freedom changes the displayed statistics.

Worked Example: Testing Two Predictor Effects

Worked Example: Testing Two Predictor Effects

Consider six observations, two responses, an intercept, and two predictors. Let the predictor columns be \(x_1=(1,-1,0,0,0,0)'\) and \(x_2=(0,0,1,-1,0,0)'\), and let the design matrix have columns \(1,x_1,x_2\). The observed responses are specified by \(Y_1=10+2x_1+0.5x_2+c\) and \(Y_2=20-x_1+3x_2+d\), where \(c=(1,1,-1,-1,0,0)'\) and \(d=(0,0,0,0,1,-1)'\).

The vectors \(1,x_1,x_2\) are mutually orthogonal, with squared lengths \(6,2,2\). Both \(c\) and \(d\) are orthogonal to all three design columns, and \(c'd=0\). Therefore the unrestricted coefficient estimate and residual cross-product are

$$ \hat B= \begin{pmatrix} 10&20\\ 2&-1\\ 0.5&3 \end{pmatrix}, \qquad E= \begin{pmatrix} c'c&c'd\\ d'c&d'd \end{pmatrix} = \begin{pmatrix} 4&0\\ 0&2 \end{pmatrix}. $$

Test the joint null that both predictor rows are zero. Here \(C\) selects the last two rows of \(B\). Since their design columns have squared length 2, the restriction covariance factor is \(C(X'X)^{-1}C'=\tfrac12 I_2\). With \(D\) denoting the two estimated predictor rows, the hypothesis matrix is

$$ D= \begin{pmatrix} 2&-1\\ 0.5&3 \end{pmatrix}, \qquad H=D'(\tfrac12 I_2)^{-1}D =2D'D = \begin{pmatrix} 8.5&-1\\ -1&20 \end{pmatrix}. $$

Thus the restricted residual cross-product matrix is \(E_0=E+H=\begin{pmatrix}12.5&-1\\-1&22\end{pmatrix}\). Its determinant is \(12.5(22)-1=274\), while \(|E|=4(2)=8\). The likelihood ratio and deviance are

$$ \Lambda=\frac{8}{274}=\frac{4}{137}, \qquad -2\log(\Lambda^{n/2}) =6\log\left(\frac{137}{4}\right) \approx 21.20. $$

For the score statistic, \(E_0^{-1}=\frac1{274}\begin{pmatrix}22&1\\1&12.5\end{pmatrix}\). Taking the trace of the product gives

$$ \operatorname{tr}(E_0^{-1}H) =\frac{22(8.5)+1(-1)+1(-1)+12.5(20)}{274} =\frac{435}{274}, \qquad S_{\mathrm{score}}=6\frac{435}{274}\approx9.53. $$

For the Wald statistic, \(E^{-1}=\operatorname{diag}(1/4,1/2)\), so \(S_{\mathrm{Wald}}=6\{(1/4)(8.5)+(1/2)(20)\}=72.75\). The three statistics are quite different numerically in this small sample. The calculations are internally consistent: the LR uses the ratio of determinants, the score uses \(E_0^{-1}\), and the Wald statistic uses \(E^{-1}\).

Worked Example: Testing One Predictor While Retaining the Other

Worked Example: Testing One Predictor While Retaining the Other

Use the same fitted model, but now test \(H_0:\beta_1=0\), where \(\beta_1=(2,-1)\) is the coefficient row for \(x_1\). The coefficient of \(x_2\) remains unrestricted and acts as a nuisance coefficient. Since \(x_1\) is orthogonal to \(x_2\) and the intercept, the additional hypothesis matrix is

$$ H_1=2\beta_1'\beta_1 = \begin{pmatrix} 8&-4\\ -4&2 \end{pmatrix}, \qquad E_{01}=E+H_1= \begin{pmatrix} 12&-4\\ -4&4 \end{pmatrix}. $$

The determinants are \(|E|=8\) and \(|E_{01}|=12(4)-16=32\). Hence \(\Lambda_1=8/32=1/4\), and the LR deviance is \(6\log(4)\approx8.32\).

The inverse restricted matrix is \(E_{01}^{-1}=\frac1{32}\begin{pmatrix}4&4\\4&12\end{pmatrix}\). Therefore \(\operatorname{tr}(E_{01}^{-1}H_1)=\{4(8)+4(-4)+4(-4)+12(2)\}/32=3/4\), giving \(S_{\mathrm{score}}=6(3/4)=4.5\).

For the Wald statistic, \(\operatorname{tr}(E^{-1}H_1)=(1/4)(8)+(1/2)(2)=3\), so \(S_{\mathrm{Wald}}=18\). This is a partial test: the restricted model removes \(x_1\) but still fits the intercept and \(x_2\). Omitting \(x_2\) as well would answer a different scientific question and produce a different \(E_0\).

Worked Example: Testing a Predictor Contrast

Worked Example: Testing a Predictor Contrast

Still using the same model, test whether the sum of the two predictor coefficient rows is zero: \(H_0:\beta_1+\beta_2=0\). This is \(CB=0\) with \(C=(0,1,1)\). The estimated contrast is \((2.5,2)\), and \(C(X'X)^{-1}C'=1/2+1/2=1\). Consequently,

$$ H_2= \begin{pmatrix} 2.5\\2 \end{pmatrix} \begin{pmatrix} 2.5&2 \end{pmatrix} = \begin{pmatrix} 6.25&5\\ 5&4 \end{pmatrix}, \qquad E_{02}=E+H_2= \begin{pmatrix} 10.25&5\\ 5&6 \end{pmatrix}. $$

The determinant of \(E_{02}\) is \(10.25(6)-25=36.5\). Thus \(\Lambda_2=8/36.5=16/73\), and the LR deviance is \(6\log(73/16)\approx9.10\).

The inverse is \(E_{02}^{-1}=\frac1{36.5}\begin{pmatrix}6&-5\\-5&10.25\end{pmatrix}\). Hence \(\operatorname{tr}(E_{02}^{-1}H_2)=\{6(6.25)-5(5)-5(5)+10.25(4)\}/36.5=28.5/36.5=57/73\), and \(S_{\mathrm{score}}=342/73\approx4.68\).

Finally, \(\operatorname{tr}(E^{-1}H_2)=(1/4)(6.25)+(1/2)(4)=57/16\), so \(S_{\mathrm{Wald}}=6(57/16)=21.375\). The restriction is on the sum of the coefficient rows, not on either predictor separately; the contrast matrix specifies exactly which effect is being tested.

Reading the Results Correctly

These calculations demonstrate the roles of the fitted models. The unrestricted fit supplies \(\hat B\) and \(E\); the restricted fit supplies \(E_0\); their difference \(H=E_0-E\) records the residual variation attributable to imposing the tested restriction. The likelihood-ratio statistic summarizes the change through a determinant ratio, while score and Wald summarize it through covariance-weighted traces.

The statistics do not have the same finite-sample numerical value, and a statistic alone is not a p-value. Inference requires the appropriate reference distribution and its regularity conditions for the model and test being used. For the computations here, the key discipline is to state the restriction, fit the correct restricted model, and use the corresponding covariance estimate in each formula.

Common Mistakes

  • Using the wrong fit for the covariance. The score formula uses \(E_0/n\); the Wald formula uses \(E/n\). Swapping them reverses the defining distinction between the tests.
  • Testing the wrong set of coefficients. A partial test retains nuisance predictors in the restricted model. Removing them changes both \(E_0\) and the hypothesis being tested.
  • Replacing a determinant with a trace. The LR statistic depends on \(|E_0|/|E|\), not on \(\operatorname{tr}(E_0)/\operatorname{tr}(E)\).
  • Forgetting the sample-size scaling. With maximum-likelihood covariance estimates, \(S_{\mathrm{score}}\) and \(S_{\mathrm{Wald}}\) include the factor \(n\). Residual-degree-of-freedom covariance estimates require a different formula.
Key takeaway: For a specified restriction, construct \(H\) and \(E_0=E+H\) from the restricted and unrestricted fits. Then use \(n\log(|E_0|/|E|)\) for the LR deviance, \(n\operatorname{tr}(E_0^{-1}H)\) for the score statistic, and \(n\operatorname{tr}(E^{-1}H)\) for the Wald statistic.

Check Your Understanding

Use the definitions and worked calculations above to answer these questions.

  1. Which fitted model supplies the covariance estimate in the score statistic, and which supplies it in the Wald statistic?
  2. For a partial test of one predictor, why must the other predictors remain in the restricted model?
  3. In the first example, how is \(E_0\) obtained from \(E\) and \(H\)?
  4. Why does the likelihood-ratio statistic use a determinant ratio rather than a trace?
  5. What scientific restriction is represented by \(C=(0,1,1)\) in the third example?

Previous

Next

Back to Tutorials
Clinical Biostats