OPTIONS NOCENTER ;
TITLE 'Mixed-effects Analyses of San Diego Homeless Data';
DATA one;
INFILE 'C:\data\sandiego\sdhouse.DAT';
INPUT Id Housing Int Section8 Time1 Time2 Time3 Sect8T1 Sect8T2 Sect8T3
      NonSect8 LinTime S8LinT ;

      IF Housing = 999 then Housing = .;

/* fixed-effects proportional odds model */
PROC LOGISTIC DESCENDING DATA=one;
MODEL Housing = Time1 Time2 Time3 Section8 Sect8T1 Sect8T2 Sect8T3;
RUN;


/* Random-intercept proportional odds model via GLIMMIX */
PROC GLIMMIX DATA=one METHOD=QUAD(QPOINTS=21) NOCLPRINT;
CLASS id;
MODEL Housing(DESC) = Time1 Time2 Time3 Section8 Sect8T1 Sect8T2 Sect8T3 /
      SOLUTION DIST=MULTINOMIAL LINK=CUMLOGIT;
RANDOM INTERCEPT / SUBJECT=id;
RUN;


/* Random-intercept proportional odds model via NLMIXED */
PROC NLMIXED DATA=one QPOINTS=21;
PARMS bTime1=1.4 bTIme2=1.8 bTIme3=2.0 bSection8=.4 bSect8T1=.9 bSect8T2=.9 bSect8T3=.4
g1=.2 g2=2.3 varu=1;
z = bTime1*Time1 + bTime2*Time2 + bTime3*Time3 + bSection8*Section8
       + bSect8T1*Sect8T1 + bSect8T2*Sect8T2 + bSect8T3*Sect8T3+ u;
IF (Housing=0) THEN
  p = 1 / (1 + EXP(-(g1-z)));
ELSE IF (Housing=1) THEN
  p = (1/(1 + EXP(-(g2-z)))) - (1/(1 + EXP(-(g1-z))));
ELSE IF (Housing=2) THEN
  p = 1 - (1 / (1 + EXP(-(g2-z))));
LogLike = LOG(p);
MODEL Housing ~ GENERAL(LogLike);
RANDOM u ~ NORMAL(0,varu) SUBJECT=Id;
ESTIMATE 'icc' varu/((((ATAN(1)*4)**2)/3) + varu);
RUN;


/* Random-intercept NON-proportional odds model via NLMIXED */
PROC NLMIXED DATA=one QPOINTS=21;
PARMS g1_Time1=1.4 g1_TIme2=1.8 g1_TIme3=2.0 g1_Section8=.4 g1_Sect8T1=.9 g1_Sect8T2=.9 g1_Sect8T3=.4
      g2_Time1=1.4 g2_TIme2=1.8 g2_TIme3=2.0 g2_Section8=.4 g2_Sect8T1=.9 g2_Sect8T2=.9 g2_Sect8T3=.4
      g1=.2 g2=2.3 varu=1;
z1 = g1_Time1*Time1 + g1_Time2*Time2 + g1_Time3*Time3 + g1_Section8*Section8
       + g1_Sect8T1*Sect8T1 + g1_Sect8T2*Sect8T2 + g1_Sect8T3*Sect8T3+ u;
z2 = g2_Time1*Time1 + g2_Time2*Time2 + g2_Time3*Time3 + g2_Section8*Section8
       + g2_Sect8T1*Sect8T1 + g2_Sect8T2*Sect8T2 + g2_Sect8T3*Sect8T3+ u;
IF (HOUSING=0) THEN
  p = 1 / (1 + EXP(-(g1-z1)));
ELSE IF (HOUSING=1) THEN
  p = (1/(1 + EXP(-(g2-z2)))) - (1/(1 + EXP(-(g1-z1))));
ELSE IF (HOUSING=2) THEN
  p = 1 - (1 / (1 + EXP(-(g2-z2))));
LogLike = LOG(p);
MODEL Housing ~ GENERAL(LogLike);
RANDOM u ~ NORMAL(0,varu) SUBJECT=Id;
ESTIMATE 'icc' varu/((((ATAN(1)*4)**2)/3) + varu);
RUN;