Extracting the "p adj" values from a TukeyHSD test

Question

With a dummy data set:

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

I have multiple Species and about 50 Variables (Var50). I would like to perform a One-way Anova on the paired grouping variable (Species) for each response variable and get the output of frequencies that are statistically significant at the 95% CI, for example. I began writing a function to do this as follows:

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

Ultimately, I would like the final table (dummy answers) to look like

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

I know I can use [[]] to access data within the results for the Tukey test. I have tried to use tukey.results[[1]][,1] up to tukey.results[[1]][,3] to no avail. tukey.results[[1]] returns all the columns from the Tukey test.

Also, I am thinking I might have to use cbind somewhere in the function to get the data in their respective columns. OR I was thinking it would be possible to use the apply command but I don't know how to keep the grouping variable constant while varying the response variable at every iteration.

Any suggestions would be deeply appreciated.

Answer 1

Try this if you are also looking for variables:

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

Answer 2

The answer did not work on my system. Here is my solution starting with the code Metrics provided.

 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