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

/* Results from MIXOR analysis */;

x0 = { 0 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};
x1 = { 1 0 0 0 0 0 0,
       1 1 0 0 1 0 0,
       1 0 1 0 0 1 0,
       1 0 0 1 0 0 1};

int   = {-0.220};
sd    = { 1.459};

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

thresh = {2.744};


/* 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= 0      - (int + x0*beta + sd*bq[q]);
   z0b= thresh - (int + x0*beta + sd*bq[q]);
   z1a= 0      - (int + x1*beta + sd*bq[q]);
   z1b= thresh - (int + x1*beta + sd*bq[q]);

   mprb0a = mprb0a + (  1.0   / (1.0 + EXP(0 - z0a)))*wq[q];
   mprb0b = mprb0b + (  1.0   / (1.0 + EXP(0 - z0b)))*wq[q];
   mprb1a = mprb1a + (  1.0   / (1.0 + EXP(0 - z1a)))*wq[q];
   mprb1b = mprb1b + (  1.0   / (1.0 + EXP(0 - 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    = 3.141592654;
nt    = 4;
ivec  = 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 = ivec * sd * T(sd) * T(ivec) + emat;

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

z0a= (0      - (int + x0*beta)) / sdy;
z0b= (thresh - (int + x0*beta)) / sdy;
z1a= (0      - (int + x1*beta)) / sdy;
z1b= (thresh - (int + x1*beta)) / sdy;

mprb0a =   1.0   / (1.0 + EXP(0 - z0a));
mprb0b =   1.0   / (1.0 + EXP(0 - z0b));
mprb1a =   1.0   / (1.0 + EXP(0 - z1a));
mprb1b =   1.0   / (1.0 + EXP(0 - 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 */


int   = {-0.322};
sd    = { 1.457};

beta1 = { 0.592,
          2.297,
          3.345,
          2.821,
          0.566,
         -0.958,
         -0.366};
beta2 = { 0.323,
          1.079,
          1.645,
          2.145,
          2.023,
          2.016,
          1.073};

thresh = {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= 0      - (int + x0*beta1 + sd*bq[q]);
   z0b= thresh - (int + x0*beta2 + sd*bq[q]);
   z1a= 0      - (int + x1*beta1 + sd*bq[q]);
   z1b= thresh - (int + x1*beta2 + sd*bq[q]);

   mprb0a = mprb0a + (  1.0   / (1.0 + EXP(0 - z0a)))*wq[q];
   mprb0b = mprb0b + (  1.0   / (1.0 + EXP(0 - z0b)))*wq[q];
   mprb1a = mprb1a + (  1.0   / (1.0 + EXP(0 - z1a)))*wq[q];
   mprb1b = mprb1b + (  1.0   / (1.0 + EXP(0 - 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];