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 = .;


/* Nominal logistic regression via LOGISTIC */
PROC LOGISTIC DATA=one;
MODEL Housing(REF=FIRST) = Time1 Time2 Time3 Section8 Sect8T1 Sect8T2 Sect8T3 / LINK=GLOGIT;
RUN;

/* Nominal random intercept logistic regression via GLIMMIX */
PROC GLIMMIX DATA=one METHOD=QUAD(QPOINTS=11) NOCLPRINT;
CLASS Id Housing;
MODEL Housing(REF=FIRST) = Time1 Time2 Time3 Section8 Sect8T1 Sect8T2 Sect8T3 / SOLUTION DIST=MULTINOMIAL LINK=GLOGIT;
RANDOM INTERCEPT / SUBJECT=Id GROUP=Housing;
RUN;

/* Nominal random intercept logistic regression via NLMIXED */
PROC NLMIXED DATA=one QPOINTS=11;
PARMS g1_0=-.49 g1_Time1=1.63 g1_Time2=2.37 g1_Time3=1.8 g1_Section8=.43 g1_Sect8T1=-.46 g1_Sect8T2=-2.12 g1_Sect8T3=-1.1 var1=1
      g2_0=-1.7 g2_Time1=1.90 g2_Time2=2.97 g2_Time3=2.9 g2_Section8=.54 g2_Sect8T1=1.13 g2_Sect8T2=-.04  g2_Sect8T3=-.07 var2=1;
z1 = g1_0 + g1_Time1*Time1 + g1_Time2*Time2 + g1_Time3*Time3 + g1_Section8*Section8
          + g1_Sect8T1*Sect8T1 + g1_Sect8T2*Sect8T2 + g1_Sect8T3*Sect8T3 + u1;
z2 = g2_0 + g2_Time1*Time1 + g2_Time2*Time2 + g2_Time3*Time3 + g2_Section8*Section8
          + g2_Sect8T1*Sect8T1 + g2_Sect8T2*Sect8T2 + g2_Sect8T3*Sect8T3 + u2;
IF (Housing=0) THEN
  p = 1 / (1 + EXP(z1) + EXP(z2));
ELSE IF (Housing=1) THEN
  p = EXP(z1) / (1 + EXP(z1) + EXP(z2));
ELSE IF (Housing=2) THEN
  p = EXP(z2) / (1 + EXP(z1) + EXP(z2));
LogLike = LOG(p);
MODEL Housing ~ GENERAL(LogLike);
RANDOM u1 u2 ~ NORMAL([0,0],[var1,0,var2]) SUBJECT=Id;
ESTIMATE 'ICC1' var1 /((((ATAN(1)*4)**2)/3) + var1);
ESTIMATE 'ICC2' var2 /((((ATAN(1)*4)**2)/3) + var2);
RUN;


/* Nominal random intercept logistic regression via NLMIXED - HELMERT CONTRASTS */
PROC NLMIXED DATA=one QPOINTS=11;
PARMS h1_0=-.49 h1_Time1=1.63 h1_Time2=2.37 h1_Time3=1.8 h1_Section8=.43 h1_Sect8T1=-.46 h1_Sect8T2=-2.12 h1_Sect8T3=-1.1 varh1=1
      h2_0=-1.7 h2_Time1=1.90 h2_Time2=2.97 h2_Time3=2.9 h2_Section8=.54 h2_Sect8T1=1.13 h2_Sect8T2=-.04  h2_Sect8T3=-.07 varh2=1;
z1 = h1_0 + h1_Time1*Time1 + h1_Time2*Time2 + h1_Time3*Time3 + h1_Section8*Section8
          + h1_Sect8T1*Sect8T1 + h1_Sect8T2*Sect8T2 + h1_Sect8T3*Sect8T3 + u1;
z2 = h2_0 + h2_Time1*Time1 + h2_Time2*Time2 + h2_Time3*Time3 + h2_Section8*Section8
          + h2_Sect8T1*Sect8T1 + h2_Sect8T2*Sect8T2 + h2_Sect8T3*Sect8T3 + u2;
IF (Housing=0) THEN
  p = EXP(-2/3*z1) /
      (EXP(-2/3*z1) + EXP(1/3*z1 - 1/2*z2) + EXP(1/3*z1 + 1/2*z2));
ELSE IF (Housing=1) THEN
  p = EXP(1/3*z1 - 1/2*z2) /
      (EXP(-2/3*z1) + EXP(1/3*z1 - 1/2*z2) + EXP(1/3*z1 + 1/2*z2));
ELSE IF (Housing=2) THEN
  p = EXP(1/3*z1 + 1/2*z2) /
      (EXP(-2/3*z1) + EXP(1/3*z1 - 1/2*z2) + EXP(1/3*z1 + 1/2*z2));
LogLike = LOG(p);
MODEL Housing ~ GENERAL(LogLike);
RANDOM u1 u2 ~ NORMAL([0,0],[varh1,0,varh2]) SUBJECT=Id;
ESTIMATE 'ICCh1' varh1 /((((ATAN(1)*4)**2)/3) + varh1);
ESTIMATE 'ICCh2' varh2 /((((ATAN(1)*4)**2)/3) + varh2);
RUN;


/* Nominal random intercept logistic regression via NLMIXED - correlated RE */
PROC NLMIXED DATA=one QPOINTS=11;
PARMS g1_0=-.49 g1_Time1=1.63 g1_Time2=2.37 g1_Time3=1.8 g1_Section8=.43 g1_Sect8T1=-.46 g1_Sect8T2=-2.12 g1_Sect8T3=-1.1 var1=1
      g2_0=-1.7 g2_Time1=1.90 g2_Time2=2.97 g2_Time3=2.9 g2_Section8=.54 g2_Sect8T1=1.13 g2_Sect8T2=-.04  g2_Sect8T3=-.07 var2=1
      cov12=0;
z1 = g1_0 + g1_Time1*Time1 + g1_Time2*Time2 + g1_Time3*Time3 + g1_Section8*Section8
          + g1_Sect8T1*Sect8T1 + g1_Sect8T2*Sect8T2 + g1_Sect8T3*Sect8T3 + u1;
z2 = g2_0 + g2_Time1*Time1 + g2_Time2*Time2 + g2_Time3*Time3 + g2_Section8*Section8
          + g2_Sect8T1*Sect8T1 + g2_Sect8T2*Sect8T2 + g2_Sect8T3*Sect8T3 + u2;
IF (Housing=0) THEN
  p = 1 / (1 + EXP(z1) + EXP(z2));
ELSE IF (Housing=1) THEN
  p = EXP(z1) / (1 + EXP(z1) + EXP(z2));
ELSE IF (Housing=2) THEN
  p = EXP(z2) / (1 + EXP(z1) + EXP(z2));
LogLike = LOG(p);
MODEL Housing ~ GENERAL(LogLike);
RANDOM u1 u2 ~ NORMAL([0,0],[var1,cov12,var2]) SUBJECT=Id;
ESTIMATE 'RE Corr' cov12 / SQRT(var1*var2);
RUN;


/* Nominal random intercept logistic regression via NLMIXED - HELMERT CONTRASTS - correlated RE */
PROC NLMIXED DATA=one QPOINTS=11;
PARMS h1_0=-.49 h1_Time1=1.63 h1_Time2=2.37 h1_Time3=1.8 h1_Section8=.43 h1_Sect8T1=-.46 h1_Sect8T2=-2.12 h1_Sect8T3=-1.1 varh1=1
      h2_0=-1.7 h2_Time1=1.90 h2_Time2=2.97 h2_Time3=2.9 h2_Section8=.54 h2_Sect8T1=1.13 h2_Sect8T2=-.04  h2_Sect8T3=-.07 varh2=1
      covh12=0;
z1 = h1_0 + h1_Time1*Time1 + h1_Time2*Time2 + h1_Time3*Time3 + h1_Section8*Section8
          + h1_Sect8T1*Sect8T1 + h1_Sect8T2*Sect8T2 + h1_Sect8T3*Sect8T3 + u1;
z2 = h2_0 + h2_Time1*Time1 + h2_Time2*Time2 + h2_Time3*Time3 + h2_Section8*Section8
          + h2_Sect8T1*Sect8T1 + h2_Sect8T2*Sect8T2 + h2_Sect8T3*Sect8T3 + u2;
IF (Housing=0) THEN
  p = EXP(-2/3*z1) /
      (EXP(-2/3*z1) + EXP(1/3*z1 - 1/2*z2) + EXP(1/3*z1 + 1/2*z2));
ELSE IF (Housing=1) THEN
  p = EXP(1/3*z1 - 1/2*z2) /
      (EXP(-2/3*z1) + EXP(1/3*z1 - 1/2*z2) + EXP(1/3*z1 + 1/2*z2));
ELSE IF (Housing=2) THEN
  p = EXP(1/3*z1 + 1/2*z2) /
      (EXP(-2/3*z1) + EXP(1/3*z1 - 1/2*z2) + EXP(1/3*z1 + 1/2*z2));
LogLike = LOG(p);
MODEL Housing ~ GENERAL(LogLike);
RANDOM u1 u2 ~ NORMAL([0,0],[varh1,covh12,varh2]) SUBJECT=Id;
ESTIMATE 'RE Corr' covh12 / SQRT(varh1*varh2);
RUN;