/*********************************************************************************** Happiness Data from Agresti (2013) http://www.stat.ufl.edu/~aa/cda/data.html which is from GSS. Table 8.5, example of ordinal logistic regression ************************************************************************************/ data gss; input race trauma happy ; /*where race is 0=white, 1=black */ datalines; 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 3 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 3 0 2 1 0 2 1 0 2 1 0 2 1 0 2 1 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 2 0 2 3 0 3 1 0 3 2 0 3 2 0 3 2 0 3 2 0 3 2 0 3 2 0 3 2 0 3 2 0 3 2 0 3 3 0 4 1 0 4 2 0 4 2 0 4 2 0 4 2 0 5 3 0 5 3 1 0 2 1 0 3 1 1 2 1 1 2 1 1 2 1 1 3 1 2 2 1 2 2 1 2 2 1 2 3 1 3 2 1 3 2 1 3 3 ; run; /*************************************************************** Proportional Odds model ***************************************************************/ ods graphics on; ods html; title 'Different ways to fit proportional odds models'; proc genmod; class race; model happy = trauma race / dist=multinomial link=clogit lrci type3; effectplot slicefit; run; proc logistic data=gss plots=effect(at(race=all)); class race; model happy = trauma race ; output out=predict p=hi_hat lower=LCL upper=UCL; proc print data=predict; run; proc logistic data=gss plots=all; class race; model happy = trauma race; run; ods html close; ods graphics off; run; /********************************************************** Adjacent Categories ***********************************************************/ title 'Check for zeros'; proc freq data=gss; tables race*trauma*happy / nopercent norow nocol sparse out=table; run; data fillin; set table; count2=count+ .01; urace = (3-happy)*race; utrauma=(3-happy)*trauma; run; title 'Adjacent Categories (WLS)--- had to add small number due to 0s'; proc catmod data=fillin; weight count2; response alogits; population race trauma; direct trauma race ; model happy = _response_ race trauma ; run; title 'Adjacent Categories (WLS)--- had to add small number due to 0s'; proc catmod data=fillin; weight count2; response alogits; population race trauma; direct trauma race ; model happy = race trauma ; run; quit; title 'Adjacent Categories (MLE)-- tabled data'; proc nlmixed data=fillin; parms a1-a2=.1 br=.1 bt=.1; /* Linear predictors */ eta1 = a1 + br*(3-1)*race + bt*(3-1)*trauma; eta2 = a2 + br*(3-2)*race + bt*(2-1)*trauma; /* Define likelihood */ if happy=1 then prob= exp(eta1)/(1 + exp(eta1) + exp(eta2)); if happy=2 then prob= exp(eta2)/(1 + exp(eta1) + exp(eta2)); if happy=3 then prob= 1 /(1 + exp(eta1) + exp(eta2)); /* To make sure that probabilities are valid ones */ p = (prob>0 and prob<=1)*prob + (prob<=0)*1e-8 + (prob>1); mle = log(p); /* Specify distribution for response variable and random effect */ model happy ~ general(mle); replicate count; run; title 'Adjacent Categories (MLE)-- un-collapsed data'; proc nlmixed data=gss; parms a1-a2=.1 br=.1 bt=.1; /* Linear predictors */ eta1 = a1 + br*(3-1)*race + bt*(3-1)*trauma; eta2 = a2 + br*(3-2)*race + bt*(2-1)*trauma; /* Define likelihood */ if happy=1 then prob= exp(eta1)/(1 + exp(eta1) + exp(eta2)); if happy=2 then prob= exp(eta2)/(1 + exp(eta1) + exp(eta2)); if happy=3 then prob= 1 /(1 + exp(eta1) + exp(eta2)); /* To make sure that probabilities are valid ones */ p = (prob>0 and prob<=1)*prob + (prob<=0)*1e-8 + (prob>1); mle = log(p); /* Specify distribution for response variable and random effect */ model happy ~ general(mle); run;