Cannot run lmer from within a function

Question

I am running into a problem trying to embed lmer within a function. Here is a reproducible example using data from lexdec. If I run lmer on the data frame directly, there is no problem. For example, say that I want to see whether reading times in a lexical decision task differed as a function of Trial. There were 2 types of word stimuli, "animal" (e.g. "dog") and "plant" (e.g. "cherry"). I can compute a mixed-effects model for animal words:

library(languageR)       #load lexdec data
library(lme4)            #load lmer()
s <- summary(lmer(RT ~ Trial + (1|Subject) + (1|Word), data = lexdec[lexdec$Class== "animal", ]))
s                        #this works well

However, if I embed the lmer model inside a function (say to not type the same command for each level of class) I get an error message. Do you know why? Any suggestions will be much appreciated!

#lmer() is now embedded in a function
compute.lmer <- function(df,class) {
  m <- lmer(RT ~ Trial + (1|Subject) + (1|Word),data = df[df$Class== class, ])
  m <- summary(m)
  return(m)
}

#Now I can use this function to iterate over the 2 levels of the **Class** factor
for (c in levels(lexdec$Class)){
 s <- compute.lmer(lexdec,c)
 print(c)
 print(s)
}

#But this gives an error message
Error in `colnames<-`(`*tmp*`, value = c("Estimate", "Std. Error", "df",  : 
  length of 'dimnames' [2] not equal to array extent 

Answer 1

I don't know what the problem is, your code runs just fine for me. (Are your packages up to date? What R version are you running? Have you cleaned your workspace and tried your code from scratch?)

That said, this is a great use case for plyr::dlply. I would do it like this:

library(languageR) 
library(lme4)
library(plyr)

stats <- dlply(lexdec,
      .variables = c("Class"),
      .fun=function(x) return(summary(lmer(RT ~ Trial + (1 | Subject) +
                                                (1 | Word), data = x))))

names(stats) <- levels(lexdec$Class)

Which then yields

> stats[["plant"]]
Linear mixed model fit by REML ['lmerMod']
Formula: RT ~ Trial + (1 | Subject) + (1 | Word)
   Data: x

REML criterion at convergence: -389.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2647 -0.6082 -0.1155  0.4502  6.0593 

Random effects:
 Groups   Name        Variance Std.Dev.
 Word     (Intercept) 0.003718 0.06097 
 Subject  (Intercept) 0.023293 0.15262 
 Residual             0.028697 0.16940 
Number of obs: 735, groups: Word, 35; Subject, 21

Fixed effects:
              Estimate Std. Error t value
(Intercept)  6.3999245  0.0382700  167.23
Trial       -0.0001702  0.0001357   -1.25

Correlation of Fixed Effects:
      (Intr)
Trial -0.379

When I run your code (copied and pasted without modification), I get similar output. It's identical except for the Data: line.

stats = list()  

compute.lmer <- function(df,class) {
    m <- lmer(RT ~ Trial + (1|Subject) + (1|Word),data = df[df$Class== class, ])
    m <- summary(m)
    return(m)
}

for (c in levels(lexdec$Class)){
    s <- compute.lmer(lexdec,c)
    stats[[c]] <- s
}

> stats[["plant"]]
Linear mixed model fit by REML ['lmerMod']
Formula: RT ~ Trial + (1 | Subject) + (1 | Word)
   Data: df[df$Class == class, ]

REML criterion at convergence: -389.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2647 -0.6082 -0.1155  0.4502  6.0593 

Random effects:
 Groups   Name        Variance Std.Dev.
 Word     (Intercept) 0.003718 0.06097 
 Subject  (Intercept) 0.023293 0.15262 
 Residual             0.028697 0.16940 
Number of obs: 735, groups: Word, 35; Subject, 21

Fixed effects:
              Estimate Std. Error t value
(Intercept)  6.3999245  0.0382700  167.23
Trial       -0.0001702  0.0001357   -1.25

Correlation of Fixed Effects:
      (Intr)
Trial -0.379