*include goptions; goptions vsize=7.5in horigin=.5in; title h=1.6 'Logistic Regression: Proportional Odds Model'; data arthrit; length treat $7. sex $6. ; input id treat $ sex $ age improve @@ ; better = (improve > 0); _treat_ = (treat ='Treated') ; _sex_ = (sex = 'Female'); agesex = age*_sex_ ; agetrt = age*_treat_; sextrt = _sex_*_treat_; age2 = age*age ; label age='Age' sex='Sex'; cards ; 57 Treated Male 27 1 9 Placebo Male 37 0 46 Treated Male 29 0 14 Placebo Male 44 0 77 Treated Male 30 0 73 Placebo Male 50 0 17 Treated Male 32 2 74 Placebo Male 51 0 36 Treated Male 46 2 25 Placebo Male 52 0 23 Treated Male 58 2 18 Placebo Male 53 0 75 Treated Male 59 0 21 Placebo Male 59 0 39 Treated Male 59 2 52 Placebo Male 59 0 33 Treated Male 63 0 45 Placebo Male 62 0 55 Treated Male 63 0 41 Placebo Male 62 0 30 Treated Male 64 0 8 Placebo Male 63 2 5 Treated Male 64 1 80 Placebo Female 23 0 63 Treated Male 69 0 12 Placebo Female 30 0 83 Treated Male 70 2 29 Placebo Female 30 0 66 Treated Female 23 0 50 Placebo Female 31 1 40 Treated Female 32 0 38 Placebo Female 32 0 6 Treated Female 37 1 35 Placebo Female 33 2 7 Treated Female 41 0 51 Placebo Female 37 0 72 Treated Female 41 2 54 Placebo Female 44 0 37 Treated Female 48 0 76 Placebo Female 45 0 82 Treated Female 48 2 16 Placebo Female 46 0 53 Treated Female 55 2 69 Placebo Female 48 0 79 Treated Female 55 2 31 Placebo Female 49 0 26 Treated Female 56 2 20 Placebo Female 51 0 28 Treated Female 57 2 68 Placebo Female 53 0 60 Treated Female 57 2 81 Placebo Female 54 0 22 Treated Female 57 2 4 Placebo Female 54 0 27 Treated Female 58 0 78 Placebo Female 54 2 2 Treated Female 59 2 70 Placebo Female 55 2 59 Treated Female 59 2 49 Placebo Female 57 0 62 Treated Female 60 2 10 Placebo Female 57 1 84 Treated Female 61 2 47 Placebo Female 58 1 64 Treated Female 62 1 44 Placebo Female 59 1 34 Treated Female 62 2 24 Placebo Female 59 2 58 Treated Female 66 2 48 Placebo Female 61 0 13 Treated Female 67 2 19 Placebo Female 63 1 61 Treated Female 68 1 3 Placebo Female 64 0 65 Treated Female 68 2 67 Placebo Female 65 2 11 Treated Female 69 0 32 Placebo Female 66 0 56 Treated Female 69 1 42 Placebo Female 66 0 43 Treated Female 70 1 15 Placebo Female 66 1 71 Placebo Female 68 1 1 Placebo Female 74 2 ; *-- proportional odds model; proc logistic nosimple data=arthrit; model improve = _sex_ _treat_ ; output out=results p=predict l=lower u=upper xbeta=logit / alpha=.33; data results; set results; if _level_=0 then better = (improve > 0); else better = (improve > 1); *-- Change direction of probabilities & logit; predict = 1 - predict; lower = 1 - lower; upper = 1 - upper; logit = -logit; proc print data=results(obs=50); id id treat sex; var improve _level_ predict lower upper logit; *-- Collapse results data set to one obs per group, find observed probability for each group (mean of improve); proc summary data=results nway; class treat sex _level_; var better; id predict lower upper logit; output out=res2 mean=obs; data res2; set res2; treatl = treat||put(_level_,1.0); proc sort data=res2; by treatl; proc print; id treatl sex _level_; *-- Annotate data set for confidence intervals and obs. prob.; data bars; set res2; by treatl; retain xsys ysys '2'; length text $21; if treat='Placebo' then color='BLACK'; else color='RED'; xc = sex; line=33; text='circle'; y = obs; function='SYMBOL '; output; y = upper; function='MOVE '; output; text='-'; function='LABEL '; output; y = lower; function='DRAW '; output; text='-'; function='LABEL '; output; if last.treatl then do; y = predict; position='C'; size=1.3; text = ' '|| treat; function='LABEL'; output; position='E'; size=1; if _level_ = 0 then text=' None vs. Some/Marked'; else text=' None/Some vs. Marked'; output; end; proc gplot data=res2; plot predict * sex = treatl / vaxis=axis1 haxis=axis2 frame nolegend anno=bars; axis1 label=(h=1.4 a=90 'Prob. Improvement (67% CI)') value=(h=1.2) order=(0 to 1 by .2); axis2 label=(h=1.5) value=(h=1.4) offset=(6,22); symbol1 v=+ h=1.5 i=join l=3 c=black; symbol2 v=+ h=1.5 i=join l=3 c=black; symbol3 v=$ h=1.5 i=join l=1 c=red; symbol4 v=$ h=1.5 i=join l=1 c=red; %gskip; *-- Proportional Odds Model: Effects of treat, sex and age; proc logistic data=arthrit nosimple; model improve = _sex_ _treat_ age ; output out=results p=predict l=lower u=upper xbeta=logit stdxbeta=selogit / alpha=.33; data results; set results; treatl = treat||put(_level_,1.0); if _level_=0 then better = (improve > 0); else better = (improve > 1); *-- Change direction of probabilities & logit; predict = 1 - predict; lower = 1 - lower; upper = 1 - upper; logit = -logit; proc print data=results(obs=50); id id treat sex; var improve _level_ predict lower upper logit; format predict lower upper logit selogit 7.3; proc sort data=results; by sex treatl age; data bars; set results; by sex treatl; length text $8; xsys='2'; ysys='2'; if treat='Placebo' then color='BLACK'; else color='RED'; x = age; line=33; *-- plot confidence limits ; y = upper; function='MOVE '; output; text='-'; function='LABEL '; output; y = lower; function='DRAW '; output; text='-'; function='LABEL '; output; if last.treatl then do; y = predict; x = age+1; position='C'; size=1.4; text = treat; function='LABEL'; output; position='F'; if _level_ = 0 then text='> None'; else text='> Some';; output; end; if first.sex then do; ysys ='1'; y=90; xsys ='1'; x=10; size=1.5; text = sex; function='LABEL'; output; end; /*: :*/ goptions hby=0; proc gplot; plot predict * age = treatl / vaxis=axis1 haxis=axis2 nolegend anno=bars ; by sex; axis1 label=(h=1.4 a=90 'Prob. Improvement (67% CI)') value=(h=1.2) order=(0 to 1 by .2); axis2 label=(h=1.4) value=(h=1.2) order=(20 to 80 by 10) offset=(2,5); symbol1 v=+ h=1.4 i=join l=3 c=black; symbol2 v=+ h=1.4 i=join l=3 c=black; symbol3 v=$ h=1.4 i=join l=1 c=red; symbol4 v=$ h=1.4 i=join l=1 c=red; %gfinish;