margins accounting for indirect effects after GSEM mediation model? (pred prob's at specific values after GSEM with complete mediation)

[Rajan Sonik]

Hi all,

I am running a GSEM model in Stata/MP 16.0 with both a binary mediator and binary final outcome of interest. I am interested in an interaction between membership in a particular group (yes/no) and exposure to a particular intervention (also yes/no), yielding 4 groups of interest: (i) not group, no intervention, (ii) not group, yes intervention, (iii) yes group, no intervention, (iv) yes group, yes intervention. More precise code below, but broadly/simply, the path looks like: (group intervention groupXintervention covariates) --> binaryoutcome1 --> binaryoutcome2. For theoretical reasons, I am examining binaryoutcome1 as a complete mediator of the relationship between group##intervention and binaryoutcome2 (i.e., no direct path from group##intervention to binaryoutcome 2), though I think my question would be relevant even if that weren't the case.

Ultimately, I would like to know the predicted probabilities of being "yes" on each of the binary outcomes for each of groups (i) through (iv). I believe I am able to produce this for binaryoutcome1, but I am struggling about how to obtain similar predicted probabilities for each of these groups for binaryoutcome2.

To be more precise, a simplified version of the original model looks like the following:

gsem (binaryoutcome1 <- group intervention groupXintervention covariates, family(binomial) link(probit)) (binaryoutcome2 <- binaryoutcome1 covariates, family(binomial) link(probit)), nocapslatent

To get the values I am interested in for binaryoutcome1, I use:
margins, at(group=0 intervention=0 groupXintervention=0) predict(outcome(binaryoutcome1)) force
margins, at(group=0 intervention=1 groupXintervention=0) predict(outcome(binaryoutcome1)) force
margins, at(group=1 intervention=0 groupXintervention=0) predict(outcome(binaryoutcome1)) force
margins, at(group=1 intervention=1 groupXintervention=1) predict(outcome(binaryoutcome1)) force

However, when I use similar coding for binaryoutcome2--i.e.:
margins, at(group=0 intervention=0 groupXintervention=0) predict(outcome(binaryoutcome2)) force
margins, at(group=0 intervention=1 groupXintervention=0) predict(outcome(binaryoutcome2)) force
margins, at(group=1 intervention=0 groupXintervention=0) predict(outcome(binaryoutcome2)) force
margins, at(group=1 intervention=1 groupXintervention=1) predict(outcome(binaryoutcome2)) force

--I get the same value for all 4 groups. In other words, Stata appears to not recognize any effects of the values on group, intervention, or groupXintervention on the predicted probability of being a yes on binaryoutcome2.

I am presuming that this because is because the margins commands I am using are not accounting for indirect effects. This also makes me concerned in general that, even if I didn't have a full mediation model, the outputs from the margins commands above would be wrong, as they would not be accounting for both the direct and indirect effects (i.e., they would only be using the direct effects). For example, when I do the above where I don't model the mediation as complete--i.e., using this code:

gsem (binaryoutcome1 <- group intervention groupXintervention covariates, family(binomial) link(probit)) (binaryoutcome2 <- group intervention groupXintervention binaryoutcome1 covariates, family(binomial) link(probit)), nocapslatent

--and then the margins commands above, I get different values after each of the 4 margins commands, but now I think that is simply based on the direct effects of group##intervention on binaryoutcome2 rather than the total effects, which is what I would be interested in. Either way, the problem is the same: getting these predicted probabilities at specific values of predictors that accounts for the indirect effects.

I hope I have provided the necessary/relevant information; my apologies if not. I would be happy to provide any more details if needed to provide thoughts on this.

Thanks so much in advance for any thoughts!
Rajan

[William Lisowski]

I am not a user of gsem, so the following comes from general principles of Stata estimation commands. It may be inapplicable to your situation; if so, my apologies in advance.

I am surprised to see that you seem not to be using Stata's factor variable notation, which in general is considered a requirement for postestimation commands like margins to perform correctly.

I would have expected you to have described your model with something more like the following

gsem (binaryoutcome1 <- i.group##i.intervention covariates, family(binomial) link(probit)) (binaryoutcome2 <- binaryoutcome1 covariates, family(binomial) link(probit)), nocapslatent

and a margins command more like the following

margins i.group##i.intervention, predict(outcome(binaryoutcome1))

[Rajan Sonik]

Hi William,
Thanks for your message. The reason for not using factor notation in this instance is a bit complicated and outside of the parameters of the main question here. However, the basic reason is that the relationship between the intervention and binaryoutcome1 has a problem with confounding that I am addressing via instrumental variables. In order to interact the group variable with the intervention variable after doing the equivalent of a biprobit model with instrumental variables prevents me from using the factor notation.
All best,
Rajan

[William Lisowski]

Thank you for confirming that the choice to ignore factor variable notation was a deliberate one.

That knowledge is in fact an important parameter of the main question here, which I now see is more fully described as one of using margins in a situation where you find it necessary to overlook the explicit advice in the margins documentation to use factor variable notation to fully communicate the relationships between variables in your model.

Since as I indicated I'm not a user of gsem, I can't comment further on whether it introduces additional complications beyond the need to set the values of your interaction variable along with those of the two interacted variables, as the documentation suggests. Perhaps someone else will be able to advise you.

[Rajan Sonik]

Hi all,

I think I was able to figure it out. Using path diagram I wasn't able to use factor notation to interact the terms I referenced above, but using direct coding I was. So the ultimate code for the gsem command was (simplified as above to get at the heart of the question):

gsem (binaryoutcome1 <- i.group##i.intervention covariates, family(binomial) link(probit)) (binaryoutcome2 <- binaryoutcome1 covariates, family(binomial) link(probit)), nocapslatent

To get the adjusted predicted probabilities for binaryoutcome1, then, as above, I just ran:

margins group#intervention, predict(outcome(binaryoutcome2)) force

To get the adjusted predicted probabilities for binaryoutcome2 (the one I was struggling with), I ran:

margins group#intervention, predict(exp(normal(eta(binaryoutcome1)*_b[binaryoutcome2:1.binaryoutcome1] + eta(binaryoutcome2)))) force

indirect effects on binaryoutcome2 though binaryoutcome1 are accounted for by "eta(binaryoutcome1)*_b[binaryoutcome2:1.binaryoutcome1]"

direct effects (i.e., of the covariates) on binaryoutcome2 are accounted for by eta(binaryoutcome2)

I believe this solves it, though if any of this appears off to anyone would be curious to hear your thoughts.

Thanks,
Rajan