从统计测试中提取输出

Question

我在R(在RStudio中运行)中运行统计测试。我使用变量名保存结果。我想提取结果的一部分。我该怎么做呢？下面是一个例子，最后是R代码。我设置了一个有四个处理的实验，并收集了数据。接下来，我运行ANOVA并执行Tukey HSD测试。结果存储在一个名为“posthoc”的变量中。

我查看并注意到RStudio是一个1的列表。在posthoc中，我看到名称左侧有一个蓝色小箭头，单击它会显示更多信息。我不确定如何解释它，以便我可以用来回答我自己的问题。

我可以打印(Posthoc)并得到以下结果。

#  Tukey multiple comparisons of means
#    95% family-wise confidence level
#
#Fit: aov(formula = Expt1$Treat1 ~ Expt1$Trt)
#
#$`Expt1$Trt`
#          diff        lwr          upr     p adj
#B-A   6.523841   2.664755  10.38292569 0.0001372
#C-A  18.584160  14.725075  22.44324507 0.0000000
#D-A   2.643719  -1.215367   6.50280370 0.2854076
#C-B  12.060319   8.201234  15.91940456 0.0000000
#D-B  -3.880122  -7.739207  -0.02103681 0.0482260
#D-C -15.940441 -19.799527 -12.08135619 0.0000000

我也可以输入class(posthoc)，结果是:1 "TukeyHSD“"multicomp”

在本例中，我需要的是一个新变量中的所有p值。一般的问题是R给了我输出，我需要知道如何从输出中提取特定的元素。我可能在使用aov，lm，nlme，或者其他的东西。

Mean1=3.2
Sd1=3.2
Mean2=9.4
Sd2=2.4
Mean3=21.4
Sd3=6.4
Mean4=3.9
Sd4=10.7

Size1=30

Treat1=rnorm(Size1,mean=Mean1, sd=Sd1)
Trt="A"
Treat1M <- data.frame(Treat1, Trt)
Treat1=rnorm(Size1,mean=Mean2, sd=Sd2)
Trt="B"
Treat2M <- data.frame(Treat1, Trt)
Treat1=rnorm(Size1,mean=Mean3, sd=Sd3)
Trt="C"
Treat3M <- data.frame(Treat1, Trt)
Treat1=rnorm(Size1,mean=Mean4, sd=Sd4)
Trt="D"
Treat4M <- data.frame(Treat1, Trt)

Expt1=rbind(Treat1M, Treat2M, Treat3M, Treat4M)

Expt1R<-aov(Expt1$Treat1 ~ Expt1$Trt)
posthoc <-TukeyHSD(x=Expt1R, 'Expt1$Trt', conf.level=.95)

Answer 1

只需使用posthoc变量的子集。

posthoc$`Expt1$Trt`[,4]

或者，您可以尝试使用broom包。

library(broom)
res <- tidy(posthoc)
res
       term comparison   estimate    conf.low conf.high  adj.p.value
1 Expt1$Trt        B-A   5.904138   1.3639293 10.444346 5.223255e-03
2 Expt1$Trt        C-A  16.886340  12.3461316 21.426548 3.919087e-14
3 Expt1$Trt        D-A   4.283597  -0.2566111  8.823805 7.189220e-02
4 Expt1$Trt        C-B  10.982202   6.4419940 15.522410 3.226398e-08
5 Expt1$Trt        D-B  -1.620540  -6.1607487  2.919668 7.886097e-01
6 Expt1$Trt        D-C -12.602743 -17.1429509 -8.062534 3.235351e-10

tidy函数的输出是一个data.frame。因此，您可以使用res$adj.p.value访问p值。

class(res)
[1] "data.frame"

Answer 2

金宝已经提供了一个很好的解决方案。我也会选择broom，特别是在用ggplot2做绘图的时候。我只想阐述一下John Coleman的评论。

您可以使用str检查对象。在posthoc的情况下，

str(posthoc)

给出

List of 1
 $ Expt1$Trt: num [1:6, 1:4] 6.46 18.19 -0.76 11.74 -7.22 ...

打字

posthoc$`Expt1$Trt`

给出

     diff     lwr     upr     p adj
B-A   6.4562   2.130  10.782 9.530e-04
C-A  18.1922  13.866  22.519 2.454e-14
D-A  -0.7598  -5.086   3.566 9.680e-01
C-B  11.7360   7.410  16.062 7.427e-10
D-B  -7.2160 -11.542  -2.890 1.725e-04
D-C -18.9521 -23.278 -14.626 1.588e-14

因此您可以通过键入以下命令来访问第四列

posthoc$`Expt1$Trt`[,4]

或

posthoc$`Expt1$Trt`[,'p adj']

对于像Expt1R这样的一些对象，str()的输出有时可能是压倒性的。使用names()查看其中的不同对象也很有帮助。

names(Expt1R)

 [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values" "assign"       
 [7] "qr"            "df.residual"   "contrasts"     "xlevels"       "call"          "terms"        
[13] "model" 

所以

Expt1R$df.residual

会给你余量的自由度。