如何正确使用 stack() 函数在 R 中提取 ANCOVA 的边际均值?

How to use stack() function correctly to extract marginal means of ANCOVA in R?

我是 R 的初学者,我想从对 200 多个结果变量执行的 ANCOVA 测试中提取边际均值。当我只对一个结果变量使用 stack() 时效果很好,但是当我同时使用 stack()lapply().

时出现错误

这里我使用内置数据集“iris”来展示问题。 数据集“iris”在 Species 有三个级别,我使用 Petal.Width 作为协变量,Species 作为预测变量,前三列变量作为结果变量。

我的目的是同时提取对应结果变量的多个边际均值,而不是一个一个提取。

#load data and packages
data("iris")
library(car); library(compute.es); library(effects); library(ggplot2);
library(multcomp); library(pastecs); library(WRS)

#set contrasts for the following ANCOVA tests
contrasts(iris$Species) <- contr.poly(3)

#perform 
list2 <- lapply(colnames(iris)[1:3], function(x){
anova_fit = aov(reformulate(c("Petal.Width","Species"),x), data = iris)
summary(effect("Species",anova_fit, se=TRUE))
})

在我提出前一个问题 () 之后,在@StupidWolf 的帮助下,上面的代码运行良好。 然后当我执行以下代码时出现错误:

means.all <- stack(lapply(colnames(iris)[1:3], function(x){
anova_fit = aov(reformulate(c("Petal.Width","Species"),x), data = iris)
summary(effect("Species",anova_fit, se=TRUE))[[5]][1]
}))[2:1]

错误是 Error in rep.int(factor(names(x), unique(names(x))), lengths(x)) : invalid 'times' value.

但是当我只提取一个结果变量的边际均值时(以Sepal.Length为例),我可以用下面的代码提取边际均值:

anova_fit = aov(reformulate(c("Petal.Width","Species"),"Sepal.Length"), data = iris)
means1 <- summary(effect("Species",anova_fit, se=TRUE))[[5]][1]

我不知道如何正确使用 stack()lapply() 来提取边际均值。

非常感谢!

艾拉

我不确定您希望最终的预期输出如何。

或许,你可以试试这个方法:

do.call(rbind, lapply(list2, function(x) 
  data.frame(prop = c('effect', 'lower', 'upper'), 
            rbind(x$effect, x$lower, x$upper))))


#    prop setosa versicolor virginica
#1 effect   5.88       5.82      5.83
#2  lower   5.49       5.68      5.49
#3  upper   6.27       5.96      6.17
#4 effect   4.17       2.67      2.33
#5  lower   3.93       2.58      2.11
#6  upper   4.42       2.76      2.54
#7 effect   2.43       4.13      4.71
#8  lower   2.13       4.02      4.44
#9  upper   2.74       4.24      4.98

您还可以通过将 do.call + rbind 替换为 purrrmap_df 来简化此操作:

purrr::map_df(list2, function(x) data.frame(prop = c('effect', 'lower', 'upper'), 
                                 rbind(x$effect, x$lower, x$upper)))