Logistic regression with PROC GLIMMIX - specification of random statement w/ multiple random effects

[bretthouston]

Hi, 

I'm running a logistic regression model, hoping to evaluate the impact of a patient's surgical team on their likelihood of receiving a specific medication. I would like to incorporate the patient's anesthesiologist & surgeon as separate random effects. I've structured the model as below, but I just wanted to confirm that this is a reasonable approach, as most of the glimmix documents I've read specify the random statement with "random int / subject = X". In my case, I haven't structured it that way because I have two random effects and I would rather not study the composite variability of (anesthesiologist*surgeon) if possible (I'd like to evaluate the variability in their likelihood of using a medication separately if possible). Is this a reasonable way to structure the code given those objectives?

proc glimmix data=TXAHipA method=Laplace;
class anesth surgeon sex urgency year;
effect age_spline = spline(age / details naturalcubic basis=tpf(noint)
knotmethod=percentiles(5));
effect Hb_spline = spline(Hb / details naturalcubic basis=tpf(noint)
knotmethod=percentiles(5));
model TXAstatus = age_spline Hb_spline sex urgency year / solution dist=binomial link=logit;
random anesth surgeon;
COVTEST / wald;
output out=glimmixout pred( blup ilink) = PredProb
pred(noblup ilink) = PredProb_PA;
run;

Another way I've considered is as follows, but it creates issues with identifying the separate impacts of the two individuals:

proc glimmix data=TXAHipA method=Laplace;
class anesth surgeon sex urgency year;
effect age_spline = spline(age / details naturalcubic basis=tpf(noint)
knotmethod=percentiles(5));
effect Hb_spline = spline(Hb / details naturalcubic basis=tpf(noint)
knotmethod=percentiles(5));
model TXAstatus = age_spline Hb_spline sex urgency year / solution dist=binomial link=logit;
random intercept / subject=anesth*surgeon;
COVTEST / wald;
output out=glimmixout pred( blup ilink) = PredProb
pred(noblup ilink) = PredProb_PA;
run;

Any advice would be much appreciated. I've had a difficult time figuring out the differences between these two syntaxes and want to make my code reflects my intentions. 

Thanks,

Brett

[SteveDenham]

Nothing wrong with the first batch of code.  The main reason you see a lot of RANDOM intercept/subject=X examples is that by subject processing is much faster.  You might consider having two RANDOM statements, each modeling a random intercept, such as:

RANDOM intercept/subject=anesth;

RANDOM intercept/subject=surgeon;

Be sure the data set is sorted by anesth and surgeon before trying this, however.  I have some other ideas in case this doesn't work out, but again, your first batch of code ought to give what you want, at the price of computation time.

Regarding the Wald option in COVTEST, check out this Wikipedia link (down at the part on Alternatives to the Wald Test);https://en.wikipedia.org/wiki/Wald_test .  I would go with the classical option, which does a likelihood ratio test.  At least that doesn't assume that you have a good estimate of the variability of the variance components.

SteveDenham

[bretthouston]

Hi Steve, 

Thank-you for this feedback - this is incredibly helpful. I've run the code the second way you suggested (with two random statements), and the results are essentially identical. 

I appreciate your thoughts about the likelihood ratio test. I had explicitly used the covtest / wald statement because otherwise SAS doesn't provide p-values (maybe appropriately so!) for the covariance parameters in this setting. Is there an alternative way you request the likelihood ratio test in proc glimmix (ie, the classical option)?

Thanks again,

Brett

[SteveDenham]

COVTEST classical;

That should give a chi^2 test for each individual component, and

COVTEST zerog;;

Should give a joint test that all of the variance components are zero.  If either of these leads to problems (like can't find a feasible starting set of parameters) try adding the RESTART option, so that the null model is fit from scratch, rather than from the fully fitted model.

SteveDenham

[bretthouston]

Thanks Steve. 

I can get the covtest zerog working, although the covtest classical doesn't produce a covariance test result (the table remains limited to estimate and standard estimate). Could this be because the method is specified as laplace? (I get non-convergence if I leave the method unspecified).

proc glimmix data=TXAHipA method=laplace;
class anesth surgeon sex admitCategory year;
effect age_spline = spline(age / details naturalcubic basis=tpf(noint)
knotmethod=percentiles(5));
effect Hb_spline = spline(Hb / details naturalcubic basis=tpf(noint)
knotmethod=percentiles(5));
model TXAstatus = age_spline Hb_spline sex admitCategory CS1 CS2 CS3 year / solution dist=binomial link=logit;
random anesth surgeon;
covtest / classical;
output out=glimmixout pred( blup ilink) = PredProb
pred(noblup ilink) = PredProb_PA;
run;

Thanks,

Brett

[SteveDenham]

Remove the / before classical and see if that helps.

SteveDenham