library(conflicted)
library(tidyverse)
library(vegan)
Note
This post and its twin were conceived together. In my mind, the two topics are tightly linked, but a tangent about different software implementations of a diversity index grew so large that it deserved its own post. If you’re interested in both, you can read that one after this.
There is a tradition in ecology to represent community data in a “community matrix” (Legendre and Legendre 2012). This what we can call a “wide” format, with species on columns and samples on rows. I prefer this convention to having them the other way around, since we can treat species as “variables” and samples as “observations”.
The popular package vegan interprets this literally as matrix objects. Matrices in R have nice uses, for example if you want to exploit matrix algebra for the fastest calculations, but they are raw and ugly to print and view, definitely not what you would want for general data analysis. The package tidyvegan seems to share this sentiment, adapting the most useful functions in the package to the use of tibbles and the pipe-friendly “data as first argument” convention.
The thing is, community data is often awfully sparse. Wouldn’t it be nice to work in a long format? That has often been my preference, so I’ll show you how to reconcile the two worlds without too much pain. Spoiler: this mostly involves the use of pivot_wider() and pivot_longer(), but I think the practical examples and some added details make this worth sharing.
Setup
Just loading some packages.
Example 1: pairwise dissimilarities
The following example uses a small, cleaned-up subset of the data in Cervellini et al. (2026). Several observers have recorded species cover in the same sampling unit (this is just one sampling unit and one exercise):
training <- read_rds("tidyveg_ex1.rds")
training# A tibble: 175 × 3
observer species cover
<chr> <chr> <dbl>
1 a Ruscus aculeatus 10
2 a Hedera helix 2
3 a Melica uniflora 0.1
4 a Rubia peregrina 0.5
5 a Cyclamen repandum 1
6 a Quercus cerris 0.1
7 B Lonicera etrusca 0.1
8 B Cyclamen repandum 0.5
9 B Melica uniflora 0.5
10 B Hedera helix 1
# ℹ 165 more rowsThis long tibble is nice: you can easily calculate observer species richness, filter out rows, transform the percent cover; it’s not the format vegdist() wants to see though. You could just pivot the tibble and feed it to the function, this actually works because vegdist() silently converts it to a matrix:
plain1 <- training |>
pivot_wider(names_from = species, values_from = cover, values_fill = 0) |>
column_to_rownames("observer") |>
# default Bray-Curtis (Steinhaus, Motyka, Odum, percentage difference or however you want to call it)
vegdist()
plain1 a B c d E f g
B 0.37628866
c 0.50035386 0.23443019
d 0.59124088 0.34743202 0.20819602
E 0.35370685 0.16987591 0.21955096 0.34746573
f 0.59479222 0.36860827 0.26540951 0.11345394 0.36791126
g 0.30291698 0.51232302 0.53798227 0.58219431 0.34083045 0.63436203
h 0.42393205 0.24168773 0.21192294 0.32967522 0.12320200 0.41688501 0.42443485
i 0.43859649 0.11675741 0.16246709 0.29105047 0.15499675 0.37092432 0.47615241
l 0.49898692 0.31496423 0.17263279 0.27286082 0.30997911 0.34016393 0.55771024
m 0.63031423 0.44891717 0.39170183 0.23830963 0.46952856 0.14808081 0.67762308
n 0.44652206 0.18074160 0.07150635 0.23863131 0.20249740 0.24658331 0.52344740
o 0.31958128 0.22515953 0.25679376 0.38076152 0.08411677 0.41259199 0.34318039
p 0.40981938 0.12837705 0.24051565 0.33743067 0.16105820 0.36128273 0.47115136
q 0.32626326 0.23101703 0.29649309 0.40094181 0.13450000 0.43125309 0.36942675
h i l m n o p
B
c
d
E
f
g
h
i 0.17394726
l 0.26001195 0.29563614
m 0.49541948 0.46770217 0.32832750
n 0.28706864 0.19449902 0.21701493 0.36657033
o 0.16958546 0.18881442 0.37072596 0.52446625 0.26178358
p 0.20639199 0.15928144 0.28549849 0.41123541 0.18454441 0.21694075
q 0.21011412 0.22333722 0.39847328 0.52607014 0.26021890 0.05654281 0.21946782There are some issues with this. First of all, vegan really likes to use row and column names to identify the pairs. There’s not much choice with matrices, but maybe in the original data you have a separate identifier column. If you want that to be transferred, you need to transform it into row names. Tibbles hate row names so much that as soon as I call that second function the data stops being a tibble and becomes a plain dataframe. The documentation in ?tibble::rownames explains this stance convincingly: row names are just a character variable, but much more volatile (you can lose it from one operation to the next) and with a special treatment. If you want a character variable, make a column for it.
The second issue is with the output format. I really think the dist class is an ugly and error-prone convention started by base R. What you get from the above looks like a matrix (more specifically, the lower triangle), but is in fact a weird dressed-up vector. The same is true even if you set diag = TRUE and upper = TRUE in vegdist(), even if the printed output looks identical to a matrix.
plain2 <- training |>
pivot_wider(names_from = species, values_from = cover, values_fill = 0) |>
column_to_rownames("observer") |>
vegdist(diag = T, upper = T)
plain2 a B c d E f g
a 0.00000000 0.37628866 0.50035386 0.59124088 0.35370685 0.59479222 0.30291698
B 0.37628866 0.00000000 0.23443019 0.34743202 0.16987591 0.36860827 0.51232302
c 0.50035386 0.23443019 0.00000000 0.20819602 0.21955096 0.26540951 0.53798227
d 0.59124088 0.34743202 0.20819602 0.00000000 0.34746573 0.11345394 0.58219431
E 0.35370685 0.16987591 0.21955096 0.34746573 0.00000000 0.36791126 0.34083045
f 0.59479222 0.36860827 0.26540951 0.11345394 0.36791126 0.00000000 0.63436203
g 0.30291698 0.51232302 0.53798227 0.58219431 0.34083045 0.63436203 0.00000000
h 0.42393205 0.24168773 0.21192294 0.32967522 0.12320200 0.41688501 0.42443485
i 0.43859649 0.11675741 0.16246709 0.29105047 0.15499675 0.37092432 0.47615241
l 0.49898692 0.31496423 0.17263279 0.27286082 0.30997911 0.34016393 0.55771024
m 0.63031423 0.44891717 0.39170183 0.23830963 0.46952856 0.14808081 0.67762308
n 0.44652206 0.18074160 0.07150635 0.23863131 0.20249740 0.24658331 0.52344740
o 0.31958128 0.22515953 0.25679376 0.38076152 0.08411677 0.41259199 0.34318039
p 0.40981938 0.12837705 0.24051565 0.33743067 0.16105820 0.36128273 0.47115136
q 0.32626326 0.23101703 0.29649309 0.40094181 0.13450000 0.43125309 0.36942675
h i l m n o p
a 0.42393205 0.43859649 0.49898692 0.63031423 0.44652206 0.31958128 0.40981938
B 0.24168773 0.11675741 0.31496423 0.44891717 0.18074160 0.22515953 0.12837705
c 0.21192294 0.16246709 0.17263279 0.39170183 0.07150635 0.25679376 0.24051565
d 0.32967522 0.29105047 0.27286082 0.23830963 0.23863131 0.38076152 0.33743067
E 0.12320200 0.15499675 0.30997911 0.46952856 0.20249740 0.08411677 0.16105820
f 0.41688501 0.37092432 0.34016393 0.14808081 0.24658331 0.41259199 0.36128273
g 0.42443485 0.47615241 0.55771024 0.67762308 0.52344740 0.34318039 0.47115136
h 0.00000000 0.17394726 0.26001195 0.49541948 0.28706864 0.16958546 0.20639199
i 0.17394726 0.00000000 0.29563614 0.46770217 0.19449902 0.18881442 0.15928144
l 0.26001195 0.29563614 0.00000000 0.32832750 0.21701493 0.37072596 0.28549849
m 0.49541948 0.46770217 0.32832750 0.00000000 0.36657033 0.52446625 0.41123541
n 0.28706864 0.19449902 0.21701493 0.36657033 0.00000000 0.26178358 0.18454441
o 0.16958546 0.18881442 0.37072596 0.52446625 0.26178358 0.00000000 0.21694075
p 0.20639199 0.15928144 0.28549849 0.41123541 0.18454441 0.21694075 0.00000000
q 0.21011412 0.22333722 0.39847328 0.52607014 0.26021890 0.05654281 0.21946782
q
a 0.32626326
B 0.23101703
c 0.29649309
d 0.40094181
E 0.13450000
f 0.43125309
g 0.36942675
h 0.21011412
i 0.22333722
l 0.39847328
m 0.52607014
n 0.26021890
o 0.05654281
p 0.21946782
q 0.00000000This gets especially confusing when you try to do normal matrix operations like subsetting, since the underlying object is one-dimensional:
try(plain2[1, 1])Error in plain2[1, 1] : incorrect number of dimensionsI think it’s best to call as.matrix() right away and work with the full symmetrical matrix, then handle it as needed.
Finally, maybe I would like to let stats and vegan enjoy their matrices and look at the output in a tibble instead. Let’s say we want to use those distance values as data in a model, as we did in the paper above. To do all this behind the scenes, I made two helper functions, which I have now adapted for more general use:
vegan_cook <- function(x, species_col = "species", value_col = "cover", sample_name = NULL){
if(is.null(sample_name)) sample_name <- setdiff(names(x), c(species_col, value_col))
if(length(sample_name) > 1) stop("Please provide a unique identifier column for samples")
x |>
pivot_wider(names_from = all_of(species_col), values_from = all_of(value_col), values_fill = 0) |>
column_to_rownames(sample_name) |>
as.matrix()
}
tidyvegdist <- function(x, sample_names = "sample", diag = F, tri = F, ...){
if(length(sample_names) == 1) sample_names <- str_c(sample_names, 1:2, sep = "_")
vegdist(x, ...) |>
as.matrix() |>
(\(m) {if (tri) m[upper.tri(m, diag = F)] <- NA; m})() |>
(\(m) {if (!diag) diag(m) <- NA; m})() |>
as_tibble(rownames = sample_names[1]) |>
pivot_longer(!all_of(sample_names[1]), names_to = sample_names[2]) |>
drop_na()
}vegan_cook() handles the forward transformation into a vegan matrix. You can specify the names of the columns to pivot from with species_col and value_col. The row (sample) name columns is gathered automatically from the remaining column, or can be specified manually. tidyvegdist() calls vegdist() and returns a long tibble. You can choose how to call the columns with the pair names, either with an atomic string, which gets suffixes, or with a character vector. It defaults to "sample"; you can choose to include the diagonal values (usually 0) with diag, or to get a triangular matrix instead of the full symmetric matrix with tri; finally, you can pass other arguments to vegdist(), such as method.
training |>
vegan_cook() |>
tidyvegdist()# A tibble: 210 × 3
sample_1 sample_2 value
<chr> <chr> <dbl>
1 a B 0.376
2 a c 0.500
3 a d 0.591
4 a E 0.354
5 a f 0.595
6 a g 0.303
7 a h 0.424
8 a i 0.439
9 a l 0.499
10 a m 0.630
# ℹ 200 more rowsExample 2: transect data
The second example is a snippet from the data of De Benedictis et al. (2026). This is a peculiar case, because it comes from a 1000-unit looping transect which has been resampled with overlap at several scales, so there are many sampling units, but each has few records because they come from beech forest understory. This is just one scale of one transect. For each sampling unit X, row represent observations of a species. When the sampling units are aggregates as in this case, these can be used to calculate frequencies.
transect <- read_rds("tidyveg_ex2.rds")
transect# A tibble: 7,200 × 2
X species
<int> <chr>
1 152 Rosa arvensis
2 153 Rosa arvensis
3 154 Rosa arvensis
4 155 Rosa arvensis
5 156 Rosa arvensis
6 157 Rosa arvensis
7 157 Brachypodium sylvaticum
8 158 Rosa arvensis
9 158 Brachypodium sylvaticum
10 158 Brachypodium sylvaticum
# ℹ 7,190 more rowsThere are 7200 observations, and so rows. What would happen if we had this data in a wide format?
too_wide <- transect |>
mutate(presence = T) |>
pivot_wider(names_from = species, values_from = presence, values_fn = length)
glimpse(too_wide)Rows: 564
Columns: 11
$ X <int> 152, 153, 154, 155, 156, 157, 158, 159, 160,…
$ `Rosa arvensis` <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 4, 4, 4, 5,…
$ `Brachypodium sylvaticum` <int> NA, NA, NA, NA, NA, 1, 2, 3, 4, 5, 6, 7, 8, …
$ `Daphne laureola` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Festuca heterophylla` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Aremonia agrimonoides` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Clinopodium nepeta` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Dactylis glomerata` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Luzula forsteri` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Lactuca muralis` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ `Hieracium murorum` <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …I warned you the data were sparse! Out of 5640 values, 4039 are NA. And we’re lucky because only the 10 species recorded in this transect are listed.
And that’s not even fair to the long format, because this has already aggregated frequency counts in the rows. We can do the same with the long tibble:
transect <- transect |>
count(X, species)
transect# A tibble: 1,601 × 3
X species n
<int> <chr> <int>
1 152 Rosa arvensis 1
2 153 Rosa arvensis 1
3 154 Rosa arvensis 1
4 155 Rosa arvensis 1
5 156 Rosa arvensis 1
6 157 Brachypodium sylvaticum 1
7 157 Rosa arvensis 1
8 158 Brachypodium sylvaticum 2
9 158 Rosa arvensis 1
10 159 Brachypodium sylvaticum 3
# ℹ 1,591 more rowsAside from this cosmetic difference, what can we do with the long format? We can join some traits or other information for example, and keep it neatly along the records.
traits <- read_rds("traits.rds")
transect |>
left_join(traits)# A tibble: 1,601 × 13
X species n leaf_area_mm2 nmass_mg_g lma_g_m2 plant_height_m
<int> <chr> <int> <dbl> <dbl> <dbl> <dbl>
1 152 Rosa arvensis 1 302. 24.3 45.4 1.16
2 153 Rosa arvensis 1 302. 24.3 45.4 1.16
3 154 Rosa arvensis 1 302. 24.3 45.4 1.16
4 155 Rosa arvensis 1 302. 24.3 45.4 1.16
5 156 Rosa arvensis 1 302. 24.3 45.4 1.16
6 157 Brachypodium sy… 1 2195. 21.6 23.1 0.662
7 157 Rosa arvensis 1 302. 24.3 45.4 1.16
8 158 Brachypodium sy… 2 2195. 21.6 23.1 0.662
9 158 Rosa arvensis 1 302. 24.3 45.4 1.16
10 159 Brachypodium sy… 3 2195. 21.6 23.1 0.662
# ℹ 1,591 more rows
# ℹ 6 more variables: diaspore_mass_mg <dbl>, ssd_combined_mg_mm3 <dbl>,
# BBRsize <dbl>, persistence <dbl>, offspring <dbl>, spread <dbl>Or we can easily convert frequencies to proportions:
transect |>
mutate(prop = n / sum(n), .by = X)# A tibble: 1,601 × 4
X species n prop
<int> <chr> <int> <dbl>
1 152 Rosa arvensis 1 1
2 153 Rosa arvensis 1 1
3 154 Rosa arvensis 1 1
4 155 Rosa arvensis 1 1
5 156 Rosa arvensis 1 1
6 157 Brachypodium sylvaticum 1 0.5
7 157 Rosa arvensis 1 0.5
8 158 Brachypodium sylvaticum 2 0.667
9 158 Rosa arvensis 1 0.333
10 159 Brachypodium sylvaticum 3 0.75
# ℹ 1,591 more rowsHow would you do that with the wide format? Well, for the traits I’m not sure where I would fit them, but for the proportions you need a rowwise operation, which isn’t usually as easy:
too_wide <- too_wide |>
# you can avoid this step by using values_fill in pivot_wider, but I wanted to count the NAs
mutate(across(everything(), \(x) replace_na(x, 0)))
# rowwise is usually a sign that you're doing something wrong. Also, this is slow
too_wide |>
rowwise() |>
mutate(across(!X, \(x) x / sum(c_across(!X))))# A tibble: 564 × 11
X `Rosa arvensis` `Brachypodium sylvaticum` `Daphne laureola`
<int> <dbl> <dbl> <dbl>
1 152 1 0 0
2 153 1 0 0
3 154 1 0 0
4 155 1 0 0
5 156 1 0 0
6 157 0.5 0.5 0
7 158 0.333 0.667 0
8 159 0.25 0.75 0
9 160 0.333 0.667 0
10 161 0.375 0.625 0
# ℹ 554 more rows
# ℹ 7 more variables: `Festuca heterophylla` <dbl>,
# `Aremonia agrimonoides` <dbl>, `Clinopodium nepeta` <dbl>,
# `Dactylis glomerata` <dbl>, `Luzula forsteri` <dbl>,
# `Lactuca muralis` <dbl>, `Hieracium murorum` <dbl># since we're doing this, might as well give up and convert to matrix
too_wide |>
column_to_rownames("X") |>
as.matrix() |>
proportions(margin = 1) |>
head() Rosa arvensis Brachypodium sylvaticum Daphne laureola Festuca heterophylla
152 1.0 0.0 0 0
153 1.0 0.0 0 0
154 1.0 0.0 0 0
155 1.0 0.0 0 0
156 1.0 0.0 0 0
157 0.5 0.5 0 0
Aremonia agrimonoides Clinopodium nepeta Dactylis glomerata Luzula forsteri
152 0 0 0 0
153 0 0 0 0
154 0 0 0 0
155 0 0 0 0
156 0 0 0 0
157 0 0 0 0
Lactuca muralis Hieracium murorum
152 0 0
153 0 0
154 0 0
155 0 0
156 0 0
157 0 0Or we could calculate Rao’s quadratic entropy!
References
Cervellini, Marco, Luciano Ludovico Maria De Benedictis, Leonardo Salvatori, et al. 2026. “Training and Intercalibration Reduce Observer-Induced Variability in Forest Vegetation Surveys.” Vegetation Ecology and Diversity 63 (May): e189819. https://doi.org/10.3897/ved.189819.
De Benedictis, Luciano Ludovico Maria, Stefano Chelli, Zhengxue Zhu, et al. 2026. “Forest Structure and Understory Functional Diversity at Multiple Scales: The Importance of Median Tree Height.” Ecological Indicators 183 (February): 114687. https://doi.org/10.1016/j.ecolind.2026.114687.
Legendre, Pierre, and Louis Legendre. 2012. Numerical Ecology. Elsevier.