4

一个常见的任务是必须对数据集的不同子集执行某种统计分析（如 anova、glm 或混合模型），并将输出表与汇总系数和 p 值组合在单个数据框中。我正在寻找一个通用函数，该函数将采用模型类型（例如aov(...)orlm(...)或glm(...)or glmer(...)）和特定输出项，根据某些分组变量，必须为每个复制分析返回系数和 p 值（ s) 在一个人的数据集中。

假设我有一个数据帧，我想在其中对数据帧中不同级别的因子“复制”进行某种分析data：

data(iris)
library(car)
data=data.frame()
for (i in 1:10) {data=rbind(data,cbind(replicate=i,iris))}

使用broom+dplyr，我可以例如对此数据帧的每个子集进行方差分析（按复制分组）并使用“物种”一词保留 p 值

library(devtools)
install_github("dgrtwo/broom")
library(broom)
library(dplyr)

group_by(data, replicate) %>% do(tidy(Anova(aov(Sepal.Length ~ Species, data = .),type="III"))) %>% filter(term=="Species")

Source: local data frame [10 x 6]
Groups: replicate [10]

   replicate    term    sumsq    df statistic      p.value
       (int)   (chr)    (dbl) (dbl)     (dbl)        (dbl)
1          1 Species 189.6364     2  362.6614 2.580311e-94
2          2 Species 189.6364     2  362.6614 2.580311e-94
3          3 Species 189.6364     2  362.6614 2.580311e-94
4          4 Species 189.6364     2  362.6614 2.580311e-94
5          5 Species 189.6364     2  362.6614 2.580311e-94
6          6 Species 189.6364     2  362.6614 2.580311e-94
7          7 Species 189.6364     2  362.6614 2.580311e-94
8          8 Species 189.6364     2  362.6614 2.580311e-94
9          9 Species 189.6364     2  362.6614 2.580311e-94
10        10 Species 189.6364     2  362.6614 2.580311e-94

（我在这里使用了 10 个相同的数据子集作为示例）

我正在寻找一个更通用的函数“ Anovabygroup”，它将采用数据框、分组变量（这里replicate，但它也可以是几个分组变量的组合）、要运行的模型类型（例如在这个case 'aov(Sepal.Length ~ Species, data = .)'，但它也可以是 lm、glm、lme、lmer 或 glmer 模型或由Anova()) 以及返回系数和 p 值的因子（可能使用选项“all”返回所有内容）作为参数（给定的任何其他选项都可以传递给对 Anova 的调用）。任何人都知道如何使用与上面使用的代码类似的代码来执行此操作，但可以概括为采用这些参数？我不知道该怎么做的主要事情是将模型（例如在这种情况下为'aov（Sepal.Length ~ Species，data = .）'）作为参数传递并对其进行评估。或者它可能已经存在于某个包中？我认为这可能很有用，因为我总是发现自己一遍又一遍地编写这个任务......

PS 我使用了 github 版本的 broom 包，因为当前的 CRAN 版本似乎不能很好地处理 Anova 输出

4

1 回答 1

4

回答：

您可以通过创建一个解析文本输入的包装函数来解决这个问题。在这里，我使用parseand eval。

理想情况下，该函数还将检查传递的表达式（、、等）的有效性，但这里有一个示例函数可以帮助您入门glm。lm

它还允许您Anova根据要求将其他选项传递给：

anova_wrapper <- function(data, model_expression_as_string, grouping_variable,...) {
  f_wrap <- paste0('function(.) {',model_expression_as_string,'}') %>%    parse(text=.) %>% eval
  data %>% group_by_(grouping_variable) %>% 
     do(f_wrap(.) %>% Anova(...=...) %>% tidy) %>% return
}

例子：

使用您的示例代码（感谢您发布好的代码！）：

data=data.frame()
for (i in 1:10) {data=rbind(data,cbind(replicate=i,iris))}

aov_model_expression_as_string = 'aov(Sepal.Length ~ Species, data = .)'
lm_model_expression_as_string = 'lm(Sepal.Length ~ Sepal.Width + Petal.Length , data = .)'
grouping_variable = 'replicate'


data %>% 
    anova_wrapper(model_expression_as_string = aov_model_expression_as_string,
                  grouping_variable = grouping_variable,type="III")



    Source: local data frame [30 x 6]
    Groups: replicate [10]

    replicate        term      sumsq    df statistic       p.value
    (int)       (chr)      (dbl) (dbl)     (dbl)         (dbl)
    1          1 (Intercept) 1253.00180     1 4728.1630 1.134286e-113
    2          1     Species   63.21213     2  119.2645  1.669669e-31
    3          1   Residuals   38.95620   147        NA            NA
    4          2 (Intercept) 1253.00180     1 4728.1630 1.134286e-113
    5          2     Species   63.21213     2  119.2645  1.669669e-31
    6          2   Residuals   38.95620   147        NA            NA
    7          3 (Intercept) 1253.00180     1 4728.1630 1.134286e-113
    8          3     Species   63.21213     2  119.2645  1.669669e-31
    9          3   Residuals   38.95620   147        NA            NA
    10         4 (Intercept) 1253.00180     1 4728.1630 1.134286e-113
    ..       ...         ...        ...   ...       ...           ...

并使用 alm而不是aov和不同的参数Anova：

data %>% 
    anova_wrapper(model_expression_as_string = lm_model_expression_as_string,
                  grouping_variable = grouping_variable,type="III")



    replicate         term    sumsq    df statistic      p.value
    (int)        (chr)    (dbl) (dbl)     (dbl)        (dbl)
    1          1  Sepal.Width  8.19627     1  73.78707 1.163254e-14
    2          1 Petal.Length 84.42733     1 760.05861 5.847914e-60
    3          1    Residuals 16.32876   147        NA           NA
    4          2  Sepal.Width  8.19627     1  73.78707 1.163254e-14
    5          2 Petal.Length 84.42733     1 760.05861 5.847914e-60
    6          2    Residuals 16.32876   147        NA           NA
    7          3  Sepal.Width  8.19627     1  73.78707 1.163254e-14
    8          3 Petal.Length 84.42733     1 760.05861 5.847914e-60
    9          3    Residuals 16.32876   147        NA           NA
    10         4  Sepal.Width  8.19627     1  73.78707 1.163254e-14
    ..       ...          ...      ...   ...       ...          ...

我经常在蒙特卡罗模拟中重复这些类型的回归/方差分析，但我通常为每种分析创建一个单独的函数。R 社区可能会对重复执行这些分析的包感兴趣！

于 2015-11-20T19:10:51.553 回答