[R] How to run Hutcheson t-test on R?

[Karl Schilling]

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