Below, you’ll find several examples of experimental designs or other multivariate data sets that we wish to test using constrained ordination methods and restricted permutation (Monte Carlo) tests.
In this example we will study the effect of sod-cutting and sod-addition on the number of sporocarps of ectomycorrhizal fungi in Pine stands differing in age and soil type using data from an experiment reported in Baar & Ter Braak (1996). The experimental layout is as follows. Six stands (St1–St6) with Scots pine of different age and soil type were selected. Within each stand, plots were laid out, with a size of 15 m x 15 m. Treatments applied to stands were:
S — sod cutting, in which the litter and humus
layers and the herbaceous vegetation were removed,
A — sod addition, in which removed litter and humus
were added onto the existing litter and humus layers, and
C — control, in which the litter and humus layers
were left unchanged.
Treatment A (sod addition) was not done in the oldest
two stands (St5 and St6) as these were the
oldest stands and already had thick layers of humus and litter. The
treatments were replicated 4 times each, hence in the younger stands we
have 12 plots (4 replicates by 3 treatments), while in the oldest stands
we have 8 plots each (4 replicates by 2 treatments). Within stands,
treatments were randomized.
library("vegan")## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-7library("permute")
library("readr")
library("tidyr")
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, unionlibrary("tibble")
library("forcats")
library("stringr")
fungi <- read_csv(url("https://bit.ly/fungi-spp"), skip = 1L,
col_types = "-cddddddddddddddddddddddddddddddddddd") %>%
mutate(across(Amangemm:Xerosubt, ~ replace_na(.x, 0))) %>%
select(-(Ncarp:Nspec)) %>%
rename(label = "...2")## New names:
## • `` -> `...1`
## • `` -> `...2`design <- read_csv(url("https://bit.ly/fungi-design"), skip = 1L,
col_types = "-cccdc") %>%
rename(label = "...2", Treatment = `LH-treat`) %>%
mutate(Treatment = fct_relevel(Treatment, "C", after = 0L),
Treatment = fct_recode(Treatment, control = "C",
sod_cutting = "S", sod_addition = "A"),
Stand = factor(Stand),
SoilType = factor(SoilType))## New names:
## • `` -> `...1`
## • `` -> `...2`Using restricted permutations for split-plot designs, we can test the effects of variables that vary between whole-plots (but which are constant at the split-lot level). However, this is only possible if we make the design balanced by
St5 and
St6, orsod_addition plots.For this particular hypothesis we are testing the effect of a
variable that variables at the whole-plot level (the stand level), it is
better to keep as many stands as possible, so the
sod_addition plots will be deleted.
fungi_bal <- fungi %>%
left_join(design) %>%
filter(Treatment != "sod_addition") %>%
select(label:"Xerosubt") %>%
mutate(across(-label, ~ log1p(.x))) %>%
column_to_rownames("label")## Joining with `by = join_by(label)`design_bal <- design %>%
filter(Treatment != "sod_addition") %>%
column_to_rownames("label")m1 <- rda(fungi_bal ~ Age + SoilType, data = design_bal)
m1## Call: rda(formula = fungi_bal ~ Age + SoilType, data = design_bal)
##
## Inertia Proportion Rank
## Total 27.4042 1.0000
## Constrained 7.8451 0.2863 2
## Unconstrained 19.5591 0.7137 31
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2
## 6.393 1.453
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 7.208 4.524 2.071 1.490 0.908 0.820 0.747 0.479
## (Showing 8 of 31 unconstrained eigenvalues)RsquareAdj(m1)## $r.squared
## [1] 0.2862729
##
## $adj.r.squared
## [1] 0.2545517how() to test the effects of the whole plot variables that
only vary at the stand level (whole-plots).whole_h <- with(design_bal,
how(plots = Plots(Stand, type = "free"),
within = Within(type = "none", constant = FALSE),
nperm = 999))set.seed(42)
anova(m1, permutations = whole_h)## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## Permutation test for rda under reduced model
## Plots: Stand, plot permutation: free
## Permutation: none
## Number of permutations: 719
##
## Model: rda(formula = fungi_bal ~ Age + SoilType, data = design_bal)
## Df Variance F Pr(>F)
## Model 2 7.8451 9.0247 0.1097
## Residual 45 19.5591We can test the effect of Age while controlling (or
adjusting) for the effect of SoilType by partialling out
the effect of SoilType using Condition().
Age while
adjusting for the effect of SoilTypem2 <- rda(fungi_bal ~ Age + Condition(SoilType), data = design_bal)
m2## Call: rda(formula = fungi_bal ~ Age + Condition(SoilType), data =
## design_bal)
##
## Inertia Proportion Rank
## Total 27.40423 1.00000
## Conditional 1.60791 0.05867 1
## Constrained 6.23718 0.22760 1
## Unconstrained 19.55914 0.71373 31
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1
## 6.237
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 7.208 4.524 2.071 1.490 0.908 0.820 0.747 0.479
## (Showing 8 of 31 unconstrained eigenvalues)Age adjusted for SoilTypeset.seed(42)
anova(m2, permutations = whole_h)## 'nperm' >= set of all permutations: complete enumeration.## Set of permutations < 'minperm'. Generating entire set.## Permutation test for rda under reduced model
## Plots: Stand, plot permutation: free
## Permutation: none
## Number of permutations: 719
##
## Model: rda(formula = fungi_bal ~ Age + Condition(SoilType), data = design_bal)
## Df Variance F Pr(>F)
## Model 1 6.2372 14.35 0.07639 .
## Residual 45 19.5591
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The treatment varies at the split-plot (within) level. We now want to
test the effects of the treatment on the abundance of the fungi. This is
best done by randomising samples within blocks (stands), while keeping
the blocks fixed. To do this correctly, we also need to include the
block variable as a covariate to be partialled out using
Condition().
We can use all the data for this test, so we prepare the full data set for analysis
fungi_all <- fungi %>%
mutate(across(-label, ~ log1p(.x))) %>%
column_to_rownames("label")
design_all <- design %>%
column_to_rownames("label")m3 <- rda(fungi_all ~ Treatment + Condition(Stand), data = design_all)
m3## Call: rda(formula = fungi_all ~ Treatment + Condition(Stand), data =
## design_all)
##
## Inertia Proportion Rank
## Total 25.1883 1.0000
## Conditional 11.1065 0.4409 5
## Constrained 5.0201 0.1993 2
## Unconstrained 9.0617 0.3598 31
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2
## 4.871 0.149
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 1.8502 1.5930 1.3938 0.9987 0.7120 0.5380 0.5085 0.3620
## (Showing 8 of 31 unconstrained eigenvalues)split_h <- with(design_all,
how(blocks = Stand,
within = Within(type = "free"),
nperm = 999))set.seed(42)
anova(m3, permutations = split_h)## Permutation test for rda under reduced model
## Blocks: Stand
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = fungi_all ~ Treatment + Condition(Stand), data = design_all)
## Df Variance F Pr(>F)
## Model 2 5.0201 15.512 0.001 ***
## Residual 56 9.0617
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Additionally, we might be interested in testing the effect of the treatment differs among stands of different age, adjusting for a possible treatment by soil type interaction.
m4 <- rda(fungi_all ~ Treatment:Age +
Condition(Stand + Treatment + Treatment * SoilType),
data = design_all)
m4## Call: rda(formula = fungi_all ~ Treatment:Age + Condition(Stand +
## Treatment + Treatment * SoilType), data = design_all)
##
## Inertia Proportion Rank
## Total 25.18831 1.00000
## Conditional 16.71618 0.66365 9
## Constrained 0.96202 0.03819 2
## Unconstrained 7.51011 0.29816 31
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1 RDA2
## 0.8105 0.1515
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 1.5265 1.3627 0.9535 0.8659 0.5603 0.5044 0.3872 0.3433
## (Showing 8 of 31 unconstrained eigenvalues)anova(m4, permutations = split_h)## Permutation test for rda under reduced model
## Blocks: Stand
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = fungi_all ~ Treatment:Age + Condition(Stand + Treatment + Treatment * SoilType), data = design_all)
## Df Variance F Pr(>F)
## Model 2 0.9620 3.3305 0.001 ***
## Residual 52 7.5101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m4 <- rda(fungi_all ~ Age + Treatment + Treatment:Age +
Condition(Stand + SoilType + Treatment:SoilType),
data = design_all)
anova(m4, permutations = split_h, by = "margin")## Permutation test for rda under reduced model
## Marginal effects of terms
## Blocks: Stand
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = fungi_all ~ Age + Treatment + Treatment:Age + Condition(Stand + SoilType + Treatment:SoilType), data = design_all)
## Df Variance F Pr(>F)
## Age:Treatment 2 0.9620 3.3305 0.001 ***
## Residual 52 7.5101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(42)
anova(m4, permutations = split_h)## Permutation test for rda under reduced model
## Blocks: Stand
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = fungi_all ~ Age + Treatment + Treatment:Age + Condition(Stand + SoilType + Treatment:SoilType), data = design_all)
## Df Variance F Pr(>F)
## Model 2 0.9620 3.3305 0.001 ***
## Residual 52 7.5101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In this example you’ll look at the effect of pH on plant species in a wetland. The samples were laid out in grid of 20 samples arranged in five rows and four columns over the wetland.
The response data are abundances of 100 plant species at each of the 20 plots. We also have some associated data to go with the species abundances. The plots were sampled in 1977 and 1988, and pH and water depth were recorded at each location in the grid for each year.
For this example we’ll look only at the effect of pH on wetland plant abundance in 1988, so we will remove the older samples from the data before proceeding. We also need to ensure that the data are arranged in the data set in grid-column order. The code below loads the data, filters to retain the 1988 samples only, and then orders the data in grid-column order.
grid_veg <- read_csv(url("https://bit.ly/grid-veg"),
skip = 1L) %>%
mutate(across(ACHILMIL:`POPUL-SP`, ~ replace_na(.x, 0))) %>%
select(-("...1")) %>%
rename(label = "...2")## New names:
## Rows: 40 Columns: 102
## ── Column specification
## ──────────────────────────────────────────────────────────────── Delimiter: "," chr
## (2): ...1, ...2 dbl (100): ACHILMIL, ACHILPTA, AGROS-SP, ANGELSYL, ANTHOODO, CALAMCAN,
## CAPSE...
## ℹ 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`grid <- read_csv(url("https://bit.ly/grid-predictors"), skip = 1L,
col_types = "-cccdddd") %>%
rename(label = "...2") %>%
mutate(Column = factor(str_sub(Plot, 2L, 2L)),
Plot = factor(Plot),
Row = factor(Row)) %>%
relocate(Column, .after = Row) %>%
arrange(Year, Column, Row)## New names:
## • `` -> `...1`
## • `` -> `...2`# grid 1988 data
g88_veg <- grid %>%
left_join(grid_veg) %>%
filter(Year == 1988) %>%
select(c(label, ACHILMIL:`POPUL-SP`)) %>%
column_to_rownames("label")## Joining with `by = join_by(label)`g88 <- grid %>%
filter(Year == 1988) %>%
column_to_rownames("label")
g88_h <- how(within = Within(type = "grid", nrow = 5, ncol = 4,
mirror = TRUE))
# test
perm <- shuffle(nrow(g88), control = g88_h)
tmp <- matrix(0, nrow = 5, ncol = 4)
tmp[] <- as.character(g88$Plot)
tmp## [,1] [,2] [,3] [,4]
## [1,] "A1" "A2" "A3" "A4"
## [2,] "B1" "B2" "B3" "B4"
## [3,] "C1" "C2" "C3" "C4"
## [4,] "D1" "D2" "D3" "D4"
## [5,] "E1" "E2" "E3" "E4"tmp[] <- as.character(g88$Plot)[perm]
tmp # works!!## [,1] [,2] [,3] [,4]
## [1,] "D2" "D1" "D4" "D3"
## [2,] "C2" "C1" "C4" "C3"
## [3,] "B2" "B1" "B4" "B3"
## [4,] "A2" "A1" "A4" "A3"
## [5,] "E2" "E1" "E4" "E3"g1 <- rda(g88_veg ~ pH88, data = g88)
g1## Call: rda(formula = g88_veg ~ pH88, data = g88)
##
## Inertia Proportion Rank
## Total 23.4000 1.0000
## Constrained 5.8478 0.2499 1
## Unconstrained 17.5522 0.7501 18
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1
## 5.848
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 3.851 2.287 1.737 1.435 1.192 1.051 1.013 0.908
## (Showing 8 of 18 unconstrained eigenvalues)set.seed(88)
anova(g1, permutations = g88_h)## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.## Permutation test for rda under reduced model
## Permutation: grid mirrored 5 rows 4 columns
## Number of permutations: 79
##
## Model: rda(formula = g88_veg ~ pH88, data = g88)
## Df Variance F Pr(>F)
## Model 1 5.8478 5.997 0.0125 *
## Residual 18 17.5522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# or
set.seed(42)
perms <- shuffleSet(nrow(g88), 499, control = g88_h, check = FALSE)
anova(g1, permutations = perms)## Permutation test for rda under reduced model
## Permutation: grid mirrored 5 rows 4 columns
## Number of permutations: 499
##
## Model: rda(formula = g88_veg ~ pH88, data = g88)
## Df Variance F Pr(>F)
## Model 1 5.8478 5.997 0.014 *
## Residual 18 17.5522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This leaves us with 20 samples arranged in a grid.
We want to look at how the species abundances vary with
pH, while controlling for potential spatial autocorrelation
across the wetland. Fit an appropriate constrained ordination to
estimate the effect of pH on species abundances, and test
the estimated pH effect using an appropriate restricted
permutation test.
The code below loads the Ohraz data and the experimental design data using slightly more modern code than was used in the slides. Repeated observations of composition from an experiment
Results in 2 x 2 x 2 = 8 combinations in three replicates leading to 24 2m x 2m plots.
The data a repeated observations that include the baseline (before treatment), the interaction of treatment and time is of particular interest.
Test 1 of the hypotheses
There are no directional changes in species composition in time that are common to all treatments or specific treatments
library("vegan")
library("readr")
library("dplyr")
spp <- read_csv(url("https://bit.ly/ohraz-spp")) %>%
rename(label = "...1") %>%
janitor::clean_names()## New names:
## 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.
## • `` -> `...1`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`When we test for the interactions mentioned above, we’ll use
plotid as a covariate (conditioning variable) when setting
up the model.
We fit an RDA
a1 <- rda(spp ~ year + year:mowing + year:fertilizer + year:removal + Condition(plotid),
data = env)
a1## Call: rda(formula = spp ~ year + year:mowing + year:fertilizer +
## year:removal + Condition(plotid), data = env)
##
## Inertia Proportion Rank
## Total 748.2536 1.0000
## Conditional 345.2618 0.4614 23
## Constrained 96.8457 0.1294 4
## Unconstrained 306.1461 0.4091 68
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 64.30 21.52 7.60 3.42
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 67.29 26.65 23.15 20.00 16.99 15.68 14.56 12.52
## (Showing 8 of 68 unconstrained eigenvalues)To set up the appropriate permutation design we need to restrict
permutations to be within the 2m x 2m plots. We can do this
using the blocks argument, or via the plots
argument. As some of the other hypotheses we might address with these
data requires permuting under a split plot design, we will use the
plots formulation to test H1
h1 <- with(env,
how(nperm = 999,
within = Within(type = "free"),
plots = Plots(strata = plotid, type = "none")))set.seed(42)
anova(a1, permutations = h1, model = "reduced")## Permutation test for rda under reduced model
## Plots: plotid, plot permutation: none
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = spp ~ year + year:mowing + year:fertilizer + year:removal + Condition(plotid), data = env)
## Df Variance F Pr(>F)
## Model 4 96.846 5.3777 0.001 ***
## Residual 68 306.146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hypothesis 2 is
The temporal trend in species is independent of the treatment. I.e. the individuals treatments do not differ in their temporal dynamics
To address this, we need to remove the average effect of
year and focus on the year effects
within each treatment.
a2 <- rda(spp ~ year:mowing + year:fertilizer + year:removal + Condition(year + plotid),
data = env)
a2## Call: rda(formula = spp ~ year:mowing + year:fertilizer + year:removal
## + Condition(year + plotid), data = env)
##
## Inertia Proportion Rank
## Total 748.25364 1.00000
## Conditional 404.87144 0.54109 24
## Constrained 37.23609 0.04976 3
## Unconstrained 306.14610 0.40915 68
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3
## 24.016 8.362 4.858
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 67.29 26.65 23.15 20.00 16.99 15.68 14.56 12.52
## (Showing 8 of 68 unconstrained eigenvalues)For this hypothesis, under the null hypothesis (no effect) the
temporal changes in the individual plots are independent of the
treatments and are thus exchangeable. However, we should keep the
samples for each plot together (i.e. not permute within each
plot). This design requires the use of plots;
blocks are never permuted, so we can’t shuffle entire
blocks of samples using that mechanism. To create the correct
design-based permutation test then we use
h2 <- with(env,
how(nperm = 999,
within = Within(type = "none"),
plots = Plots(strata = plotid, type = "free")))Now we can use anova() to test hypothesis 2
set.seed(42)
anova(a2, permutations = h2, model = "reduced")## Permutation test for rda under reduced model
## Plots: plotid, plot permutation: free
## Permutation: none
## Number of permutations: 999
##
## Model: rda(formula = spp ~ year:mowing + year:fertilizer + year:removal + Condition(year + plotid), data = env)
## Df Variance F Pr(>F)
## Model 3 37.236 2.7569 0.001 ***
## Residual 68 306.146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The remaining 3 hypotheses are
Fertilisation has no effect on the temporal changes in species composition
Mowing has no effect on the temporal changes in species composition
Removal has no effect on the temporal changes in species composition
To answer each of these hypotheses, would require us to modify
a2 such that all but one of the year -
treatment interactions is removed from the model and added to the
Condition(), e.g.
a3 <- rda(spp ~ year:mowing + Condition(year + year:fertilizer + year:removal + plotid),
data = env)
a3## Call: rda(formula = spp ~ year:mowing + Condition(year +
## year:fertilizer + year:removal + plotid), data = env)
##
## Inertia Proportion Rank
## Total 748.25364 1.00000
## Conditional 430.84672 0.57580 26
## Constrained 11.26081 0.01505 1
## Unconstrained 306.14610 0.40915 68
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1
## 11.261
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 67.29 26.65 23.15 20.00 16.99 15.68 14.56 12.52
## (Showing 8 of 68 unconstrained eigenvalues)We can avoid having to do this by hand for each hypothesis by
recognising that the by = "margin" option of
anova() is going to do this for us. We can replicate the
results generated by Canoco presented in the text book using:
set.seed(42)
anova(a2, by = "margin", permutations = h2, model = "reduced", parallel = 4)## Permutation test for rda under reduced model
## Marginal effects of terms
## Plots: plotid, plot permutation: free
## Permutation: none
## Number of permutations: 999
##
## Model: rda(formula = spp ~ year:mowing + year:fertilizer + year:removal + Condition(year + plotid), data = env)
## Df Variance F Pr(>F)
## year:mowing 1 11.261 2.5012 0.006 **
## year:fertilizer 1 19.790 4.3957 0.001 ***
## year:removal 1 6.185 1.3739 0.205
## Residual 68 306.146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The code below loads the crayfish data and the scale
predictor variables (design) using slightly more modern
code than was used in the slides.
crayfish <- read_csv(url("https://bit.ly/crayfish-spp")) %>%
rename(label = "...1") %>%
janitor::clean_names() %>%
filter(!is.na(label)) %>%
column_to_rownames("label")## New names:
## Rows: 568 Columns: 11
## ── Column specification
## ──────────────────────────────────────────────────────────────── Delimiter: "," chr
## (1): ...1 dbl (10): punct-L, punct-S, march-L, march-S, hubbs-L, hubbs-S, oz-L, oz-S,
## ...
## ℹ 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`design <- read_csv(url("https://bit.ly/crayfish-design")) %>%
rename(label = "...1") %>%
filter(!is.na(label)) %>%
mutate(across(Watershed:`Run Nested`, ~ factor(.x))) %>%
column_to_rownames("label")## New names:
## Rows: 568 Columns: 8
## ── Column specification
## ──────────────────────────────────────────────────────────────── Delimiter: "," chr
## (8): ...1, Watershed, Stream, Reach, Run, Stream Nested, Reach Nested, R...
## ℹ 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`Repeat the exercise from the slides, testing the different levels of the hierarchical scales of variation in crayfish species abundance.
A number of samples have 0 crayfish, which excludes unimodal methods, hence we will use RDA for this analysis.
We first start at the watershed scale
m.ws <- rda(crayfish ~ Watershed, data = design)
m.ws## Call: rda(formula = crayfish ~ Watershed, data = design)
##
## Inertia Proportion Rank
## Total 9.3580 1.0000
## Constrained 1.7669 0.1888 6
## Unconstrained 7.5911 0.8112 10
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## 0.7011 0.5540 0.3660 0.1064 0.0381 0.0013
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 3.0957 1.2109 0.9717 0.7219 0.5333 0.3838 0.2772 0.2040 0.1879 0.0048summary(eigenvals(m.ws, constrained = TRUE))## Argument `constrained` is deprecated; use `model` instead.## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Eigenvalue 0.7011 0.5540 0.3660 0.1064 0.03814 0.0012791
## Proportion Explained 0.3968 0.3135 0.2072 0.0602 0.02159 0.0007239
## Cumulative Proportion 0.3968 0.7103 0.9175 0.9777 0.99928 1.0000000set.seed(1)
ctrl <- how(nperm = 499, within = Within(type = "none"),
plots = with(design, Plots(strata = Stream, type = "free")))
anova(m.ws, permutations = ctrl)## Permutation test for rda under reduced model
## Plots: Stream, plot permutation: free
## Permutation: none
## Number of permutations: 499
##
## Model: rda(formula = crayfish ~ Watershed, data = design)
## Df Variance F Pr(>F)
## Model 6 1.7669 21.724 0.002 **
## Residual 560 7.5911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m.str <- rda(crayfish ~ Stream + Condition(Watershed), data = design)
m.str## Call: rda(formula = crayfish ~ Stream + Condition(Watershed), data =
## design)
##
## Inertia Proportion Rank
## Total 9.3580 1.0000
## Conditional 1.7669 0.1888 6
## Constrained 1.1478 0.1227 10
## Unconstrained 6.4433 0.6885 10
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10
## 0.4928 0.2990 0.2058 0.0782 0.0372 0.0224 0.0063 0.0030 0.0029 0.0002
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 2.7853 0.8528 0.7737 0.6317 0.5144 0.2808 0.2517 0.1923 0.1559 0.0046summary(eigenvals(m.str, constrained = TRUE))## Argument `constrained` is deprecated; use `model` instead.## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 0.4928 0.2990 0.2058 0.07824 0.03719 0.02235 0.006326
## Proportion Explained 0.4293 0.2605 0.1793 0.06816 0.03240 0.01947 0.005511
## Cumulative Proportion 0.4293 0.6898 0.8691 0.93731 0.96971 0.98918 0.994694
## RDA8 RDA9 RDA10
## Eigenvalue 0.003042 0.002894 0.0001546
## Proportion Explained 0.002651 0.002521 0.0001347
## Cumulative Proportion 0.997344 0.999865 1.0000000set.seed(1)
ctrl <- how(nperm = 499, within = Within(type = "none"),
plots = with(design, Plots(strata = Reach, type = "free")),
blocks = with(design, Watershed))
anova(m.str, permutations = ctrl)## Permutation test for rda under reduced model
## Blocks: with(design, Watershed)
## Plots: Reach, plot permutation: free
## Permutation: none
## Number of permutations: 499
##
## Model: rda(formula = crayfish ~ Stream + Condition(Watershed), data = design)
## Df Variance F Pr(>F)
## Model 14 1.1478 6.9477 0.004 **
## Residual 546 6.4433
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m.re <- rda(crayfish ~ Reach + Condition(Stream), data = design)
m.re## Call: rda(formula = crayfish ~ Reach + Condition(Stream), data =
## design)
##
## Inertia Proportion Rank
## Total 9.3580 1.0000
## Conditional 2.9148 0.3115 20
## Constrained 1.4829 0.1585 10
## Unconstrained 4.9603 0.5301 10
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10
## 0.6292 0.2706 0.2146 0.1414 0.1123 0.0467 0.0344 0.0270 0.0064 0.0003
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 2.1635 0.6080 0.5605 0.5166 0.3749 0.2212 0.2052 0.1588 0.1477 0.0040summary(eigenvals(m.re, constrained = TRUE))## Argument `constrained` is deprecated; use `model` instead.## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 0.6292 0.2706 0.2146 0.14135 0.11229 0.04669 0.0344
## Proportion Explained 0.4243 0.1825 0.1447 0.09532 0.07572 0.03148 0.0232
## Cumulative Proportion 0.4243 0.6068 0.7515 0.84682 0.92254 0.95403 0.9772
## RDA8 RDA9 RDA10
## Eigenvalue 0.02704 0.006381 0.0003486
## Proportion Explained 0.01824 0.004303 0.0002351
## Cumulative Proportion 0.99546 0.999765 1.0000000set.seed(1)
ctrl <- how(nperm = 499, within = Within(type = "none"),
plots = with(design, Plots(strata = Run, type = "free")),
blocks = with(design, Stream))
anova(m.re, permutations = ctrl)## Permutation test for rda under reduced model
## Blocks: with(design, Stream)
## Plots: Run, plot permutation: free
## Permutation: none
## Number of permutations: 499
##
## Model: rda(formula = crayfish ~ Reach + Condition(Stream), data = design)
## Df Variance F Pr(>F)
## Model 42 1.4829 3.5875 0.002 **
## Residual 504 4.9603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m.run <- rda(crayfish ~ Run + Condition(Reach), data = design)
m.run## Call: rda(formula = crayfish ~ Run + Condition(Reach), data = design)
##
## Inertia Proportion Rank
## Total 9.3580 1.0000
## Conditional 4.3977 0.4699 62
## Constrained 1.8225 0.1948 10
## Unconstrained 3.1378 0.3353 10
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10
## 0.8541 0.3141 0.1679 0.1393 0.1328 0.0835 0.0474 0.0429 0.0390 0.0016
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## 1.3137 0.4165 0.3832 0.2759 0.2378 0.1725 0.1215 0.1130 0.1016 0.0021summary(eigenvals(m.run, constrained = TRUE))## Argument `constrained` is deprecated; use `model` instead.## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## Eigenvalue 0.8541 0.3141 0.16791 0.1393 0.13278 0.08345 0.04738
## Proportion Explained 0.4687 0.1723 0.09213 0.0764 0.07286 0.04579 0.02600
## Cumulative Proportion 0.4687 0.6410 0.73311 0.8095 0.88237 0.92815 0.95415
## RDA8 RDA9 RDA10
## Eigenvalue 0.04289 0.03904 0.0016298
## Proportion Explained 0.02353 0.02142 0.0008943
## Cumulative Proportion 0.97768 0.99911 1.0000000set.seed(1)
ctrl <- how(nperm = 499, within = Within(type = "free"),
blocks = with(design, Reach))
anova(m.run, permutations = ctrl)## Permutation test for rda under reduced model
## Blocks: with(design, Reach)
## Permutation: free
## Number of permutations: 499
##
## Model: rda(formula = crayfish ~ Run + Condition(Reach), data = design)
## Df Variance F Pr(>F)
## Model 126 1.8225 1.7425 0.002 **
## Residual 378 3.1378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Post (1986) carried out an experiment to investigate the effect of time of ploughing on the subsequent weed composition in summer barley fields. Three different ploughin 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.
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:
## 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.
## • `` -> `...1`
## • `` -> `...2`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()). As the
data are counts, a log-transformation is suggested (because we can’t do
a multivariate Poisson GLM to analyses these data.) Some zero abundances
are recorded as NA values, so we must replace those with 0s
before we proceed.
plough <- plough %>%
mutate(across(che_alb:ape_spi, .fns = tidyr::replace_na, replace = 0))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(che_alb:ape_spi, .fns = tidyr::replace_na, replace = 0)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))log_plough <- log1p(plough)The model is specified as
ord <- rda(log_plough ~ ploughing_time + Condition(block), data = plough_design)
ord## Call: rda(formula = log_plough ~ ploughing_time + Condition(block),
## data = plough_design)
##
## Inertia Proportion Rank
## Total 11.3151 1.0000
## Conditional 5.2128 0.4607 3
## Constrained 2.5286 0.2235 2
## Unconstrained 3.5738 0.3158 6
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2
## 1.8360 0.6926
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 1.6459 0.7674 0.5594 0.3276 0.1645 0.1089We see that almost half the variation in the species composition across the samples is explained by block-to-block differences. Ploughing time explains ~22% of the total variation in log species composition, and ~41% (2.5286 / (11.3151 - 5.2128)) of the variation after removing the block-to-block differences.
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))Using this design we can test the effect of ploughing time
anova(ord, permutations = h)## Set of permutations < 'minperm'. Generating entire set.## Permutation test for rda under reduced model
## Blocks: block
## Permutation: free
## Number of permutations: 1295
##
## Model: rda(formula = log_plough ~ ploughing_time + Condition(block), data = plough_design)
## Df Variance F Pr(>F)
## Model 2 2.5286 2.1226 0.027 *
## Residual 6 3.5738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Notice 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. In such circumstances, if we randomly generated the permutations we would get many repeated permutations. Instead, it is better to generate the entire set of permutations and do an exact test. Hence, even though we asked for 999 permutations, the results indicate that 1295 permutations were used. This corresponds to the number of possible permutations of the data
# this includes the observed ordering
numPerms(log_plough, control = h)## [1] 1296Another option for margin that is new in recent versions
of vegan, is "onedf". This will perform tests of
one degree of freedom contrasts. For these data it means testing the
difference between ploughing time 1 and 2, and between
ploughing time 1 and 3
anova(ord, permutations = h, by = "onedf")## Set of permutations < 'minperm'. Generating entire set.## Permutation test for rda under reduced model
## Sequential test for contrasts
## Blocks: block
## Permutation: free
## Number of permutations: 1295
##
## Model: rda(formula = log_plough ~ ploughing_time + Condition(block), data = plough_design)
## Df Variance F Pr(>F)
## ploughing_timepltime 2 1 1.0407 1.7472 0.081 .
## ploughing_timepltime 3 1 1.4879 2.4980 0.014 *
## Residual 6 3.5738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This isn’t so useful here, but if you had a control level and it was coded as the reference level of the factor, this would give you tests of control versus treatment level contrasts.
This example uses a data set of Swedish pine forest undergrowth
assemblages that were experimentally manipulated through the addition of
nitrogen N and phosphorus P. There are 32
samples in 4 complete, randomized bloacks of 8 plots each. In each
block, all combinations of 4 levels of nitrogen N (applied
as ammonium nitrate) and 2 levels of phosphorus P (applied
as compound PK fertilizer) were applied.
library("tidyr")
library("dplyr")
library("readr")
e40_spp <- read_csv(url("https://bit.ly/e40-spp"), skip = 1L) %>%
rename(label = "...1") %>%
select(-`...2`) %>%
mutate(across(-label, ~ replace_na(.x, 0)),
across(-label, ~ log1p(.x))) %>%
tibble::column_to_rownames("label")## New names:
## Rows: 32 Columns: 105
## ── Column specification
## ──────────────────────────────────────────────────────────────── Delimiter: "," dbl
## (105): ...1, ...2, agroscap, anthoodo, arctouva, aulacpal, betulpen, bet...
## ℹ 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`e40_design <- read_csv(url("https://bit.ly/e40-design"), skip = 1L) %>%
rename(label = "...1") %>%
select(-`...2`) %>%
mutate(across(-label, ~ factor(.x))) %>%
tibble::column_to_rownames("label")## New names:
## Rows: 32 Columns: 5
## ── Column specification
## ──────────────────────────────────────────────────────────────── Delimiter: "," chr
## (3): N, P, Block dbl (2): ...1, ...2
## ℹ 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`The levels of N are
with(e40_design, levels(N))## [1] "N0" "N1" "N2" "N3"where N0 means no N addition
The levels of P are
with(e40_design, levels(P))## [1] "P0" "P1"where P0 means no P addition.
The variable Block codes for the 4 blocks of 8
plots.
In 1987 the undergrowth of the plots was surveyed. The percentage
cover of plant species (including tree saplings, bryophytes, and
lichens) was estimated and recorded on a ten-point scale. There are 103
species recorded. The ten-point scale values were re-expressed as the
mid-points of the ten-point percentage scale in e40_spp.
The species data are log(x + 1) transformed.
The aim of the analysis is to test the effects of
using a combination of models and appropriate permutation tests.
m_n <- rda(e40_spp ~ N + Condition(Block + P), data = e40_design)
m_n## Call: rda(formula = e40_spp ~ N + Condition(Block + P), data =
## e40_design)
##
## Inertia Proportion Rank
## Total 9.0194 1.0000
## Conditional 1.2056 0.1337 4
## Constrained 4.9100 0.5444 3
## Unconstrained 2.9038 0.3219 24
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3
## 4.471 0.397 0.041
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.8460 0.5398 0.3780 0.1977 0.1840 0.1575 0.1329 0.0977
## (Showing 8 of 24 unconstrained eigenvalues)h_blk <- with(e40_design, how(blocks = Block, nperm = 999))
set.seed(245)
anova(m_n, permutations = h_blk)## Permutation test for rda under reduced model
## Blocks: Block
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp ~ N + Condition(Block + P), data = e40_design)
## Df Variance F Pr(>F)
## Model 3 4.9100 13.527 0.001 ***
## Residual 24 2.9038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_p <- rda(e40_spp ~ P + Condition(Block + N), data = e40_design)
m_p## Call: rda(formula = e40_spp ~ P + Condition(Block + N), data =
## e40_design)
##
## Inertia Proportion Rank
## Total 9.01938 1.00000
## Conditional 5.65417 0.62689 6
## Constrained 0.46143 0.05116 1
## Unconstrained 2.90378 0.32195 24
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1
## 0.4614
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.8460 0.5398 0.3780 0.1977 0.1840 0.1575 0.1329 0.0977
## (Showing 8 of 24 unconstrained eigenvalues)set.seed(48)
anova(m_p, permutations = h_blk)## Permutation test for rda under reduced model
## Blocks: Block
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp ~ P + Condition(Block + N), data = e40_design)
## Df Variance F Pr(>F)
## Model 1 0.46143 3.8138 0.002 **
## Residual 24 2.90378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1by = "margin"# or
m_margin <- rda(e40_spp ~ N + P + Condition(Block), data = e40_design)
m_margin## Call: rda(formula = e40_spp ~ N + P + Condition(Block), data =
## e40_design)
##
## Inertia Proportion Rank
## Total 9.01938 1.00000
## Conditional 0.74421 0.08251 3
## Constrained 5.37139 0.59554 4
## Unconstrained 2.90378 0.32195 24
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 4.712 0.492 0.127 0.040
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.8460 0.5398 0.3780 0.1977 0.1840 0.1575 0.1329 0.0977
## (Showing 8 of 24 unconstrained eigenvalues)set.seed(23)
anova(m_margin, permutations = h_blk, by = "margin")## Permutation test for rda under reduced model
## Marginal effects of terms
## Blocks: Block
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp ~ N + P + Condition(Block), data = e40_design)
## Df Variance F Pr(>F)
## N 3 4.9100 13.5271 0.001 ***
## P 1 0.4614 3.8138 0.029 *
## Residual 24 2.9038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_np <- rda(e40_spp ~ N + P + Condition(Block), data = e40_design)
m_np## Call: rda(formula = e40_spp ~ N + P + Condition(Block), data =
## e40_design)
##
## Inertia Proportion Rank
## Total 9.01938 1.00000
## Conditional 0.74421 0.08251 3
## Constrained 5.37139 0.59554 4
## Unconstrained 2.90378 0.32195 24
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 4.712 0.492 0.127 0.040
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.8460 0.5398 0.3780 0.1977 0.1840 0.1575 0.1329 0.0977
## (Showing 8 of 24 unconstrained eigenvalues)set.seed(23)
anova(m_np, permutations = h_blk)## Permutation test for rda under reduced model
## Blocks: Block
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp ~ N + P + Condition(Block), data = e40_design)
## Df Variance F Pr(>F)
## Model 4 5.3714 11.099 0.001 ***
## Residual 24 2.9038
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_int <- rda(e40_spp ~ N:P + Condition(Block + N + P), data = e40_design)
m_int## Call: rda(formula = e40_spp ~ N:P + Condition(Block + N + P), data =
## e40_design)
##
## Inertia Proportion Rank
## Total 9.01938 1.00000
## Conditional 6.11560 0.67805 7
## Constrained 0.31669 0.03511 3
## Unconstrained 2.58709 0.28684 21
## Inertia is variance
## Some constraints or conditions were aliased because they were redundant
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3
## 0.15960 0.11883 0.03825
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.8178 0.5116 0.2798 0.1967 0.1543 0.1409 0.1073 0.0961
## (Showing 8 of 21 unconstrained eigenvalues)h_blk <- with(e40_design, how(blocks = Block, nperm = 999))
set.seed(42)
anova(m_int, permutations = h_blk)## Permutation test for rda under reduced model
## Blocks: Block
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp ~ N:P + Condition(Block + N + P), data = e40_design)
## Df Variance F Pr(>F)
## Model 3 0.31669 0.8569 0.646
## Residual 21 2.58709m_int2 <- rda(e40_spp ~ N + P + N:P + Condition(Block), data = e40_design)
m_int2## Call: rda(formula = e40_spp ~ N + P + N:P + Condition(Block), data =
## e40_design)
##
## Inertia Proportion Rank
## Total 9.01938 1.00000
## Conditional 0.74421 0.08251 3
## Constrained 5.68808 0.63065 7
## Unconstrained 2.58709 0.28684 21
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## 4.727 0.548 0.185 0.110 0.065 0.040 0.013
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.8178 0.5116 0.2798 0.1967 0.1543 0.1409 0.1073 0.0961
## (Showing 8 of 21 unconstrained eigenvalues)set.seed(42)
anova(m_int2, permutations = h_blk, by = "margin")## Permutation test for rda under reduced model
## Marginal effects of terms
## Blocks: Block
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp ~ N + P + N:P + Condition(Block), data = e40_design)
## Df Variance F Pr(>F)
## N:P 3 0.31669 0.8569 0.511
## Residual 21 2.58709In general, design-based tests of interaction terms are not possible. However, for simple two-way designs such as the one we have here, it is possible to test for an interaction with a design-based test. For this to work we need to think of the data are being split into 4 blocks each with
NPAnd the ordering of the levels of N and P
must be the same within each block of samples. If the data are arranged
precisely like this, then we can permute the rows and the columns of
these blocks of data to test the interaction effect.
Arrange the data in the correct order
e40 <- e40_spp %>%
bind_cols(e40_design) %>%
group_by(Block) %>%
arrange(N, P, .by_group = TRUE) %>%
ungroup()
e40_spp2 <- e40 %>%
select(agroscap:veronoff)
e40_design2 <- e40 %>%
select(N:Block)
ord <- rda(e40_spp2 ~ N:P + Condition(N + P + Block), data = e40_design2)
h_int <- with(e40_design2,
how(nperm = 999, plots = Plots(strata = N, type = "free"),
blocks = Block,
within = Within(type = "free")))
anova(ord, permutations = h_int)## Permutation test for rda under reduced model
## Blocks: Block
## Plots: N, plot permutation: free
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = e40_spp2 ~ N:P + Condition(N + P + Block), data = e40_design2)
## Df Variance F Pr(>F)
## Model 3 0.31669 0.8569 0.746
## Residual 21 2.58709