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

/* Results from GLIMMIX analysis - uncorrelated random effects */;

u0 = { 1 0 0 0 0 0 0 0,
       1 1 0 0 0 0 0 0,
       1 0 1 0 0 0 0 0,
       1 0 0 1 0 0 0 0};
u1 = { 1 0 0 0 1 0 0 0,
       1 1 0 0 1 1 0 0,
       1 0 1 0 1 0 1 0,
       1 0 0 1 1 0 0 1};

varu = { 1.099, 3.249};

gam1 = { -0.635,
          1.914,
          2.686,
          1.982,
          0.518,
         -0.639,
         -2.499,
         -1.234};

gam2 = { -2.600,
          2.332,
          3.677,
          3.735,
          0.675,
          1.886,
          0.556,
          0.256};


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


/* initialize to zero */
grp0a = J(4,1,0);
grp0b = J(4,1,0);
grp0c = J(4,1,0);
grp1a = J(4,1,0);
grp1b = J(4,1,0);
grp1c = J(4,1,0);

DO q1 = 1 to nq;
DO q2 = 1 to nq;

z1_0 = u0*gam1 + SQRT(varu[1])*bq[q1];
z2_0 = u0*gam2 + SQRT(varu[2])*bq[q2];
z1_1 = u1*gam1 + SQRT(varu[1])*bq[q1];
z2_1 = u1*gam2 + SQRT(varu[2])*bq[q2];

denom0 = 1 + exp(z1_0) + exp(z2_0);
denom1 = 1 + exp(z1_1) + exp(z2_1);

grp0a=grp0a + (1/denom0)*wq[q1]*wq[q2];
grp0b=grp0b + (exp(z1_0)/denom0)*wq[q1]*wq[q2];
grp0c=grp0c + (exp(z2_0)/denom0)*wq[q1]*wq[q2];

grp1a=grp1a + (1/denom1)*wq[q1]*wq[q2];
grp1b=grp1b + (exp(z1_1)/denom1)*wq[q1]*wq[q2];
grp1c=grp1c + (exp(z2_1)/denom1)*wq[q1]*wq[q2];

END;
END;

print 'Uncorrelated random effects model fit - Quadrature method - 10 points';

print 'marginal probability for group 0 - category 0' grp0a [FORMAT=8.4];
print 'marginal probability for group 0 - category 1' grp0b [FORMAT=8.4];
print 'marginal probability for group 0 - category 2' grp0c [FORMAT=8.4];

print 'marginal probability for group 1 - category 0' grp1a [FORMAT=8.4];
print 'marginal probability for group 1 - category 1' grp1b [FORMAT=8.4];
print 'marginal probability for group 1 - category 2' grp1c [FORMAT=8.4];




/* Results from NLMIXED analysis - correlated random effects */;


varu = { 2.702 2.884,
         2.884 5.821};

/* get the Cholesky in lower triangular form */
chol = T(ROOT(varu));

/* select the two rows of the Cholesky for the two category contrasts */
chol1 = chol[1,];
chol2 = chol[2,];

gam1 = { -0.622,
          2.374,
          3.345,
          2.589,
          0.652,
         -0.331,
         -2.475,
         -1.160};

gam2 = { -2.599,
          2.888,
          4.394,
          4.310,
          0.831,
          1.969,
          0.287,
          0.202};


/* initialize to zero */
grp0a = J(4,1,0);
grp0b = J(4,1,0);
grp0c = J(4,1,0);
grp1a = J(4,1,0);
grp1b = J(4,1,0);
grp1c = J(4,1,0);

DO q1 = 1 to nq;
DO q2 = 1 to nq;

bqq = bq[q1] // bq[q2];

z1_0 = u0*gam1 + chol1*bqq;
z2_0 = u0*gam2 + chol2*bqq;
z1_1 = u1*gam1 + chol1*bqq;
z2_1 = u1*gam2 + chol2*bqq;

denom0 = 1 + exp(z1_0) + exp(z2_0);
denom1 = 1 + exp(z1_1) + exp(z2_1);

grp0a=grp0a + (1/denom0)*wq[q1]*wq[q2];
grp0b=grp0b + (exp(z1_0)/denom0)*wq[q1]*wq[q2];
grp0c=grp0c + (exp(z2_0)/denom0)*wq[q1]*wq[q2];

grp1a=grp1a + (1/denom1)*wq[q1]*wq[q2];
grp1b=grp1b + (exp(z1_1)/denom1)*wq[q1]*wq[q2];
grp1c=grp1c + (exp(z2_1)/denom1)*wq[q1]*wq[q2];


END;
END;

print 'Correlated random effects model fit - Quadrature method - 10 points';

print 'marginal probability for group 0 - category 0' grp0a [FORMAT=8.4];
print 'marginal probability for group 0 - category 1' grp0b [FORMAT=8.4];
print 'marginal probability for group 0 - category 2' grp0c [FORMAT=8.4];

print 'marginal probability for group 1 - category 0' grp1a [FORMAT=8.4];
print 'marginal probability for group 1 - category 1' grp1b [FORMAT=8.4];
print 'marginal probability for group 1 - category 2' grp1c [FORMAT=8.4];