Message-ID: <b3b12e54-0ce8-f9c4-e665-1555d8577834@uni-bonn.de>
Date: 2020-09-08T18:01:55Z
From: Karl Schilling
Subject: How to run Hutcheson t-test on R?
In-Reply-To: <88ef05c2-dda5-7e9e-e8f0-cb7f168cb7ad@sapo.pt>
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