TITLE1 'San Diego Homeless Data - Estimated Marginal Probabilities';
PROC IML;

/* Marginalization using Quadrature */;

/* Results from NLMIXED analysis */;

x0 = { 0 0 0 0 0 0 0,
       1 0 0 0 0 0 0,
       0 1 0 0 0 0 0,
       0 0 1 0 0 0 0};
x1 = { 0 0 0 1 0 0 0,
       1 0 0 1 1 0 0,
       0 1 0 1 0 1 0,
       0 0 1 1 0 0 1};

varu  = { 2.128};

beta  = { 1.736,
          2.315,
          2.499,
          0.497,
          1.408,
          1.173,
          0.638};

thresh = {.220, 2.964};


/* 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};
wq = {   0.00000431065265,   0.00075807095698,   0.01911158107317,
         0.13548370704150,   0.34464234526294,   0.34464234526294,
         0.13548370704150,   0.01911158107317,   0.00075807095698,
         0.00000431065265};

mprb0a = J(4,1,0);
mprb0b = J(4,1,0);
mprb1a = J(4,1,0);
mprb1b = J(4,1,0);

DO q = 1 to nq;

   z0a= thresh[1] - (x0*beta + SQRT(varu)*bq[q]);
   z0b= thresh[2] - (x0*beta + SQRT(varu)*bq[q]);
   z1a= thresh[1] - (x1*beta + SQRT(varu)*bq[q]);
   z1b= thresh[2] - (x1*beta + SQRT(varu)*bq[q]);

   mprb0a = mprb0a + (  1.0   / (1.0 + EXP(- z0a)))*wq[q];
   mprb0b = mprb0b + (  1.0   / (1.0 + EXP(- z0b)))*wq[q];
   mprb1a = mprb1a + (  1.0   / (1.0 + EXP(- z1a)))*wq[q];
   mprb1b = mprb1b + (  1.0   / (1.0 + EXP(- z1b)))*wq[q];

END;

print 'Random-intercepts model: Quadrature method - 10 points';

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',
       (1-mprb0b) [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',
       (1-mprb1b) [FORMAT=8.4];



/* Approximate Marginalization Method - no quadrature */;

pi    = ATAN(1)*4;
nt    = 4;
ivec  = j(nt,1,1);
evec  = (15/16)**2 * (pi**2)/3 * ivec;
zmat  = 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;

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

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

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

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',
       (1-mprb0b) [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',
       (1-mprb1b) [FORMAT=8.4];


/* non-proportional odds model */


varu = { 2.124};

gam1 = { 2.297,
         3.346,
         2.823,
         0.592,
         0.567,
        -0.959,
        -0.369};
gam2 = { 1.080,
         1.646,
         2.147,
         0.324,
         2.022,
         2.016,
         1.071};

thresh = {0.322, 2.377};

mprb0a = J(4,1,0);
mprb0b = J(4,1,0);
mprb1a = J(4,1,0);
mprb1b = J(4,1,0);

DO q = 1 to nq;

   z0a= thresh[1] - (x0*gam1 + SQRT(varu)*bq[q]);
   z0b= thresh[2] - (x0*gam2 + SQRT(varu)*bq[q]);
   z1a= thresh[1] - (x1*gam1 + SQRT(varu)*bq[q]);
   z1b= thresh[2] - (x1*gam2 + SQRT(varu)*bq[q]);

   mprb0a = mprb0a + (  1.0   / (1.0 + EXP(- z0a)))*wq[q];
   mprb0b = mprb0b + (  1.0   / (1.0 + EXP(- z0b)))*wq[q];
   mprb1a = mprb1a + (  1.0   / (1.0 + EXP(- z1a)))*wq[q];
   mprb1b = mprb1b + (  1.0   / (1.0 + EXP(- z1b)))*wq[q];

END;

print 'NON-PROPORTIONAL ODDS MODEL';
print 'Random-intercepts model: Quadrature method - 10 points';

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',
       (1-mprb0b) [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',
       (1-mprb1b) [FORMAT=8.4];