数据子集的成对方差分析答案

【问题标题】：pairwise ANOVA of subset of data数据子集的成对方差分析
【发布时间】：2017-12-19 19:01:48
【问题描述】：

我需要在 R 中执行多个成对方差分析，并使用 bonferroni 更正 p 值。但是，我不需要将每个 CLASS 相互比较。下面是我的数据格式和selcontrasts：我需要对比log10relquant的那对。你们中有人知道我该如何执行吗？我使用dplyr、lsmeans 和broom 包。

SEX      EXPERIENCED    AGE  CLASS compound    relquant log10relquant

1 FEMALE          NO     1D     1F      C14 0.004012910     -2.396541
2 FEMALE          NO     1D     1F      C14 0.003759812     -2.424834
3 FEMALE          NO     1D     1F      C14 0.003838553     -2.415832
4 FEMALE          NO     1D     1F      C14 0.003582754     -2.445783
5   MALE          NO     1D     1M      C14 0.005099237     -2.292495
6   MALE          NO     1D     1M      C14 0.005379093     -2.269291

selcontrasts <- c("1F - 1M", "4F - 4M", "4EF - 4EM", 
                  "7F - 7M", "7EF - 7EM", # sex differences
              "1M - 4M", "4M - 7M", "1M - 7M", "1F - 4F", 
              "4F - 7F", "1F - 7F", # age differences
              "4M - 4EM", "7M - 7EM", "4F - 4EF", 
              "7F - 7EF" # social experience)

x=list(selcontrasts)

目前我正在使用它来配对整个数据集（以便比较每个类）而不是选定的对比：

pvalsage=data.frame(datagr %>% 
    do( data.frame(summary(contrast(lsmeans(
          aov(log10relquant ~ CLASS, data = .), ~ CLASS ),               
          method="pairwise",adjust="none"))) ))

为了只在列表 x 中选择对比，我尝试了：

pvalsage=data.frame(datagr %>%  
    do(data.frame(summary(contrast(lsmeans(
        aov(log10relquant ~ CLASS, data = .),~ CLASS),
        method = x, adjust="none"))) ))

但我得到了错误：

 error in contrast.ref.grid(lsmeans(aov(log10relquant ~ CLASS, data = .),  : 
 Nonconforming number of contrast coefficients

【问题讨论】：

不要使用”pairwise” in the contrast()` 函数——给出一个具有所需对比度系数的命名欲望。请参阅文档。此外，您可能应该升级到继续 lsmeans 的 emmeans 包
什么是对比度系数？
列表看起来如何？
我不断收到错误消息：'不合格的对比度系数数'
你必须给出一个数字向量的命名列表，每个向量的长度等于被对比的平均值的数量。查看contrast 帮助页面中的示例。例如，要比较六个均值中的前两个，该列表元素可能看起来像 con12 = c(-1,1,0,9,0,0)

标签： r subset anova contrast pairwise

【解决方案1】：

如果我正确理解了这个问题（我很可能不会），实际上涉及三个因素：SEX（两个级别）、EXPERIENCED（两个级别）和AGE（三个级别，即1、4 和 7)。并且需要对其他两个因素的每个组合分别比较每个因素的水平。

如果是这种情况，那么将这三个因素合并为一个名为 CLASS 的因素只会让事情变得更加困难，因为单独跟踪因素的水平会变得更加困难。更简单的是拟合考虑所有三个因素的模型，估计每个因素组合的均值，然后使用by 变量进行所需的比较。因此，对于每个数据集 dat，您可以：

require(emmeans)
mod = aov(log10relquant ~ SEX * EXPERIENCED * AGE, data = dat)
emm = emmeans(mod, ~ SEX * EXPERIENCED * AGE)
rbind(pairs(emm, by = c("EXPERIENCED", "AGE")),
      pairs(emm, by = c("SEX", "EXPERIENCED")),
      pairs(emm, by = c("SEX", "AGE")),
      adjust = "bonferroni")

我没有尝试将其嵌入到函数式编程范式中；我会把它留给 OP 来弄清楚这些细节。

注意：emmeans 包（估计边际均值）是 lsmeans 的延续，将来会被弃用。它的工作方式相同。

PS -- 查看问题中的代码，我担心最终结果不会显示正在比较的实际估计（EMM），只有比较；以及命名中的进一步暗示，实际上，只有 P 值被寻求。这让我很恼火。我不喜欢看到人们不看被测试的数量就直接进行统计测试。

【讨论】：

【解决方案2】：

无论如何，您都可以进行成对对比，然后将仅包含您的 selcontrasts 的行过滤到一个新的数据框中，然后将 p.adjust= bonferroni 与您感兴趣的对比。

或者您可以编写一个 mycontr.lmsc 函数并定义 selcontrasts 并将其用作方法 =

(Y)

【讨论】：