goptions vsize=7.5in htext=1.8;
title 'Logistic Regression: Proportional Odds Model';
%include catdata(arthrit);
 
proc logistic data=arthrit descending;
   class sex (ref=last) treat (ref=first) / param=ref;
   model  improve = sex  treat  age ;
   output out=results p=prob l=lower u=upper
          xbeta=logit stdxbeta=selogit / alpha=.33;

proc print data=results(obs=6);
   id id treat sex;
   var improve _level_ prob lower upper logit;
   format prob lower upper logit selogit 6.3;
run;

*-- combine treatment and _level_, set bar color;
data results;
   set results;
   treatl = trim(treat)||put(_level_,1.0); 
   if treat='Placebo' then col='BLACK';
                      else col='RED';
proc sort data=results;
   by sex treatl age;

*-- Error bars, on prob scale;
%bars(data=results, var=prob, 
	class=age, cvar=treatl, by=age,
	lower=lower, upper=upper, 
    color=col, out=bars);
proc sort data=bars;
   by sex treatl age;

*-- Custom legends, for treat-level and sex;
%label(data=results, y=prob, x=age, xoff=1, cvar=treatl,
	by=sex, subset=last.treatl, out=label1, pos=6, text=treatl);
%label(data=results, y=0.9, x=20, size=2,
	by=sex, subset=first.sex, out=label2, pos=6, text=sex);

*-- Combine the annotate data sets;
data bars;
    set label1 label2 bars;
	by sex;

goptions hby=0;
proc gplot data=results;
   plot prob * age = treatl / 
       vaxis=axis1 haxis=axis2 hminor=1 vminor=1
       nolegend anno=bars name='glogist2a';
   by sex;
   axis1 label=(a=90 'Prob. Improvement (67% CI)')
         order=(0 to 1 by .2);
   axis2 order=(20 to 80 by 10)
         offset=(2,5);
   symbol1 v=circle  i=join l=3 c=black;
   symbol2 v=circle  i=join l=3 c=black;
   symbol3 v=dot     i=join l=1 c=red;
   symbol4 v=dot     i=join l=1 c=red;
run;