sysuse auto, clear qenvnormal weight, gen(lower upper) reps(20000) qenvnormal weight, overall reps(20000) gen(lower2 upper2) su weight scalar mean = r(mean) scalar sd = r(sd) qplot weight lower2 upper2 lower upper, ms(oh none ..) /// c(. l l l l) lc(gs8 gs8 gs8 gs12 gs12) /// legend(order(2 "Overall" 4 "Pointwise") /// subtitle("95% reference" "region") ) /// ytitle("Weight, lb") xtitle(Normal quantiles) /// trscale(`=mean' + `=sd' * invnormal(@))