Glimmix and predictive margins

[AJJ]

Hello all

My current work involves profiling of hospital performance regarding a specific binary quality indicator (Patients are nested within hospitals, i.e. 2-level data)

Performance will be adjusted for various risk factors, e.g. age, diagnosis etc.

My intention is then to use a random hospital intercept binary logit model (a fixed effect model  - estimating intercepts for each hospital using dummies - could perhaps  be an alternative ) , 

As a risk adjusted measure of hospital performance I would use so called "predictive margins"; i.e.:

1. Estimate the random intercept model including the empirical bayesian estimates (BLUPs) of the hospital specific intercepts

2. Assume all patients belong to hospital 1 - values on risk factors remain unchanged

3. calculate the predicted probability using the estimated fixed parameters and BLUP for hospital 1 for each patient and average over patients

Do 2 and 3 for every hospital.

Questions:

1. This is a build in command in stata (the margins command), but I guess i could obtain the same results using "store" in glimmix and proc plm?

2. the margins command in Stata uses the delta method to obtain standard errors of the predictive margins. Is there some way I could get the standard errors in SAS without having to bootstrap (time consuming since > 100 hospitals)

Thanks in advance.

Anders

[ballardw]

A general SAS observation when you say you want to "do the same thing for each hospita"l: BY processing. Sort the data by the variable(s) of interest and then add a BY statement to the procedure. You will get output for each level of the variable, or combination of variables, specified.

There are few exceptions but generally the BY group processsing is the preferred approach.

[SteveDenham]

Does this example in the PROC HPLMIXED documentation help in any way (Example 9.1 Computing BLUPs for a Large Number of Subjects)?  You would need to "know" the covariance parameters, but approximate values might be obtained from a preliminary analysis using PROC MIXED/GLIMMIX.  I realize the HPLMIXED is for a normal variate, and yours is binomial, but it sure looks like the right starting point to me.

Steve Denham

[lvm · Accepted solution]

You can get the EBLUPs easy enough with a solution option on a random statement:

random hospital / s;

You could store these with the approrpriate ODS OUTPUT statement (look up the name for the random effect predictions). You can get the inverse link for each subject using ESTIMATE statements. For instance, assume you have one predictor variable X, the EBLUP for hospital one and three are :

class hospital .....<other terms> ...;

model y = X / dist=binomial link=logit;

random hospital / s;

estimate 'x 1 hosp 1' int 1 X 1 | hospital 1 / ilink;

estimate 'x 1 hosp 2' int 1 X 1 | hospital 0 0 1 / ilink;

My code gets the EBLUP for X=1 (if X is continuous) or the first level of X if it is a factor.

The ILINK option uses the delta method for the SEs. I have never heard of predictive margins, and have no idea why one assumes that all observations belong, in turn, to each subject. I assume there are good reasons. But based on what you wrote, you would need to write code in a data step (or multiple data steps) to do what you want. It would take some time.

[AJJ]

Thank you all for your advice and help! I guess there is no easy way to do this - but I did not expect that either.

Anders