Skip to content

How to run Hutcheson t-test on R?

8 messages · Luigi Marongiu, Rolf Turner, Bernard McGarvey +3 more

Rolf Turner
# 
On Mon, 7 Sep 2020 11:17:36 +0200
Luigi Marongiu <marongiu.luigi at gmail.com> wrote:

            

      
      
      
    

        
Almost surely.  With R, all things are possible. :-)

            

      
      
      
    

        
Program it up?

cheers,

Rolf Turner

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Bernard McGarvey
    

    

      
      #
    

  


  

      
This website has an example calculation shown in Excel Which might help in programming it in R.

https://www.dataanalytics.org.uk/comparing-diversity/


Bernard
Sent from my iPhone so please excuse the spelling!"

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      David Winsemius
    

    

      
      #
    

  


  

      
On 9/7/20 3:17 PM, Rolf Turner wrote:

            

      
      
      
    

        
To Luigi;


Citing a 50 year-old paper that sits behind a paywall seems a bit 
ineffective in getting coding support.

Seems this might be a more appropriate question on the R-SIG-ecology or 
R-SIG-phylo mailing lists. (It would also have been appropriate to 
indicate what sort of searching has been done. My efforts at searching 
led me to the vegan package and this tutorial: 
https://cran.r-project.org/web/packages/vegan/vignettes/diversity-vegan.pdf 
. It doesn't appear to have a Hutcheson t-test, but I'm guessing that is 
because there are more modern and more sophisticated tests currently in 
use.)


See: https://www.r-project.org/mail.html

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Luigi Marongiu
    

    

      
      #
    

  


  

      
I cited that paper to show what test I was referring to, I was hoping it
was already coded...
Thanks, I will look into vegan
Best regards

        
On Tue, 8 Sep 2020, 01:23 David Winsemius, <dwinsemius at comcast.net> wrote:

        

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Karl Schilling
    

    

      
      #
    

  


  

      
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

  
  


  


      

    

    

      
      

        


  

        

      Rui Barradas
    

    

      
      #
    

  


  

      
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:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      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: