options pagesize=150 linesize=68;
filename grafout 'g09.ps';
goptions device=pslmono gsfname=grafout 
   rotate=landscape gsfmode=append;
/*****************************************************/
/* Housing data from Cox and Snell (1980), Example W,*/
/* on satisfaction with housing in Copenhagen.  Data */
/* was extracted from R package MASS.  Data may also */
/* be found at                                       */
/* www.stat.ucla.edu/data/cox-and-snell/exampleW.data*/
/*****************************************************/
data housing; infile 'housing.dat'; 
   input nn sat $ infl $ type $ contact $ freq ;
   nsat=0;if sat="Medium" then nsat=1; 
   if sat="High" then nsat=2;
   ninf=0;if infl="Medium" then ninf=1; 
   if infl="High" then ninf=2;
   ncont=0;if contact="Medium" then ncont=1; 
   if contact="High" then ncont=2;
   run;
data housing; set housing; rsat=sat;
   if sat="Medium" then rsat="Intermediate"; run;
/**************************************************************/
/* Mark C: Polytomous Regression.  Baseline logit model first.*/
/**************************************************************/
*We have seen logistic regression, with two categories, via    ;
*proc genmod.  Proc logistic provides an alternative tool.     ;
*We now will model categorical responses with multiple catego- ;
*ries.  Proc logistic (to be revisited) does something related.;
*proc catmod parameterizes the model in such a way as to make  ;
*the sum of parameters equal zero.  In comparison with usual   ;
*dummy coding, parameter estimates are half of what they were. ;
title 'Baseline logit model';
proc catmod data=housing; weight freq;
   model rsat=contact; run;
/**************************************************************/
/* Mark D: Cumulative logit model.  First reuse a known tech- */
/* nique.  Simple model exhibits different parameterizations. */
/**************************************************************/
title 'Binary response for comparison purpose';
data housing; set housing; oksat=min(1,nsat); run;
* For catmod represent multiple observations using weight.;
* For genmod and logistic use freq.;
proc logistic data=housing; class contact ;
   model oksat=contact ; freq freq;  run;
proc genmod data=housing; class contact ; freq freq;
   model oksat=contact /dist=binomial; run;
/**************************************************************/
/* The following data set will not give us what we want, be-  */
/* cause it will order sat incorrectly.  Recode data to put   */
/* medium between high and low, and add descending to make    */
/* low the baseline.  First fit the simplest model.  This     */
/* should give logits of raw proportions.  Check this.        */
/**************************************************************/
proc freq data=housing; table  rsat; weight freq; run;
title 'Ordinal Logistic Regression with wrong order';
proc logistic data=housing; class sat ;
   model sat=; freq freq; run;
title 'Ordinal Logistic Reg. with right order, empty model';
proc logistic data=housing descending;
   class rsat ; model rsat=; freq freq; run;
title 'Ordinal Logist. Reg. with added covariates';
proc logistic data=housing descending;
   class rsat infl type contact; freq freq;
   model rsat=infl type contact; run;
/**************************************************************/
/* Mark E: Cf. the complimentary log log.  Odds ratios are no */
/* longer given, since they don't come naturally out of the   */
/* model.  Intercepts have a completely different meaning and */
/* won't be compared.  Linear parameters are multiplied by a  */
/* constant, approx., and z statistics are similar.           */
/**************************************************************/
title 'Complimentary log log';
proc logistic data=housing descending;
   class rsat infl type contact; freq freq;
   model rsat=infl type contact/link=cloglog; run;