I am summarizing percentage estimates across cities graphically, using margins and marginsplot. However, for particular cities I wish to make the marker a different symbol. (And I want the legend to include both the primary symbol and that special symbol.) While I can change markers by margin using plot1opt, plot2opt, etc., I don't see how to change marker according to categorical value --- that is, according to the categories listed along the x-axis.

Here is reproducible example of my problem, below. The first graph works fine, but imagine that I wanted to change the SECOND marker to a square rather than a circle. I can't figure out how to do it; neither of the two attempts at the bottom are successful.

I'm open to code that gets at the same figure without margins(plot). I tried unsuccessful hacks with both combomarginsplot and coefplot.


Code:
sysuse auto.dta, clear
gen heavy=weight>3500
gen cat = substr(make, 1, strpos(make, " ") - 1) 
bysort cat: gen x=_n
bysort cat: egen y=max(x)
drop if x==y
encode(cat), gen(emake)
drop if emake==.

reg heavy i.emake
margins emake

* this works fine, but let's say I want the SECOND market to be a square
marginsplot, recast(scatter) ///
    plotopts(msymbol(O))

* this doesn't work; it's for various margins not various values
marginsplot, recast(scatter) ///
    plotopts(msymbol(O)) plot2opts(msymbol(S))
    
* this felt like it might work, but it doesn't
marginsplot, recast(scatter) ///
    plotopts(msymbol(O S O O O O O O O O O O O O O O O O O O))
​​​