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
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
