Using SUEST to compare marginal effects of categorical interactions across two models?

[Meagan Rainock]

I am analyzing the interaction of two categorical independent variables (x1 (insurance) and x2 (ethnicity)) on two dichotomous dependent variables y1 (mental health care) and y2 (physical health care). My current code is:

logistic y1 i.x1##i.x2 ///
m1 m2 m3 if sample1==1, or nolog
margins x1##x2if e(sample), atmeans
margins if e(sample), dydx(x1) at(x2=(1 2 3 4 5 6)) atmeans post

logistic y2 i.x1##i.x2 ///
m1 m2 m3 if sample1==1, or nolog
margins x1##x2if e(sample), atmeans
margins if e(sample), dydx(x1) at(x2=(1 2 3 4 5 6)) atmeans post

And this provides me with the average marginal effects of insurance (x1) for each of the six categories of ethnicity (x2). Then, I am using suest to see if these AME's are statistically significantly different between the model predicting mental health care and the model predicting physical health care:

logistic y1 i.x1##i.x2 ///
m1 m2 m3 if sample1==1, or nolog
margins x1##x2if e(sample), atmeans
est store MODEL1

logistic y2 i.x1##i.x2 ///
m1 m2 m3 if sample1==1, or nolog
margins x1##x2if e(sample), atmeans
est store MODEL2

suest MODEL1 MODEL2

The trouble is, since I use an interaction effect in the logistic regression model, I am having trouble with the syntax for statistical significance tests that would compare the average marginal effects (AME's) of insurance (x1) for each category of ethnicity (x2) across model 1 and model 2. I can compare different coefficients across the two models using the test command in STATA, but suest just reports the coefficients for x1, x2, and x1##x2 and that isn't really what I want to compare. I've tried using margins after suest to no avail:

margins, dydx(x1) predict(MODEL1) atmeans

I get the following error:
Warning: cannot perform check for estimable functions.
option MODEL1 not allowed

Any suggestions for moving forward from here would be much appreciated!

[Andrew Musau]

margins x1##x2if e(sample), atmeans

These are not AMEs. They are adjusted predictions. There is no such thing as a marginal effect for an interaction term, the reason being that you cannot change the interaction term while holding its constituent variables constant. See more on definitions and this point in the following presentation by Richard Williams: https://www3.nd.edu/~rwilliam/stats/Margins01.pdf. For the comparison that you want, why not just test the coefficients?

[Meagan Rainock]

Thank you for your speedy response. You are correct that the highlighted line of syntax will report adjusted predictions - but it is my understanding that the next line of syntax ( "margins if e(sample), dydx(x1) at(x2=(1 2 3 4 5 6)) atmeans post" ) provides the average marginal effects of x1 for each of the six categories of x2. It is these 6 AME's that I would like to compare across the two different models.

I want to determine if the association between insurance status and mental health or physical health is the same or different for each of the 6 ethnic categories. I use the interaction term in order to determine this, the average marginal effects, as well as a test of second differences (which I have not depicted in the code here). But now I would like to determine if the AME's of x1 for each of the 6 categories of x2 are statistically significantly different between the model that uses mental health as the dependent variable, and the model that uses physical health as the dependent variable. In essence, I want to perform cross-model difference testing (see section 6.4 of the following article by Mize et al https://journals.sagepub.com/doi/pdf...81175019852763)

The example that is used in the Mize et al article is one where there are no interaction effects; you can compare the coefficients directly using wald tests after running suest. But I am not sure how to apply this to a model where an interaction effect is used, since an interaction effect coefficient cannot be interpreted in such a straightforward way. Is there a way to compare the results of "margins if e(sample), dydx(x1) at(x2=(1 2 3 4 5 6)" across the two models using suest?

Thank you very much for your help so far.

[Meagan Rainock]

Aha! I have found a way to do this through another route. I will post the instructions I found here for anyone else who finds themselves in a similar situation in the future.

https://www.trentonmize.com/software/mecompare
I discovered that Trenton Mize has stata do files here on his website that complement the article I mentioned above. This file uses gsem and vce(robust) to achieve the same results as the suest command would. I can then run the syntax "margins if e(sample), dydx(x1) at(x2=(1 2 3 4 5 6)) atmeans post" and it will give me the results that I desire, and then I will use mlincom to compare these AME's across the two models.

Thank you very much for your help.

[Andrew Musau]

Thank you for your speedy response. You are correct that the highlighted line of syntax will report adjusted predictions - but it is my understanding that the next line of syntax ( "margins if e(sample), dydx(x1) at(x2=(1 2 3 4 5 6)) atmeans post" ) provides the average marginal effects of x1 for each of the six categories of x2. It is these 6 AME's that I would like to compare across the two different models.

You could also check out cmp from SSC which can allow you to estimate both equations simultaneously. However, you would need to switch from logit to probit.

[Andrew Musau]

Empirically, it shouldn't matter whether you use logit or probit as they give very similar transformations from \(Xb\) to \(Pr(Xb)\). Here is an example:

webuse lbw, clear
set seed 07112024
gen outcome2= runiformint(0,1)
cmp (low= age i.race smoke) (outcome2= age i.race smoke), ind($cmp_probit $cmp_probit)
margins, dydx(age) at(race=(1 2 3)) predict(pr eq(#1)) predict(pr eq(#2))

Res.:

. cmp (low= age i.race smoke) (outcome2= age i.race smoke), ind($cmp_probit $cmp_probit)

Fitting individual models as starting point for full model fit.
Note: For programming reasons, these initial estimates may deviate from your specification.
      For exact fits of each equation alone, run cmp separately on each.

Iteration 0:  Log likelihood =   -117.336  
Iteration 1:  Log likelihood = -109.26071  
Iteration 2:  Log likelihood = -109.21408  
Iteration 3:  Log likelihood = -109.21408  

Probit regression                                       Number of obs =    189
                                                        LR chi2(4)    =  16.24
                                                        Prob > chi2   = 0.0027
Log likelihood = -109.21408                             Pseudo R2     = 0.0692

------------------------------------------------------------------------------
         low | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |  -.0218726   .0201108    -1.09   0.277    -.0612891    .0175439
             |
        race |
      Black  |   .6183845   .2994882     2.06   0.039     .0313984    1.205371
      Other  |   .6455825   .2401702     2.69   0.007     .1748575    1.116308
             |
       smoke |   .6758233   .2210635     3.06   0.002     .2425469      1.1091
       _cons |  -.6062843   .5187651    -1.17   0.243    -1.623045    .4104765
------------------------------------------------------------------------------

Iteration 0:  Log likelihood = -130.23923  
Iteration 1:  Log likelihood = -130.01781  
Iteration 2:  Log likelihood =  -130.0178  

Probit regression                                       Number of obs =    189
                                                        LR chi2(4)    =   0.44
                                                        Prob > chi2   = 0.9788
Log likelihood = -130.0178                              Pseudo R2     = 0.0017

------------------------------------------------------------------------------
    outcome2 | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |  -.0114787   .0177833    -0.65   0.519    -.0463334    .0233759
             |
        race |
      Black  |  -.0245494   .2849971    -0.09   0.931    -.5831336    .5340347
      Other  |  -.0494902   .2166297    -0.23   0.819    -.4740767    .3750962
             |
       smoke |  -.0028934   .2001845    -0.01   0.988    -.3952478     .389461
       _cons |   .1756254   .4744082     0.37   0.711    -.7541976    1.105448
------------------------------------------------------------------------------

Fitting constant-only model for LR test of overall model fit.

Fitting full model.

Iteration 0:  Log likelihood = -239.06234  
Iteration 1:  Log likelihood = -239.02267  
Iteration 2:  Log likelihood = -239.02267  

Mixed-process regression                                Number of obs =    189
                                                        LR chi2(8)    =  16.65
Log likelihood = -239.02267                             Prob > chi2   = 0.0340

------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
low          |
         age |   -.021752   .0200635    -1.08   0.278    -.0610757    .0175717
             |
        race |
      Black  |   .6185183   .2992926     2.07   0.039     .0319156    1.205121
      Other  |   .6436665   .2398411     2.68   0.007     .1735865    1.113747
             |
       smoke |   .6737483   .2207688     3.05   0.002     .2410495    1.106447
       _cons |  -.6077316   .5176542    -1.17   0.240    -1.622315    .4068519
-------------+----------------------------------------------------------------
outcome2     |
         age |  -.0114729     .01778    -0.65   0.519     -.046321    .0233753
             |
        race |
      Black  |  -.0242417   .2848028    -0.09   0.932     -.582445    .5339616
      Other  |   -.049205   .2165518    -0.23   0.820    -.4736386    .3752287
             |
       smoke |  -.0027503   .2000944    -0.01   0.989    -.3949281    .3894275
       _cons |   .1753136   .4744906     0.37   0.712    -.7546708    1.105298
-------------+----------------------------------------------------------------
/atanhrho_12 |    .079936   .1237268     0.65   0.518    -.1625641     .322436
-------------+----------------------------------------------------------------
      rho_12 |   .0797661   .1229395                      -.161147    .3117079
------------------------------------------------------------------------------

. 
. margins, dydx(age) at(race=(1 2 3)) predict(pr eq(#1)) predict(pr eq(#2))

Average marginal effects                                   Number of obs = 189
Model VCE: OIM

dy/dx wrt: age

1._predict: Pr(low), predict(pr eq(#1))
2._predict: Pr(outcome2), predict(pr eq(#2))

1._at: race = 1
2._at: race = 2
3._at: race = 3

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
age          |
_predict#_at |
        1 1  |  -.0059135   .0055168    -1.07   0.284    -.0167261    .0048992
        1 2  |    -.00796   .0072249    -1.10   0.271    -.0221206    .0062005
        1 3  |  -.0079999   .0072721    -1.10   0.271     -.022253    .0062531
        2 1  |  -.0045493   .0070343    -0.65   0.518    -.0183362    .0092376
        2 2  |  -.0045378   .0069916    -0.65   0.516    -.0182411    .0091654
        2 3  |  -.0045233   .0069731    -0.65   0.517    -.0181902    .0091437
------------------------------------------------------------------------------

.