Legend marker with both the dstat line and its confidence interval

[ericmelse]

Dear Ben,
Is it possible to have a legend marker with both the dstat line and its confidence interval?
Using this code:

local var pir
local vart PIR
local binary clip_in
qui sum `var'
local N = r(N)
sum `var' if `binary'==1
local Nr = r(N)
local Nrp = round((r(N)/`N')*100,.1)
sum `var' if `binary'==0
local Nx = r(N)
local Nxp = round((r(N)/`N')*100,.1)

mylabels 0(10)100, local(myxlab) myscale(@/100) suffix("%")
dstat density `var' , over(`binary') n(200) k(e) ll(0) ul(1) unconditional /// if `constraint'
	gr(aspect(1) xsize(8) ysize(8) graphreg(fc(white) ilc(white) lc(white) m(l-2 r-0 b-0 t-0)) ///
	xsc(noex) xlab(`myxlab', labs(*.8) ) ylab(0(.1)2, labs(*.7) angle(none) nogrid) ///
	p1(lc(stred) ciopts(color(stred%20))) ///
	p2(lc(stblue) ciopts(color(stblue%20))) ///
	xtitle("Person isotonic R² (PIR)", m(l+1 r-1 t+1 b-1)) ytitle("Density (Epanechnikov kernel function)", m(t+1 b-1)) ///
	legend(order(2 "95% CI" 4 "Retained" "n=`Nr', `Nrp'%" 1 "95% CI" 3 "Excluded" "n=`Nx', `Nxp'%") ///
	region(c(none)) size(*.7) col(1) pos(11) rowgap(0) ring(0)) ///
	ysc(range(0 2) noex) ylabel(0(.1)2) flip name(HHdstatPIRclip, replace)) 

The result plot is:

But, what I am trying to get is:

Perhaps this is not possible, but if it is, time permitting, I really appreciate your helping hand,

[benjann]

Hi Eric, this is difficult to do in Stata because the legend() option is not flexible enough (cannot print multiple keys on top of each other). The only approach I see is to create a custom legend using using twoway commands in addplot(). The difficulty with this approach is that you have to use the plot region's coordinates to position the elements. Here's an example.

sysuse nlsw88, clear
dstat pdf wage, over(union) notable

// colors
local c1 stred
local c2 stblue

// legend
local l 30   // left x-coordinate of legend keys
local r 35   // right x-coordinate of legend keys
local t 36   // x-coordinate of legend labels
local y1 .14 // vertical position of first key
local y2 .12 // vertical position of second key
local bopts recast(bar) horizontal barwidth(.01) lc(%0) base(`l')
local lopts recast(line)
local topts place(e) just(left)
local legend ///
    /// first key
       scatteri `r' `y1', color(`c1'%20) `bopts' ///
    || scatteri `y1' `l' `y1' `r', lc(`c1') `lopts' ///
       text(`y1' `t' "Nonunion" "n = 1417", `topts') ///
    /// second key
    || scatteri `r' `y2', color(`c2'%20) `bopts' ///
    || scatteri `y2' `l' `y2' `r', lc(`c2') `lopts' ///
       text(`y2' `t' "Union" "n = 461", `topts')

// plot
dstat graph, overlay legend(off) ///
    p1(lc(`c1') ciopts(color(`c1'%20))) ///
    p2(lc(`c2') ciopts(color(`c2'%20))) ///
    addplot(`legend')

ben

[ericmelse]

Thank you Ben. Certainly this is a practical solution. I can replicate your example just fine. But, when I applied your solution to my own data for some reason the recast(bar) is refused to process a graphical result although there is no error message and the plot is generated. So, I replaced that part with recast(line) and by tweaking the line width etc. all is as I want to have it. Here is my code:

// colors
local c1 stred
local c2 stblue
// statistics
local var pir
local binary clip_in
qui sum `var'
local N = r(N)
qui sum `var' if `binary'==1
local Nr = r(N)
local Nrp = round((r(N)/`N')*100,.1)
qui sum `var' if `binary'==0
local Nx = r(N)
local Nxp = round((r(N)/`N')*100,.1)
// legend
local l1 0    // left x-coordinate of first key
local r1 .1   // right x-coordinate of first key
local l2 .9   // left x-coordinate of second key
local r2 1    // right x-coordinate of second key
local t1 .12  // x-coordinate of of first key label
local t2 .88  // x-coordinate of of first key label
local y1 2    // vertical position of first key
local y2 2    // vertical position of second key
local lopts recast(line)
local topts1 place(e) just(left) size(*.6)
local topts2 place(w) just(right) size(*.6)
local legend ///
    /// first key
       scatteri `y1' `l1' `y1' `r1', lc(`c2'%20) lw(20pt) `lopts' ///
    || scatteri `y1' `l1' `y1' `r1', lc(`c2') `lopts' ///
       text(`y1' `t1' "Retained" "n=`Nr', `Nrp'%", `topts1') ///
    /// second key
    || scatteri `y2' `l2' `y2' `r2', lc(`c1'%20) lw(20pt) `lopts' ///
    || scatteri `y2' `l2' `y2' `r2', lc(`c1') `lopts' ///
       text(`y2' `t2' "Excluded" "n=`Nx', `Nxp'%", `topts2')
// calculate the distribution density   
dstat density `var' , over(`binary') n(200) k(e) ll(0) ul(1) unconditional notable
// set x axis labels with % character
mylabels 0(10)100, local(myxlab) myscale(@/100) suffix("%")
// plot
dstat graph, overlay legend(off) ///
    p1(lc(`c1') ciopts(color(`c1'%20))) ///
    p2(lc(`c2') ciopts(color(`c2'%20))) ///
	aspect(1) xsize(8) ysize(8) graphreg(fc(white) ilc(white) lc(white) m(l-2)) ///
	xsc(noex) xlab(`myxlab', labs(*.8) )  ///
	ysc(range(0 2) noex) ylab(0(.1)2, labs(*.7) angle(none) nogrid) ///
	xtitle("Person isotonic R² (PIR)", m(l+1 r-1 t+1 b-1)) ///
	ytitle("Density (Epanechnikov kernel function)", m(t+1 b-1)) ///
    addplot(`legend')

which results in:

What I like about this particular positioning of the two legend keys is the symmetry while the whole plot is sort of balanced as well. Each plot area corner now has 'something going', i.e. the two legend keys at the top and the excluded cases at the low and high end of the distribution at the bottom. The distribution of the retained cases now stands on its own (visually), which is the objective of this plot: communicating about the distribution difference between thrustworthy respondents and those we should suspect not to be that.

[benjann]

But, when I applied your solution to my own data for some reason the recast(bar) is refused to process a graphical result although there is no error message and the plot is generated.

Usage of recast(bar) is a bit tricky because you need to specify the coordinates in reverse order and then apply horizontal, and provide one of the x-coordinates in base(). I guess something got mixed up so that invisible zero-length bars were plotted.

[benjann]

Hi Eric, it has been on my list since long to put together a do-it-yourself legend generator, mostly because Stata's legend option cannot combine multiple symbols in a single key. I now released such a command; see github.com/benjann/mklegend. You may find the command useful in applications such as the one above. Ben

[benjann]

Hi Erich, note that I discontinued mklegend and relaunched the package as addlegend. ben

[benjann]

Using addlegend you could create the legend in your example as follows:

dstat graph ...
addlegend, Y(2) lskip(0) text(size(*.6)): ///
       (line, lc(%20) lw(20pt)) (line) "Retained" "n=`Nr', `Nrp'%", lc(`c2') X(1) W(-.1) ///
    || (line, lc(%20) lw(20pt)) (line) "Excluded" "n=`Nx', `Nxp'%", lc(`c1') X(0) W(.1)

Or with slighly less typing:

dstat graph ...
local symbl (line, lc(%20) lw(20pt)) (line)
addlegend, Y(2) lskip(0) text(size(*.6)): ///
       `symbl' "Retained" "n=`Nr', `Nrp'%", lc(`c2') X(1) W(-.1) ///
    || `symbl' "Excluded" "n=`Nx', `Nxp'%", lc(`c1') X(0) W(.1)

ben

[ericmelse]

Dear Ben,
This a most welcome package. Thank you so much!
Time permitting, would you consider to publish about addlegend in TSJ?
Mundane as legends might be considered to be, it is one of those 'little things' that are so necessary to include in graphs, and are so much discussed on the Statalist because it is complicated if not impossible to control.

[benjann]

Yes, I'm tempted :)