1
$\begingroup$

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)
df1 <- read.dta13("https://www.stata-press.com/data/r8/regress.dta")

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
readreg mathreg
readreg 20.7636 20.7163
mathreg 20.7163 20.7876

The correlations of the residuals
readreg  mathreg
readreg 1.000000 0.997146
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)
restriction <- "readreg_x1- mathreg_x1"
test <- linearHypothesis(fitsur, restriction, test = "Chisq")

Yielding

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

Hypothesis:
readreg_x1 - mathreg_x1 = 0

Model 1: restricted model
Model 2: fitsur

  Res.Df Df  Chisq Pr(>Chisq)
1    292                     
2    291  1 0.2768     0.5988
$\endgroup$
1
$\begingroup$

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)

nls <- read.dta("nls_saveold.dta")

### 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 
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.