使用 purr 和 map (tidyverse) 迭代因变量

iterate dependent variable using purr and map (tidyverse)

我有一个简单的数据集,我想使用 aovtidyverse 迭代因变量。然后我想从这些输出中计算 Tukey HSD 测试。我在 for 循环结构中工作,但我正在尽最大努力摆脱这种心态。我看到 this post 关于使用自变量迭代 aov 函数。试图将此逻辑合并到我的工作流程中,但效果不佳。任何 tidyverse 爱好者可以在这里引导我朝着正确的方向前进吗?

library(tidyverse)
library(data.table)

pfuel <- fread("data/CFL.csv") %>%
  mutate(AFCL = AFCL*10,
         LCW = LCW*10,
         DCW = DCW*10,
         LiDe = ifelse(Status == "Li", "Live", "Dead")) %>%
  filter(S.F == "S") %>%
  group_by(Site, Year, Age, Plot) %>%
  select(LiFol, DeFol, Li.1hr, De.1hr, Li.10hr, De.10hr, Li.100hr, De.100hr) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  mutate(sb_age = paste0(Year, Age))

aov.models = pfuel %>%
  select (-c(Year, Age)) %>%
  select(LiFol, DeFol, Li.1hr, De.1hr, Li.10hr, De.10hr, Li.100hr, De.100hr, Site, Plot, sb_age) %>%
  map(~ aov(.x ~ sb_age + Site/Plot, data = pfuel))

aov.models 运行时,我生成此错误:

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
  NA/NaN/Inf in 'y'
In addition: Warning message:
In model.response(mf, "numeric") : NAs introduced by coercion

我还没有参加 Tukey 测试,因为我无法通过 aov 函数。任何建议将不胜感激!

您可以在这里找到数据:https://www.dropbox.com/s/yb8rh860fc7fff2/CFL.csv?dl=0

谢谢!

@Z.Lin 在您的指导下,我找到了解决问题第一部分的方法。可能不是最优雅的,但它至少现在可以工作了!欢迎任何改进,但谢谢。

pfuel_var <- pfuel %>%
  select(Site, Plot, sb_age) %>%
  mutate(Site = as.factor(Site),
         Plot = as.factor( Plot),
         sb_age = as.factor(sb_age))

aov.models <- pfuel %>%
  select(LiFol, DeFol, Li.1hr, De.1hr, Li.10hr, De.10hr, Li.100hr, De.100hr) %>%
  map(~ aov(.x ~ pfuel_var$sb_age + pfuel_var$Site/pfuel_var$Plot, data = pfuel))

我问题的第二部分是如何将此输出从 agricolae 包提供给 HSD.test。有人对此有想法吗?

我当时的想法是:

t <- aov.models %>% 
   map(~ HSD.test(.x, "pfuel_var$sb_age", alpha=0.1))

但这不能正常工作。非常感谢。

将数据转换为长格式可能更容易,按响应拆分,然后拟合模型并将输出提供给 HSD.test 函数,例如

aov.models <- pfuel %>%
  select(-Year, -Age) %>%
  gather(variable, value, -sb_age, -Site, -Plot) %>%
  split(.$variable) %>%
  map(~ aov(value ~ sb_age + Site/Plot, data = .x)) %>%
  map(HSD.test, trt = 'sb_age')

我还删除了其中一个 select() 语句,因为它选择了所有列。