TITLE1 'NIMH Schizophrenia Study - Estimated Marginal Probabilities';
PROC IML;

/* Results from NLMIXED analysis: random-intercepts model */;

x0 = { 0 0.00000 0,
       0 1.00000 0,
       0 1.73205 0,
       0 2.44949 0};
x1 = { 1 0.00000 0.00000,
       1 1.00000 1.00000,
       1 1.73205 1.73205,
       1 2.44949 2.44949};

varu  = { 3.774};

beta  = {-0.058,
         -0.766,
         -1.206};

thresh = {-5.859, -2.826, -.709};

/* Approximate Marginalization Method */;

pi    = ATAN(1)*4;
nt    = 4;
ivec  = j(nt,1,1);
zvec  = j(nt,1,1);
evec  = (15/16)**2 * (pi**2)/3 * ivec;

/* nt by nt matrix with evec on the diagonal and zeros elsewhere */;
emat = diag(evec);

/* variance variance-covariance matrix of underlying latent variable */;
vary = zvec * varu * T(zvec) + emat;

sdy  = sqrt(vecdiag(vary) / vecdiag(emat));

z0a= (thresh[1] - x0*beta) / sdy;
z0b= (thresh[2] - x0*beta) / sdy;
z0c= (thresh[3] - x0*beta) / sdy;
z1a= (thresh[1] - x1*beta) / sdy;
z1b= (thresh[2] - x1*beta) / sdy;
z1c= (thresh[3] - x1*beta) / sdy;

mprb0a =   1.0   / (1.0 + EXP(- z0a));
mprb0b =   1.0   / (1.0 + EXP(- z0b));
mprb0c =   1.0   / (1.0 + EXP(- z0c));
mprb1a =   1.0   / (1.0 + EXP(- z1a));
mprb1b =   1.0   / (1.0 + EXP(- z1b));
mprb1c =   1.0   / (1.0 + EXP(- z1c));

print 'Approximate Marginalization Method', (1 / sdy) [FORMAT=8.4];

print 'marginal probability for group 0 - 1st category',
       mprb0a [FORMAT=8.4];
print 'marginal probability for group 0 - 2nd category',
       (mprb0b-mprb0a) [FORMAT=8.4];
print 'marginal probability for group 0 - 3rd category',
       (mprb0c-mprb0b) [FORMAT=8.4];
print 'marginal probability for group 0 - 4th category',
       (1-mprb0c) [FORMAT=8.4];
print 'marginal probability for group 1 - 1st category',
       mprb1a [FORMAT=8.4];
print 'marginal probability for group 1 - 2nd category',
       (mprb1b-mprb1a) [FORMAT=8.4];
print 'marginal probability for group 1 - 3rd category',
       (mprb1c-mprb1b) [FORMAT=8.4];
print 'marginal probability for group 1 - 4th category',
       (1-mprb1c) [FORMAT=8.4];


/* Results from NLMIXED analysis: random-intercept & trend model */;

x0 = { 0 0.00000 0,
       0 1.00000 0,
       0 1.73205 0,
       0 2.44949 0};
x1 = { 1 0.00000 0.00000,
       1 1.00000 1.00000,
       1 1.73205 1.73205,
       1 2.44949 2.44949};

varu  = { 7.001 -1.510,
         -1.510  2.010};

beta  = { 0.057,
         -0.883,
         -1.695};

thresh = {-7.321, -3.419, -.813};

zmat  = {1 0.00000,
         1 1.00000,
         1 1.73205,
         1 2.44949};

/* variance variance-covariance matrix of underlying latent variable */;
vary = zmat * varu * T(zmat) + emat;

sdy  = sqrt(vecdiag(vary) / vecdiag(emat));

z0a= (thresh[1] - x0*beta) / sdy;
z0b= (thresh[2] - x0*beta) / sdy;
z0c= (thresh[3] - x0*beta) / sdy;
z1a= (thresh[1] - x1*beta) / sdy;
z1b= (thresh[2] - x1*beta) / sdy;
z1c= (thresh[3] - x1*beta) / sdy;

mprb0a =   1.0   / (1.0 + EXP(0 - z0a));
mprb0b =   1.0   / (1.0 + EXP(0 - z0b));
mprb0c =   1.0   / (1.0 + EXP(0 - z0c));
mprb1a =   1.0   / (1.0 + EXP(0 - z1a));
mprb1b =   1.0   / (1.0 + EXP(0 - z1b));
mprb1c =   1.0   / (1.0 + EXP(0 - z1c));

print 'Approximate Marginalization Method', (1 / sdy) [FORMAT=8.4];

print 'marginal probability for group 0 - 1st category',
       mprb0a [FORMAT=8.4];
print 'marginal probability for group 0 - 2nd category',
       (mprb0b-mprb0a) [FORMAT=8.4];
print 'marginal probability for group 0 - 3rd category',
       (mprb0c-mprb0b) [FORMAT=8.4];
print 'marginal probability for group 0 - 4th category',
       (1-mprb0c) [FORMAT=8.4];
print 'marginal probability for group 1 - 1st category',
       mprb1a [FORMAT=8.4];
print 'marginal probability for group 1 - 2nd category',
       (mprb1b-mprb1a) [FORMAT=8.4];
print 'marginal probability for group 1 - 3rd category',
       (mprb1c-mprb1b) [FORMAT=8.4];
print 'marginal probability for group 1 - 4th category',
       (1-mprb1c) [FORMAT=8.4];