/* From Logistic Regression Examples A55201 */
/* Example 12, pp. 100-109 */
%version(6.10);
*include goptions;
title 'Diabetes data: logistic regression';
goptions cback=white colors=(black) htext=1.3 vsize=8in hsize=7in;

proc format;
   value gp 3='(3) Overt Diabetic ' 2='(2) Chem. Diabetic' 1='(1) Normal';
run;

data diabetes;
   infile cards eof=endfile;
   input patient relwt glufast glutest instest sspg group;
   label relwt   = 'Relative weight'
         glufast = 'Fasting Plasma Glucose'
         glutest = 'Test Plasma Glucose'
         instest = 'Plasma Insulin during Test'
         sspg    = 'Steady State Plasma Glucose'
         group   = 'Clinical Group';
   output;
   return;
   endfile:  do glutest=250 to 1600 by 25;
                output;
             end;
   cards;
   1  0.81  80  356 124   55  1
   2  0.95  97  289 117   76  1
   3  0.94 105  319 143  105  1
   4  1.04  90  356 199  108  1
   5  1.00  90  323 240  143  1
   6  0.76  86  381 157  165  1
   7  0.91 100  350 221  119  1
   8  1.10  85  301 186  105  1
   9  0.99  97  379 142   98  1
  10  0.78  97  296 131   94  1
  11  0.90  91  353 221   53  1
  12  0.73  87  306 178   66  1
  13  0.96  78  290 136  142  1
  14  0.84  90  371 200   93  1
  15  0.74  86  312 208   68  1
  16  0.98  80  393 202  102  1
  17  1.10  90  364 152   76  1
  18  0.85  99  359 185   37  1
  19  0.83  85  296 116   60  1
  20  0.93  90  345 123   50  1
  21  0.95  90  378 136   47  1
  22  0.74  88  304 134   50  1
  23  0.95  95  347 184   91  1
  24  0.97  90  327 192  124  1
  25  0.72  92  386 279   74  1
  26  1.11  74  365 228  235  1
  27  1.20  98  365 145  158  1
  28  1.13 100  352 172  140  1
  29  1.00  86  325 179  145  1
  30  0.78  98  321 222   99  1
  31  1.00  70  360 134   90  1
  32  1.00  99  336 143  105  1
  33  0.71  75  352 169   32  1
  34  0.76  90  353 263  165  1
  35  0.89  85  373 174   78  1
  36  0.88  99  376 134   80  1
  37  1.17 100  367 182   54  1
  38  0.85  78  335 241  175  1
  39  0.97 106  396 128   80  1
  40  1.00  98  277 222  186  1
  41  1.00 102  378 165  117  1
  42  0.89  90  360 282  160  1
  43  0.98  94  291  94   71  1
  44  0.78  80  269 121   29  1
  45  0.74  93  318  73   42  1
  46  0.91  86  328 106   56  1
  47  0.95  85  334 118  122  1
  48  0.95  96  356 112   73  1
  49  1.03  88  291 157  122  1
  50  0.87  87  360 292  128  1
  51  0.87  94  313 200  233  1
  52  1.17  93  306 220  132  1
  53  0.83  86  319 144  138  1
  54  0.82  86  349 109   83  1
  55  0.86  96  332 151  109  1
  56  1.01  86  323 158   96  1
  57  0.88  89  323  73   52  1
  58  0.75  83  351  81   42  1
  59  0.99  98  478 151  122  2
  60  1.12 100  398 122  176  1
  61  1.09 110  426 117  118  1
  62  1.02  88  439 208  244  2
  63  1.19 100  429 201  194  2
  64  1.06  80  333 131  136  1
  65  1.20  89  472 162  257  2
  66  1.05  91  436 148  167  2
  67  1.18  96  418 130  153  1
  68  1.01  95  391 137  248  1
  69  0.91  82  390 375  273  1
  70  0.81  84  416 146   80  1
  71  1.10  90  413 344  270  2
  72  1.03 100  385 192  180  1
  73  0.97  86  393 115   85  1
  74  0.96  93  376 195  106  1
  75  1.10 107  403 267  254  1
  76  1.07 112  414 281  119  1
  77  1.08  94  426 213  177  2
  78  0.95  93  364 156  159  1
  79  0.74  93  391 221  103  1
  80  0.84  90  356 199   59  1
  81  0.89  99  398  76  108  1
  82  1.11  93  393 490  259  1
  83  1.19  85  425 143  204  2
  84  1.18  89  318  73  220  1
  85  1.06  96  465 237  111  2
  86  0.95 111  558 748  122  2
  87  1.06 107  503 320  253  2
  88  0.98 114  540 188  211  2
  89  1.16 101  469 607  271  2
  90  1.18 108  486 297  220  2
  91  1.20 112  568 232  276  2
  92  1.08 105  527 480  233  2
  93  0.91 103  537 622  264  2
  94  1.03  99  466 287  231  2
  95  1.09 102  599 266  268  2
  96  1.05 110  477 124   60  2
  97  1.20 102  472 297  272  2
  98  1.05  96  456 326  235  2
  99  1.10  95  517 564  206  2
 100  1.12 112  503 408  300  2
 101  0.96 110  522 325  286  2
 102  1.13  92  476 433  226  2
 103  1.07 104  472 180  239  2
 104  1.10  75  455 392  242  2
 105  0.94  92  442 109  157  2
 106  1.12  92  541 313  267  2
 107  0.88  92  580 132  155  2
 108  0.93  93  472 285  194  2
 109  1.16 112  562 139  198  2
 110  0.94  88  423 212  156  2
 111  0.91 114  643 155  100  2
 112  0.83 103  533 120  135  2
 113  0.92 300 1468  28  455  3
 114  0.86 303 1487  23  327  3
 115  0.85 125  714 232  279  3
 116  0.83 280 1470  54  382  3
 117  0.85 216 1113  81  378  3
 118  1.06 190  972  87  374  3
 119  1.06 151  854  76  260  3
 120  0.92 303 1364  42  346  3
 121  1.20 173  832 102  319  3
 122  1.04 203  967 138  351  3
 123  1.16 195  920 160  357  3
 124  1.08 140  613 131  248  3
 125  0.95 151  857 145  324  3
 126  0.86 275 1373  45  300  3
 127  0.90 260 1133 118  300  3
 128  0.97 149  849 159  310  3
 129  1.16 233 1183  73  458  3
 130  1.12 146  847 103  339  3
 131  1.07 124  538 460  320  3
 132  0.93 213 1001  42  297  3
 133  0.85 330 1520  13  303  3
 134  0.81 123  557 130  152  3
 135  0.98 130  670  44  167  3
 136  1.01 120  636 314  220  3
 137  1.19 138  741 219  209  3
 138  1.04 188  958 100  351  3
 139  1.06 339 1354  10  450  3
 140  1.03 265 1263  83  413  3
 141  1.05 353 1428  41  480  3
 142  0.91 180  923  77  150  3
 143  0.90 213 1025  29  209  3
 144  1.11 328 1246 124  442  3
 145  0.74 346 1568  15  253  3
;

proc logistic data=diabetes descending;
   model group=glutest;
   output out=probs predicted=prob xbeta=logit;
   format group gp.;
run;

proc print data=probs (obs=10);
   var patient group glutest _level_ prob;
   title1 'Printout of first 10 observations of PROBS';
run;

title1;

proc sort data=probs;
   by glutest;
run;

proc gplot data=probs;
   plot logit*glutest=_level_ / frame;
   label logit='Logit (Diabetic vs. Normal)';
   title2 'Plot of Parallel Logits';
run;
   plot prob*glutest=_level_ / frame;
   title2 'Plot of Parallel Probability Lines';
run;

data overt(rename=(prob=p3)) chem(rename=(prob=p2));
   set probs;
   if mod(_n_,2)=1 then output overt;
   else output chem;
run;

data probs2;
   merge overt chem;
   keep glutest prob_3 prob_2 prob_1;
   prob_3=p3;
   prob_2=p2-p3;
   prob_1=1-p2;
run;

proc sort data=probs2;
   by glutest;
run;

symbol1 i=join line=1 value=x height=.75;
symbol2 i=join line=2 value=plus height=.75;
symbol3 i=join line=3 value=circle height=.75;

axis1 label=(angle=-90 rotate=90 'Probability');
axis2 label=('Test Plasma Glucose');

title2 'Individual Density Plots';

footnote1
   height=.75 font=zapf 'x x x ' h=1.0 font=zapf ' Overt '
   height=.75 font=zapf '+ + + ' h=1.0 font=zapf ' Chemical '
   height=.75 font=special 'H H H ' h=1.0 font=zapf ' Normal';

proc gplot data=probs2;
   plot prob_3*glutest=1
        prob_2*glutest=2
        prob_1*glutest=3
        / overlay vaxis=axis1 haxis=axis2 frame;
run;
footnote1;

/* pp. 107 */
proc logistic data=diabetes descending outest=parms;
   model group=glutest;
   format group gp.;
run;

data parms;
   set parms;
   rename glutest=b1;
run;

data parms2;
   set parms;
   do glutest=250 to 1600 by 25;
      logit1=intercp1 + b1*glutest;
      logit2=intercp2 + b1*glutest;
      p1=exp(logit1)/(1+exp(logit1));
      p2=exp(logit2)/(1+exp(logit2));
      prob_1=p1;
      prob_2=p2-p1;
      prob_3=1-p2;
      output;
   end;
run;

proc gplot data=parms2;
   plot prob_3*glutest
        prob_2*glutest
        prob_1*glutest
        / overlay frame;
run;

/* pp. 108-109 */
proc catmod data=diabetes;
   direct glutest;
   response logits / out=cat_prob;
   model group=glutest;
run;