title 'Logistic Regression of Piceance Mule Deer Fawns';
proc format;
   value status 0='Died' 1='Lived';
data fawns;
   *change the input file statement to the e drive (e:\fawns.dat);
   infile 'w:\courses\wbio595\solutions\lab11_knownfate\fawns.dat';
   input area $ 1-9 freq 11-17 sex $ 19-24
         weight 26-29 length 31-35 time 37-39;
   if time = . then status=1;
   else status=0;
   one=1;
   format status status.;
   label area='Treatment or Control area'
         freq='Individual id -- radio frequency'
         sex='Sex of fawn'
         weight='Weight in kg'
         length='Length in cm'
         time='Days until death'
         status='Lived or Died';

proc sort;
   by sex length;

proc genmod;
   class sex area;
   model status/one=sex weight length area /
      link=logit
      dist=binomial
      type3
      lrci
      dscale;
   title2 'Logistic regression with all variables in the model, SCALE est.';

proc genmod;
   class sex;
   model status/one=sex length /
      link=logit
      dist=binomial
      type3
      lrci
      obstats          ;
   title2 'Reduced model with only sex and length';
   make 'obstats' out=obstats;

data obstats;
   merge fawns obstats;

proc print data=obstats;
   run;

*goptions device=hpljs2;
goptions device=win;
goptions ftext=swissb;
symbol1 color=blue v=dot i=none;
symbol2 color=yellow i=join line=1 width=3;
symbol3 color=yellow i=join line=4 width=2;

proc gplot data=obstats;
   by sex;
   plot pred*length=2 status*length=1 lower*length=3 upper*length=3 / overlay;
   label pred='Pr(Lived)';
   run;