0

我需要在 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
4

2 回答 2

0

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

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

(Y)

于 2018-01-20T20:49:47.050 回答
0

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

如果是这样的话，那么将这三个因素组合成一个命名CLASS只会让事情变得更加困难，因为这使得单独跟踪因素的水平变得更加困难。更简单的是拟合考虑所有三个因素的模型，估计每个因素组合的均值，然后使用by变量进行所需的比较。因此，对于每个 dataset 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 值被寻求。这让我很恼火。我不喜欢看到人们不看被测试的数量就直接进行统计测试。

于 2018-01-14T22:23:14.743 回答