SUEST vs Stacking Regressions

[Marie Gri]

I have been running a test using suest, but recently realized that I get a very different p-value for the same test when I stack the regressions. In my specific context (the effect for a policy by race, gender, SES and ethnicity), it makes a big difference. The p-value with suest is 0.03 and stacking the regressions 0.22 or 0.34 depending if I only stack the demographic groups or stack all the triple interaction variables. I have clustering in all my regressions at the household level. That suest provides an asymptotic test vs finite-sample F-test for the stacked regression does not seem to be driving the difference in my setting (I have a few thousand observations in total).

I do not have a great reproduceable example of the problem (unfortunately I have not figured out yet what is causing this sizable difference), but the below nevertheless highlights that the p-values across the specifications are not (exactly) the same in contrast to what seems to be indicated in the stata base reference manual pg2975-2976 for suest.

sysuse auto, clear


summ mpg, det 
gen highmpg=mpg>`r(p50)' 


reg price foreign 
estimates store m1 
reg price highmpg 
estimates store m2 
suest m1 m2 
test [m1_mean]foreign=[m2_mean]highmpg //0.0382 

//calculating the finite-sample F-statistic 
local chi2=r(chi2)
local q=r(df)
local N=e(N)
local k=e(rank)
local F=`chi2'/`q'
dis "F: " `F', "p: " Ftail(`q',`N'-`k',`F') //0.042 

*stacking the regression 
gen id=_n 
gen zero=0 
gen one=1 
gen two=2 

stack id price foreign zero one   id price zero highmpg two, into(id y z1 z2 model)
gen model2=model==2 
reg y model2 z1 z2, vce(cluster id)
test _b[z1]=_b[z2] //0.0438 F-test

Any advice as to what suest does differently relative to the stacked regression would be most welcome! I am trying to figure out which method I should rely on for my test.

[Andrew Musau]

The difference usually arises from the distribution used for the test statistic. regress reports test statistics based on the Student's t-distribution, whereas suest reports z-statistics. The chi-square distribution is the limiting distribution of the Wald statistic as the sample size becomes large. Likewise, the z-statistic is the large-sample (asymptotic) limit of the t-statistic. Therefore, you may see differences when the sample size is small—for example, when you have very few observations. When clustering, what matters is the number of clusters, because the asymptotics are based on the number of clusters rather than the number of observations. With a small number of clusters, cluster-robust standard errors may be unreliable, and small-sample corrections may be needed. Otherwise, if you use an estimator that reports z-statistics (such as glm), then the stacking method and suest will yield identical results because both rely on asymptotic distributions.

sysuse auto, clear
summ mpg, det 
gen highmpg=mpg>`r(p50)' 

reg price foreign 
estimates store m1 
reg price highmpg 
estimates store m2 
suest m1 m2 
test [m1_mean]foreign=[m2_mean]highmpg 

*stacking the regression 
gen id=_n 
gen zero=0 
gen one=1 
gen two=2 

stack id price foreign zero one   id price zero highmpg two, into(id y z1 z2 model)
gen model2=model==2 
glm y model2 z1 z2, vce(cluster id)
test _b[z1]=_b[z2]

Res.:

. suest m1 m2 

Simultaneous results for m1, m2                             Number of obs = 74

------------------------------------------------------------------------------
             |               Robust
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
m1_mean      |
     foreign |   312.2587   696.9581     0.45   0.654    -1053.754    1678.271
       _cons |   6072.423   428.2447    14.18   0.000     5233.079    6911.767
-------------+----------------------------------------------------------------
m1_lnvar     |
       _cons |    15.9902   .2260545    70.74   0.000     15.54714    16.43325
-------------+----------------------------------------------------------------
m2_mean      |
     highmpg |   -1587.01   654.7856    -2.42   0.015    -2870.366    -303.654
       _cons |   6937.316   525.7846    13.19   0.000     5906.797    7967.835
-------------+----------------------------------------------------------------
m2_lnvar     |
       _cons |   15.91643   .2317151    68.69   0.000     15.46228    16.37058
------------------------------------------------------------------------------

. test [m1_mean]foreign=[m2_mean]highmpg 

 ( 1)  [m1_mean]foreign - [m2_mean]highmpg = 0

           chi2(  1) =    4.29
         Prob > chi2 =    0.0382

. 
. *stacking the regression 
. gen id=_n 

. gen zero=0 

. gen one=1 

. gen two=2 

. 
. stack id price foreign zero one   id price zero highmpg two, into(id y z1 z2 model)
Warning: data in memory will be lost.
         Press any key to continue, Ctrl-Break to abort.

. gen model2=model==2 

. glm y model2 z1 z2, vce(cluster id)

Iteration 0:  Log pseudolikelihood = -1388.5709  

Generalized linear models                         Number of obs   =        148
Optimization     : ML                             Residual df     =        144
                                                  Scale parameter =    8486551
Deviance         =   1222063317                   (1/df) Deviance =    8486551
Pearson          =   1222063317                   (1/df) Pearson  =    8486551

Variance function: V(u) = 1                       [Gaussian]
Link function    : g(u) = u                       [Identity]

                                                  AIC             =   18.81853
Log pseudolikelihood = -1388.570876               BIC             =   1.22e+09

                                    (Std. err. adjusted for 74 clusters in id)
------------------------------------------------------------------------------
             |               Robust
           y | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      model2 |   864.8927   374.6168     2.31   0.021     130.6574    1599.128
          z1 |   312.2587   696.9581     0.45   0.654    -1053.754    1678.271
          z2 |   -1587.01   654.7856    -2.42   0.015    -2870.366    -303.654
       _cons |   6072.423   428.2447    14.18   0.000     5233.079    6911.767
------------------------------------------------------------------------------

. test _b[z1]=_b[z2]

 ( 1)  [y]z1 - [y]z2 = 0

           chi2(  1) =    4.29
         Prob > chi2 =    0.0382

[Marie Gri]

Thanks for the reply! I have 2.8k clusters, so that was not the issue.

However, I actually have figured out the issue. SUEST seems to not be able to cope with variables that carry the same name across specifications but are different in each specification. Maybe worth flagging somewhere in the helpfile of suest?

rename foreign variable 
reg price variable 
estimates store m1 
rename variable foreign 
rename highmpg variable 
reg price variable 
estimates store m2 
suest m1 m2 
test [m1_mean]variable=[m2_mean]variable

This gives a p-value of <0.0001, in contrast to the original (and correct) p-value of 0.0382.

[Andrew Musau]

SUEST seems to not be able to cope with variables that carry the same name across specifications but are different in each specification.
.

This does not make much sense when you think about it. The only way Stata can identify whether a variable is the same or different is through its name. When calculating the Wald statistic, the formula is

\[W = \frac{(\hat{\beta}_1 - \hat{\beta}_2)^2}
{\hat{\sigma}_{11} + \hat{\sigma}_{22} - 2\hat{\sigma}_{12}}\]

The only component of the Wald test that changes when the underlying data differ is the cross-equation covariance term, \(\hat{\sigma}_{12}\). So when you rename the variable the second time, you have variable= highmpg and you lose the fact that variable= foreign in the first model refers to a different data vector. As a result, the cross-equation covariance changes, and therefore the Wald statistic and its significance also change. Note that test requires all variables in the combined models to be in their correct original form in order to properly calculate the Wald test.

[Marie Gri]

Great, thanks for confirming that the variables need to be in their original state!

Yes, agreed that the above example does not make much sense, it is only there for illustrative purposes.
My real-life setting made more sense: I was running the specifications in a loop, and I did not care for keeping a new variable that I only needed for running that specification. So I was generating it under the same name and dropping it after running the specification in each loop. It would have been easy to generate the variables under different names, but not knowing this would be a concern that was not what I was doing. I was thinking that stata had stored all the info it needed with eststo and was not going to go back to the variables.

Thanks again, I really appreciate the responses!