如何使用 effect() 函数提取多个变量的边际均值?

How to extract marginal means of multiple variables with effect() function?

我是 R 的初学者,我想在包含超过 200 列结果变量的数据集中执行 ANCOVA。 对我来说最重要的是提取这些变量的 p 值和边际均值。 我在 lapply() 函数的帮助下成功提取了 p 值,但是当我提取边际意味着我得到了这样的错误 Error in eval(predvars, data, env) : object 'x' not found.

这里我以内置数据集“iris”为例来展示我的问题

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)

#同时对多个结果变量进行ANCOVA(这里我比较了不同Specie水平下的多个结果变量,以Petal.Width作为协变量)

list1 <- lapply(iris[, 1:3], function(x) Anova(aov(x ~ Petal.Width  + Species, data = iris), type="III"))
str(list1)

#提取主要测试的 p 值

pvalues <- stack(lapply(iris[, 1:3], function(x) Anova(aov(x ~ Petal.Width  + Species, data = iris), type="III")[3, 4]))[2:1]

上面的代码运行良好,但是当我使用 effect() 函数提取边际均值时出现错误: #提取边际均值

list2 <- lapply(iris[, 1:3], function(x) summary(effect("Species", aov(x ~ Petal.Width + Species, data = iris)), se=TRUE))

Error in eval(predvars, data, env) : object 'x' not found

marginal.means <- stack(lapply(iris[, 1:3], function(x) summary(effect("Species", aov(x ~ Petal.Width + Species, data = iris)), se=TRUE)[[5]][[1]][1]))[2:1]

Error in eval(predvars, data, env) : object 'x' not found

#当我提取某个变量的边际均值时(例如Sepal.Length),不使用

marginal.mean1 <- summary(effect("Species", aov(Sepal.Length ~ Petal.Width + Species, data = iris)), se=TRUE)
marginal.mean1

输出:

 Species
    setosa versicolor  virginica 
  5.880113   5.819859   5.830028 

 Lower 95 Percent Confidence Limits
Species
    setosa versicolor  virginica 
  5.490905   5.676927   5.485953 

 Upper 95 Percent Confidence Limits
Species
    setosa versicolor  virginica 
  6.269322   5.962791   6.174102

由于结果变量超过200列,我想一次提取边际均值而不是一个一个地提取。

非常感谢您的帮助,

艾拉

您收到该错误是因为函数 effect() 调用 update() 并尝试重新拟合模型,此时,它无法再访问您的 x。 (好吧,也许我没有解释得太清楚)你可以阅读这篇文章 book chapter 来了解函数是如何工作的。

尝试将所有内容保留在 data.frame 内,而不是提供适合不同变量的公式:

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

如您所见,这也可以应用于您的其他功能。