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

/* Results from MIXOR 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};

int   = { 5.858};
sd    = { 1.944};

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

thresh = {3.033, 5.152};

/* Approximate Marginalization Method */;

pi    = 3.141592654;
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 * sd * T(sd) * T(zvec) + emat;

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

z0a= (0         - (int + x0*beta)) / sdy;
z0b= (thresh[1] - (int + x0*beta)) / sdy;
z0c= (thresh[2] - (int + x0*beta)) / sdy;
z1a= (0         - (int + x1*beta)) / sdy;
z1b= (thresh[1] - (int + x1*beta)) / sdy;
z1c= (thresh[2] - (int + 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];

/* Results from MIXOR 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};

int   = { 7.309};
chol  = { 2.669 0,
         -0.588 1.308};

beta  = { 0.111,
         -0.875,
         -1.724};

thresh = {3.912, 6.528};

/* Approximate Marginalization Method */;

pi    = 3.141592654;
nt    = 4;
ivec  = j(nt,1,1);
zmat  = {1 0.00000,
         1 1.00000,
         1 1.73205,
         1 2.44949};
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 = zmat * chol * T(chol) * T(zmat) + emat;

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

z0a= (0         - (int + x0*beta)) / sdy;
z0b= (thresh[1] - (int + x0*beta)) / sdy;
z0c= (thresh[2] - (int + x0*beta)) / sdy;
z1a= (0         - (int + x1*beta)) / sdy;
z1b= (thresh[1] - (int + x1*beta)) / sdy;
z1c= (thresh[2] - (int + 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];