library("vegan")
## Loading required package: permute
library("readr")
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("forcats")

spp <- read_csv(url("https://bit.ly/ohraz-spp")) %>%
    rename(label = "...1") %>%
    janitor::clean_names()
## New names:
## • `` -> `...1`
## Rows: 96 Columns: 87
## ── Column specification ──────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): ...1
## dbl (86): molicaer, nardstri, festrubr, brizmed, anthodor, siegdecu, holclan...
## 
## ℹ 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.
molinia <- spp %>%
    select(label:molicaer)

spp <- spp %>%
    select(-molicaer) %>%
    tibble::column_to_rownames("label")

env <- read_csv(url("https://bit.ly/ohraz-env")) %>%
    rename(label = "...1") %>%
    mutate(across(c(mowing:removal, plotid), ~ factor(.x))) %>%
    tibble::column_to_rownames("label")
## New names:
## Rows: 96 Columns: 6
## ── Column specification
## ────────────────────────────────────────────────────────────────────── Delimiter: "," chr
## (5): ...1, mowing, fertilizer, removal, plotid dbl (1): year
## ℹ 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`

We need to recode the variables to make this easier to understand for PRC

env2 <- env %>%
mutate(mowed = fct_recode(mowing, mowed = "Yes", unmowed = "No"),
removed = fct_recode(removal, control = "No", removed = "Yes"),
fertilized = fct_recode(fertilizer, fertilized = "Yes", unfertilized = "No"),
removed = fct_relevel(removed, "control"),
treatment = fct_cross(removed, mowed, fertilized),
fyear = factor(year))

And we should log transform the species data as this is an RDA:

log_spp <- log1p(spp)

The PRC can then be fitted using

ohraz_prc <- prc(log_spp, env2$treatment, env2$fyear)

The permutation deisgn to use is

h <- with(env2, how(plots = Plots(strata = plotid), nperm = 999))

And with this we can now test the first PRC (RDA) axis

anova(ohraz_prc, permutations = h, first = TRUE)
## Permutation test for rda under reduced model
## Plots: plotid, plot permutation: none
## Permutation: free
## Number of permutations: 999
## 
## Model: prc(response = log_spp, treatment = env2$treatment, time = env2$fyear)
##          Df Variance      F Pr(>F)    
## RDA1      1   2.8841 11.151  0.001 ***
## Residual 64  16.5534                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Are the time structured effects of the treatment significant?

We might also restrict the permutations within plots for time series

h_time <- with(env2, how(plots = Plots(strata = plotid),
    within = Within(type = "series"), nperm = 999))

This gives a stricter test, especially if there is residual autocorrelation

anova(ohraz_prc, permutations = h_time, first = TRUE)
## Permutation test for rda under reduced model
## Plots: plotid, plot permutation: none
## Permutation: series
## Number of permutations: 999
## 
## Model: prc(response = log_spp, treatment = env2$treatment, time = env2$fyear)
##          Df Variance      F Pr(>F)    
## RDA1      1   2.8841 11.151  0.001 ***
## Residual 64  16.5534                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Does the effect remain after performing the stricter test?

Plot the PRC, without the species scores

plot(ohraz_prc, species = FALSE, legpos = "bottomleft", col = 1:8, lty = 1, lwd = 2)

Which treatment levels/combinations lead to the largest changes in the species composition by the end of the experiment?

We can plot the species scores on a separate plot, retaining only those species that have a fit greater than 5%

fit <- goodness(ohraz_prc, display = "species", choices = 1, scaling = "symmetric")
scrs <- scores(ohraz_prc, display = "species", choices = 1, scaling = "symmetric")
linestack(scrs[fit > 0.05, , drop = FALSE])
axis(side = 2)

Suggest which species were affected most in the treatments that affected the composition the most.