How to run Hutcheson t-test on R?
Maybe the following is a solution: # load package needed # QSutils is on Bioconductor library(QSutils) # here some exemplary data - these are the data from Pilou 1966 that are used # in the second example of Hutcheson, J theor Biol 129:151-154 (1970) earlier <- c(0,0,146,0,5,46,0,1,295,0,0,3,0,0,0,0,0) later <- c(0,0,142,0,5,46,0,1,246,0,0,3,0,0,0,0,0) # here starts the code ; you may replace the variables "earlier" and "later" # by your own numbers. # calculate h, var(h) etc h1 <- Shannon(earlier) varh1 <- ShannonVar(earlier) n1 <- sum (earlier) h2 <- Shannon(later) varh2 <- ShannonVar(later) n2 <- sum (later) degfree <- (varh1 + varh2)**2 /(varh1**2/n1 + varh2**2 /n2) # compare numbers with those in the paper h1 h2 varh1 varh2 Here var(h) is calculated as in ref 1 cited by Rui Barradas - I guess that explains the minor numerical differences obtained with the code above and the published variances. # this is the actual t-test t <- (h1-h2) /sqrt(varh1 + varh2) p <- 2*pt(-abs(t),df= degfree) p that's it Best Karl
On 08.09.2020 16:55, Rui Barradas wrote:
Hello,
No, it's not. That's the Shannon diversity index, the test the OP is
looking for is a t-test for Shannon diversity index equality. The index
itself is easy to code. A very simple example, based on ?vegan::diversity:
library(vegan)
data(BCI)
H <- diversity(BCI[1,])??? # just first row
divers <- function(n){
? p <- n/sum(n)
? log_p <- numeric(length(n))
? log_p[n != 0] <- log(p[n != 0])
? -sum(p * log_p)
}
HRui <- divers(BCI[1,])
identical(H, HRui)
#[1] TRUE
The vegan function is more general, it applies this and other indices
calculations to a matrix or array.
The t-test doesn't seem difficult to code.
The variance formula in the paper and in the OP's posted link [1] are
not the same, the original has one more term, but the degrees of freedom
formula are the same. It all seems straightforward coding.
Luigi: Maybe later today I will have time but I am not making promises.
[1] https://www.dataanalytics.org.uk/comparing-diversity/
Hope this helps,
Rui Barradas
?s 12:35 de 08/09/20, Karl Schilling escreveu:
Could it be that the test you are looking for is implemented in the vegan package (function diversity(... index = "shannon" ...), and/or the BiodiversityR package, function "diversityresult (..., index = "Shannon",...) best, Karl Schilling