
LOG TO COMPARE PRELIMINARY OUTPUTS FROM LOGISTIC REGRESSION AND USUAL
LINEAR REGRESSION FOR POLYNESIAN ISLANDS DEFORESTATION DATA
========================================================== 12/8/08

libname SAS "SASdat";
options nocenter linesize=75 nodate ;

proc contents data=SAS.polynest; run;


Alphabetic List of Variables and Attributes

 #    Variable    Type    Len

13    Age         Num       8
 2    Archip      Char      8
 6    Area        Num       8
 5    Deforstn    Num       8
15    Dust        Num       8
10    Elev        Num       8
 4    ForstRep    Num       8
 1    ID          Num       8
 3    Islnd       Char      8
 8    Isol_25     Num       8
 9    Isol_75     Num       8
11    LatitS      Num       8
 7    NrArea      Num       8
14    PctMakat    Num       8
12    Rain        Num       8
16    Tephra      Num       8


* Recall response variable Deforstn takes values 2..5 , 
  38% with values 3.5, 4, 5 with one missing which 
  we delete;

data tmpPoly;
   set SAS.polynest;
   if Deforstn ne . ;
   HIdefrst = int(Deforstn/3);
   abslat = abs(LatitS);
run;                   /* 79 obs, 18 varbls */

proc reg data=tmppoly;
   model deforstn = age area dust elev isol_25 abslat
              pctmakat rain tephra / selection=forward
              sle = .05 aic r ;
   run;

* Generated LOTS of output. ; 

                     Summary of Forward Selection

      Variable   Number   Partial    Model
Step  Entered    Vars In  R-Square  R-Square   C(p)    F Value  Pr > F

  1   Rain           1     0.3940    0.3940   54.0527    47.45  <.0001
  2   abslat         2     0.1013    0.4952   35.1533    14.45  0.0003
  3   Area           3     0.1068    0.6020   15.1229    19.05  <.0001
  4   Age            4     0.0350    0.6370    9.9063     6.74  0.0115
  5   Isol_25        5     0.0269    0.6639    6.3535     5.52  0.0216

             Parameter     Standard
Variable      Estimate        Error   Type II SS  F Value  Pr > F

Intercept      1.20003      0.51592      1.39680     5.41  0.0230
Age            0.48832      0.14511      2.92364    11.32  0.0013
Area       -0.00001318   0.00000308      4.74569    18.38  <.0001
Isol_25     0.00029317   0.00012473      1.42626     5.52  0.0216
abslat         0.06910      0.01086     10.45083    40.48  <.0001
Rain       -0.00020982   0.00005244      4.13380    16.01  0.0002


proc reg data=tmppoly;
    model deforstn = elev isol_25 abslat rain tephra;
    plot residual.*predicted.;
    plot obs. * pred. ;
   run;

/* Now try for some interactions */

data tmppoly2;
   set tmppoly;
    latElv = abslat*elev;  latIsol = abslat*isol_25;
    latRain= abslat*rain;  latTeph = abslat*tephra;
    elvIso = elev*isol_25;  elvRain = elev*rain;
    elvTeph = elev*tephra; IsoRain = isol_25*rain;
    rainTeph = rain*tephra; IsoTep = isol_25*tephra;

/* some were missing ! */

proc reg data=tmppoly2;
    model deforstn = elev isol_25 abslat rain tephra
     latElv latIsol latRain latTeph elvIso elvRain 
     elvTeph IsoRain rainTeph IsoTep / 
       selection = forward sle = .05 aic r include=5;
   run;
                Parameter     Standard
  Variable      Estimate        Error   Type II SS  F Value  Pr > F

  Intercept      2.13987      0.52294      5.49853    16.74  0.0001
* Elev       -0.00011512   0.00007756      0.72347     2.20  0.1422
* Isol_25     0.00012537   0.00014696      0.23899     0.73  0.3965
* abslat         0.11586      0.02730      5.91416    18.01  <.0001
* Rain       -0.00025990   0.00006320      5.55294    16.91  0.0001
* Tephra         0.39122      0.22040      1.03458     3.15  0.0802
  latTeph       -0.03393      0.01136      2.92652     8.91  0.0039

           * Forced into the model by the INCLUDE= option


The REG Procedure
Model: MODEL1
Dependent Variable: Deforstn

Forward Selection: Step 1

Bounds on condition number: 13.756, 220.06
---------------------------------------------------------------------------


No other variable met the 0.0500 significance level for entry into the
model.


                     Summary of Forward Selection

      Variable   Number   Partial    Model
Step  Entered    Vars In  R-Square  R-Square   C(p)    F Value  Pr > F

  1   latTeph        6     0.0512    0.5921    1.4495     8.91  0.0039


/* Now try the same with logistic regression, for HIdefrst */

proc logistic data=tmppoly descending;
   model HIdefrst = age area dust elev isol_25 abslat
              pctmakat rain tephra / selection=forward
              sle = .05  ;  run;


          Response Profile

 Ordered                      Total
   Value     HIdefrst     Frequency

       1            1            58
       2            0            17

Probability modeled is HIdefrst=1.


                    Summary of Forward Selection

            Effect                Number         Score
    Step    Entered       DF          In    Chi-Square    Pr > ChiSq

       1    Rain           1           1       25.3367        <.0001
       2    Tephra         1           2       13.0101        0.0003


             Analysis of Maximum Likelihood Estimates

                               Standard          Wald
Parameter    DF    Estimate       Error    Chi-Square    Pr > ChiSq

Intercept     1      9.3746      2.8857       10.5537        0.0012
Rain          1    -0.00118    0.000380        9.6518        0.0019
Tephra        1     -2.0578      0.7917        6.7551        0.0093


            Odds Ratio Estimates

               Point          95% Wald
Effect      Estimate      Confidence Limits

Rain           0.999       0.998       1.000
Tephra         0.128       0.027       0.603

             Partition for the Hosmer and Lemeshow Test

                            HIdefrst = 1            HIdefrst = 0
   Group       Total    Observed    Expected    Observed    Expected

       1           8           2        1.40           6        6.60
       2           8           1        2.91           7        5.09
       3           8           7        4.65           1        3.35
       4           8           5        6.46           3        1.54
       5           8           8        7.70           0        0.30
       6           8           8        7.95           0        0.05
       7           8           8        7.96           0        0.04
       8           8           8        7.98           0        0.02
       9          11          11       10.99           0        0.01


Hosmer and Lemeshow Goodness-of-Fit Test

Chi-Square       DF     Pr > ChiSq

    7.2906        7         0.3993

================================================

ILLUSTRATION OF THE OUTPUT STATEMENT AND KEYWORDS 
FOR USE IN PRODUCING RESIDUALS PLOTS IN LOGISTIC 
REGRESSION.


proc logistic data=tmppoly descending;
   model HIdefrst = rain tephra ; 
   OUTPUT out= tmpout  pred = defrPred reschi = defrResd;
  run;

/* The syntax here allows us to specify keywords for the 
LOGISTIC REGRESSION predictors probabilities, and also 
for some kinds of residuals.

The "reschi" keyword specifies the "Pearson" residual which
is a lot like the studentized residual (keyword "stud") 
in PROC REG output statements. However, to produce "raw" 
residuals OBS - PRED, we have to use a separate DATA step. */

data tmpout ;
   set tmpout;
   RAWRES = Hidefrst - defrpred; run;

symbol v=circle color=black;

proc gplot data=tmpout;
   plot RAWRES * defrpred; 
   plot defrResd * defrpred; run;

/* The residual plot based on Pearson residuals are better
  to look at than the plot for raw residuals, but still not really
  interpretable ! */
