I am trying to plot estimates from multiple models using coefplot. I need to group all the estimates (for each independent variable) in one model, and then plot the models side by side.

Here's my code:

local indnumber "23 31 42 52 54 61 62 71 92"

foreach x of local indnumber{

eststo r`x'_2: reg ln_hourwage i.sex i.treat_pub b1.educ i.treat_pub#b1.educ i.treat_pub#i.sex b1.educ#i.sex i.treat_pub#b1.educ#i.sex ///
i.statecensus i.statecensus#b2010.year b1.month Age AgeSq b4.race b0.ptft b0.sector b43.occupation if industry == `x', vce(cluster clustervar)

eststo m`x'_2: nlcom ///
(GPG_1: _b[1.sex]) (GPG_2: _b[1.sex]+_b[2.educ#1.sex]) (GPG_3: _b[1.sex]+_b[3.educ#1.sex]) (GPG_4: _b[1.sex]+_b[4.educ#1.sex]) ///
(Male_1: _b[1.treat_pub]) (Male_2: _b[1.treat_pub]+_b[1.treat_pub#2.educ]) ///
(Male_3: _b[1.treat_pub]+_b[1.treat_pub#3.educ]) (Male_4: _b[1.treat_pub]+_b[1.treat_pub#4.educ]) ///
(TGPG_1: _b[1.treat_pub#1.sex]) (TGPG_2: _b[1.treat_pub#1.sex]+_b[1.treat_pub#1.sex#2.educ]) ///
(TGPG_3: _b[1.treat_pub#1.sex]+_b[1.treat_pub#1.sex#3.educ]) (TGPG_4: _b[1.treat_pub#1.sex]+_b[1.treat_pub#1.sex#4.educ]), post

}

I need to plot the estimates Male_1 Male_2 Male_3 and Male_4 side-by-side (with some offset) for each of the industries (indicated by the NAICS category in indnumber). The four estimates would be grouped for each industry and all 9 industry estimates would be in the same graph. I also need to rename Male_* to different education groups, and label the industries.

Any help would be great.