Distance-based RDA

Effect of ploughing time on weeds

Post (1986) carried out an experiment to investigate the effect of time of ploughing on the subsequent weed composition in summer barley fields. Three different ploughing times were used. There are 13 weed species observed in 12 plots in a randomized block experiment of four complete blocks of three plots each.

Traditionally, a multivariate ANOVA could have been used to analyses data like this. However, as there are more species (response variables; 13) than there are samples (12), MANOVA cannot be used. Instead, RDA is a useful option.

library("vegan")

## Loading required package: permute

library("dplyr")

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library("readr")

plough <- read_csv(url("https://bit.ly/plough-spp"), skip = 1) |>
    rename(label = "...1") |>
    select(-`...2`) |>
    janitor::clean_names() |>
    filter(!is.na(label)) |>
    tibble::column_to_rownames("label")

## New names:
## • `` -> `...1`
## • `` -> `...2`

## Rows: 12 Columns: 15
## ── Column specification ──────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): ...1, ...2
## dbl (13): Che alb, Pol per, Pol avi, Spe arv, Mat rec, Gna uli, Cap bur, Ste...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

plough_design <- read_csv(url("https://bit.ly/plough-design"), skip = 1) |>
    rename(label = "...1", ploughing_time = "Treatmnt") |>
    select(-`...2`) |>
    janitor::clean_names() |>
    filter(!is.na(label)) |>
    tibble::column_to_rownames("label")

## New names:
## Rows: 12 Columns: 4
## ── Column specification
## ────────────────────────────────────────────────────────────────────── Delimiter: "," chr
## (4): ...1, ...2, Block, Treatmnt
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ Specify the column
## types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
## • `` -> `...2`

We are interested in the effect of ploughing time and wish to remove any possible effect of the blocks. We achieve this by using block as a covariate (in Condition()). Some zero abundances are recorded as NA values, so we must replace those with 0s before we proceed. Instead of log-transforming the species counts, we might use a non-Euclidean dissimilarity coefficient with dbRDA instead of RDA.

plough <- plough |>
  mutate(across(che_alb:ape_spi,
    .fns = \(x) tidyr::replace_na(x, replace = 0)))

A correct permutation test for these data should keep samples contained at the block level and only permute within, not between, blocks.

h <- with(plough_design,
    how(blocks = block, nperm = 999))

The model and tests are specified as

ord <- dbrda(plough ~ ploughing_time + Condition(block), data = plough_design,
  distance = "bray", add = "lingoes")
ord

## 
## Call: dbrda(formula = plough ~ ploughing_time + Condition(block), data =
## plough_design, distance = "bray", add = "lingoes")
## 
##               Inertia Proportion Rank
## Total          1.2211     1.0000     
## Conditional    0.6843     0.5604    3
## Constrained    0.2934     0.2403    2
## Unconstrained  0.2434     0.1993    6
## 
## Inertia is Lingoes adjusted squared Bray distance
## 
## Eigenvalues for constrained axes:
##  dbRDA1  dbRDA2 
## 0.22484 0.06857 
## 
## Eigenvalues for unconstrained axes:
##    MDS1    MDS2    MDS3    MDS4    MDS5    MDS6 
## 0.09176 0.06705 0.04044 0.03076 0.01020 0.00320 
## 
## Constant added to distances: 0.01433247

anova(ord, permutations = h)

## Set of permutations < 'minperm'. Generating entire set.

## Permutation test for dbrda under reduced model
## Blocks:  block 
## Permutation: free
## Number of permutations: 1295
## 
## Model: dbrda(formula = plough ~ ploughing_time + Condition(block), data = plough_design, distance = "bray", add = "lingoes")
##          Df SumOfSqs      F Pr(>F)   
## Model     2  0.29341 3.6163  0.004 **
## Residual  6  0.24341                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notice how we get a warning/message from permute

Set of permutations < ‘minperm’. Generating entire set.

This means that there are not a large number of permutations of these data owing to the block restriction.

plot(ord)