options ls=80 nodate pageno=1 NOFMTERR; 

 *  Class18_10.sas    ;

 
 **********************  log binomial *****************;

 * from Class20_08;

 data greenland; * from Am J Epidemiol (2004) 160:301-305;
  input receptor$ stage deaths survivors;
  total = deaths + survivors;
  stage2= (stage=2);
  stage3= (stage=3);
  cards;
 low 1 2 10
 high 1 5 50
 low 2 9 13
 high 2 17 57
 low 3 12 2
 high 3 9 6
 ;

 * proc print;
 title3;
  proc logistic data=greenland;
   class receptor(ref="high") stage(ref="1");
   model deaths/total = receptor stage receptor*stage / CLodds=PL;

  proc logistic data=greenland;
   class receptor(ref="high") stage(ref="1");
   model deaths/total = receptor stage  / CLodds=PL;* scale=deviance;

 proc genmod  data=greenland;
   title3;
   class receptor(ref=first) stage(ref=first);  * ref command doesn't work;
   model deaths/total = receptor stage / type3 dist=bin link=log;*scale=deviance;
   estimate "receptor low vs. high" receptor -1 1 / exp E;
   estimate "stage2 vs stage1" stage -1 1 0 / exp E;
   estimate "stage 3 vs stage1" stage -1 0 1/exp;
   ODS output ParameterEstimates = reg_coef;

 proc print data=reg_coef;

   run; quit;
*******************ordinal regression******************;

 
    * from  Lect21B_08.sas ;

 * sets up NHANES age - BMI dataset for logistic regression;
 * logistic regression and ordinal regression;


/* Proc SGplot data=ph6470.age_bmi;*/
/*   loess y = obese x = age ;*/
/*  run; */
/* Proc SGplot data=ph6470.age_bmi;*/
/*   loess y = overweight x = age ;*/
/*  run; */

 /*
    data ph6470.age_bmi_sample;
   set ph6470.age_bmi;
   age_category = "young";
   if (age > 39) then age_category = "old";
  bmi_cat = "1_normal";
  if (25.0 LE bmi < 30.0) then bmi_cat = "2_overwt";
  if (bmi GE 30.0) then bmi_cat = "3_obese";
  gender = "female";
  if (female=0) then gender="male";


proc freq data=ph6470.age_bmi_sample;
  tables bmi_cat ;
  tables age_category*(bmi_cat obese overweight) /nocol nopercent relrisk;


proc logistic  data=ph6470.age_bmi_sample;
  class age_category gender  /param=glm;
*  model obese(event='1') = age_category  gender age_category * gender /  CLodds = PL;
*  model obese(event='1') = age_category  gender;
  model obese(event='1') = age_category;
*  output out=z predicted = phat xbeta = logit stdxbeta = logitse;

 proc logistic  data=ph6470.age_bmi_sample;
  class age_category gender  /param=glm;
 * model overweight(event='1') = age_category  gender;
  model overweight(event='1') = age_category  ;


 proc logistic descending data=ph6470.age_bmi_sample;
  class age_category gender  /param=glm;
  model bmi_cat = age_category ;


 proc logistic descending data=ph6470.age_bmi_sample;
  class age_category gender  /param=glm;
  model bmi_cat = age_category  gender ;

 proc logistic descending data=ph6470.age_bmi_sample;
  class age_category gender  /param=glm;
  model bmi_cat = age_category  gender  / link=glogit  CLodds = PL;

proc logistic descending data=ph6470.age_bmi_sample;
  class age_category gender  /param=glm;
  model bmi_cat = age_category | gender / link=glogit  CLodds = PL;
  run; quit;

run; quit;
 *  another working of this example with decade instead of old/young
    but left as is in the notes;
 Proc Freq data=ph6470.age_bmi;
   tables bmi_cat * decade / nopercent norow;

 Proc Logistic descending data=ph6470.age_bmi;
  class decade (ref=first);*  /param=glm;
  model bmi_cat = decade ;


 proc logistic descending data=ph6470.age_bmi;
  class decade (ref=first) female (ref=first);* /param=glm;
  model bmi_cat = decade  female ;

 run; quit;/*