Getting missing values for the random effects in PRIC GLIMMIX

[DocMartin]

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

[Rick_SAS]

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;

[jiltao]

What if you add METHOD=LAPLACE or METHOD=QUAD option to the PROC GLIMMIX statement?

Jill

[SteveDenham]

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

[DocMartin]

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!!!!

[Rick_SAS]

What does PROC LOGISTIC report? What are the parameter estimates? Any convergence issues or notes in the log?

[DocMartin]

It runs fine in logistic regression. I included hospital as a covariate in the model.

[Rick_SAS]

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?

[DocMartin]

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 

[Rick_SAS]

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."

[jiltao · Accepted solution]

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

[DocMartin]

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!

[jiltao]

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