TITLE1 'Reisby Data - Estimated Marginal Probabilities';
PROC IML;

/* covariate matrix for low drug level observations */
x0 = {1 0};
/* covariate matrix for high drug level observations */
x1 = {1 1};

/* GEE analysis - exchangeable working correlation structure */;
beta   = {-0.397, 1.092};
z0     = x0*beta;
z1     = x1*beta;
mprb0 =  1.0   / (1.0 + EXP(0 - z0));
mprb1 =  1.0   / (1.0 + EXP(0 - z1));

PRINT 'gee analysis - exchangeable working correlation';
PRINT 'marginal prob for low drug'  mprb0  [FORMAT=8.4];
PRINT 'marginal prob for high drug' mprb1  [FORMAT=8.4];



/* random-intercepts model */
beta   = {-0.661,  1.842};
sd     = 2.004;
pi     = 3.141592654;

betastd = beta / SQRT(((16/15)**2) * (3/pi**2) * sd**2 + 1);

z0     = x0*betastd;
z1     = x1*betastd;
mprb0 =  1.0   / (1.0 + EXP(0 - z0));
mprb1 =  1.0   / (1.0 + EXP(0 - z1));

PRINT 'random-intercepts model - marginalized results';
PRINT 'marginal prob for low drug'  mprb0   [FORMAT=8.4];
PRINT 'marginal prob for high drug' mprb1   [FORMAT=8.4];


/* get the estimated marginal probabilities using quadrature */;
/* number of quadrature points, quadrature nodes & weights   */;

nq = 10;
bq = {  -4.85946282833231,  -3.58182348355193,  -2.48432584163895,
        -1.46598909439116,  -0.48493570751550,   0.48493570751550,
         1.46598909439116,   2.48432584163895,   3.58182348355193,
         4.85946282833231};
aq = {   0.00000431065265,   0.00075807095698,   0.01911158107317,
         0.13548370704150,   0.34464234526294,   0.34464234526294,
         0.13548370704150,   0.01911158107317,   0.00075807095698,
         0.00000431065265};

mprb0  = 0;
mprb1  = 0;

DO q=1 TO nq;
   z0 = x0*beta + sd*bq[q];
   z1 = x1*beta + sd*bq[q];
   mprb0 = mprb0 + (1.0   / (1.0 + EXP(0 - z0)))*aq[q];
   mprb1 = mprb1 + (1.0   / (1.0 + EXP(0 - z1)))*aq[q];
END;

PRINT 'random-int model: Quad method - 10 points';
PRINT 'marginal prob for low drug'  mprb0 [FORMAT=8.4];
PRINT 'marginal prob for high drug' mprb1 [FORMAT=8.4];


--------
 Output
--------

                         Reisby Data - Estimated Marginal Probabilities                        


                        gee analysis - exchangeable working correlation


                                                            MPRB0

                              marginal prob for low drug   0.4020


                                                             MPRB1

                              marginal prob for high drug   0.6671


                         random-intercepts model - marginalized results


                                                            MPRB0

                              marginal prob for low drug   0.3947


                                                             MPRB1

                              marginal prob for high drug   0.6822


                           random-int model: Quad method - 10 points


                                                            MPRB0

                              marginal prob for low drug   0.4014


                                                             MPRB1

                              marginal prob for high drug   0.6726