[R] lm, gee and lme

[Bliese, Paul D MAJ WRAIR-Wash DC]

Behavioral science data is often collected from nested structures (students
in schools, in districts, etc.). This can produce nonindependence among
responses from individuals in the same groups.  Consequently, researchers
are advised to model the nested nature of the data to avoid biases in SE
estimates.

Failing to account for nonindependence can lead to SE estimates that are too
large or too small depending ones model.  In the literature on linear mixed
effect models (lme) it is well documented that ignoring nonindependence
leads to SE values that are too small when the predictor varies only across
groups (is a level-2 predictor).  So, for instance, using linear regression
(lm) to regress individual performance (level-1 outcome) on group size
(level-2 predictor) using data where a lot of individuals come from the same
groups will lead to an SE value that is too small.  Using lme in this case
leads to a larger (and generally agreed upon as being less biased) SE
estimate.

What is trickier is what happens when both the predictor and outcome vary
among individuals in the same group.  An example of this type of model would
be one where individual performance is regressed upon individual age, but
where individual peformance is partially influenced by group membership.  My
understanding here is that ignoring nonindependence (i.e., using lm)
actually results in SE estimates that are too large, while modeling the
nonindependence reduces SE and increases power.

Here is an example:
# lme model
> mod.lme<-lme(GWB.ADD4~HOR,random=~1|GRP,data=TBH)
> VarCorr(mod.lme)
GRP = pdLogChol(1) 
            Variance  StdDev   
(Intercept) 0.3160445 0.5621783
Residual    0.7449425 0.8631005
> 0.3160445/(0.3160445+0.7449425)
[1] 0.2978778  #Note the large ICC (high nonindependence)

> summary(mod.lme)$tTable
                Value  Std.Error   DF  t-value       p-value
(Intercept) 1.9846214 0.06819493 7282 29.10219 3.041808e-176
HOR         0.2493643 0.01189157 7282 20.96984  7.533401e-95

#lm model
> summary(mod.lm)$coef
             Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 1.9404202 0.04437889 43.72395 0.000000e+00
HOR         0.2537195 0.01391557 18.23278 1.098509e-72

Notice the SE of .012 (lme) versus .014 (lm) and the higher t-value in lme.

These results make sense to me in that lme basically "sets aside" the
level-2 variance (tau) which is substantial in my example (see the high ICC
value).  As a result, the level-1 variable (HOR) only has to "explain"
level-1 variance (sigma-squared).  This means results in an increase in
power.

Now my question....why don't I see the same power increase if I use gee?
Notice below that the gee model SE values are basically identical to the lm
model values:

#gee model

>
mod.gee<-geese(GWB.ADD4~HOR,id=GRP,family="gaussian",corstr="exch",data=TBH)
> summary(mod.gee)
 Coefficients:
             estimate     san.se     wald p
(Intercept) 1.9847217 0.06992705 805.5802 0
HOR         0.2493582 0.01379557 326.7142 0
....
 Estimated Correlation Parameters:
       estimate     san.se     wald            p
alpha 0.3176938 0.04224159 56.56357 5.440093e-14

The SE value in the gee model (0.0138) is virtually identical to the SE
value in the lm model ignoring the nonindependence (0.0139).  Note that the
alpha estimate of .32 in gee is close to the ICC estimate in the lme model
(are these basically the same?).

I would have expected lme and gee to provide more similar answers.  Any
thoughts on why gee and lme are not more similar would be appreciated as
would any clarification about when ignoring nonindepedence leads to too
small versus too large SE values....

Paul

MAJ Paul Bliese, Ph.D.
Walter Reed Army Institute of Research
Phone: (301) 319-9873
Fax: (301) 319-9484
paul.bliese at na.amedd.army.mil