我如何编写一个函数来泛化具有不同响应变量的多个 lmer 模型？

Question

我正在尝试概括一个函数，该函数将在不同的响应变量上运行所有这些不同的模型组合。

`
model_suite <- function(response, df){
mod1.2 <- lmer(log(response) ~ traffic + (1 | site), data = df, na.action = na.exclude, REML = FALSE)}
#all eelgrass dry weight by site type with site as a random effect #
mod2.2 <- lmer(log(response) ~ type + (1 | site), data = df, na.action = na.exclude, REML = FALSE)
#all data by season with site as a random effect
mod3.2 <- lmer(log(response) ~ season + (1 | site), data = df, na.action = na.exclude, REML = FALSE)
#all data by traffic class and season with site as a random effect
mod4.2 <- lmer(log(response) ~ traffic * season + (1 | site), data = df, na.action = na.exclude, REML = FALSE)
mod5.2<- lmer(log(response) ~ traffic + season + (1 | site), data = df, na.action = na.exclude, REML = FALSE)
#all data by site type and season with site as a random effect
mod6.2 <- lmer(log(response) ~ type * season + (1 | site), data = df, na.action = na.exclude, REML = FALSE)
mod7.2 <- lmer(log(response) ~ type + season + (1 | site), data = df, na.action = na.exclude, REML = FALSE)
mod.sel <- model.sel(mod1.2, mod2.2, mod3.2, mod4.2, mod5.2, mod6.2, mod7.2)
print(mod.sel)
}`

这样当我运行这个函数时：

`model_suite(eelgrass_length, all_eelgrass) `

它返回此错误： eval(predvars, data, env) 中的错误：找不到对象“eelgrass_length”

我如何概括它以循环这些？

Answer 1

在 R 中使用公式时必须小心。公式中的名称不会按照您期望的方式进行替换。例如，如果我们这样做：

my_formula <- function(variable) {
  formula(variable ~ other_variable)
}

然后运行：

my_formula(banana)
#> variable ~ other_variable
#> <environment: 0x000001de2e0d8900>

我们可以看到公式的左侧仍然是

variable

- 它并未如您所期望的那样替换为

banana

。

有多种方法可以解决这个问题。一种是使用

deparse(substitute(...))

从输入名称创建字符串，然后使用

as.formula

将其转换为公式。

所以我们可以将上面的函数改成如下：

my_formula <- function(variable) {
  
  as.formula(paste(deparse(substitute(variable)), '~ other_variable'))
}

现在函数的输入将被代入公式中

my_formula(banana)
#> banana ~ other_variable
#> <environment: 0x000001de2e379048>

这正是您的功能的问题。

response

没有在您的公式中被替换，因此您的函数实际上是在查找名为

response

的列，该列不存在。

使用

as.formula

创建公式的机制还使我们有机会删除函数中的一些重复项。我们可以创建一个公式

list

并在

lmer

内循环调用

lapply

。然后我们可以直接在此列表上调用

model.sel

。

library(lme4)
library(MuMIn)

model_suite <- function(response, df){
  lhs <- paste0('log(', deparse(substitute(response)), ') ~')
  rhs <- paste(c('traffic', 'type', 'season', 'traffic * season', 
                 'traffic + season', 'type * season', 'type + season'),
               '+ (1 | site)')
  formulas <- lapply(paste(lhs, rhs), as.formula)
  models <- lapply(formulas, function(x) {
    lmer(x, data = df, na.action = na.exclude, REML = FALSE)
  })
  do.call('model.sel', models)
}

这允许：

model_suite(eelgrass_length, all_eelgrass)
#> Model selection table 
#>     (Int) trf typ ssn ssn:trf ssn:typ df   logLik  AICc delta weight
#> 3 -0.7877           +                  4 -138.417 285.3  0.00  0.427
#> 6 -0.3129       +   +               +  8 -134.480 286.7  1.43  0.209
#> 2 -0.8611       +                      5 -138.462 287.6  2.33  0.133
#> 7 -0.7211       +   +                  6 -137.754 288.5  3.22  0.086
#> 1 -0.8698   +                          5 -138.923 288.6  3.25  0.084
#> 5 -0.7298   +       +                  6 -138.223 289.5  4.15  0.054
#> 4 -0.7336   +       +       +          8 -137.701 293.2  7.87  0.008
#> Models ranked by AICc(x) 
#> Random terms (all models): 
#>   1 | site

---

使用的数据

我们没有您用于运行此函数的数据，并且该函数似乎依赖于传递的数据帧中的特定列。我创建了一个与您的数据同名的小数据框，以便上面的代码也适用于您的真实数据：

set.seed(123)
all_eelgrass <- expand.grid(traffic = c('light', 'medium', 'heavy'),
                            type = c('type1', 'type2', 'type3'),
                            season = c('dry', 'rainy'),
                            site = c('1', '2', '3', '4', '5'))

all_eelgrass$eelgrass_length <- rexp(90, c(0.8, 1, 1.2, 1, 1.2, 1.4, 1.2, 1.4,
    1.6, 1, 1.2, 1.4, 1.2,  1.4, 1.6, 1.4, 1.6, 1.8, 1, 1.2, 1.4, 1.2, 1.4, 
    1.6, 1.4, 1.6, 1.8, 1.2, 1.4, 1.6, 1.4, 1.6, 1.8, 1.6, 1.8, 2, 1.2, 1.4, 
    1.6, 1.4, 1.6, 1.8, 1.6, 1.8, 2, 1.4, 1.6, 1.8, 1.6, 1.8, 2, 1.8, 
    2, 2.2, 1.4, 1.6, 1.8, 1.6, 1.8, 2, 1.8, 2, 2.2, 1.6, 1.8, 2, 1.8, 2, 2.2, 
    2, 2.2, 2.4, 1.6, 1.8, 2, 1.8, 2, 2.2, 2, 2.2, 2.4, 1.8, 2, 2.2, 2, 2.2, 
    2.4, 2.2, 2.4, 2.6))