Coefplot for several estimation models

[Noemi Seng]

Dear community,

I want to do a plot that shows the coefficients for the same variable but estimated with 5 different models. I include code for a reproducible example with 3 different models here:

use "http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA" , clear
egen imp = group(isoimp)
egen exp = group(isoexp)

preserve
keep if category == "MANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model1, replace
restore

preserve
keep if category == "NONMANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model2, replace
restore

preserve
keep if category == "TOTAL"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model3, replace
restore

estimates use model1.ster
estimates store model1

estimates use model2.ster
estimates store model2

estimates use model3.ster
estimates store model3


coefplot model1 model2 model3, ///
         keep(fta) ///
         xline(0) ///
         ciopts(recast(rcap) color(blue)) ///
         msymbol(O) ///
         ytitle("Coefficient of lndistw") ///
         xtitle("Model") ///
         xlabel(1 "sector1" 2 "sector2" 3 "sector3" , angle(0) labsize(medium)) ///
         vertical ///

This yields:

However, I want the coefficients and confidence intervals nicely spread over the whole range of the x-axis and I want each of them labelled with the name of the estimation model on the x-axis.

Is this possible with coefplot?



Best
Noemi

---

[Andrew Musau]

coefplot is from SSC, as you are asked to explain (FAQ Advice #12).

use "http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA" , clear
egen imp = group(isoimp)
egen exp = group(isoexp)

preserve
keep if category == "MANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model1, replace
restore

preserve
keep if category == "NONMANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model2, replace
restore

preserve
keep if category == "TOTAL"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model3, replace
restore

estimates use model1.ster
estimates store model1

estimates use model2.ster
estimates store model2

estimates use model3.ster
estimates store model3


coefplot (model1, aseq("Sector 1")) (model2, aseq("Sector 2")) (model3, aseq("Sector 3")),  ///
         keep(fta) ///
         xline(0) ///
         ciopts(recast(rcap) color(blue)) ///
         msymbol(O) mcolor(blue) ///
         ytitle("Coefficient of lndistw") ///
         xtitle("") ///
         xlabel( , angle(0) labsize(medium)) ///
         vertical aseq swapnames nokey

[Noemi Seng]

Dear Andrew Musau

thanks a lot. The only thing: When I execute the exact same code, the "Sector 1" labels are missing on the x-axis.

use "http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA" , clear
egen imp = group(isoimp)
egen exp = group(isoexp)

preserve
keep if category == "MANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model1, replace
restore

preserve
keep if category == "NONMANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model2, replace
restore

preserve
keep if category == "TOTAL"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model3, replace
restore

estimates use model1.ster
estimates store model1

estimates use model2.ster
estimates store model2

estimates use model3.ster
estimates store model3



coefplot (model1, aseq("Sector 1")) (model2, aseq("Sector 2")) (model3, aseq("Sector 3")),  ///
         keep(fta) ///
         xline(0) ///
         ciopts(recast(rcap) color(blue)) ///
         msymbol(O) mcolor(blue) ///
         ytitle("Coefficient of lndistw") ///
         xtitle("") ///
         xlabel( , angle(0) labsize(medium)) ///
         vertical aseq swapnames nokey

[Noemi Seng]

Edit:

Actually, it might be a problem of the font size as I could identify tiny tiny tiny dots below the 3 dashes on the x-axis. They might be the labels, but I'm not sure.

[Andrew Musau]

ssc install coefplot, replace

and try again.

[Noemi Seng]

Dear Andrew Musau

it worked after quitting Stata and opening it again!

However, as I have a lot of models (25), the labels ("Sector 1", "Sector 2", etc.) do all overlap such that they can't be read anymore. Do you maybe know if I could also just label them with "1", "2", etc. "25"?

Best,
Noemi

[Andrew Musau]

xlab(1/25) xtitle("Sector")

[Andrew Musau]

If the sector names are meaningful, use horizontal CIs.

[Noemi Seng]

Dear Andrew Musau

thanks a lot. I was also thinking about marking individual coefficients (& confidence intervals) by displaying them in a different colour. Can I also code this manually somehow?

[Andrew Musau]

Do you know 25 colors? You could do it, but that many colors don’t add anything here. In any case, each sector has its own label.

[Noemi Seng]

Ah no, I was just thinking about marking two or three of the coefficients that e.g. deviate a lot from the overall pattern in another color.

[Andrew Musau]

use "http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA" , clear
egen imp = group(isoimp)
egen exp = group(isoexp)

preserve
keep if category == "MANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model1, replace
restore

preserve
keep if category == "NONMANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model2, replace
restore

preserve
keep if category == "TOTAL"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model3, replace
restore

estimates use model1.ster
estimates store model1

estimates use model2.ster
estimates store model2

estimates use model3.ster
estimates store model3


coefplot (model1) (model2) ///
         (model3, ciopts(color(red) recast(rcap)) mcolor(red) ),  ///
         keep(fta) ///
         xline(0) ///
         ciopts(recast(rcap) color(blue))  ///
         msymbol(O) mcolor(blue) ///
         ytitle("Coefficient of lndistw") ///
         xtitle("") ///
         xlabel( , angle(0) labsize(medium)) ///
         vertical aseq swapnames nokey

[Noemi Seng]

Dear Andrew Musau

thank you very much, that was very helpful. For my plot, I want to compare the coefficients of all of the separate 25 models with the coefficient of a benchmark model. At the moment, I added the benchmark coefficient with a dashed line in red as well as dashed lines in green for the confidence interval of the benchmark estimate:

coefplot (sector1, aseq("Sector 1")) (sector2, aseq("Sector 2")) (sector3, aseq("Sector 3")) (sector4, aseq("Sector 4")) ///
        (sector5, aseq("Sector 5")) (sector6, aseq("Sector 6")) (sector7, aseq("Sector 7")) (sector8, aseq("Sector 8")) ///
        (sector9, aseq("Sector 9")) (sector10, aseq("Sector 10")) (sector11, aseq("Sector 11")) (sector12, aseq("Sector 12")) ///
        (sector13, aseq("Sector 13")) (sector14, aseq("Sector 14")) (sector15, aseq("Sector 15")) (sector16, aseq("Sector 16")) ///
        (sector17, aseq("Sector 17")) (sector18, aseq("Sector 18")) (sector19, aseq("Sector 19")) (sector20, aseq("Sector 20")) ///
        (sector21, aseq("Sector 21")) (sector22, aseq("Sector 22")) (sector23, aseq("Sector 23")), ///
         keep(lndistw) ///
         ciopts(recast(rcap) color(blue)) ///
         msymbol(O) mcolor(blue) ///
         ytitle("Coefficient of lndistw") ///
         xtitle(" ") ///
         yline (-0.405, lcolor(red)) ///
         yline (0, lcolor(black)) ///
          yline (-0.597, lcolor(green)) ///
         yline (-0.221, lcolor(green)) ///
         xlabel( , angle(0) labsize(medium)) ///
         vertical aseq swapnames nokey ///
         level(90) ///
         xlab(1/23) xtitle("Sector")

Can you maybe think of a visually more appealing way of including the benchmark estimate, especially the boundaries of the confidence intervals?



Best
Noemi

[Andrew Musau]

Maybe area shading. Below the CIs are specified manually.

use "http://fmwww.bc.edu/RePEc/bocode/e/EXAMPLE_TRADE_FTA_DATA" , clear
egen imp = group(isoimp)
egen exp = group(isoexp)

preserve
keep if category == "MANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model1, replace
restore

preserve
keep if category == "NONMANUF"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model2, replace
restore

preserve
keep if category == "TOTAL"
ppmlhdfe trade fta, a(imp#year exp#year imp#exp) cluster(imp#exp)
estimates save model3, replace
restore

estimates use model1.ster
estimates store model1

estimates use model2.ster
estimates store model2

estimates use model3.ster
estimates store model3

gen ll=.1102
gen ul=.2747
gen range = cond(_n==1, _n-.5, _n+.5)


coefplot (model1) (model2) ///
         (model3, ciopts(color(red) recast(rcap)) mcolor(red) ),  ///
         keep(fta) ///
         xline(0) ///
         ciopts(recast(rcap) color(blue))  ///
         msymbol(O) mcolor(blue) ///
         ytitle("Coefficient of lndistw") ///
         xtitle("") ///
         xlabel( , angle(0) labsize(medium)) ///
         vertical aseq swapnames nokey ///
         addplot(rarea ll ul range, color(gs8%40))

[Noemi Seng]

Dear Andrew Musau

I have worked a lot on my plot and it looks like this now, done with the following code:

gen ll=-0.799
gen ul=-0.34
gen range = cond(_n==1, _n-.5, _n+.5)    
coefplot (sector1, aseq("Sector 1")) (sector2, aseq("Sector 2")) ///
        (sector3, aseq("Sector 3") ciopts(color(orange) recast(rcap)) mcolor(orange)) ///
        (sector4, aseq("Sector 4")) ///
        (sector5, aseq("Sector 5")) (sector6, aseq("Sector 6")) (sector7, aseq("Sector 7")) (sector8, aseq("Sector 8")) ///
        (sector9, aseq("Sector 9")) (sector10, aseq("Sector 10")) (sector11, aseq("Sector 11")) ///
                (sector12, aseq("Sector 12") ciopts(color(orange) recast(rcap)) mcolor(orange)) ///
        (sector13, aseq("Sector 13")) (sector14, aseq("Sector 14")) (sector15, aseq("Sector 15")) (sector16, aseq("Sector 16")) ///
        (sector17, aseq("Sector 17")) ///
        (sector18, aseq("Sector 18") ciopts(color(orange) recast(rcap)) mcolor(orange)) ///
        (sector19, aseq("Sector 19")) (sector20, aseq("Sector 20")) ///
        (sector21, aseq("Sector 21")) (sector22, aseq("Sector 22")) ///
            (sector23, aseq("Sector 23") ciopts(color(orange) recast(rcap)) mcolor(orange)), ///
         keep(lndistw) ///
         ciopts(recast(rcap) color(blue)) ///
         msymbol(O) mcolor(blue) ///
         ytitle("Coefficient of distance (log)") ///
         xtitle(" ") ///
         yline (-0.57, lcolor(red)) ///
         yline (0, lcolor(black)) ///
         xlabel( , angle(0) labsize(medium)) ///
         vertical aseq swapnames nokey ///
         level(90) ///
         xlab(1/23) xtitle("Sector") ///
         addplot(rarea ll ul range, color(gs8%40)) ///
              xlabel(1 "11" 2 "21" 3 "22" 4 "23" ///
        5 "31" 6 "32" 7 "33" 8 "42" ///
        9 "44" 10 "45" 11 "48" ///
        12 "49" 13 "51" 14 "52" ///
        15 "53" 16 "54" 17 "55" ///
        18 "56" 19 "61" 20 "62" ///
        21 "71" 22 "72" 23 "81", labsize(*0.64))
drop ll ul range

Do you maybe know how I can put the numbers of the sector, whose coefficient & CI are in orange, in orange as well? I.e. the numbers 22, 49, 56, 81 at the x-axis should appear in orange. I tried it with adding labcolor(orange), but this changes the color for all numbers at the x-axis. Also I tried to add several xlabel(...) commands one after the other, but it seems that the last command overwrites the preceding ones.

I appreciate any hint!

Best
Noemi