Dear community,

I'd like to use the "mrobust" command by Young & Holsteen (2017) for estimating model robustness. I investigate the effects of duration in preschool on school success. My logistic regression contains a squared term, the independent variable, and looks as follows:
Code:
logit success c.preschool##c.preschool c.loginc c.logwealth i.educ i.mig i.aspiration, vce(cluster region)
Unfortunately, the ado is not usable with interactions and other ploynomial terms. Therefore, my idea was to generate the squared term manually:
Code:
gen preschool_squared = preschool*preschool
This way, the mrobust command works, and with the grouping option I can make sure that the variables preschool and preschool_squared always are included in the robustness regressions.
Code:
mrobust logit success (preschool preschool_squared)  c.loginc c.logwealth i.educ i.mig i.aspiration, vce(cluster region)
The remaining problem is, that stata thinks of preschool_squared as a control variable, not as second independent variable. When estimating the model robustness, only the main term preschool coefficients are considered. As the squared term is negative, values of the mrobust results are higher than in my actual regressions.

Does anyone here have experience with the mrobust command and the possibilities for solving the polynomial problem? Or knows a similiar way to check on model robustness with another command than mrobust but similiar content?

Thank you very much,
Susanne

[References: Young, C. & K. Holsteen, 2017: Model Uncertainty and Robustness. Sociological Methods & Research 46: 3–40.]