Fixed vs Random Effects

Question

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:

http://www.ats.ucla.edu/stat/mult_pkg/glmm.htm

http://www.ats.ucla.edu/stat/r/dae/melogit.htm

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)

Call:
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  

Coefficients:
              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?

Answer 1

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)

Answer 2

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)

with

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

Answer 3

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.