/* altitude versus east/north .. */

%let which=%str( ABS(east - north) <= 1 );

*%let mean = %str(85 +2*east - 5*(north-4)**2  );

%let mean = %str(10 +2*east +8*north );

OPTIONS LS=75 PS=45 NODATE ; RUN:

DATA a;
Do east = 0 to 8;
  Do north =  0 to 8;
    altitude = &mean + 0.001* rannor(5434567);
    alt_10 = round(altitude,10);
   IF &which then output;
  end;
end;

TITLE Raw data... ;
PROC PRINT DATA = a;  RUN;

TITLE Raw data: altitude rounded to nearest 10;
TITLE2 i.e., a 3 means 25-35,  5 means 45-55,  etc.. ;
TITLE3  ;
TITLE4 ;
TITLE5 Correct equation: 10 + 2*east + 8*north + (a little noise) ;

PROC PLOT data=a; Plot north*east=alt_10 /
      vpos= 20 hpos = 25 ; RUN;
* ---------------------------------------------------------- ;

TITLE mean altitude at each distance east (ignoring how far north);
PROC MEANS DATA = a MEAN MAXDEC=1 ; VAR altitude ;  CLASS east  ;  RUN;

TITLE mean altitude at each distance north (ignoring how far east);
PROC MEANS DATA = a MEAN MAXDEC=1; VAR altitude ;  CLASS north ; RUN;

* ---------------------------------------------------------- ;

TITLE Regress altitude on east ...;
TITLE2   ;
TITLE3 and save residuals as alt_e ;

PROC REG DATA = a ; MODEL altitude = east  ;
 OUTPUT out = a1 residual = alt_e ; RUN;

* ---------------------------------------------------------- ;

TITLE Then regressing alt_e on north ...;
TITLE2   ;
TITLE3 gives us INCORRECT picture of process.. ;

PROC REG DATA = a1 ; MODEL alt_e = north  ; RUN;

* ---------------------------------------------------------- ;

TITLE The REASON is that we are re-using north ;
TITLE2 But in the dataset north is also partly east    ;
TITLE3 So.. calculate  north, NET of east;
TITLE4 i.e., regress north on east, and using what remains of it ;
TITLE5  ..... denote this residual (of north on east) by north_e ;

PROC REG DATA = a1 ; MODEL north = east  ;
 OUTPUT out = a2 residual = north_e ; RUN;

* ---------------------------------------------------------- ;

TITLE NOW , regress alt_e on NET north (after removing east);
TITLE2 and rewrite in fitted model in terms of original variables;

PROC REG DATA = a2 ; MODEL alt_e = north_e  ; RUN;

* ---------------------------------------------------------- ;

options ls=75 ps=40;  RUN;

TITLE CHECK against model where east and north are fitted simultaneously..;
TITLE3 and out out interest, request PARTIAL REGRESSION LEVERAGE PLOTS;
PROC REG DATA = a ; MODEL altitude = east  north / PARTIAL  ; RUN;

* ---------------------------------------------------------- ;