2
votes

Consider a dataset like this:

library(lme4)
data(sleepstudy)
head(sleepstudy, 10)

   Reaction Days Subject
1  249.5600    0     308
2  258.7047    1     308
3  250.8006    2     308
4  321.4398    3     308
5  356.8519    4     308
6  414.6901    5     308
7  382.2038    6     308
8  290.1486    7     308
9  430.5853    8     308
10 466.3535    9     308

I'm trying to run some mixed-effects linear models when stratifying based on a grouping variable. I can do:

library(dplyr)
libary(tidyr)
library(purrr)

sleepstudy %>%
 group_by(grp = gl(2, n()/2)) %>%
 nest() %>%
 mutate(fit = map(data, ~ lmer(Reaction ~ Days + (Days | Subject), data = .)))

obtaining:

  grp             data fit       
  <fct> <list<df[,3]>> <list>    
1 1           [90 × 3] <lmrMdLmT>
2 2           [90 × 3] <lmrMdLmT>

Now I would like to use the effect() function from effects library:

library(effects)

sleepstudy %>%
 group_by(grp = gl(2, n()/2)) %>%
 nest() %>%
 mutate(fit = map(data, ~ effect(term = "Days",
                                 mod = lmer(Reaction ~ Days + (Days | Subject), data = .))))

however, I'm always getting an error:

Error in is.data.frame(data) : object '.' not found

When I use it outside of tidyverse context:

sleepstudy2 <- sleepstudy %>%
 mutate(grp = gl(2, n()/2))

effect(term = "Days", mod = lmer(Reaction ~ Days + (Days | Subject), 
                                 data = sleepstudy2,
                                 subset = grp == 1))

then it works normally:

 Days effect
Days
       0        2        4        7        9 
252.2916 267.0686 281.8455 304.0110 318.7880 

Any idea how to fix this problem and do it using the tidy approach?

1

1 Answers

3
votes

So the effects() function seems to have a real problem when the data you are using is not in the global environment. It's referenced in the "Warnings and Limitations" section of the ?effects help page. They point to the "embedding" vignette from the car packages for possible work around. The workaround they suggestion is making a copy of data data in the global environment. It's not elegant but this works

poo <- sleepstudy %>%
  group_by(grp = gl(2, n()/2)) %>%
  nest() %>%
  mutate(fit = map(data, function(x) {
    assign(".dta", x, env=.GlobalEnv)
    eff <- effect(term="Days", mod=lmer(Reaction ~ Days + (Days | Subject), data=.dta))
    remove(".dta", envir=.GlobalEnv) 
    eff
  }))