*--------------- NWK Table 9.4 Life Insurance Example -----------------;
options ;
data INSURE;
     input MANAGER X1 X2  LIFEINS;
     label X1       ='Annual Income'
           X2       ='Risk Aversion score'
           X1SQ     ='Income squared'
           X2Sq     ='Risk aver. squared'
           X1X2     ='X1 x X2 Interaction'
           LIFEINS  = Life Insurance Carried;
 
     X1SQ = X1 * X1;      * X1 quadratic term ;
     X2SQ = X2 * X2;      * X2     "      "   ;
     X1X2 = X1 * X2;      * X1 x X2 interaction;
datalines;
 1 66.290   7   196
 2 40.964   5    63
 3 72.996  10   252
 4 45.010   6    84
 5 57.204   4   126
 6 26.852   5    14
 7 38.122   4    49
 8 35.840   6    49
 9 75.796   9   266
10 37.408   5    49
11 54.376   2   105
12 46.186   7    98
13 46.130   4    77
14 30.366   3    14
15 39.060   5    56
16 79.380   1   245
17 52.766   8   133
18 55.916   6   133
;
proc print;
     id MANAGER;
Title 'Life Insurance Example (NWK Table 9.4)';
 
data INSURE2;
     set INSURE;
     X1BAR=50.037; X2BAR=5.389;       * Means, calculated by hand;
     drop X1BAR X2BAR;
 
     D1 = X1 - X1BAR;        * Deviations from means;
     D2 = X2 - X2BAR;        * ... reduce multicolinearity;
 
     D1SQ = D1 * D1;         * quadratic ;
     D2SQ = D2 * D2;
     D1D2 = D1 * D2;
     label D1       = Annual Income (deviation)
           D2       = Risk Aversion (deviation)
           D1D2     = Interaction;
proc print;
     id MANAGER;
     var LIFEINS X1 X2 D1 D2 D1SQ D2SQ D1D2;
Title2 'Deviation scores added to dataset';
/*
proc reg data=INSURE2;
FULLRAW:    model LIFEINS =  X1 X2 X1SQ X2SQ X1X2 / SS1 VIF;
FULLDEV:    model LIFEINS =  D1 D2 D1SQ D2SQ D1D2 / SS1 VIF;
Title2 'Full Quadratic Model with interaction term';
*/ 
data grid;
	lifeins=.; 
     do X1 = 20 to 80 by 5;
	  x1sq = x1*x1;
     do X2 =  1 to 10 by 1;
        output;
        end; end;

data insure2;
	set insure2 grid;
proc reg data=INSURE2;
REDUCRAW:   model LIFEINS =  X1 X2 X1SQ ;
   output out=results
      p=yhat r=residual rstudent=rstudent h=hatvalue;
*REDUCDEV:   model LIFEINS =  D1 D2 D1SQ / SS1 XPX I;
Title2 'Reduced model (significant terms only)';
 
data fit;
   set results;
   if lifeins = .;         *-- select predicted values on grid;
   lifeins = yhat;
proc print data=fit(obs=30);
run; 
data needles;         /* draw residuals and label points */
   set results;
   if lifeins ^=.;               *-- select observed data;
   length caseid text $12;
   caseid = ' ';
	drop caseid rstudent hatvalue manager d1--d1d2;
   if abs(rstudent) > 2 * tinv(.975, 45-4, 0)
      |   hatvalue  > 2 * 3/45 then caseid = put(manager,2.);
 
   xsys='2'; ysys='2'; zsys='2';
   x = x1;
   y = x2;
   z = yhat;     function = 'MOVE    '; output;
   if residual > 0 then do;
      line=1;  color='BLACK'; position='2'; end;
   else do;
      line=20; color='RED';   position='8'; end;
 
   z = lifeins; function = 'DRAW   '; output;        *-- needle line;
 
   function = 'SYMBOL'; text='+';      output;        *-- point;
   function = 'LABEL' ; text=caseid;   output;        *-- label;
proc print data=needles(obs=6);

%include goptions;
goptions hsize=7.5in vsize=7in  htext=1.5;
proc g3d data=fit  ;
     plot x2 * x1 = lifeins /
           anno=needles
/*           tilt= 60, 50, 40 
           rotate= -170, -170, -170 
*/           xticknum=3
           yticknum=3
           caxis=black
           ;
     format x1 x2 lifeins  5.0;
	  label x1='Inc' x2='Risk'
	  	lifeins='Ins';
%gfinish;