I have fitted a two - level cross-classified logistic models using runmlwin (mcmc). I would like to plot/graph an interaction between two variables but I have not been able to find an alternative to the margins command. This is the regression with factor notation that I run:

Code:
quietly runmlwin def1 cons c.Age##c.Age  0.employment i.Educ#i.employment dmarital1  i.practice i.Educ#i.practice  i.dev1 i.Educ#i.dev1 divorce i.class income i.Nchildren , level2(Period:cons) level1(Cohort:) discrete(dist(binomial) link(logit) denom(cons)) nopause

runmlwin def1 cons c.Age##c.Age  0.employment i.Educ#i.employment dmarital1  i.practice i.Educ#i.practice  i.dev1 i.Educ#i.dev1 divorce i.class income i.Nchildren, level2(Period:cons) level1(Cohort:) discrete(dist(binomial) link(logit) denom(cons)) mcmc(cc) initsprevious nopause
runmlwin, noheader noretable or
I then removed factor notation and added all the dummies, then set all variables to zero except the two variables of interest (education x employment interactions) and attempted to predict the outcome + SE :


Code:
quietly runmlwin def1 cons Age agesq  demploy1 deduc2#demploy1 deduc3#demploy1 deduc4#demploy1 dmarital1  dprac2 dprac3 demploy1 deduc2#dprac1 deduc2#dprac2  deduc2#dprac3 deduc3#dprac1 deduc3#dprac2  deduc3#dprac3     deduc4#dprac1 deduc4#dprac2  deduc4#dprac3  devd1 deduc2#devd1 deduc3#devd1 deduc4#devd1 divorce devd1 dclass2 dclass3 income dchild2 dchild3 dchild4 dchild5, level2(Period:cons) level1(Cohort:) discrete(dist(binomial) link(logit) denom(cons)) nopause

runmlwin def1 cons Age agesq  demploy1 deduc2#demploy1 deduc3#demploy1 deduc4#demploy1 dmarital1  dprac2 dprac3 demploy1 deduc2#dprac1 deduc2#dprac2  deduc2#dprac3 deduc3#dprac1 deduc3#dprac2  deduc3#dprac3     deduc4#dprac1 deduc4#dprac2  deduc4#dprac3  devd1 deduc2#devd1 deduc3#devd1 deduc4#devd1 divorce devd1 dclass2 dclass3 income dchild2 dchild3 dchild4 dchild5, level2(Period:cons) level1(Cohort:) discrete(dist(binomial) link(logit) denom(cons)) mcmc(cc) initsprevious nopause
runmlwin, noheader noretable or

local unwanted = "Age agesq dmarital1 dprac2 dprac3 devd1 divorce dclass2 income dchild2 dchild3 dchild4 dchild5 "
foreach var in `unwanted' {
    replace `var' = 0
}

predict yhat
predict yhat_se, stdp
gen yhat_lb = yhat - 1.96 * yhat_se
gen yhat_ub = yhat + 1.96 * yhat_se
The predict yhat returns the following error:
variable _0b_deduc2_1_demploy1 not found


For reference, I am attempting to achieve this:

Code:
margins, at(Educ=(0(1)3)) over(employment)
marginsplot