/*
   In evaluating the fit of your model there are three related
     residual plots you can examine.

   Plot the residuals against:
     1. EACH X variable
     2. the predicted values (Y-hat)
     3. and plot Y against Y-hat

   Remember there are three kinds of residuals:
     1. Measured in the orgiinal metric
     2. Standardized
     3. Jacknifed or studentized

   See p. 223 of Kleinbaum et al. for more information.
*/

regress y  x
predict y_resid_org, residuals
predict y_resid_std, rstandard
predict y_resid_stu, rstudent
predict yhat1      , xb

regress y1 x
predict y1_resid_org, residuals
predict y1_resid_std, rstandard
predict y1_resid_stu, rstudent
predict yhat2      , xb

regress y2 x
predict y2_resid_org, residuals
predict y2_resid_std, rstandard
predict y2_resid_stu, rstudent
predict yhat3      , xb


/* Show Residuals, Standardized Residuals, and Studentized Residuals */
quietly {
  sort x
  local n=_N
  noisily di _n _col(3) "{text}" "{hline 74}" ///
             _n _skip(12) "Original Data" _skip(12) "Scenario #1"  _skip(13) "Scenario #2"  ///
             _n _skip(4)  "X"             _col(12)  "Res."         _col(20)  "Std. Res."    ///
                                          _col(36)  "Res."         _col(44)  "Std. Res."    ///
                                          _col(60)  "Res."         _col(68)  "Std. Res." _n

  forvalues i=1(1)`n' {
  noisily di "{result}" %5.0g x[`i' ] ///
             %7.2f  _skip(3) y_resid_org[`i'] _skip(3) %7.2f  y_resid_std[`i'] _skip(7) ///
             %7.2f y1_resid_org[`i'] _skip(3) %7.2f y1_resid_std[`i'] _skip(7)          ///
             %7.2f y2_resid_org[`i'] _skip(3) %7.2f y2_resid_std[`i']
  }
   noisily di _col(3) "{text}" "{hline 74}"
}


/* Plot Residuals in Y Metric Against X */
#delimit ;
graph twoway (scatter y_resid_org  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) 
   title("Original Data")
   xtitle("X")
   ytitle("Y Residuals")
   name(resid1, replace) nodraw
;

graph twoway (scatter y1_resid_org x, yline(0, lwidth(thick))),
   xlabel(0(1)10) 
   title("Scenario #1")
   xtitle("X")
   ytitle("Y Residuals")
   name(resid2, replace) nodraw
;

graph twoway (scatter y2_resid_org x, yline(0, lwidth(thick))),
   xlabel(0(1)10) 
   title("Scenario #2")
   xtitle("X")
   ytitle("Y Residuals")
   name(resid3, replace) nodraw
;

graph combine resid1 resid2 resid3,
   title("Residuals Plotted Against X")
   ycommon
   name(resid_x, replace) nocopies
;
#delimit cr


/* Plot Standardized Residuals Against X */
#delimit ;
local tplus =invttail(10-2,.05/2);
local tminus=invttail(10-2,.05/2)*-1;

graph twoway (scatter y_resid_std  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) ylabel(-3(1)3)
   title("Original Data")
   xtitle("X")
   ytitle("Standardized Residuals")
   yline(`tplus' `tminus', lp(dash) lw(thick))
   name(residstd1, replace) nodraw
;

graph twoway (scatter y1_resid_std  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) ylabel(-3(1)3)
   title("Scenario #1")
   xtitle("X")
   ytitle("Standardized Residuals")
   yline(`tplus' `tminus', lp(dash) lw(thick))
   name(residstd2, replace) nodraw
;

graph twoway (scatter y2_resid_std  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) ylabel(-3(1)3)
   title("Scenario #2")
   xtitle("X")
   ytitle("Standardized Residuals")
   yline(`tplus' `tminus', lp(dash) lw(thick))
   name(residstd3, replace) nodraw
;

graph combine residstd1 residstd2 residstd3,
   title("Standardized Residuals Plotted Against X")
   name(resid_std, replace) nocopies
;
#delimit cr


/* Plot Studentized Residuals Against X */
#delimit ;
local tplus =invttail(10-2,.05/2);
local tminus=invttail(10-2,.05/2)*-1;

graph twoway (scatter y_resid_stu  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) 
   title("Original Data")
   xtitle("X")
   ytitle("Jacknifed Residuals")
   yline(`tplus' `tminus', lp(dash) lw(thick))
   name(residstu1, replace) nodraw
;

graph twoway (scatter y1_resid_stu  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) 
   title("Scenario #1")
   xtitle("X")
   ytitle("Jacknifed Residuals")
   yline(`tplus' `tminus', lp(dash) lw(thick))
   name(residstu2, replace) nodraw
;

graph twoway (scatter y2_resid_stu  x, yline(0, lwidth(thick))),
   xlabel(0(1)10) 
   title("Scenario #2")
   xtitle("X")
   ytitle("Jacknifed Residuals")
   yline(`tplus' `tminus', lp(dash) lw(thick))
   name(residstu3, replace) nodraw
;

graph combine residstu1 residstu2 residstu3,
   title("Studentized (Jacknifed) Residuals Plotted Against X")
   ycommon
   name(resid_stu, replace) nocopies
;
#delimit cr


/*
   Test of significance for Studentized Outliers
     You need to use the Bonferroni correction
     for making multiple comparisons.

   t-critical is equal to t with (1-alpha/n*j) and
     n-j-1 df
*/
quietly {
  local studcrit1=invttail(8,.025/20)
  local studcrit2=`studcrit1'*-1
  noisily di _n _col(2) "{text}" "|t-critical|=" %7.3f `studcrit1'
  noisily list x y  y_resid_stu if  y_resid_stu > `studcrit1' |  y_resid_stu < `studcrit2'
  noisily list x y y1_resid_stu if y1_resid_stu > `studcrit1' | y1_resid_stu < `studcrit2'
  noisily list x y y2_resid_stu if y2_resid_stu > `studcrit1' | y2_resid_stu < `studcrit2'
}


/* Plot Residuals Against The Fitted Values (p. 223) */
#delimit ;
graph twoway (scatter  y_resid_org yhat1),
   yline(0, lwidth(thick))
   title("Original Data")
   ytitle("") xtitle("")
   name(yhat1, replace) nodraw
;

graph twoway (scatter y1_resid_org yhat2),
   yline(0, lwidth(thick))
   title("Scenario #1")
   ytitle("") xtitle("")
   name(yhat2, replace) nodraw
;

graph twoway (scatter y2_resid_org yhat3),
   yline(0, lwidth(thick))
   title("Scenario #2")
   ytitle("") xtitle("")
   name(yhat3, replace) nodraw
;

graph combine yhat1 yhat2 yhat3,
   title("Residuals Plotted Against Predicted Values")
   ycommon
   name(residyhat, replace) nocopies
;
#delimit cr


/* Plot Standardized Residuals Against The Fitted Values (p. 223) */
#delimit ;
local tplus =invttail(10-2,.05/2);
local tminus=invttail(10-2,.05/2)*-1;

graph twoway (scatter  y_resid_std yhat1),
   yline(0, lwidth(thick))
   yline(`tplus' `tminus', lp(dash) lw(thick))
   ylabel(-3(1)3)
   title("Original Data")
   ytitle("") xtitle("")
   name(yhat4, replace) nodraw
;

graph twoway (scatter y1_resid_std yhat2),
   yline(0, lwidth(thick))
   yline(`tplus' `tminus', lp(dash) lw(thick))
   ylabel(-3(1)3)
   title("Scenario #1")
   ytitle("") xtitle("")
   name(yhat5, replace) nodraw
;

graph twoway (scatter y2_resid_std yhat3),
   yline(0, lwidth(thick))
   yline(`tplus' `tminus', lp(dash) lw(thick))
   ylabel(-3(1)3)
   title("Scenario #2")
   ytitle("") xtitle("")
   name(yhat6, replace) nodraw
;

graph combine yhat4 yhat5 yhat6,
   title("Standardized Residuals Plotted Against Predicted Values")
   ycommon
   name(residyhatstd, replace) nocopies
;
#delimit cr


/* Plot Jacknifed Residuals Against The Fitted Values (p. 223) */
#delimit ;
local tplus =invttail(10-2,.05/2);
local tminus=invttail(10-2,.05/2)*-1;

graph twoway (scatter  y_resid_stu yhat1),
   yline(0, lwidth(thick))
   yline(`tplus' `tminus', lp(dash) lw(thick))
   title("Original Data")
   ytitle("") xtitle("")
   name(yhat7, replace) nodraw
;

graph twoway (scatter y1_resid_stu yhat2),
   yline(0, lwidth(thick))
   yline(`tplus' `tminus', lp(dash) lw(thick))
   title("Scenario #1")
   ytitle("") xtitle("")
   name(yhat8, replace) nodraw
;

graph twoway (scatter y2_resid_stu yhat3),
   yline(0, lwidth(thick))
   yline(`tplus' `tminus', lp(dash) lw(thick))
   title("Scenario #2")
   ytitle("") xtitle("")
   name(yhat9, replace) nodraw
;

graph combine yhat7 yhat8 yhat9,
   title("Studentized Residuals Plotted Against Predicted Values")
   ycommon
   name(yhatjck, replace) nocopies
;
#delimit cr


/* Plot Y Against the Predicted Values (Y-hat) */
#delimit ;
graph twoway (scatter  y yhat1),
   title("Original Data")
   ytitle("") xtitle("")
   name(yyhat1, replace) nodraw
;

graph twoway (scatter y1 yhat2),
   title("Scenario #1")
   ytitle("") xtitle("")
   name(yyhat2, replace) nodraw
;

graph twoway (scatter y2 yhat3),
   title("Scenario #2")
   ytitle("") xtitle("")
   name(yyhat3, replace) nodraw
;

graph combine yyhat1 yyhat2 yyhat3,
   title("Y Plotted Against Predicted Values")
   ycommon xcommon
   name(yyhat, replace) nocopies
;
#delimit cr


drop y_resid_org-yhat3