How to plot logistic regression probability predictions together with the cumulative logistic distribution

[ericmelse]

Following the 'not so serious' example of Nick Cox, how to plot 'the line of perfect agreement or the line of equality' in a scatter plot, I worked on a similar example of how to plot logistic regression probability predictions together with the cumulative logistic distribution:

sysuse citytemp, clear
gen b_jul=.
qui sum tempjul
replace b_jul=0 if tempjul<=r(mean)
replace b_jul=1 if tempjul>r(mean)

logistic b_jul tempjan
predict p_jul

twoway scatter p_jul tempjan, ms(Oh) mc(%50) jitterseed(2951413) jitter(12)  ///
|| fpfit p_jul tempjan    /// fit current data
|| function logistic(x), range(-4 4) lp(dash) lc(blue*.7) xaxis(2) /// 
graphregion(fcolor(white) margin(l-2 r-1 b-2 t-5) ) xsize(10) ysize(10)    ///
legend(order(2 "Fitted values" 3 "Logistic function") col(1) ring(0) pos(5) symx(*.4)) ///
ytitle("Below/above mean July temp (Pr.)") ysc(noextend) yla(0(.1)1, ang(h)) /// 
xtitle("`: var label tempjan' ({&degree}F)") xsc(titlegap(*5) noextend) ///
xtitle(" " , axis(2)) xla(-4(1)4, axis(2) labc(blue*.7)) xsc(lc(blue*.7) noextend axis(2))

which should provide the following figure:

This possibly is of educational interest.