# How to test significance in shift of a variable taken an other variable into the model (suest) in R?

I want to test in R whether an additional variable $x_2$ taken into a model influences an other variable significantly.

(Note: This a revised repost of one in Stack Overflow where there was no answer and which I thus have deleted.)

1. $$y = \beta_1x_1 + \epsilon$$
2. $$y = \beta_1x_1 + \beta_2x_2 + \epsilon$$

In Stata I achieved this with suest (see example below). I am looking for an appropriate method for R and @carl_pch's comment to this question suggested this method in R which appears to be related. Adapting it to my problem I got quite different results, though. As my statistical knowledge is just basic I'm not sure if it's the same thing.

Stata's help suest says:

suest combines the estimation results -- parameter estimates and associated (co)variance matrices -- stored under namelist into one parameter vector and simultaneous (co)variance matrix of the sandwich/robust type. This (co)variance matrix is appropriate even if the estimates were obtained on the same or on overlapping data.

How can we test in R whether the shift of $\beta_1x_1$ in the example, taking $x_2$ into the model, is statistically significant mimicking Stata's suest somehow?

The best now would be I show you in conclusion what I've done in Stata and my probably failed attempt in R:

Stata example

use http://www.stata-press.com/data/r8/regress.dta, clear

reg y x1
eststo fit1

reg y x1 x2
eststo fit2


$\beta_1x_1$ of model 1 is $7.812835$ whereas $\beta_1x_1$ of model 2 is $8.566015$. To test now we can use the postestimation command suest, which is related to "Seemingly unrelated estimation", and later apply a Wald test whether the shift in $\beta_1x_1$ is statistically significant:

suest fit1 fit2


This yields following output:

Simultaneous results for fit1, fit2

Number of obs     =        148

------------------------------------------------------------------------------
|               Robust
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
fit1_mean    |
x1 |   7.812835   .8834308     8.84   0.000     6.081343    9.544328
_cons |  -2.257346   2.532697    -0.89   0.373     -7.22134    2.706649
-------------+----------------------------------------------------------------
fit1_lnvar   |
_cons |     3.0332   .1356275    22.36   0.000     2.767375    3.299025
-------------+----------------------------------------------------------------
fit2_mean    |
x1 |   8.566015   .9999644     8.57   0.000      6.60612    10.52591
x2 |   1.056696   1.115586     0.95   0.344    -1.129813    3.243204
_cons |  -4.213927    2.76533    -1.52   0.128    -9.633873     1.20602
-------------+----------------------------------------------------------------
fit2_lnvar   |
_cons |   3.034356   .1345578    22.55   0.000     2.770627    3.298084
------------------------------------------------------------------------------


Wald test:

test _b[fit1_mean:x1] = _b[fit2_mean:x1]


The test yields:

 ( 1)  [fit1_mean]x1 - [fit2_mean]x1 = 0

chi2(  1) =    0.89
Prob > chi2 =    0.3465


Thus, in this case the shift is not significant.

R example

library(readstata13)

fit1 <- y ~ x1
fit2 <- y ~ x1 + x2

library(systemfit)
fitsur <- systemfit(list(readreg = fit1, mathreg = fit2), data=df1)


Yielding the output:

> summary(fitsur)

systemfit results
method: OLS

N  DF     SSR detRCov   OLS-R2 McElroy-R2
system 296 291 6045.68 2.46044 0.381443  -0.236551

N  DF     SSR     MSE    RMSE       R2   Adj R2
readreg 148 146 3031.48 20.7636 4.55671 0.379675 0.375426
mathreg 148 145 3014.20 20.7876 4.55934 0.383211 0.374703

The covariance matrix of the residuals
mathreg 20.7163 20.7876

The correlations of the residuals
mathreg 0.997146 1.000000

OLS estimates for 'readreg' (equation 1)
Model Formula: y ~ x1

Estimate Std. Error  t value Pr(>|t|)
(Intercept) -2.257346   2.519743 -0.89586   0.3718
x1           7.812835   0.826488  9.45306   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.556705 on 146 degrees of freedom
Number of observations: 148 Degrees of Freedom: 146
SSR: 3031.480349 MSE: 20.763564 Root MSE: 4.556705
Multiple R-Squared: 0.379675 Adjusted R-Squared: 0.375426

OLS estimates for 'mathreg' (equation 2)
Model Formula: y ~ x1 + x2

Estimate Std. Error  t value   Pr(>|t|)
(Intercept) -4.21393    3.31082 -1.27277    0.20514
x1           8.56601    1.16888  7.32842 1.4993e-11 ***
x2           1.05670    1.15897  0.91176    0.36341
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.55934 on 145 degrees of freedom
Number of observations: 148 Degrees of Freedom: 145
SSR: 3014.199626 MSE: 20.787584 Root MSE: 4.55934
Multiple R-Squared: 0.383211 Adjusted R-Squared: 0.374703


Hypothesis test

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


Yielding

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

Hypothesis:

Model 1: restricted model
Model 2: fitsur

Res.Df Df  Chisq Pr(>Chisq)
1    292
2    291  1 0.2768     0.5988


I was able to replicate Stata's suest using geepack in R. I based my replication on:

Zhuan Pei, Jörn-Steffen Pischke & Hannes Schwandt (2019) Poorly Measured Confounders are More Useful on the Left than on the Right, Journal of Business & Economic Statistics, 37:2, 205-216

The data and replication code for Stata is available under the article's supplemental material here

library(foreign)
library(geepack)
library(sandwich)

### GET SUBSET OF DATA USED BY AUTHORS IN REPLICATION FILE ###

lm.fit <- lm(lwage76 ~ ed76 + exp + exp2 + black + reg76r + smsa76r + momed + libcrd14 + height73, data = nls)
sum(!is.element(rownames(nls), names(lm.fit$$residuals))) nls <- nls[is.element(rownames(nls), names(lm.fit$$residuals)),]
#nls <- nls[!is.na(nls$$iq) & !is.na(nls$$kww),] # THIS IS THE DATA USED IN THE PUBLISHED ARTICLE


We want to estimate the effect of years of education on log weekly earnings and want to compare the coefficient on education when adding mother's education as a covariate, which is a potential confounder

f01 <- lwage76 ~ ed76 + exp + exp2 + black + smsa76r + smsa66r + reg76r + reg661 + reg662 + reg663 + reg664 + reg665 + reg667 + reg668 + reg669
f02 <- lwage76 ~ ed76 + exp + exp2 + black + smsa76r + smsa66r + reg76r + reg661 + reg662 + reg663 + reg664 + reg665 + reg667 + reg668 + reg669 + momed

gee.fit01 <- geese(f01, data = nls, corstr = "independence")
gee.fit02 <- geese(f02, data = nls, corstr = "independence")

comp0102 <- compCoef(gee.fit01, gee.fit02)
t0102 <- comp0102$$delta[2]/sqrt(comp0102$$variance[2,2])


Squaring the t-statistic and comparing it to $$\chi^2$$-distribution

> t0102^2
ed76
4.048876
> pchisq(t0102^2, 1, lower.tail = F)
ed76
0.04420076


yields the same result as suest in Stata (check the relication data in the link above).The test statistic for comparing the coefficient ($$\beta$$) between model 1 and 2 is

$$t = \frac{\beta_1 - \beta_2}{\sqrt{\text{Var}\left(\beta_1 - \beta_2\right)}}$$

The tricky part of this equation is that

$$\text{Var}\left(\beta_1 - \beta_2\right) = \text{Var}\left(\beta_1\right) + \text{Var}\left(\beta_2\right) - 2\text{Cov}\left(\beta_1, \beta_2\right)$$

If we use systemfit and linearHypothesis

library(systemfit)
library(sandwich)
library(car)

system <- list(f01 = f01, f02 = f02)
sys.fit0102 <- systemfit(system, data = nls)
linearHypothesis(sys.fit0102, "f01_ed76-f02_ed76=0", vcov. = function(x) vcovHC(x, type = "HC1"), test = "Chisq")


we assume that $$\text{Cov}\left(\beta_1, \beta_2\right) = 0$$. We will get the same result if we fit two lm models and calculate the statistic "by hand":

lm.fit01 <- lm(f01, data = nls)
lm.fit02 <- lm(f02, data = nls)

d <- lm.fit01$$coef[2]-lm.fit02$$coef[2]
se <- sqrt(vcovHC(lm.fit01, type = "HC1")[2,2]+vcovHC(lm.fit02, type = "HC1")[2,2])
> (d/se)^2
ed76
0.1637189
> pchisq((d/se)^2, 1, lower.tail = F)
ed76
0.6857554