*---- NWK 2nd Ed. Table 9.4 Life Insurance Example -----;
options ls=90 ps=62;
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';
 
proc reg data=INSURE2;
*REDUCRAW:   model LIFEINS =  X1 X2 X1SQ ;
REDUCDEV:   model LIFEINS =  D1 D2 D1SQ / SS1;
Title2 'Reduced model (significant terms only)';
 
* Proc GLM gives sequential (Type I) sums of squares
           for testing extra sum of squares;
 
proc GLM data=INSURE2;
         model LIFEINS =  D1 D2 D1*D1 D2*D2 D1*D2/ SS1 ;
Title2 'Reduced mode fit by GLM';
 
proc RSREG data=INSURE;
     model LIFEINS = X1 X2 / noopt;
Title2 'Proc RSREG output';
run;
* The fitted response needs to be displayed in a 3-dimensional
  plot. The steps below show two ways to do this;
 
* Calcualte fitted response (Y) values over a rectangular grid;
data FIT;
     do X1 = 20 to 80 by 2;
     do X2 =  1 to 10 by .5;
        * Calculate Yhat;
        Y  = -74.583 + .8383 * X1 + 5.685 * X2  + .0371 * X1**2;
        output;
        end; end;
Title2 'Response surface';
proc plot;
     plot X2 * X1 = Y / contour=8;
     run;
*--- Response surface, using proc RSP;
/*
proc RSP;
     model Y = -74.583 + .8383 * X1 + 5.685 * X2 + .0371 * X1**2;
     ranges X1 = 20 to 80 by 2
            X2 =  1 to 10 by .5
            Y  =  -40 to 360 by 20;
     plot X1 * X2;
     label X1=Annual Income
           X2=Risk Aversion;
*/