Graph regression line with overlaid Gaussian

[paulvonhippel]

In teaching I often see regression lines with overlaid Gaussian curves showing the residual distribution. How can I make one of these in Stata?

So you can see what I mean, here's a blog post showing how to do it in SAS:
https://blogs.sas.com/content/iml/20...reg-model.html

And here's an example in R:
https://stackoverflow.com/questions/...egression-line

There's a cooler version that's 3D, with the Gaussians seeming to rise off the line plot. I'll see if I can find it.

[Alan Neustadtl]

Here is an example that I have used in my classes that show one way to produce that kind of graphic. Note the use of the horizontal option:

#d ;
tw (function y=normalden(x,  0, 1)+0.0, range(-4 4) hori lw(medthick))
   (function y=normalden(x,  0, 1)+0.5, range(-4 4) hori lw(medthick))
   (function y=normalden(x,  0, 1)+1.0, range(-4 4) hori lw(medthick))
   (function y=normalden(x,  0, 1)-.75, range(-4 4) hori lw(thick) lc(maroon))

   (scatteri -3.50 -0.75  3.50 -0.75, ms(i) c(l) lw(thin) lc(black))
   (scatteri -3.50  0.00  3.50  0.00, ms(i) c(l) lw(thin) lc(black))
   (scatteri -3.50  0.50  3.50  0.50, ms(i) c(l) lw(thin) lc(black))
   (scatteri -3.50  1.00  3.50  1.00, ms(i) c(l) lw(thin) lc(black))

   (scatteri -7.5 -.25 -7.5 1.50, ms(i) c(l) lw(medthick) lc(black))
   (scatteri -7.5 -.25 10.0 -.25, ms(i) c(l) lw(medthick) lc(black)),
     ytitle("Y")
     note("SS{sub:between} =   0" "SS{sub:within}   =100" "SS{sub:total}    =100", j(right))
     xscale(off) yscale(off)
     ylabel(, nogrid)
     legend(off)
     yline(0, lp(dash))
     text(8 0.20 "{it:j}=1", size(medlarge))
     text(8 0.65 "{it:j}=2", size(medlarge))
     text(8 1.20 "{it:j}=3", size(medlarge))
     xsize(16) ysize(16) aspect(1)
     graphregion(color(white))
     name(g1, replace) nodraw;
#d cr

#d ;
tw (function y=normalden(x, -2.5,  .33)+0.00, range(-4.00 -1.00) hori lw(medthick))
   (function y=normalden(x, -0.0,  .33)+0.50, range(-1.50  1.50) hori lw(medthick))
   (function y=normalden(x,  2.5,  .33)+1.00, range( 1.18  3.82) hori lw(medthick))
   (function y=normalden(x,  0.0, 1.00)-0.75, range(-4.00  4.00) hori lw(thick) lc(maroon))

   (scatteri -4.00  0.00 -1.00  0.00, ms(i) c(l) lw(thin) lc(black))
   (scatteri -1.50  0.50  1.50  0.50, ms(i) c(l) lw(thin) lc(black))
   (scatteri  1.15  1.00  3.75  1.00, ms(i) c(l) lw(thin) lc(black))  /* Base purple */
   (scatteri -4.00 -0.76  4.00 -0.76, ms(i) c(l) lw(thin) lc(black))  /* Base SSTO   */

   (scatteri -2.50  0.00 -2.50  1.23, ms(i) c(l) lw(thin) lc(cranberry) lp(dash))
   (scatteri  2.50  1.00  2.50  2.23, ms(i) c(l) lw(thin) lc(cranberry) lp(dash))

   (scatteri -7.50 -0.25 -7.50  1.50, ms(i) c(l) lw(medthick) lc(black))
   (scatteri -7.50 -0.25 10.00 -0.25, ms(i) c(l) lw(medthick) lc(black)),
     ytitle("Y")
     note("SS{sub:between}=  59" "SS{sub:within}   =  41" "SS{sub:total}    =100", j(right))
     xscale(off) yscale(off)
     ylabel(, nogrid)
     legend(off)
     yline(0, lp(dash))
     text(8 0.55 "{it:j}=1", size(medlarge))
     text(8 1.00 "{it:j}=2", size(medlarge))
     text(8 1.55 "{it:j}=3", size(medlarge))
     xsize(16) ysize(16) aspect(1)
     graphregion(color(white))
     name(g2, replace) nodraw;
#d cr

graph combine g1 g2, nocopies graphregion(color(white)) ysize(2.5) name(anovaexplanation, replace)