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

/* BINARY OUTCOME */
/* 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.387};
sd    = { 2.116};

beta  = {-0.025,
         -1.500,
         -1.015};

/* 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));

z0= (int + x0*beta) / sdy;
z1= (int + x1*beta) / sdy;

mprb0 =   1.0   / (1.0 + EXP(0 - z0));
mprb1 =   1.0   / (1.0 + EXP(0 - z1));

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

print 'marginal response probability for group 0',
       mprb0 [FORMAT=8.4];
print 'marginal response probability for group 1',
       mprb1 [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   = { 6.025};
chol  = { 2.726 0,
         -0.829 1.561};

beta  = { 0.281,
         -1.477,
         -1.587};

/* 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));

z0= (int + x0*beta) / sdy;
z1= (int + x1*beta) / sdy;

mprb0 =   1.0   / (1.0 + EXP(0 - z0));
mprb1 =   1.0   / (1.0 + EXP(0 - z1));

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

print 'marginal response probability for group 0',
       mprb0 [FORMAT=8.4];
print 'marginal response probability for group 1',
       mprb1 [FORMAT=8.4];