I have very recently started learning about Generalised Linear Mixed Models and was using R to explore what difference it makes to treat group membership as either fixed or random effect. In particular, I am looking at the example dataset discussed here:



As outlined in this tutorial, the effect of Doctor ID is appreciable and I was expecting the mixed model with a random intercept to give better results. However, comparing AIC values for the two methods suggest that this model is worse:

> require(lme4) ; hdp = read.csv("http://www.ats.ucla.edu/stat/data/hdp.csv")
> hdp$DID = factor(hdp$DID) ; hdp$Married = factor(hdp$Married)
> GLM = glm(remission~Age+Married+IL6+DID,data=hdp,family=binomial);summary(GLM)

glm(formula = remission ~ Age + Married + IL6 + DID, family = binomial, 
data = hdp)

Deviance Residuals: 
Min       1Q   Median       3Q      Max  
-2.5265  -0.6278  -0.2272   0.5492   2.7329  

              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.560e+01  1.219e+03  -0.013    0.990    
Age         -5.869e-02  5.272e-03 -11.133  < 2e-16 ***
Married1     2.688e-01  6.646e-02   4.044 5.26e-05 ***
IL6         -5.550e-02  1.153e-02  -4.815 1.47e-06 ***
DID2         1.805e+01  1.219e+03   0.015    0.988    
DID3         1.932e+01  1.219e+03   0.016    0.987   


DID405       1.566e+01  1.219e+03   0.013    0.990    
DID405       1.566e+01  1.219e+03   0.013    0.990    
DID406      -2.885e-01  3.929e+03   0.000    1.000    
DID407       2.012e+01  1.219e+03   0.017    0.987    
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 10353  on 8524  degrees of freedom
Residual deviance:  6436  on 8115  degrees of freedom
AIC: 7256

Number of Fisher Scoring iterations: 17

> GLMM = glmer(remission~Age+Married+IL6+(1|DID),data=hdp,family=binomial) ; m

Generalized linear mixed model fit by the Laplace approximation 
Formula: remission ~ Age + Married + IL6 + (1 | DID) 
Data: hdp 
AIC  BIC logLik deviance
7743 7778  -3867     7733
Random effects:
Groups Name        Variance Std.Dev.
DID    (Intercept) 3.8401   1.9596  
Number of obs: 8525, groups: DID, 407

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.461438   0.272709   5.359 8.37e-08 ***
Age         -0.055969   0.005038 -11.109  < 2e-16 ***
Married1     0.260065   0.063736   4.080 4.50e-05 ***
IL6         -0.053288   0.011058  -4.819 1.44e-06 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
         (Intr) Age    Marrd1
Age      -0.898              
Married1  0.070 -0.224       
IL6      -0.162  0.012 -0.033

> extractAIC(GLM) ; extractAIC(GLMM)

[1]  410.000 7255.962
[1]    5.000 7743.188

Thus, my questions are:

(1) Is it appropriate to compare the AIC values provided by the two functions? If so, why does the fixed effect model do better?

(2) What is the best way to identify if fixed or random effects are more important (ie to quantify that the variability due to the doctor is more important than patient characteristics?


Fixed effects models and random effects models ask different questions of the data. Specifying a set of group-level dummy variables essentially controls for all group-level unobserved heterogeneity in the average response, leaving your estimates to reflect only variability within units. Random effects models start with the assumption that there is a meta-population of (whatever effect), and that your sample reflects many draws from that population. So rather than anchoring your results around heterogeneous intercepts, your data will be used to elucidate the parameters of that (usually normal) distribution from which your data were supposedly drawn.

It is often said that fixed effects models are good for conducting inference on the data that you have, and that random effects models are good for trying to conduct inference on some larger population from which your data is a random sample.

When I learned about fixed effects models, they were motivated using error components and panel data. Take multiple observations of a given unit, and a random treatment in time $t$.

$$y_{it} = \alpha_i + \beta T_{it} + \epsilon_{it}$$

You can break your error term out into that component of your error term that varies in time, and one that doesn't:

$$y_{it} = \alpha_i + \beta T_{it} + e_i + u_{it}$$

Now subtract the groupwise mean from both sides:

$$y_{it} - \bar y_i = \alpha_i - \bar \alpha_i + \beta \left(T_{it}- \bar T_i\right) + e_i - \bar e_i+ u_{it}- \bar u_it$$

Things that aren't subscripted by $t$ come out of the equation by basic subtraction -- which is to say that the average over time is the same as it is at any time if it never changes. This includes your non-time-varying component of your error term. Thus your estimates are unconfounded by time-invariant heterogeneity.

This doesn't quite work for a random effects model -- your non-$t$-indexed variables won't be sopped up by that transformation (the "within" transformation). As such, you can draw inference on the effects of things that don't vary within group. In the real world, such things have importance. Thus, random effects are good for "modeling the data", while fixed effects models are good for getting closer to unbiased estimates of particular terms. With a random effects model, you can't make the claim to have removed that $e_i$ entirely.

In this example, time is the grouping variable. In your example, it is DID. (i.e.: it generalizes)


1) It is appropriate to make the comparison, just not with those two models. You would want to compare:

GLM <- glm(remission~Age+Married+IL6, data=hdp, family=binomial)


GLMM <- glmer(remission~Age+Married+IL6+(1|DID), data=hdp, family=binomial)

and you can do this with an anova:

anova(GLM, GLMM)

(Not sure if this will work with the glm and glmer results, as they might be different R objects. You might have to use two functions that have comparable return objects, like lme and gls, or do the anova yourself.)

The anova will do a log-likelihood ratio test to see if the addition of the random doctor effect is significant. You would need to divide that p-value by 2 before declaring significance because you are testing the null hypothesis that the random doctor effect is 0, and 0 is on the boundary of the parameter space for a variance (the actual distribution you are using in the test is a mixture of the $\chi^2_0$ and $\chi^2_1$ distribution -- but I'm near the boundary of my own ignorance at this point).

For me, the best book for understanding the process of nested model building and hypothesis testing has been West, Welsh, and Galecki (2007) Linear Mixed Models: A practical guide. They go through everything step by step.

2) If you have multiple observations per patient you would also add a random effect for patient. Then to test the relative importance of patience vs. doctor you could look at the predictive effects of patient vs. the predictive effects for doctor. The random effects terms for each will quantify the amount of variance between patients and between doctors, if that is a question you are interested in.

(Someone please correct me if I'm wrong!)

  • $\begingroup$ I'm not sure it makes sense to have DID as both a fixed effect, and a random intercept in the 2nd model. Moreover, having it as a fixed effect in the 1st model means that the choice b/t these 2 would be about which way to think about the effect of DID, not whether it needs to be included. On a different note, I notice you have an item (2); did you mean to have an item (1) somewhere? $\endgroup$ – gung Jun 18 '13 at 15:48
  • $\begingroup$ You're absolutely right; I was going from the OP's original glm formula which should not have had DID as a fixed effect in the 1st place. Now the choice is between whether treating DID as a random effect adds any value to the model. $\endgroup$ – Christopher Poile Jun 18 '13 at 15:55

The models are very different. The glm model is addressing the overall reduction in deviance (from a null-model) when all of the doctorID effects are being estimated and and are being assigned parameter estimates. You notice, of course, that Age, Married, and IL6 all have the same Wald statistics in the two model, right? My understanding (not a highly refined one I will admit) is that mixed model is treating the doctorIDs as nuisance factors or strata, namely "effects" that cannot be assumed to be drawn from any particular parent distribution. I see no reason to think that using a mixed model would improve your understanding of the "doctor-effect", quite the opposite in fact.

If your interest were in the effects of Age, Married or IL6 I would have imagined that you would not be comparing AIC across those two models but rather across differences in AIC with removal of covariates of interest within the same modeling structure.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.