options linesize=80 pagesize=54;
DATA  one; INFILE 'c:\mix\schizx1.dat';
INPUT id imps79 imps79b imps79o int tx week sweek txswk ;

/* get rid of observations with missing values */
IF imps79 > -9;

PROC FORMAT;
     VALUE tx      0 = 'placebo' 1 = 'drug';

/* fixed effects ordinal logistic regression */
PROC LOGISTIC DESCENDING;
MODEL imps79o = tx sweek tx*sweek;
RUN;


/* random intercept logistic regression via GLIMMIX */
PROC GLIMMIX DATA=one METHOD=QUAD(QPOINTS=21) NOCLPRINT;
CLASS id;
MODEL imps79o(DESC) = tx sweek tx*sweek / SOLUTION DIST=MULTINOMIAL LINK=CUMLOGIT;
RANDOM INTERCEPT / SUBJECT=id;
RUN;


/* random trend logistic regression via GLIMMIX */
PROC GLIMMIX DATA=one METHOD=QUAD(QPOINTS=11) NOCLPRINT;
CLASS id;
MODEL imps79o(DESC) = tx sweek tx*sweek / SOLUTION DIST=MULTINOMIAL LINK=CUMLOGIT;
RANDOM INTERCEPT sweek / SUBJECT=id TYPE=UN GCORR SOLUTION;
ODS LISTING EXCLUDE SOLUTIONR; ODS OUTPUT SOLUTIONR=ebest2;
RUN;

/* print out the estimated random effects */
PROC PRINT DATA=ebest2;
RUN;

/* ordinal logistic random intercept model via NLMIXED */
/* GLIMMIX descending parameterization */
PROC NLMIXED DATA=one QPOINTS=21;
PARMS b_tx=0 b_sweek=-.54 b_txswk=-.75 varu=1 int4=.42 int3=1.8 int2=3.8 ;
z = b_tx*tx + b_sweek*sweek + b_txswk*tx*sweek + u;
IF (imps79o=4) THEN
  p = 1 / (1 + EXP(-(int4+z)));
ELSE IF (imps79o=3) THEN
  p = (1/(1 + EXP(-(int3+z)))) - (1/(1 + EXP(-(int4+z))));
ELSE IF (imps79o=2) THEN
  p = (1/(1 + EXP(-(int2+z)))) - (1/(1 + EXP(-(int3+z))));
ELSE IF (imps79o=1) THEN
  p = 1 - (1 / (1 + EXP(-(int2+z))));
loglik = LOG(p);
MODEL imps79o ~ GENERAL(loglik);
RANDOM u ~ NORMAL(0,varu) SUBJECT=id;
ESTIMATE 'icc' varu/((((ATAN(1)*4)**2)/3) + varu);
RUN;


/* ordinal logistic random intercept model via NLMIXED */
PROC NLMIXED DATA=one QPOINTS=21;
PARMS b_tx=0 b_sweek=-.54 b_txswk=-.75 varu=1 g1=-3.8 g2=-1.8 g3=-.42;
z = b_tx*tx + b_sweek*sweek + b_txswk*tx*sweek + u;
IF (imps79o=1) THEN
  p = 1 / (1 + EXP(-(g1-z)));
ELSE IF (imps79o=2) THEN
  p = (1/(1 + EXP(-(g2-z)))) - (1/(1 + EXP(-(g1-z))));
ELSE IF (imps79o=3) THEN
  p = (1/(1 + EXP(-(g3-z)))) - (1/(1 + EXP(-(g2-z))));
ELSE IF (imps79o=4) THEN
  p = 1 - (1 / (1 + EXP(-(g3-z))));
loglik = LOG(p);
MODEL imps79o ~ GENERAL(loglik);
RANDOM u ~ NORMAL(0,varu) SUBJECT=id;
ESTIMATE 'icc' varu/((((ATAN(1)*4)**2)/3) + varu);
RUN;


/* ordinal logistic random trend model via NLMIXED */
PROC NLMIXED data=one QPOINTS=21;
PARMS b_tx=-.06 b_sweek=-.77 b_txswk=-1.21 v0=3.77 c01=0 v1=1 g1=-5.9 g2=-2.8 g3=-.71;
z = b_tx*tx + b_sweek*sweek + b_txswk*tx*sweek + u0 + u1*sweek;
IF (imps79o=1) THEN
  p = 1 / (1 + EXP(-(g1-z)));
ELSE IF (imps79o=2) THEN
  p = (1/(1 + EXP(-(g2-z)))) - (1/(1 + EXP(-(g1-z))));
ELSE IF (imps79o=3) THEN
  p = (1/(1 + EXP(-(g3-z)))) - (1/(1 + EXP(-(g2-z))));
ELSE IF (imps79o=4) THEN
  p = 1 - (1 / (1 + EXP(-(g3-z))));
loglik = LOG(p);
MODEL imps79o ~ GENERAL(loglik);
RANDOM u0 u1 ~ NORMAL([0,0], [v0,c01,v1]) SUBJECT=id  out=ebest2b;
ESTIMATE 're corr' c01/SQRT(v0*v1);
RUN;

/* print out the estimated random effects */
PROC PRINT DATA=ebest2b;
RUN;