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

Question

我是 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))
})

在我提出前一个问题 (How to extract marginal means of multiple variables with effect() function?) 之后，上面的代码在 @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() 来提取边际均值。

非常感谢！

艾拉

Answer 1

我不确定您希望最终的预期输出是什么样的。

或许，你可以试试这个方法：

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 替换为purrr 的map_df 来简化此操作：

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