Question about Nested Random Intercept in Proc Glimmix

[anujmehta]

I have a question about interpreting the estimates for a nested random intercept using proc glimmix. We are interested in understanding racial variation in mortality within a hospital and comparing across hospitals. We want estimates for each race within each hospital or how much each race deviates from the hospital estimate. the key variables are the binary outcome died, racial identity race3 (4 level categorical variable), and hospital id dshospid, and multiple clinical variables for risk-adjustment (var1-var10). the code we used is 

proc glimmix data=dataset method=laplace;

class dshospid race3; 

model died (event='1') = var1 var2 var3 . . . var10 / dist=binary link=logit solution;

random intercept / subject=dshospid type=un s;

random intercept / subject=race3(dshospid) type=un s;

run;

The random intercept output looks like:

My question is about how we should interpret the nested random intercept output. In my view, the nested random intercepts (race3(DSHOSPID)) values show how much each racial identity deviates from the hospital estimate. The key question is should these nested random intercept values sum to zero? 

[SteveDenham]

If race3 was considered as a fixed effect, I would expect the values to sum to zero. In this case, I think it will depend on the distribution of responses to each race, which will differ by hospital (as I see that not all races are observed at every hospital). So perhaps the variances and covariances for each hospital are not homogeneous, and could be fit with a GROUP= option. You would likely need two runs (one homogeneous, ignoring race3 and one heterogeneous, accommodating race3 differences). This approach has the advantage of being applicable to the fixed levels of race3, rather than assuming that the levels were selected randomly from a population. I would try fitting the following:

/* Homogeneous variance-covariance across race3 */
proc glimmix data=dataset method=laplace;
class dshospid race3;
model died (event='1') = var1 var2 var3 . . . var10 / dist=binary link=logit solution;
random intercept / subject=dshospid type=un s;
run;

/* Heterogeneous variance-covariance across race3 */
proc glimmix data=dataset method=laplace;
class dshospid race3;
model died (event='1') = var1 var2 var3 . . . var10 / dist=binary link=logit solution;
random intercept / subject=dshospid type=un s group=race3;
covtest homogeneity;
run;

SteveDenham