从TukeyHSD检验中提取“p adj”值

Question

使用虚拟数据集：

Species      Var1     Var2   Var3
   a          1         2      3
   a          4         5      6 
   b          7         8      9
   b          10       11      12

我有多个物种和大约50个变量（Var50）。 我想对每个响应变量的配对分组变量（物种）执行单向方差分析，并获取在95％ CI下具有统计显着性的频率输出。例如，我开始编写以下函数来执行此操作：

data<-read.table("example.txt", header=T, sep="\t")
function(y){
for(y in 2:50)
anova.r<-aov(y~Species, data = data)
result<-TukeyHSD(anova.r, conf.level = 0.95) 
f.result ## I cannot figure out how to extract the "p adj" from the results

f.result<-sum(prob.result>=0.05)
write.table(f.result, file = "anova95.csv", sep = ",",
        col.names = FALSE, append=TRUE)
}  

最终，我希望最终的表格（虚拟答案）看起来像这样：

                     Var1   Var2   Var3......Var50 
Frequency at 95% CI   106    200    45         246 

我知道可以使用 [[]] 来访问 Tukey 检验结果中的数据。我已经尝试了 tukey.results[[1]][,1] 到 tukey.results[[1]][,3]，但都不起作用。 tukey.results[[1]] 返回了 Tukey 检验的所有列。

此外，我认为我可能需要在函数中使用 cbind 将数据放入各自的列中。或者我想可能可以使用apply 命令，但是我不知道如何在每次迭代时保持分组变量不变而改变响应变量。

任何建议将不胜感激。

- user2507608

---

result$Species[4] 将会给你 "p adj"。 - eddi

@eddi，谢谢您的建议。我已经尝试了，但是我只得到了一个值（并不是“p adj”中的第四个值）。 - user2507608

只需查看 result$Species 的内容，您应该能够在那里找到所需的内容；如果没有找到，那么我猜测您的示例并不代表您真正遇到的问题。 - eddi

result$Species 给出了所有预期值。前几行是 diff lwr upr p adj A oil palm2-A oil palm1 -0.021441963 -0.06055280 0.01766887 0.8704602 AfricanOilPalm_1-A oil palm1 -0.001384416 -0.04049525 0.03772642 1.0000000 AfricanOilPalm_2-A oil palm1 -0.017210672 -0.05632151 0.02190016 0.9776883。问题可能是 result 的结构吗？我拥有的结构是 str(result$Species) num [1:120, 1:4]。 - user2507608

尝试使用 result$Species[,'p.adj'] 或 result$Species[, 4]。 - eddi

感谢eddi, result$Species[, 4] 有点用，但为了获得一列中成对比较的列表和另一列中相应的“p adj”值，我不得不使用来自reshape包的 melt(result$Species[, 4]). 不过还是非常感谢。现在我必须继续使用该函数，尝试获得我的最终结果。再次感谢。 - user2507608

Answer 1

如果你也在寻找变量，可以尝试以下方法：

summary(fm1 <- aov(breaks ~ wool + tension, data = warpbreaks))
kk<-TukeyHSD(fm1, "tension", ordered = TRUE)
kk$tension
result<-data.frame( kk$tension)
result["p.adj"]

            p.adj
M-H 0.447421021
L-H 0.001121788
L-M 0.033626219

- Metrics

---

感谢您提供的示例和帮助。但是，当我尝试将其应用到我的数据上时，答案却是[1] NA。与data = warpbreaks数据集不同，我的数据只有一组因子变量。我不确定这是否会有所不同。 - user2507608

你可能需要使用dput命令来提供可重现的代表性数据。在此之前，请确保您理解使用aov的原理和要求。 - Metrics

我理解aov的原理和要求。我认为我打开了太多的软件包，所以我会在另一台电脑上尝试。谢谢你的帮助。 - user2507608

@Metrics 这个方法只适用于单向方差分析吗？当我在我的数据上尝试时，它只返回了空值。我正在尝试让我的图基结果易于理解。 - jbest

Answer 2

这个答案在我的系统上不起作用。以下是我基于Metrics提供的代码提供的解决方案。

 summary(fm1 <- aov(breaks ~ wool + tension, data = warpbreaks))
 kk<-TukeyHSD(fm1, "tension", ordered = TRUE)
 kk<-kk$tension          #strips off some headers in kk
 kk<-as.data.frame(kk)   #converts to data frame
 kk<-kk$'p adj'          #selects relevant output
 print(kk)               #to check answer

- TimothyEbert