I’m having trouble getting a mixed-effects model to run properly in PROC GLIMMIX. I get the message that the G matrix is not positive definite, the covariance parameter estimate is 0, and the Solution for Random Effects are all missing.
The objective of my study is to see whether or not mortality is reduced when comparing a placebo phase to a treatment phase. To do that, I’m trying to construct a mixed-effects model of mortality. Fixed effects include age, weight, phase (placebo/treatment), and whether or not the patient had a Do Not Resuscitate order written (DNR). The random effects are the six hospitals that participated in both phases. There were 2,250 patients in the placebo phase and 1,890 patients in the treatment phase.
data main; input mortality phase $ age weight dnr hospital $;
cards; * Only first 3 lines of 4,140 patients;
0 p 65 153 0 hosp1
1 p 78 225 0 hosp1
0 t 74 185 0 hosp2
;
run;
proc glimmix data = main1;
class hospital phase(ref = first) ;
model mortality(ref = '0') = age weight phase dnr / dist = binary link = logit solution;;
random int / solution sub = hospital;
title "Mixed Effects Model of Mortality";
title2 “Random Intercepts Are Hospitals”;
lsmeans phase / diff cl;
run;
Mortality has a incidence of 12%, and there are plenty of patients within each category of the predictors. When I run an identical model on readmission, which has a much lower incidence (5%), GLIMMIX runs with no problems.
Do anyone have an idea what is happening and a potential solution?
Thanks!
Andrew
It looks to me that your data might not support the model you specified.
One possible scenario is that the correlation is negative, in which case, you might try the following similar model (but not identical) --
proc glimmix data = main1;
class hospital phase(ref = first) ;
model mortality(ref = '0') = age weight phase dnr / dist = binary link = logit solution;;
random _residual_/ type=cs sub = hospital;
title "Mixed Effects Model of Mortality";
title2 “Random Intercepts Are Hospitals”;
lsmeans phase / diff cl;
run;
This model does not give you random effect solutions though.
Thanks,
Jill
Try running the corresponding logistic model without the random effects. Does it converge? Any warnings or errors? Post the log if there are any problems with the logistic model.
title "Fixed Effects Model of Mortality";
proc logistic data = main;
class hospital phase(ref = first) / param=glm;
model mortality(ref = '0') = age weight phase dnr;
lsmeans phase / diff cl;
run;
What if you add METHOD=LAPLACE or METHOD=QUAD option to the PROC GLIMMIX statement?
Jill
I like @jiltao 's suggestion. Another one would be to try a different technique as an optimizer. i believe that the default method is being used, which is probably QUANEW. Check the log - if it is not converging the error message could say something with regards to QUANEW, and the optimization information table where the technique should be listed. Failure to converge might be the source of the missing estimates. Another would be something pathological with the data, so that the Hessian runs into issues. In either case, try:
nloptions tech=congra;
The conjugate gradient method requires only that the first derivatives exist (same as QUANEW) but is appropriate when the objective function and gradient are faster to evaluate than the Hessian. Watch out for local optima with this method - you may have to search a parameter grid to check that it doesn't get stuck in a local stationary point.
SteveDenham
I tried all of the suggestions and still came up with estimates of 0 for the random effect. Here is the warning message in the LOG:
NOTE: The GLIMMIX procedure is modeling the probability that mortality ='1'.
NOTE: Convergence criterion (GCONV=1E-8) satisfied.
NOTE: At least one element of the gradient is greater than 1e-3.
NOTE: Estimated G matrix is not positive definite.
NOTE: PROCEDURE GLIMMIX used (Total process time):
real time 0.34 seconds
cpu time 0.31 seconds
When I looked at the results, two of the tables stated:
Optimization Information
Optimization Technique Dual Quasi-Newton
Parameters in Optimization 6
Lower Boundaries 1
Upper Boundaries 0
Fixed Effects Not Profiled
Starting From GLM estimates
Model Information
Data Set WORK.MAIN1
Response Variable mortality
Response Distribution Binary
Link Function Logit
Variance Function Default
Variance Matrix Blocked By icuid
Estimation Technique Maximum Likelihood
Likelihood Approximation Laplace
Degrees of Freedom Method Containment
Taking out hospital as an effect does result in the model running o.k., but I need to include that variable since observations are correlated within hospitals.
Curiously, when I take out DNR as a variable the model runs properly. But DNR is such an important confounder that I'll get crushed for excluding it from the model.
Help!!!!
What does PROC LOGISTIC report? What are the parameter estimates? Any convergence issues or notes in the log?
OK, but that is not the model that I asked you to run. I asked you to run the fixed-effects model that excludes hospital. I also asked for the parameter estimates for the model. Can you cut/paste the output?
Here are the results of a logistic regression model that did not include hospital:
Model Fit Statistics
Criterion Intercept Only Intercept and
Covariates
AIC 2891.742 1697.565
SC 2898.072 1729.215
-2 Log L 2889.742 1687.565
Testing Global Null Hypothesis: BETA=0
Test Chi-Square DF Pr > ChiSq
Likelihood Ratio 1202.1767 4 <.0001
Score 1523.0580 4 <.0001
Wald 721.5999 4 <.0001
Type 3 Analysis of Effects
Effect DF Wald Chi-Square Pr > ChiSq
Age 1 259.5185 <.0001
weight 1 0.1552 0.6936
dnr 1 508.2104 <.0001
phase 1 0.021 0.9054
Analysis of Maximum Likelihood Estimates
Parameter DF Estimate Standard Error Wald Chi-Square Pr > ChiSq
Intercept 1 -5.3453 0.1803 878.6498 <.0001
Age 1 0.0354 0.00220 259.5185 <.0001
Weight 1 0.00246 0.00625 0.1552 0.6936
dnr 1 2.9583 0.1312 508.2104 <.0001
phase t 1 -0.00812 0.0658 0.0141 0.9054
Association of Predicted Probabilities and Observed
Responses
Percent Concordant 91.6 Somers' D 0.833
Percent Discordant 8.4 Gamma 0.833
Percent Tied 0.0 Tau-a 0.164
Pairs 1695560 c 0.916
Sometimes looking at the fixed-effect model first can guide the choice of effects in the model, which can help to simplify the model, which can affect whether the model converges.
The results from PROC LOGISTIC suggest that you can remove Weight from the model. Interestingly, Phase is not significant. So the objective of your study ("to see whether or not mortality is reduced when comparing a placebo phase to a treatment phase") seems to be "no."
It looks to me that your data might not support the model you specified.
One possible scenario is that the correlation is negative, in which case, you might try the following similar model (but not identical) --
proc glimmix data = main1;
class hospital phase(ref = first) ;
model mortality(ref = '0') = age weight phase dnr / dist = binary link = logit solution;;
random _residual_/ type=cs sub = hospital;
title "Mixed Effects Model of Mortality";
title2 “Random Intercepts Are Hospitals”;
lsmeans phase / diff cl;
run;
This model does not give you random effect solutions though.
Thanks,
Jill
The model converged! While I didn't get the hospital-specific random effects, they were not vital. It was adjusting for the fact that patients within hospitals were correlated.
One question: What does taking out "int" from the RANDOM statement and replacing that term with "_resid_" imply in terms of the model?
Thanks!
The random intercept model --
random int / sub = hospital;
uses a G-side random effect to model the correlations in the data.
The correlations in the data can also be modeled through the R-side random effects. For example,
random _residual_/ type=cs sub = hospital;
The documentation below has more information on the G-side and R-side random effects --
SAS Help Center: G-Side and R-Side Random Effects and Covariance Structures
If your data is normal, the above two models produce identical marginal models. And this is explained and illustrated in the PROC MIXED documentation below -
SAS Help Center: Mixed Models Theory
For generalized linear mixed models like yours, the two models are not going to be identical, but should be close.
Hope this helps,
Jill
Join us for SAS Innovate 2025, our biggest and most exciting global event of the year, in Orlando, FL, from May 6-9. Sign up by March 14 for just $795.
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.