# What test should I use with shannon diversity index values?

What sort of test should I use to get the p-value after getting a shannon diversity index? Also does anyone know if and how they can be done using PRIMER 5?

-
What is your null hypothesis? – onestop Mar 25 '12 at 18:00
global hypothesis. – Ste Trevors Mar 25 '12 at 18:33

It depends on how the data were obtained, on the diversity model, and on the null hypothesis.

### Theory

The usual setting for computing this index is an observational study in which a number of "species" (or types of things, more generally) are observed and counted. For $n$ distinct possible species the counts can be considered a vector $(k_1,k_2,\ldots,k_n)$ of natural numbers. The Shannon Diversity Index, aka entropy, is computed as

$$H = -\sum_{i=1}^n \frac{k_i}{n} \log(\frac{k_i}{n}).$$

Whenever any of the species are not observed ($k_i=0$), they are considered as not contributing to $H$ (which makes mathematical sense, because the limiting value of $p\log(p)$ as $p$ goes down to $0$ is indeed $0$). Therefore it is not necessary to know $n$ for this computation.

A null hypothesis will stipulate the relative frequencies of species occurrence. For instance, it might suppose all possible species occur with equal frequency and have equal and independent chances of appearing in the study. The alternative is that the chances are unequal: this would necessarily reduce $H$, implying the test is one-tailed: a sufficiently small value of $H$ should be considered significant. In practice, a more interesting null hypothesis would instead be formulated according to some underlying theory of species distribution. (This particular one is uninteresting in most applications because we have little reason to suppose all species occur with equal frequencies, so demonstrating they do not is not terribly useful.) The null would also have to account for how the observations were made. (In ecology, for instance, species can be counted by quadrats, sighted along transects, reported by a large group of dispersed people during a definite interval of time, and in many other ways: each leads to a different expected distribution of observed counts.)

When all these assumptions are made, the null distribution of $H$ is complicated but easily obtained by simulation. It depends both on the number of possible species $n$ and on the dataset size, which I will write as a multiple of $n$, $k*n$. (Fortunately, as $n$ grows, the distribution of $H$ shifts to higher values, so that one can at least obtain an upper bound for the p-value by assuming $n$ equals the number of distinct observed species.) Once the null distribution of $H$ is obtained, the p-value is computed as usual by referencing the observed value of $H$ (for the data) to the cumulative probability in the null distribution. For example, an observed $H$ of $1.81$ in a dataset of $32$ observations from a population of $n=8$ or more distinct species would be referred to the middle distribution $(k=4,n=8)$ in the figure below and we would see that it corresponds to a low p-value, because few of the simulated values of $H$ are any less than $1.81$.

### Examples

To provide a sense of what these null distributions look like, I simulated 10,000 datasets for $n= 2, 8, 32$ of sizes $k*n$ with $k=1, 4, 16$ assuming very large parent populations of individuals (so that finite-population corrections would not enter into consideration). For instance, the case $k=4,n=8$ estimates the null distribution of $H$ for $8$ possible species in a sample of $32$ (so each is expected to occur $4$ times in the sample, but the actual results will vary by chance).

Across each row, the sample size grows larger for a fixed $n$. With large sample sizes the value of $H$ will be close to but typically slightly less than the null value, $\log(n)$, and it will be approximately normally distributed. The smaller sample sizes (left two columns) and smaller species richness (top two rows) exhibit noticeable negative skewness, multimodality, and discrete behavior: that is a warning against always using a normal approximation.

### Implementation

I have no idea what "PRIMER 5" is, so I conducted these simulations in R, which is well suited to the task. First, we need to calculate the Shannon Index, here named entropy:

entropy <- function(x) {       # x is a vector of (positive) species indexes
k <- tabulate(x)           # k is a vector of species counts
p <- k[k != 0] / length(x) # p gives the associated relative frequencies
-sum(p * log(p))           # Note that this is the natural logarithm
}


Then we can conduct a simulation for any given k and n; e.g.,

kn <- 32; n <- 8  # I.e., dataset size is kn and there are n possible species
set.seed(17)      # To reproduce these results
trials <- replicate(100000, entropy(sample(1:n, kn, replace=TRUE)))


It's wise to look at the full histogram of the simulated null distribution:

hist(trials)


The p-value for an observed index, such as 1.81, is then computed with a comparison:

length(trials[trials <= 1.81]) / length(trials)


The response is 0.02453; i.e, about 2.5%: a modestly significant (small) p-value.

-
Primer 5 is apparently software for analyzing and designing PCR primers. That part of the question may be better suited for BioStar SE. An extremely well thought out and developed post, @whuber, thanks for taking the time to write all of that up! – Chris Simokat Mar 26 '12 at 16:19