Computing diversity indices in R with vegan

Compute species richness, Shannon and Simpson diversity from a site-by-species community matrix in R with vegan, with reusable code and a small worked example.

Step 11 of the Tidy Ecology course, Part three: describe the community.

Updated 29 September 2026: a paragraph of checks for your own data was added: empty sampling units, single-species units, where Pielou’s J is undefined, and cover data, which the count-based rarefaction tools do not accept.

Updated 25 September 2026: the post said Shannon values usually fall between about 1.5 and 3.5, a range with no source and no general truth. A new section, Reading a Shannon value, converts the index to an effective number of species and to evenness on the post’s own plots, and says why values from unequal sampling effort do not compare.

Alpha diversity, meaning how diverse a single site or sample is, is usually the first number an ecologist wants out of a species-abundance table. This post covers the three you’ll reach for most often: species richness, the Shannon index, and the Simpson index, all from the vegan package.

We’ll go from a raw community matrix to a tidy results table and a finished plot.

A community matrix

vegan expects a community matrix: one row per site (or sample), one column per species, cells holding abundances (counts or cover). Here’s a small worked example, with counts of five plant species across four meadow plots:

library(vegan)
library(dplyr)
library(ggplot2)

paper <- "#f5f4ee"; ink <- "#16241d"; forest <- "#275139"
gold  <- "#cda23f"; abandoned <- "#b5534e"; faint <- "#5d6b61"; line <- "#dad9ca"

theme_te <- theme_minimal(base_size = 12) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = line, linewidth = 0.3),
        plot.background = element_rect(fill = paper, colour = NA),
        panel.background = element_rect(fill = paper, colour = NA),
        axis.title = element_text(colour = ink),
        axis.text = element_text(colour = faint),
        plot.title = element_text(colour = ink, face = "bold"))

community <- data.frame(
  row.names = c("Plot_A", "Plot_B", "Plot_C", "Plot_D"),
  Festuca   = c(10,  0, 25,  4),
  Trifolium = c( 8, 12,  2,  5),
  Plantago  = c( 5,  3,  1,  6),
  Achillea  = c( 2,  9,  0,  5),
  Lotus     = c( 0, 15,  0,  4)
)

community

       Festuca Trifolium Plantago Achillea Lotus
Plot_A      10         8        5        2     0
Plot_B       0        12        3        9    15
Plot_C      25         2        1        0     0
Plot_D       4         5        6        5     4

Each row is a plot; each column a species. Plot_C is dominated by Festuca, while Plot_D spreads its individuals evenly across all five species.

The three indices

vegan gives you each of these in a single call:

richness <- specnumber(community)                 # number of species present
shannon  <- diversity(community, index = "shannon")
simpson  <- diversity(community, index = "simpson")

richness

Plot_A Plot_B Plot_C Plot_D 
     4      4      3      5 

shannon

   Plot_A    Plot_B    Plot_C    Plot_D 
1.2550811 1.2658568 0.4086977 1.5974167 

simpson

   Plot_A    Plot_B    Plot_C    Plot_D 
0.6912000 0.6982249 0.1964286 0.7951389 

A few things worth knowing:

• specnumber() simply counts the non-zero species. It’s the most basic measure, and blind to abundance.
• Shannon (H’) rises with both richness and evenness. Its maximum is not fixed: with S species it can reach at most ln S, so the range it can take grows with richness (the next section shows how to read a value).
• Simpson is returned here as 1 - D: the probability that two individuals drawn at random belong to different species. It runs from 0 to 1 and is less sensitive to rare species than Shannon.

Reading a Shannon value

To read a Shannon value, convert it. exp(H’) is the number of equally common species that would give the same value, and H’ divided by ln S, the most H’ can reach with S species, says how close the community is to perfectly even. The raw number has no scale of its own: it measures the uncertainty about which species the next individual drawn will belong to, and its size depends on the base of the logarithm, on how many species there are and on how hard the site was sampled. The chunk below computes each reading for the four plots:

shannon_bits      <- diversity(community, index = "shannon", base = 2)
shannon_ceiling   <- log(richness)       # H' if every species were equally common
pielou_j          <- shannon / shannon_ceiling
effective_species <- exp(shannon)        # Hill number of order q = 1
plot_totals       <- rowSums(community)  # individuals counted per plot

shannon_reading <- data.frame(
  n       = plot_totals,
  S       = richness,
  H_nats  = round(shannon, 2),
  H_bits  = round(shannon_bits, 2),
  ln_S    = round(shannon_ceiling, 2),
  J       = round(pielou_j, 2),
  exp_H   = round(effective_species, 2)
)
shannon_reading

        n S H_nats H_bits ln_S    J exp_H
Plot_A 25 4   1.26   1.81 1.39 0.91  3.51
Plot_B 39 4   1.27   1.83 1.39 0.91  3.55
Plot_C 28 3   0.41   0.59 1.10 0.37  1.50
Plot_D 24 5   1.60   2.30 1.61 0.99  4.94

example_h <- c(2.1, 2.4)                 # two published-looking values

Units. diversity() uses the natural log, so its H’ is in nats. Software or papers that use log base 2 report bits, which is every value multiplied by 1/ln 2 = 1.443: Plot_D reads 1.60 in nats and 2.30 in bits. Base 10, which some papers use, multiplies by 1/ln 10 = 0.434. Check the base before you set your value beside a published one.

The ceiling. With S species, H’ reaches its maximum of ln S only when all of them are equally common. Plot_D sits at 1.60 against a ceiling of ln 5 = 1.61; Plot_C reaches 0.41 against ln 3 = 1.10. The ratio H’/ln S is Pielou’s evenness J, 0.99 and 0.37 here. Because the ceiling moves with richness there is no general “normal” range for H’: fifty equally common species give ln 50 = 3.91, five give 1.61, and both communities are perfectly even.

Effective number of species. exp(H’) turns the index into the number of equally common species that would give the same value, the Hill number of order 1 (Jost 2006). Like H’ itself, exp(H’) computed from a sample is biased low, and Estimating diversity with Hill numbers shows how to correct it. Plot_D is worth 4.94 equally common species out of 5 present, Plot_C only 1.50 out of 3. This is the scale for saying how much more diverse one site is than another. Values of 2.1 and 2.4 nats look close, but they are 8.2 and 11.0 effective species, 35 per cent more at the second site.

Effort. Even in the same units, 2.1 at one site and 2.4 at another is not yet a comparison. The H’ of a sample is biased low, because rare species are missed, and on average it climbs as more individuals are counted; if the second site was sampled harder, part or all of the gap can be effort. When not to use Shannon shows the climb on a single community, and Pielou’s evenness under unequal sampling effort shows J ranking two communities by effort rather than by evenness. Even these four plots hold between 24 and 39 individuals. Equalise effort first, by rarefying to a common count or to a common coverage (Coverage-based rarefaction and extrapolation), and compare the values after that.

A tidy results table

Let’s pull everything into one data frame, rounded for reading:

diversity_summary <- tibble(
  plot     = rownames(community),
  richness = richness,
  shannon  = round(shannon, 2),
  simpson  = round(simpson, 2)
) |>
  arrange(desc(shannon))

diversity_summary

# A tibble: 4 × 4
  plot   richness shannon simpson
  <chr>     <int>   <dbl>   <dbl>
1 Plot_D        5    1.6     0.8 
2 Plot_B        4    1.27    0.7 
3 Plot_A        4    1.26    0.69
4 Plot_C        3    0.41    0.2 

Plot_D comes out most diverse. It carries all five species at similar abundances, so its evenness is high. Plot_A and Plot_B land in the middle and almost tie (1.26 vs 1.27, from 25 and 39 individuals, so read the tie as provisional), despite very different species makeups. Plot_C, dominated by a single grass, scores lowest even though three species are present. That gap between richness and Shannon is the evenness effect at work.

A plot

Finally, the bit you’ll actually paste into a report: Shannon diversity per plot.

ggplot(diversity_summary, aes(x = reorder(plot, shannon), y = shannon)) +
  geom_col(fill = "#2f6b3e", width = 0.7) +
  geom_text(aes(label = shannon), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(x = NULL, y = "Shannon diversity (H')") +
  ylim(0, max(diversity_summary$shannon) * 1.15) +
  theme_te

Figure 1: Shannon diversity (H’) across the four meadow plots.

Where to go next

This is the foundation. Realistic next steps:

• swap the toy matrix for your own data, read in with read.csv() into the same site-by-species shape;
• compare Pielou’s evenness (shannon / log(richness), the J column above) only at equal effort, and at equal coverage as well when richness differs: it moves with sampling effort, as Pielou’s evenness under unequal sampling effort shows;
• once you have many sites, move from alpha to beta diversity and ordination (vegdist(), metaMDS()).

Three checks belong between read.csv() and diversity() when the data are your own. First, find empty sampling units with rowSums(community) == 0: diversity() still returns a number for a row of zeros (a Shannon index of 0), and that number looks like data but describes nothing, so drop or report such units separately. Second, a unit with a single species has H’ = 0 and ln S = 0, so its Pielou’s J is 0/0 and comes back as NaN; evenness is undefined there, not zero. Third, cover values are fine for the indices themselves, which only use proportions, but not for the count-based tools: rarefy() accepts only integer counts, and rarefying to a common number of individuals has no meaning for percentage cover, so the effort advice above applies to count data.

Ordination is the next post. If you hit a snag adapting this to your own data, let me know.

References

Jost L 2006 Oikos 113(2):363-375 (10.1111/j.2006.0030-1299.14714.x)

Shannon CE 1948 Bell System Technical Journal 27(3):379-423 (10.1002/j.1538-7305.1948.tb01338.x)

Simpson EH 1949 Nature 163(4148):688 (10.1038/163688a0)

Related tutorials

• When not to use Shannon
• Rarefaction and accumulation curves
• Coverage-based rarefaction and extrapolation
• Pielou’s evenness under unequal sampling effort

Reproducibility

• Code last run: 2026-10-15
• R version 4.5.3 (2026-03-11), aarch64-apple-darwin20
• Packages: dplyr 1.2.1, ggplot2 4.0.3, vegan 2.7-6

PreviousStep 10, Standard errors and confidence intervals CourseAll 44 lessons NextStep 12, When not to use Shannon