A single model may contain a number of linear equations. In such a model, it is often unrealistic to expect that the equation errors would be uncorrelated. A set of equations that has contemporaneous cross-equation error correlation (i.e. the error terms in the regression equations are corrlated) is called a seemingly unrelated regression (SUR) system. At first look, the equations seem unrelated, but the equations are related through the correlation in the errors.
The systemfit R package allows a user to specify multiple equations and fit them
in an SUR. After doing so, one can perform tests on coefficients across the
equations.

We will illustrate SUR using the hsb2 dataset, predicting
read and math with the overlapping sets of coefficients and then
comparing some coefficients across the two equations. We will first define our
model equations as R formulas.

First let’s install the systemfit R package which requires installation of three other packages: nnet, mgcv, and quantreg.

install.packages("nnet")
install.packages("mgcv")
install.packages("quantreg")
install.packages("systemfit")
install.packages("foreign")
install.packages("car")
install.packages("Rcpp")
library(foreign)
library(systemfit)

hsb2 <- read.dta("https://stats.idre.ucla.edu/stat/stata/notes/hsb2.dta")

r1 <- read~female + as.numeric(ses) + socst
r2 <- math~female + as.numeric(ses) + science

Once the equations have been defined, they can be passed in a list to the
systemfit command. A summary of the systemfit first shows a summary
of the system (where N = 400), then the separate equations, and then how the
residuals of the two equations are related. These are followed by the OLS fits
of the separate equations.

fitsur <- systemfit(list(readreg = r1, mathreg = r2), data=hsb2)
summary(fitsur)

systemfit results 
method: OLS 

         N  DF     SSR detRCov   OLS-R2 McElroy-R2
system 400 392 22835.2 3227.86 0.405103   0.342707

          N  DF     SSR     MSE   RMSE       R2   Adj R2
readreg 200 196 12550.9 64.0351 8.0022 0.400037 0.390854
mathreg 200 196 10284.4 52.4712 7.2437 0.411171 0.402158

The covariance matrix of the residuals
        readreg mathreg
readreg 64.0351 11.4952
mathreg 11.4952 52.4712

The correlations of the residuals
        readreg mathreg
readreg 1.00000 0.19831
mathreg 0.19831 1.00000

OLS estimates for 'readreg' (equation 1)
Model Formula: read ~ female + as.numeric(ses) + socst

                  Estimate Std. Error  t value   Pr(>|t|)    
(Intercept)     20.6824980  2.9789550  6.94287 5.5019e-11 ***
femalefemale    -1.5111280  1.1510793 -1.31279    0.19079    
as.numeric(ses)  1.2183658  0.8399004  1.45061    0.14849    
socst            0.5699327  0.0562967 10.12373 < 2.22e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.002195 on 196 degrees of freedom 
Number of observations: 200 Degrees of Freedom: 196 
SSR: 12550.883066 MSE: 64.035118 Root MSE: 8.002195 
Multiple R-Squared: 0.400037 Adjusted R-Squared: 0.390854 

OLS estimates for 'mathreg' (equation 2) 
Model Formula: math ~ female + as.numeric(ses) + science 

                 Estimate Std. Error  t value   Pr(>|t|)    
(Intercept)     19.305181   2.998047  6.43925 9.0557e-10 ***
femalefemale     1.160903   1.041641  1.11449   0.266432    
as.numeric(ses)  1.399639   0.742390  1.88531   0.060867 .  
science          0.575330   0.054328 10.58993 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.243704 on 196 degrees of freedom
Number of observations: 200 Degrees of Freedom: 196 
SSR: 10284.364144 MSE: 52.471246 Root MSE: 7.243704 
Multiple R-Squared: 0.411171 Adjusted R-Squared: 0.402158 

We may be interested in comparing the effect of female on read,
controlling for ses and socst, to the effect of female on
math, controlling for ses and science. For this, we will
use the linear.hypothesis command from the car package. To do
this, we create a “restriction” on the model system. We will force the
coefficient of female to be the same in both equations and then compare
such a system fit to the one seen when the coefficients are not equal.

library(car)
restriction <- "readreg_femalefemale- mathreg_femalefemale"
linearHypothesis(fitsur, restriction, test = "Chisq")

Linear hypothesis test (Chi^2 statistic of a Wald test)

Hypothesis:
readreg_femalefemale - mathreg_femalefemale = 0

Model 1: restricted model
Model 2: fitsur

 Res.Df Df Chisq Pr(>Chisq) 
1 393 
2 392 1 2.9626 0.08521 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1