TITLE1 'nimh schizophrenia study - estimated marginal probabilities';
PROC IML;
/* results from nlmixed analysis: random-trend model */;

/* covariate matrics for placebo and drug groups */
x0 = {1 0 0.00000 0,
      1 0 1.00000 0,
      1 0 1.73205 0,
      1 0 2.44949 0};
x1 = {1 1 0.00000 0.00000,
      1 1 1.00000 1.00000,
      1 1 1.73205 1.73205,
      1 1 2.44949 2.44949};

/* nlmixed estimates of covariate effects and random effect variance covariance matrix */
beta  = { 5.814,  0.274, -1.390, -1.556};
varu  = { 6.504 -1.802, 
         -1.802  2.682};

/* marginalization of person-specific estimates */;
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 * varu * T(zmat) + emat;

/* marginalization factor */
sdy  = SQRT(VECDIAG(vary) / VECDIAG(emat));
z0 = (x0*beta) / sdy;
z1 = (x1*beta) / sdy;

grp0 =   1.0   / (1.0 + EXP(0 - z0));
grp1 =   1.0   / (1.0 + EXP(0 - z1));

print 'random intercept and trend model';
print 'marginalization of person-specific estimates';
print 'marginal response probability for group 0', grp0 [FORMAT=8.4];
print 'marginal response probability for group 1', grp1 [FORMAT=8.4];