Seemingly Unrelated Regression (SUR) in R: systemfit Package

Seemingly unrelated regression (SUR) fits a system of linear equations jointly when their error terms are correlated within the same observation. In R, the systemfit package estimates SUR with feasible generalised least squares in one line. SUR gives you smaller standard errors than running each equation through separate lm() calls, and it lets you test hypotheses that span both equations.

When should you use seemingly unrelated regression?

You reach for SUR when you have two or more regression equations whose predictors may differ, but whose errors share a common shock per observation. Households, firms, and plants react to the same unobserved jolt in the same year, row, or region, so the residuals end up correlated even though the equations look independent. systemfit() turns that correlation into an efficiency gain. Here is the whole workflow on iris, with two flower-width equations sharing petal predictors.

The off-diagonal entry, 0.32, is the cross-equation residual correlation. That number is the licence to use SUR. When two observations land above their predicted Sepal.Length, they also tend to land above their predicted Sepal.Width. Separate lm() calls would ignore that pattern; SUR uses it to sharpen both equations at once. If the off-diagonal had been near zero, SUR and OLS would return essentially the same numbers and you could skip the extra machinery.

How do you specify a system of equations for systemfit()?

Every SUR model is a list of standard R formulas. You write each one exactly as you would for lm(), then bundle them into a named list. systemfit() stitches them into a single block-diagonal design matrix behind the scenes, but the formula API stays familiar.

Every coefficient is prefixed by its equation label. That keeps eq_sepal_Petal.Length clearly distinct from eq_petal_Petal.Length even though they come from the same predictor column. You will need the exact prefixed names later when testing cross-equation restrictions, so name your equations something short and memorable.

How much efficiency does SUR gain over separate OLS?

SUR is only worth the extra setup when it produces tighter standard errors than running each equation separately. Let us measure the gain directly by fitting the same iris system two ways, equation-by-equation OLS and joint SUR, then comparing their coefficient standard errors side by side.

The reduction is exactly zero on this particular iris system because both equations share an identical set of regressors. That is the first of the two cases where SUR mathematically collapses to OLS. The cross-equation correlation exists (0.32, from the earlier block) but there is nothing left to squeeze out of it once the regressors match. To see a non-zero gain, you would need equations with different predictors, as the earlier "Species in equation 2 only" example had.

If you are comfortable skipping the math, jump to the next section, the practical code above is all you need. Otherwise, the SUR estimator is:

$$\hat{\beta}_{\text{SUR}} = \left(X^\top (\hat{\Sigma}^{-1} \otimes I_n) X\right)^{-1} X^\top (\hat{\Sigma}^{-1} \otimes I_n) y$$

Where:

• $\hat{\beta}_{\text{SUR}}$ is the stacked coefficient vector across all equations
• $\hat{\Sigma}$ is the estimated cross-equation residual covariance matrix
• $\otimes$ is the Kronecker product, which expands $\hat{\Sigma}$ to act on every observation
• $n$ is the number of observations per equation
• $X$ is the block-diagonal matrix of regressors

When $\hat{\Sigma}$ is diagonal (no cross-equation correlation) or every block of $X$ is identical, the sandwich collapses and the formula reduces to the usual OLS normal equations.

How do you test hypotheses across SUR equations?

A large reason to fit the two equations jointly, rather than separately, is to test whether a predictor has the same effect in both. The car::linearHypothesis() function reads a systemfit object directly, so you write the restriction as a string and R does the rest.

The chi-square statistic is huge (121.17) and the p-value is effectively zero, so we reject the null that Petal.Length has the same slope in both equations. That makes biological sense: petal length positively predicts sepal length but negatively predicts sepal width. The cross-equation test gives you one clean p-value for that difference, which is exactly the thing two separate lm() fits cannot do.

What assumptions and diagnostics should you check?

SUR inherits every OLS assumption and adds one more: errors may be correlated across equations within the same observation, but must be independent across observations. If you suspect endogenous regressors, SUR will not save you, the estimates are biased for the same reasons OLS is.

1. Errors are independent across observations (no serial correlation, no clustering).
2. Errors within an observation may be correlated across equations.
3. Regressors are exogenous and not correlated with the error term.
4. Each equation is correctly specified, with a linear functional form and no omitted confounders.

The first two assumptions you check with the residual correlation matrix; the other two you check with the same diagnostics you would run on any lm() fit (residual plots, VIF, reset tests).

An off-diagonal of 0.32 is a middling correlation, meaningful but not dominant. That is enough to make SUR worth running when your regressors differ across equations. If this value had come back near zero (say |r| < 0.05), you could safely fall back on separate lm() calls and lose nothing. If it came back near one, you should also check whether the two outcomes are really measuring different things or just renamed copies of each other.

Practice Exercises

Exercise 1: Fit a SUR system on mtcars with different regressors per equation

Build a SUR system on mtcars where equation 1 is mpg ~ wt + hp and equation 2 is qsec ~ wt + hp + cyl. Save the fit to my_fit and print the residual correlation between the two equations.

Exercise 2: Test whether wt has the same effect in both mtcars equations

Using the my_fit object from Exercise 1, test $H_0:\ \beta_{\text{wt, mpg\_eq}} = \beta_{\text{wt, qsec\_eq}}$. Save the test to my_test and print the chi-square statistic plus p-value.

Complete Example

Here is the full SUR workflow from start to finish on a fresh iris fit, so you can copy it as a template.

Residual correlation of 0.38 says there is real cross-equation dependence worth exploiting. The chi-square of 70.5 tells you the petal-length slopes differ strongly between the two equations. Together, those two numbers justify joint estimation and answer a scientific question that separate lm() calls could not reach in one step.

Summary

References

1. Zellner, A. (1962). An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias. Journal of the American Statistical Association, 57(298), 348-368. Link
2. Henningsen, A., and Hamann, J. D. (2007). systemfit: A Package for Estimating Systems of Simultaneous Equations in R. Journal of Statistical Software, 23(4). Link
3. systemfit CRAN vignette. Link
4. UCLA OARC Statistical Consulting. How can I perform seemingly unrelated regression in R? Link
5. Greene, W. H. (2018). Econometric Analysis (8th ed.), Chapter 10. Pearson.
6. Wikipedia. Seemingly unrelated regressions. Link

Continue Learning

• Multiple Regression in R, the parent post covering single-equation multiple regression, the natural starting point before moving to systems.
• Linear Regression Assumptions in R, every assumption SUR inherits from OLS, plus how to diagnose violations.
• Regression Diagnostics in R, residual plots, leverage, influence measures that apply equation-by-equation inside a SUR fit.