Testing Equality of Coefficients across different regressions.

[Tom Ford]

Hello,
I am regressing my treatment on two separate dependent variables.
I would like to formally test if the effect of my treatment on the first dependent variable is significantly different from the effect on the second dependent variable.
Does anyone know how I can formally test this? Ideally, I would use Bonferroni adjusted p-values.

Take the following example. I am estimating the impact of Brands on the weight and length of the cars.

sysuse auto,clear
gen brand = word(make, 1)
keep if inlist(brand, "Audi", "Fiat", "Olds", "Pont.")
encode brand, generate(brand_en)

* estimate regression of interest
reg length ib1.brand_en
eststo rlr
reg weight ib1.brand_en
eststo rwe

The first step of my analysis is to examine if the effect of going from "Audi" to "Fiat" is significantly different from zero. And I can do this by simply using the test command after my regressions like so.

reg length ibn.brand_en,nocons
eststo rlr
test 1.brand_en=2.brand_en, mtest(bonf)

reg weight ibn.brand_en,nocons
eststo rwe
test 1.brand_en=2.brand_en, mtest(bonf)

Now, I would like to examine if my treatment (going from Audi to Fiat) has a different effect on weight and length.
I tried the following

* joining the two models.
suest rlr rwe
* test difference in means
test [rlr_mean]2.brand_en-[rwe_mean]2.brand_en = 0

Does anyone know if this is the correct approach? Is there another way to test if my effect of brand is different on weight and length? Also, because I could potentially make 4 comparisons (Audi-Fiat,Olds-Fiat,Olds-Pont,Audi-pont), I would like my results to account for a Bonferroni adjustment for 4 hypothesis tests.


thanks a lot in advance for your help

[George Ford]

that looks correct.

help test.

you can mtest and make bonferroni (among other) adjustments.

[Bruce Weaver]

I looked at this earlier today, but did not have time to respond until now. If I understood what you want to do, why not use -mvreg- followed by -test-, or even better, followed by -lincom-? For example:

clear *
sysuse auto
gen brand = word(make, 1)
keep if inlist(brand, "Audi", "Fiat", "Olds", "Pont.")
encode brand, generate(brand_en)

* Estimate a multivariate regression model with length & weight as DVs
mvreg length weight = ib1.brand_en
*mvreg, coeflegend
* Use -test- to compare the 2.brand coefficients for the two DVs
test _b[length:2.brand_en] = _b[weight:2.brand_en]
* -lincom- will give the same result, but will report a CI too
lincom _b[length:2.brand_en] - _b[weight:2.brand_en]

Output from -test- and -lincom-:

. test _b[length:2.brand_en] = _b[weight:2.brand_en]

 ( 1)  [length]2.brand_en - [weight]2.brand_en = 0

       F(  1,    12) =    0.34
            Prob > F =    0.5681

. * -lincom- will give the same result, but will report a CI too
. lincom _b[length:2.brand_en] - _b[weight:2.brand_en]

 ( 1)  [length]2.brand_en - [weight]2.brand_en = 0

------------------------------------------------------------------------------
             | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |      299.5   510.2158     0.59   0.568    -812.1647    1411.165
------------------------------------------------------------------------------

[George Ford]

https://www.stata.com/support/faqs/statistics/multiple-comparison-tests/

[Bruce Weaver]

This page on Andrew Wheeler's website has some relevant examples, by the way.

• https://andrewpwheeler.com/2017/06/1...ent-variables/

[Tom Ford]

Thanks a lot George and Bruce for your help with this.

I have only one remaining concern.
In my real model I am using bootstrapped and clustered standard errors, but it seems that both, suest and mvreg are unable to account for both of these circumstances.
So the two coefficients I want to compare have standard errors thus estimated?

reg length ibn.brand_en,nocons cluster(foreign) vce(bootstrap, rep(99) seed(123))
reg weight ibn.brand_en,nocons cluster(foreign) vce(bootstrap, rep(99) seed(123))

mvref does not allow clustering and suest does not allow bootstrap.

[Bruce Weaver]

Hello Tom Ford. Have you considered using the -mixed- command to estimate your model? It would handle the clustered nature of the data quite nicely. And I believe you can use -bootstrap- as a prefix command. As for the multivariate part, I believe one can include a categorical explanatory variable indicating which of the multiple DVs is present on a given row. I know that several years ago, Ryan Black posted some examples of that to an SPSS discussion forum. I expect there could be similar examples in the archives of this discussion forum. I hope this helps.

[Bruce Weaver]

Here are some of the old SPSS discussions I mentioned in #7.

• http://spssx-discussion.165.s1.nabbl...td1092023.html
• http://spssx-discussion.165.s1.nabbl...td3392871.html

Unfortunately, the link for Ronald Heck's slides on how to estimate a multivariate model via the (SPSS) MIXED command is no longer working, and I don't think I saved a copy of the slides.

Oh dear! While digging around trying to find Dr. Heck's website, I discovered that he died earlier this year at the age of 71. How sad.

[Enrique Pinzon (StataCorp)]

Dear Tom,

If you want to estimate two regressions in a system of equations as you are trying to do with -suest-, you can play with any of Stata's commands that allow you to model multiple equations simultaneously and then use -test- as you find fit. The advantage here is that you can have access to the standard-error computations you want. Below is an example with -gmm- and another with -gsem-. First the code,

webuse cattaneo2, clear 

gmm (eq1: msmoke - {msmoke: ibn.mmarried c.mage#c.mage})          ///
    (eq2: alcohol -{alcohol: ibn.mmarried c.mage#c.mage}),       ///
    instruments(eq1: ibn.mmarried c.mage#c.mage, nocons)          ///
    instruments(eq2: ibn.mmarried c.mage#c.mage, nocons)          ///
    winitial(unadjusted, independent) vce(cluster birthmonth)

test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4
test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4, mtest(bonf)

gsem (msmoke <- ibn.mmarried c.mage#c.mage, nocons)      ///
     (alcohol <- ibn.mmarried c.mage#c.mage, nocons),    ///
     vce(cluster birthmonth)
     
test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4
test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4, mtest(bonf)

Next the output:

. gmm (eq1: msmoke - {msmoke: ibn.mmarried c.mage#c.mage})          ///
>     (eq2: alcohol -{alcohol: ibn.mmarried c.mage#c.mage}),       ///
>     instruments(eq1: ibn.mmarried c.mage#c.mage, nocons)          ///
>         instruments(eq2: ibn.mmarried c.mage#c.mage, nocons)          ///
>         winitial(unadjusted, independent) vce(cluster birthmonth)

Step 1
Iteration 0:  GMM criterion Q(b) =  .19726942  
Iteration 1:  GMM criterion Q(b) =  1.135e-26  
Iteration 2:  GMM criterion Q(b) =  7.418e-34  

Step 2
Iteration 0:  GMM criterion Q(b) =  1.257e-33  
Iteration 1:  GMM criterion Q(b) =  7.884e-34  

note: model is exactly identified.

GMM estimation 

Number of parameters =   6
Number of moments    =   6
Initial weight matrix: Unadjusted                 Number of obs   =      4,642
GMM weight matrix:     Cluster (birthmonth)

                             (Std. err. adjusted for 12 clusters in birthmonth)
-------------------------------------------------------------------------------
              |               Robust
              | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
--------------+----------------------------------------------------------------
msmoke        |
     mmarried |
 Not married  |   .7102409   .0403478    17.60   0.000     .6311607    .7893211
     Married  |   .3064983   .0397798     7.70   0.000     .2285313    .3844652
              |
c.mage#c.mage |  -.0000383   .0000503    -0.76   0.446     -.000137    .0000603
--------------+----------------------------------------------------------------
alcohol       |
     mmarried |
 Not married  |   .0325416   .0034994     9.30   0.000     .0256829    .0394003
     Married  |  -.0129123   .0108988    -1.18   0.236    -.0342736    .0084491
              |
c.mage#c.mage |    .000043    .000013     3.32   0.001     .0000176    .0000684
-------------------------------------------------------------------------------
Instruments for equation eq1: 0.mmarried 1.mmarried c.mage#c.mage
Instruments for equation eq2: 0.mmarried 1.mmarried c.mage#c.mage

. 
. test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4

 ( 1)  [msmoke]1.mmarried - [alcohol]1.mmarried = .4

           chi2(  1) =    3.85
         Prob > chi2 =    0.0497

. test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4, mtest(bonf)

 ( 1)  [msmoke]1.mmarried - [alcohol]1.mmarried = .4

--------------------------------------
       |        chi2     df   p > chi2
-------+------------------------------
   (1) |        3.85      1     0.0497*
-------+------------------------------
   All |        3.85      1     0.0497
--------------------------------------
* Bonferroni-adjusted p-values

. 
. gsem (msmoke <- ibn.mmarried c.mage#c.mage, nocons)      ///
>          (alcohol <- ibn.mmarried c.mage#c.mage, nocons),    ///
>          vce(cluster birthmonth)

Iteration 0:  Log pseudolikelihood =  -4500.528  
Iteration 1:  Log pseudolikelihood =  -4500.528  

Generalized structural equation model                    Number of obs = 4,642

Response: msmoke  
Family:   Gaussian
Link:     Identity

Response: alcohol 
Family:   Gaussian
Link:     Identity

Log pseudolikelihood = -4500.528

 ( 1)  [msmoke]_cons = 0
 ( 2)  [alcohol]_cons = 0
                              (Std. err. adjusted for 12 clusters in birthmonth)
--------------------------------------------------------------------------------
               |               Robust
               | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
msmoke         |
      mmarried |
  Not married  |   .7102409   .0421419    16.85   0.000     .6276443    .7928375
      Married  |   .3064983   .0415486     7.38   0.000     .2250645    .3879321
               |
 c.mage#c.mage |  -.0000383   .0000526    -0.73   0.466    -.0001413    .0000647
               |
         _cons |          0  (omitted)
---------------+----------------------------------------------------------------
alcohol        |
      mmarried |
  Not married  |   .0325416    .003655     8.90   0.000     .0253779    .0397052
      Married  |  -.0129123   .0113835    -1.13   0.257    -.0352235    .0093989
               |
 c.mage#c.mage |    .000043   .0000135     3.18   0.001     .0000165    .0000696
               |
         _cons |          0  (omitted)
---------------+----------------------------------------------------------------
  var(e.msmoke)|   .7716293    .030488                       .714129    .8337594
 var(e.alcohol)|   .0308859   .0021913                      .0268762    .0354937
--------------------------------------------------------------------------------

.          
. test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4

 ( 1)  [msmoke]1.mmarried - [alcohol]1.mmarried = .4

           chi2(  1) =    3.53
         Prob > chi2 =    0.0602

. test _b[msmoke:1.mmarried] - _b[alcohol:1.mmarried] = .4, mtest(bonf)

 ( 1)  [msmoke]1.mmarried - [alcohol]1.mmarried = .4

--------------------------------------
       |        chi2     df   p > chi2
-------+------------------------------
   (1) |        3.53      1     0.0602*
-------+------------------------------
   All |        3.53      1     0.0602
--------------------------------------
* Bonferroni-adjusted p-values

Everything above is just for illustration. I do not know if this is exactly what you want, but whenever I see -suest- I think that -gmm- and -gsem- would do the trick.

[George Ford]

Learned something.

Here's a simple variant.

sysuse auto, clear

eststo e1: reg mpg gear_ratio turn
eststo e2: reg weight gear_ratio turn 
suest e1 e2, r
lincom _b[e1_mean:gear_ratio] - _b[e2_mean:gear_ratio]

gmm (eq1: mpg - {mpg: gear_ratio turn _cons})                             ///
    (eq2: weight - {weight: gear_ratio turn _cons}),                         ///
    instruments(eq1: gear_ratio turn)                                    ///
    instruments(eq2: gear_ratio turn)                                    ///
    winitial(unadjusted, independent) vce(robust)
    
lincom _b[mpg:gear_ratio] - _b[weight:gear_ratio]