Dear Chris,

Your solution worked like a charm, many thanks for your guidance!!

May I follow up with one further question? Suppose we expect the effect of

*educ *to decline in

*d_lit*. Without predicted probability and marginal effects plots it is hard to interpret this logit interaction model credibly.

Is there any way to obtain the predicted probabilities of

*use* at various levels of

*educ* and

*d_lit* (ideally with CIs) as well as the marginal effects? I have seen the suggestion to use MCMC chains for that here in the forum, but I am not entirely sure how this would be worked out. Would you have any good case practice on how to get there, i.e. how "imitate" the margins command of Stata when the model was estimated in MLwiN?

Code: Select all

```
ar all
use http://www.bristol.ac.uk/cmm/media/runmlwin/bang, clear
generate lc1 = (lc==1)
generate lc2 = (lc==2)
generate lc3plus = (lc>=3)
foreach var of varlist hindu urban {
bysort district: egen M_`var' =mean(`var')
}
*Stata
logit use cons lc1 lc2 lc3plus age M_urban M_hindu i.educ##c.d_lit, vce(cluster district) // using "logit" rather than "melogit" for illustration to speed up estimation
est sto m
margins, noatlegend at(educ==(1(1)4) d_lit==(0(0.1)0.3)) atmeans post
marginsplot, xdimension(d_lit) recast(line)
est resto m
margins, dydx(d_lit) at(educ==(1(1)4)) atmeans post
marginsplot, xdimension(educ) recast(line) yline(0)
*MLwiN
gen double educ_X_d_lit = educ*d_lit // new moderation of interest
runmlwin use cons lc1 lc2 lc3plus age d_lit M_urban M_hindu educ educ_X_d_lit , ///
level2(district: cons) level1(woman) discrete(distribution(binomial) link(logit) denominator(cons)) ///
nopause forcesort zratio level(99)
```