proc bglimm

[carstenenevoldsen]

GLIMMIX gives this option which is very useful for studying discrimination and calibration by means of proc logistic:

output out=gmxout pred=xbeta pred(ilink)=predprob pred(ilink noblup)=fix_predprob pearson=pearson_residual;

BGLIMM - I cannot find a similar option.

[SteveDenham]

I am supposed to get an upgrade soon that will give me access to BGLIMM, so I have been studying the documentation.  It looks like BGLIMM does everything using ESTIMATE statements, or through post-processing.  This example

 https://documentation.sas.com/?docsetId=statug&docsetTarget=statug_bglimm_examples03.htm&docsetVersi...  presents how to post-process using the parameter estimates obtained to get predicted values for a Poisson regression.  With some tweaking (and by some, I mean, well, a lot), you could get the predicted probabilities (blup and noblup) by fitting in the parameter estimates for the fixed and random effects (blup) and for just the fixed effects (noblup).

I hope this helps, and I really, really hope that my upgrade shows up soon, so I could actually run some simulated data to get at what you are asking

SteveDenham

[carstenenevoldsen]

Hi Steve,

Thanks for the link, which I should have been aware of.

You can get access to bglimm if you use SAS University Edition (straight forward to install).

The solution suggested in the link is not transparent to me.

How should it be implemented when we have two (or more) random effect statements? I have an example with two cross-classified random effects.

Regards

Carsten

[SteveDenham]

Well, the random statements give rise to the estimates for the gamma vector, so for any instance you would have to do something like:

/* Calculate the equivalent of a pred(blup) */
numblup = exp (X'beta + Z'gamma);
prob_ind = numblup/(1+numblup);

/* or calculate the equivalent of a red(noblup) */
numnoblup = exp(X'beta);
prob_mean = numnoblup/(1+numnoblup);

It appears that you get the beta's and gamma's from a dataset created with the outpost= option in the PROC BLIMM statement.  Without having one of those at hand, I am sort of stuck. I think that you could then calculate residuals by merging the calculated values back against the original dataset (or better SQL, where you can also calculate the residuals in a single call).

SteveDenham