Thanks so much, Rick.
Extra background detail. There is a problem with marginal observation-level Predicted Probabilities from logistic regression when random effects are involved. Straight up linear predictor expressed as probability is systematically different from the 'correct' value. See this for more detail (https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-015-0046-6).
The correct way to calculated marginal probabilities according to the reference and the formula I pasted, is to integrate Predicted Probability for each row over the distribution of all random intervals in the model. I wish SAS had this correction built in into all PROCs that do multi-level regression. I requested that feature. In the absence of it, my starting point is a marginal NOBLUP Predicted Probability value for each row in the model from GLIMMIX, for example. From there, I back-transform it to linear predictor so that the LP could be integrated over the distribution of random effects or Random Intervals (the u term). For each row in the model (patients), the Predicted Probability needs to be corrected using the distribution of random intercepts (hospitals). Normal integral QUAD() gives me AUC for the transformation without any density weights. What you shared seems to do the trick. I tested a couple of values and it seems to correctly increase the input for values below 0.5 and decreases them for values above .5. I modified the code slightly and now just need to figure out how to apply it to an entire column of Predicted Probability values. I think I can simplify this a little bit by outputing Linear Predictors from the model rather than Predicted Probabilities -- skip the `iLink` option.
proc IML;
PredProb = .73; /* Marginal Predicted Probability from RE GLMM */
Eta = log(PP / (1-PP)); /* Eta = linear predictor */
start Integrand(u) global(eta);
L = logistic(u + eta);
fu = pdf("Normal", u, 0, 0.2); /* specify f(u) here. eg, PDF for N(0,0.20) */
return( L * fu );
finish;
call quad(iPP, "Integrand", {.M,.P});
print PP Eta iPP;
quit;
Thanks again. Happy to discuss further. Hope you can put in a good word with SAS development teams to add this correction to the SAS/STAT core PROCS. Would love to hear if you disagree with the diagnosis for marginal PPs.
... View more