data happy;
input Income Happiness$ count;
datalines;
3 Very 272
3 Pretty 294
3 Not 49
2 Very 454
2 Pretty 835
2 Not 131
1 Very 185
1 Pretty 527
1 Not 208
;

proc genmod data=happy descending;
 model Happiness = Income / dist=multinomial link=clogit lrci type3;
 freq count;
*output out=bla PREDPROBS=I;
run;

proc genmod data=happy descending;
 class Income;
 model Happiness = Income / dist=multinomial link=clogit lrci type3;
 freq count;
*output out=bla PREDPROBS=I;
contrast "LR test" Income 1 -2 1 /E;
*estimate "odds for happy vs. not happy when income below average" Income -1 0 1 /exp;
run;

proc logistic data=happy descending;
model Happiness = Income / aggregate scale=none;
freq count;
output out=probs predprobs=I C; *requests fitted individual and category probabilities;
run;

data probs; set probs;
logit1=log(CP_very/(1-CP_very));
logit2=log(CP_pretty/(1-CP_pretty));
run;

/* setting graphing parameters */
goptions htext=2;
axis1 label=("Income") order = (1 to 3 by 1) minor = none;
axis2 label=(angle=90 "Predicted Cumulative Logit") order = (-2 to 3 by 1) minor=none;
legend1 label=none value=('<= Very Happy' '<= Pretty Happy')
       position=(bottom center inside) mode=share cborder=black;
symbol interpol=join value=dot width=2;

proc gplot data = probs;
plot (logit1 logit2)*Income /overlay haxis=axis1 vaxis=axis2 legend=legend1;
run;
quit;

/* fitted cumulative probabilities */
axis3 label=(angle=90 "Predicted Cumulative Prob.") order = (0 to 1 by .2) minor=none;
proc gplot data = probs;
plot (CP_Very CP_Pretty)*Income /overlay haxis=axis1 vaxis=axis3 legend=legend1;
run;
quit;

/* fitted category probabilities */
axis4 label=(angle=90 "Predicted Probability") order = (0 to 1 by .2) minor=none;
legend2 label=none value=('Very Happy' 'Pretty Happy' 'Not Happy')
       position=(top center inside) mode=share cborder=black;

proc gplot data = probs;
plot (IP_Very IP_Pretty IP_Not)*Income /overlay haxis=axis1 vaxis=axis4 legend=legend2;
run;
quit;


/* Fit baseline category model */
proc logistic data=happy descending;
model Happiness(ref='Not') = Income / link=glogit aggregate scale=none;
freq count;
run;
/* Fit Adjacent Category Model with proportional odds */
data happy1;
input Income y1 y2 y3;
datalines;
3 272 294 49
2 454 835 131
1 185 527 208
;
proc nlmixed data=happy1;
eta1 = alpha1+alpha2+2*beta*Income;
eta2 = alpha2+beta*Income;
p3 = 1/(1+exp(eta1)+exp(eta2));
p1 = exp(eta1)*p3;
p2 = exp(eta2)*p3;
ll = y1*log(p1) + y2*log(p2) + y3*log(p3);
model y1 ~ general(ll);
*predict p1 out=pred1;
*predict p2 out=pred2;
*predict p3 out=pred3;
/* Predicted category probabilities */
estimate "P(Y = very happy|Income=3)" exp(alpha1+alpha2+2*beta*3)/
(1+exp(alpha1+alpha2+2*beta*3)+exp(alpha2+beta*3));
estimate "P(Y = very happy|Income=2)" exp(alpha1+alpha2+2*beta*2)/
(1+exp(alpha1+alpha2+2*beta*2)+exp(alpha2+beta*2));
estimate "P(Y = very happy|Income=1)" exp(alpha1+alpha2+2*beta)/
(1+exp(alpha1+alpha2+2*beta)+exp(alpha2+beta));

estimate "P(Y = pretty happy|Income=3)" exp(alpha2+beta*3)/
(1+exp(alpha1+alpha2+2*beta*3)+exp(alpha2+beta*3));
estimate "P(Y = pretty happy|Income=2)" exp(alpha2+beta*2)/
(1+exp(alpha1+alpha2+2*beta*2)+exp(alpha2+beta*2));
estimate "P(Y = pretty happy|Income=1)" exp(alpha2+beta)/
(1+exp(alpha1+alpha2+2*beta)+exp(alpha2+beta));

estimate "P(Y = pretty happy|Income=3)" 1/
(1+exp(alpha1+alpha2+2*beta*3)+exp(alpha2+beta*3));
estimate "P(Y = not happy|Income=2)" 1/
(1+exp(alpha1+alpha2+2*beta*2)+exp(alpha2+beta*2));
estimate "P(Y = not happy|Income=1)" 1/
(1+exp(alpha1+alpha2+2*beta)+exp(alpha2+beta));
run;


/* Fit Continuation-Ratio Model with proportional odds */
proc nlmixed data=happy1;
eta1 = alpha1+beta*Income;
eta2 = alpha2+beta*Income;
p1 = exp(eta1)/(1+exp(eta1));
p2 = exp(eta2)/((1+exp(eta1))*(1+exp(eta2)));
p3 = 1-p1-p2;
ll = y1*log(p1) + y2*log(p2) + y3*log(p3);
model y1 ~ general(ll);
/* Predicted category probabilities */
predict p1 out=pred1;
 run;
proc print data=pred1 label;
label Pred = 'P(Y = very happy | Income)';
var Income Pred;
run;
/* Fit Continuation-Ratio Model with proportional odds */
data happy2;
input stratum Income successes failures;
n=successes+failures;
datalines;
1 3 272 343
1 2 454 966
1 1 185 735
2 3 294 49
2 2 835 131
2 1 527 208
;
proc genmod data=happy2;
class stratum;
model successes/n = stratum Income /noint dist=binomial link=logit;
run;

/* Fit Cumulative Probit Model */
proc genmod data=happy descending;
 model Happiness = Income / dist=multinomial link=cprobit lrci type3;
 freq count;
 output out=predc pred=C; 
run;
proc print; run;
/* Plot fitted cumulative prob for logit and probit model */
data probs; set probs;
logit1=log(CP_very/(1-CP_very));
logit2=log(CP_pretty/(1-CP_pretty));
run;
/* setting graphing parameters */
goptions htext=2;
axis1 label=("Income") order = (1 to 3 by 1) minor = none;
axis2 label=(angle=90 "Predicted Cumulative Logit") order = (-2 to 3 by 1) minor=none;
legend1 label=none value=('<= Very Happy' '<= Pretty Happy')
       position=(bottom center inside) mode=share cborder=black;
symbol interpol=join value=dot width=2;
proc gplot data = probs;
plot (logit1 logit2)*Income /overlay haxis=axis1 vaxis=axis2 legend=legend1;
run;
quit;
/* fitted cumulative probabilities */
axis3 label=(angle=90 "Predicted Cumulative Prob.") order = (0 to 1 by .2) minor=none;
proc gplot data = probs;
plot (CP_Very CP_Pretty)*Income /overlay haxis=axis1 vaxis=axis3 legend=legend1;
run;
quit;