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

lxkprmvk  于 12个月前  发布在  其他
关注(0)|答案(1)|浏览(107)

我试图概括一个函数,它将运行所有这些不同的模型组合,但对不同的响应变量。

`
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) `


Error in eval(predvars,data,env):object 'eelgrass_length' not found(eval(predvars,data,env):对象'eelgrass_length'未找到)
我如何推广它来遍历这些呢?

ssgvzors

ssgvzors1#

在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公式,并在lapply中循环lmer调用。然后我们可以直接在此列表上调用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))

相关问题