3

I am currently working through Andy Field's book, Discovering Statistics Using R. Chapter 14 is on Mixed Modelling and he uses the lme function from the nlme package.

The model he creates, using speed dating data, is such:

speedDateModel <- lme(dateRating ~ looks + personality +
                 gender + looks:gender + personality:gender + 
                 looks:personality,
    random = ~1|participant/looks/personality)

I tried to recreate a similar model using the lmer function from the lme4 package; however, my results are different. I thought I had the proper syntax, but maybe not?

speedDateModel.2 <- lmer(dateRating ~ looks + personality + gender + 
              looks:gender + personality:gender + 
              (1|participant) + (1|looks) + (1|personality), 
              data = speedData, REML = FALSE)

Also, when I run the coefficients of these models I notice that it only produces random intercepts for each participant. I was trying to then create a model that produces both random intercepts and slopes. I can't seem to get the syntax correct for either function to do this. Any help would be greatly appreciated.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
user3585829
  • 945
  • 11
  • 24

1 Answers1

5

The only difference between the lme and the corresponding lmer formula should be that the random and fixed components are aggregated into a single formula:

dateRating ~ looks + personality +
             gender + looks:gender + personality:gender + 
             looks:personality+ (1|participant/looks/personality)

using (1|participant) + (1|looks) + (1|personality) is only equivalent if looks and personality have unique values at each nested level.

It's not clear what continuous variable you want to define your slopes: if you have a continuous variable x and groups g, then (x|g) or equivalently (1+x|g) will give you a random-slopes model (x should also be included in the fixed-effects part of the model, i.e. the full formula should be y~x+(x|g) ...)

update: I got the data, or rather a script file that allows one to reconstruct the data, from here. Field makes a common mistake in his book, which I have made several times in the past: since there is only a single observation in the data set for each participant/looks/personality combination, the three-way interaction has one level per observation. In a linear mixed model, this means the variance at the lowest level of nesting will be confounded with the residual variance.

You can see this in two ways:

  • lme appears to fit the model just fine, but if you try to calculate confidence intervals via intervals(), you get
 intervals(speedDateModel)
 ## Error in intervals.lme(speedDateModel) : 
 ##   cannot get confidence intervals on var-cov components: 
 ##   Non-positive definite approximate variance-covariance
  • If you try this with lmer you get:
## Error: number of levels of each grouping factor
##   must be < number of observations

In both cases, this is a clue that something's wrong. (You can overcome this in lmer if you really want to: see ?lmerControl.)

If we leave out the lowest grouping level, everything works fine:

sd2 <- lmer(dateRating ~ looks + personality +
                 gender + looks:gender + personality:gender + 
                 looks:personality+
                     (1|participant/looks),
            data=speedData)

Compare lmer and lme fixed effects:

all.equal(fixef(sd2),fixef(speedDateModel)) ## TRUE

The starling example here gives another example and further explanation of this issue.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Thank you, Ben. In my first model example using the lme() function it produces the correct result. Would you be able to help me understand how to re-write that example with the lmer() function as I can't seem to wrap my head around how to make it work without getting a bunch of errors. – user3585829 Nov 07 '15 at 20:22
  • if I'm going to help you more I need a reproducible example. I poked at the book a little bit on Google Books, but I can't figure out where the `speedData` data set comes from (also, you should always use the `data=` argument to explicitly specify the data frame for analysis ...) – Ben Bolker Nov 07 '15 at 20:57
  • Hello, Ben, sorry for not getting back to this yet. Things got really busy. Will try and update with more info later this week if possible. Appreciate it. – user3585829 Nov 10 '15 at 22:59