Some code to run Hutcheson t-test using R
Thank you very much for the code, that was very helpful. I got the article by Hutcheson -- I don't know if I can distribute it , given the possible copyrights, or if I can attach it here -- but it does not report numbers directly: it refers to a previous article counting bird death on a telegraph each year. The numbers are: bird_1956 <- c(4,4,190,135,56,3,2,2,1,12,41,201,1,0,131,3) bird_1957 <- c(4,111,53,66,146,222,126,61,0,2323,21) bird_1958 <- c(0,3,32,228,56,102,0,11,2,220,0) bird_1959 <- c(0,0,14,59,26,68,0) bird_1960 <- c(0,0,73,66,71,25,0,109,63,1) This for sake of the argument. As for my problem, I implemented the Shannon index with the package iNext, which only gives me the index itself and the 95% CI. Even when I implemented it with vegan, I only got the index. Essentially I don't have a count of species I could feed into the Hutcheson's. Is there a way to extract these data? Or to run a Hutcheson's on the final index? Thank you On Tue, Sep 8, 2020 at 7:43 PM Karl Schilling
<karl.schilling at uni-bonn.de> wrote:
Dear Luigi, below some code I cobbled together based on the Hutcheson paper you mentioned. I was lucky to find code to calculate h and, importantly, its variance in the R-package QSutils - you may find it on the Bioconductor website. here is the code, along with an example. I also attach the code as an R-file. Hope that helps. All my best Karl PS don't forget to adjust for multiple testing if you compare more than two groups. K # 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) # numbers of the first example used by Hutcheson were unfortunately not # available to me # 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 # I assume that minor numerical differences are due to differences in the # numerical precision of computers in the early seventies and today / KS # 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 -- Karl Schilling, MD Professor of Anatomy and Cell Biology Anatomisches Institut Rheinische Friedrich-Wilhelms-Universit?t Nussallee 10 D-53115 Bonn Germany phone ++49-228-73-2602
Best regards, Luigi