q-value and t-statistic from betadisper and envfit

[baileyjc]

Hello, I am still growing my statistical knowledge so I apologize in advance and please bear with me. I am working with some ecological data and to test the homogeneity of my sites I used betadisper and ran an anova based on that result as well as a post-hoc TukeyHSD test. I also used envfit to determine environmental correlations with my NMDS results. The question I have is that with some journals requiring test statistics now for all R-squared and p-values and for the sake of reproducibility, is there a way to calculate the q-value that is used to determine significance for the TukeyHSD test and the t-statistic that is used to determine whether a correlation R-squared value is significant for envfit? It seems that vegan would be the package to go for this but I have not found a solution but I am also new to the package so I apologize if I have overlooked a function. I created my own solution for betadisper and the TukeyHSD using the following code below but I have not discovered how to do this with the envfit results mainly because I have unbalanced sample sizes of site types. Any thoughts on how to compute a t-statistic from envfit results? I know I can use the equation sr <- sqrt((1-r)/(n-2)) to determine the standard error of r and then use t = r/sr for the t-statistic but again I have unbalanced sample sizes so maybe I just need to figure that out on my own. However, if you can provide any insight for this, please let me know and thank you for your time.

I surveyed 22 sites and there are three site type that I used as categories (O, M, and S). The number of O sites is 6, and the number of M and S sites is 8 each. Also I don't trust the output from qtukey or ptukey for determining the q-value or p-value from TukeyHSD because I used it to check results and the p-values didn't match which may be because they are adjusted in TukeyHSd but I am not sure.

pd_bd <- betadisper(pd_dist, category)
anova(pd_bd)
(pd_bd.HSD <- TukeyHSD(pd_bd))

Sites

N <- 22

Categories

k <- 3

Sum of squares

SS <- 0.137271 # from anova, the residual SS

Mean squares

MS <- SS/(N-k)

Balanced

BSE <- sqrt(MS / 8)

Unbalanced

USE <- sqrt((MS/2)*((1/6)+(1/8)))

Get differences

group <- pd_bd.HSD$group

q-values

(OM_q <- (group[1,1])/USE)
(MS_q <- (abs(group[2,1]))/BSE)
(SO_q <- (abs(group[3,1]))/USE)