Невозможно запустить lmer из функции

Я столкнулся с проблемой, пытаясь встроить lmer в функцию. Вот воспроизводимый пример с использованием данных из lexdec. Если я запускаю lmer напрямую во фрейме данных, проблем не возникает. Например, предположим, что я хочу узнать, различалось ли время чтения в задаче на лексическое решение в зависимости от пробной. Было 2 типа словесных стимулов: «животное» (например, «собака») и «растение» (например, «вишня»). Я могу вычислить модель смешанных эффектов для слов животных:

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

Однако, если я встрою модель lmer в функцию (скажем, чтобы не вводить одну и ту же команду для каждого уровня класса), я получаю сообщение об ошибке. Знаете почему? Мы будем очень признательны за любые предложения!

#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

r lme4
person Sol    schedule 20.03.2014    source источник
comment
Это резонный вопрос, но я думаю, что заголовок немного вводит в заблуждение - проблема не в сохранении результатов в список, а в запуске lmer из функции, которой передаются данные.. .   -  person Ben Bolker    schedule 21.03.2014
comment
спасибо Бен. Я попытался переформулировать вопрос, чтобы сделать его более точным.   -  person Sol    schedule 21.03.2014

Ответы (1)


arrow_upward
2
arrow_downward

Я не знаю, в чем проблема, ваш код работает отлично для меня. (Обновлены ли ваши пакеты? Какую версию R вы используете? Вы очистили свое рабочее пространство и попробовали свой код с нуля?)

Тем не менее, это отличный вариант использования plyr::dlply. Я бы сделал это так:

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)

Что тогда дает

> 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

Когда я запускаю ваш код (скопированный и вставленный без изменений), я получаю аналогичный результат. Он идентичен, за исключением строки Data:.

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
person Gregor Thomas    schedule 20.03.2014
comment
отлично. не знал о dlply - person rawr; 21.03.2014
comment
согласен, у меня работает. Меня не удивляет, что это могло потерпеть неудачу в более старых/предыдущих версиях; заставить lmer искать данные в нужной среде, независимо от того, из какой среды они были первоначально вызваны, сложнее, чем вы думаете... - person Ben Bolker; 21.03.2014
comment
@shujaa: спасибо! Вы были правы: проблема исчезла, когда я обновил lme4. Я не обновлял его, так как прочитал, что lme4 1.0 иногда дает модели с худшим соответствием, чем предыдущая версия ссылка. Способ dlply тоже хорош! - person Sol; 21.03.2014
comment
@rawr Это самое замечательное в plyr, в каждой комбинации (a,d,m,l,r)(a,d,l,_)ply, (input)(output)ply. - person Gregor Thomas; 21.03.2014