# Multivariate regression - test difference between coefficients

I am trying to figure out how to test the difference in slopes for the same explanatory variable in a multivariate linear regression. What type of test should be used and how can I perform this using R?

As an example, say I would use Petal.Length and Petal.Width to explain Sepal.Length and Sepal.Width from the Iris dataset.

First, fit the model and print the summary:

mvm <- lm(cbind(Sepal.Length, Sepal.Width) ~ Petal.Length + Petal.Width, data = iris)
summary(mvm)


This gives the following result:

Response Sepal.Length :

Call:
lm(formula = Sepal.Length ~ Petal.Length + Petal.Width, data = iris)

Residuals:
Min       1Q   Median       3Q      Max
-1.18534 -0.29838 -0.02763  0.28925  1.02320

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   4.19058    0.09705  43.181  < 2e-16 ***
Petal.Length  0.54178    0.06928   7.820 9.41e-13 ***
Petal.Width  -0.31955    0.16045  -1.992   0.0483 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4031 on 147 degrees of freedom
Multiple R-squared:  0.7663,    Adjusted R-squared:  0.7631
F-statistic:   241 on 2 and 147 DF,  p-value: < 2.2e-16

Response Sepal.Width :

Call:
lm(formula = Sepal.Width ~ Petal.Length + Petal.Width, data = iris)

Residuals:
Min       1Q   Median       3Q      Max
-1.06198 -0.23389  0.01982  0.20580  1.13488

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.58705    0.09373  38.272  < 2e-16 ***
Petal.Length -0.25714    0.06691  -3.843  0.00018 ***
Petal.Width   0.36404    0.15496   2.349  0.02014 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3893 on 147 degrees of freedom
Multiple R-squared:  0.2131,    Adjusted R-squared:  0.2024
F-statistic:  19.9 on 2 and 147 DF,  p-value: 2.238e-08


Now, just by looking at the estimated coefficients I can see that the coefficient for Petal.Length is quite different for Sepal.Length as the response vs Sepal.Width as the response and they even have the opposite sign. But how would I test the difference formally? What is the test called and how can I perform this using R?

I tried searching for the answer but couldn't find anything so any help is much appreciated!

Edit:

I found this article from UCLA https://stats.oarc.ucla.edu/stata/dae/multivariate-regression-analysis/ where they perform a multivariate regression using STATA. In the article they talk about performing the test I'm looking for. The following quote is from the article:

"...If you ran a separate OLS regression for each outcome variable, you would get exactly the same coefficients, standard errors, t- and p-values, and confidence intervals as shown above. So why conduct a multivariate regression? As we mentioned earlier, one of the advantages of using mvreg is that you can conduct tests of the coefficients across the different outcome variables."

After reading the article my understanding is that there should be a way of performing the test I want. I am not a user of STATA, but I was hoping for a similar solution using R. Is there a similar function/test in R that they use in the article (the one called "test")?

You need a multivariate method for a multivariate problem, i.e. not multiple regression. I would use a structural equation model approach where you first estimate the paths from Petal.Length to Sepal.Length and Sepal.Width freely, then constrain the two paths to equality and finally compare the two models with a likelihood ratio test.

Using the lavaan package, the unconstrained model can be defined and fitted like this:

library(lavaan)

unconstrained_model <- '
# regressions
Sepal.Length ~ Petal.Length + Petal.Width
Sepal.Width ~ Petal.Length + Petal.Width

# estimate (residual) covariances freely
Petal.Length ~~ Petal.Width
Sepal.Length ~~ Sepal.Width
'

# fit the model
fit_unconstrained <- sem(unconstrained_model, data=iris)
summary(fit_unconstrained)

# visualize the model
m <- matrix(nrow=2, col=3)
m[1,]<-c(3,0,1)
m[2,]<-c(0,4,2)
semPaths(fit_unconstrained, what="path", "est", edge.label.cex=1.2, style="lisrel",
sizeMan=10, label.norm="0" ,layout=m)


The parameter estimates from this model look like this:

Regressions:
Estimate  Std.Err  z-value  P(>|z|)
Sepal.Length ~
Petal.Length      0.542    0.069    7.899    0.000
Petal.Width      -0.320    0.159   -2.012    0.044
Sepal.Width ~
Petal.Length     -0.257    0.066   -3.882    0.000
Petal.Width       0.364    0.153    2.373    0.018

Covariances:
Estimate  Std.Err  z-value  P(>|z|)
Petal.Length ~~
Petal.Width       1.287    0.151    8.495    0.000
.Sepal.Length ~~
.Sepal.Width       0.097    0.015    6.518    0.000

Variances:
Estimate  Std.Err  z-value  P(>|z|)
.Sepal.Length      0.159    0.018    8.660    0.000
.Sepal.Width       0.148    0.017    8.660    0.000
Petal.Length      3.096    0.357    8.660    0.000
Petal.Width       0.577    0.067    8.660    0.000


The graphical version of the model can be viewed here.

It can be seen that the estimated partial regression coefficients and their standard errors closely match the ones you got in your two lm regressions. Note that there are also estimates for the variances and covariances of the (residual) variables. Intercepts could be added to the model easily but they are not required.

Let's do the constrained model then:

constrained_model <- '
# regressions with labels a and b added
Sepal.Length ~ a*Petal.Length + Petal.Width
Sepal.Width ~ b*Petal.Length + Petal.Width

# estimate (residual) covariances freely
Petal.Length ~~ Petal.Width
Sepal.Length ~~ Sepal.Width

# constrain paths
a == b
'

# fit the model
fit_constrained <- sem(constrained_model, data=iris)
summary(fit_constrained)


The parameter estimates look like this:

Regressions:
Estimate  Std.Err  z-value  P(>|z|)
Sepal.Length ~
Petl.Lngth (a)    0.105    0.061    1.727    0.084
Petal.Wdth        0.655    0.144    4.549    0.000
Sepal.Width ~
Petl.Lngth (b)    0.105    0.061    1.727    0.084
Petal.Wdth       -0.443    0.143   -3.102    0.002

Covariances:
Estimate  Std.Err  z-value  P(>|z|)
Petal.Length ~~
Petal.Width       1.287    0.151    8.495    0.000
.Sepal.Length ~~
.Sepal.Width       0.061    0.016    3.745    0.000

Variances:
Estimate  Std.Err  z-value  P(>|z|)
.Sepal.Length      0.202    0.023    8.660    0.000
.Sepal.Width       0.178    0.021    8.660    0.000
Petal.Length      3.096    0.357    8.660    0.000
Petal.Width       0.577    0.067    8.660    0.000


It can be seen that the paths of interest now have the same value (0.105), in accordance with the constraint statement. Several other estimates have changed to accommodate the constraint.

In a final step we do a likelihood ratio test of the constrained and unconstrained models. We can do this (in large samples at least) because the models are nested, i.e. the constrained model is a special case of the unconstrained model. The -2*log of the difference in the models' likelihoods is asymptotically chi-squared distributed under the null, with degrees of freedom equal to the difference in degrees of freedom between the models (1 in this case, given the single constraint). The null hypothesis here is that the more parsimonious constrained model fits as well as the less parsimonious unconstrained model. The anova function can be used:

anova(fit_constrained, fit_unconstrained)


The output is (in part):

                  Df diff Pr(>Chisq)
fit_unconstrained
fit_constrained         1  < 2.2e-16 ***


The p-value is very small, suggesting that the null hypothesis should be rejected. So, in the final analysis, the (partial) regressions of Sepal.Width and Sepal_Length on Petal.Length probably have different slopes in the population; the slopes cannot be constrained to equality without a significant loss of model fit.

• Thanks for your answer! I think I'll have to read it a couple of times before I will completely understand it. I found an article yesterday while searching for an answer. I have updated my question with it. Commented Feb 4, 2022 at 10:49
• I thought one of the advantages of multivariate regression was that the model estimates covariances between all variables, which I assume would be necessary in order to calculate a correct SE. I was hoping for a solution where only a simply function from some package would be required in order to perform my test :P Commented Feb 4, 2022 at 11:00