Birds Data
This data set is on bird assemblages in a montane forest in the Velka Fatra Mountains (Slovak Republic).
library("vegan")
library("readr")
library("tidyr")
library("dplyr")
library("tibble")We can load the bird data and process it into two data sets
- the species abundances, and
- the associated environmental variables.
birds_url <- "https://bit.ly/maed-birds"
birds <- read_csv(birds_url, na = "", skip = 1)## New names:
## Rows: 43 Columns: 51
## ── Column specification
## ────────────────────────────────────────────────────────────────────── Delimiter: "," chr
## (4): ...1, ...40, Rocks, Expos dbl (46): AegiCaud, AlauArve, AnthTriv, CertFami, ColuPalu,
## CorvCora, CucuCa... lgl (1): ...39
## ℹ 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`
## • `` -> `...39`
## • `` -> `...40`birds_spp <- birds |>
mutate(across(AegiCaud:OriOri, ~ replace_na(.x, 0))) |>
column_to_rownames("...1")
birds_spp <- birds_spp[, 1:36]
birds_env <- birds |>
select(c("...1", Altit:Expos)) |>
mutate(Rocks = factor(Rocks), Expos = factor(Expos)) |>
column_to_rownames("...1")The aim of the analysis is to find a good ordination of the species data (choose your ordination method, transformation, dissimilarity, etc according to what features of the species data you want to focus on). Then, having found a good ordination, use the environmental data to try to explain the ordination axes or mapping you have created.
The variables are:
- Altitude
- Forest cover
- Forest density
- Amount of broad leaf trees
- E2 layer (scrub layer cover)
- E2 layer conifers
- E1 layer (herb layer cover)
- E1 height
- Slope
- Presence of rocks
- Exposure
How the Canoco developers would analyse these data
They treat the data as “compositional” (despite them being the average number of breeding pairs averaged over 4 visits to each sampling location).
birds_pca <- pca(birds_spp)
birds_pca##
## Call: pca(X = birds_spp)
##
## Inertia Rank
## Total 4.555
## Unconstrained 4.555 36
##
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.9069 0.7557 0.4542 0.4004 0.3159 0.2353 0.2228 0.1959
## (Showing 8 of 36 unconstrained eigenvalues)summary(eigenvals(birds_pca))## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 0.9069 0.7557 0.45421 0.40040 0.31589 0.23532 0.22275
## Proportion Explained 0.1991 0.1659 0.09973 0.08791 0.06936 0.05167 0.04891
## Cumulative Proportion 0.1991 0.3651 0.46478 0.55269 0.62205 0.67372 0.72262
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Eigenvalue 0.19587 0.16030 0.13292 0.11956 0.1039 0.09980 0.08050
## Proportion Explained 0.04301 0.03519 0.02918 0.02625 0.0228 0.02191 0.01767
## Cumulative Proportion 0.76563 0.80082 0.83001 0.85626 0.8791 0.90097 0.91865
## PC15 PC16 PC17 PC18 PC19 PC20
## Eigenvalue 0.06215 0.05467 0.04062 0.035465 0.032653 0.028140
## Proportion Explained 0.01364 0.01200 0.00892 0.007787 0.007169 0.006178
## Cumulative Proportion 0.93229 0.94430 0.95322 0.961002 0.968172 0.974350
## PC21 PC22 PC23 PC24 PC25 PC26
## Eigenvalue 0.021960 0.020116 0.014740 0.011911 0.009596 0.008122
## Proportion Explained 0.004822 0.004417 0.003236 0.002615 0.002107 0.001783
## Cumulative Proportion 0.979172 0.983589 0.986825 0.989440 0.991547 0.993330
## PC27 PC28 PC29 PC30 PC31 PC32
## Eigenvalue 0.007375 0.004942 0.0045034 0.0034125 0.0031614 0.0025260
## Proportion Explained 0.001619 0.001085 0.0009888 0.0007492 0.0006941 0.0005546
## Cumulative Proportion 0.994950 0.996035 0.9970237 0.9977730 0.9984671 0.9990217
## PC33 PC34 PC35 PC36
## Eigenvalue 0.0019434 0.0012993 0.0007002 0.0005127
## Proportion Explained 0.0004267 0.0002853 0.0001537 0.0001126
## Cumulative Proportion 0.9994484 0.9997337 0.9998874 1.0000000plot(birds_pca)
ev_pca <- envfit(birds_pca ~ Altit, data = birds_env)
ev_pca##
## ***VECTORS
##
## PC1 PC2 r2 Pr(>r)
## Altit -0.96310 0.26913 0.4583 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(ev_pca, add = TRUE)
ordisurf(birds_pca ~ Altit, data = birds_env, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 5.21 total = 6.21
##
## REML score: 257.1538How I might have approached this
birds_ca <- ca(birds_spp)
birds_ca##
## Call: ca(X = birds_spp)
##
## Inertia Rank
## Total 1.476
## Unconstrained 1.476 35
##
## Inertia is scaled Chi-square
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.21320 0.17737 0.11419 0.11092 0.09765 0.09175 0.07321 0.06225
## (Showing 8 of 35 unconstrained eigenvalues)summary(eigenvals(birds_ca))## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.2132 0.1774 0.11419 0.11092 0.09765 0.09175 0.07321
## Proportion Explained 0.1444 0.1201 0.07734 0.07513 0.06614 0.06215 0.04959
## Cumulative Proportion 0.1444 0.2645 0.34188 0.41701 0.48315 0.54529 0.59488
## CA8 CA9 CA10 CA11 CA12 CA13 CA14
## Eigenvalue 0.06225 0.05758 0.05531 0.04542 0.03964 0.03699 0.03535
## Proportion Explained 0.04216 0.03900 0.03746 0.03076 0.02685 0.02505 0.02394
## Cumulative Proportion 0.63704 0.67604 0.71350 0.74426 0.77111 0.79616 0.82011
## CA15 CA16 CA17 CA18 CA19 CA20 CA21
## Eigenvalue 0.03089 0.02835 0.02686 0.02382 0.02193 0.01941 0.01854
## Proportion Explained 0.02092 0.01920 0.01819 0.01613 0.01485 0.01315 0.01256
## Cumulative Proportion 0.84103 0.86023 0.87842 0.89456 0.90941 0.92255 0.93511
## CA22 CA23 CA24 CA25 CA26 CA27
## Eigenvalue 0.01745 0.013675 0.012507 0.011456 0.009717 0.007770
## Proportion Explained 0.01182 0.009262 0.008471 0.007759 0.006581 0.005263
## Cumulative Proportion 0.94693 0.956191 0.964661 0.972421 0.979002 0.984265
## CA28 CA29 CA30 CA31 CA32 CA33
## Eigenvalue 0.006807 0.005525 0.003228 0.002477 0.001760 0.001593
## Proportion Explained 0.004610 0.003742 0.002187 0.001677 0.001192 0.001079
## Cumulative Proportion 0.988875 0.992617 0.994804 0.996481 0.997673 0.998752
## CA34 CA35
## Eigenvalue 0.0012725 0.0005700
## Proportion Explained 0.0008619 0.0003861
## Cumulative Proportion 0.9996139 1.0000000plot(birds_ca)
ev_ca <- envfit(birds_ca ~ Altit, data = birds_env)
ev_ca##
## ***VECTORS
##
## CA1 CA2 r2 Pr(>r)
## Altit 0.47349 -0.88080 0.5684 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(ev_ca, add = TRUE)
ordisurf(birds_ca ~ Altit, data = birds_env, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 6.76 total = 7.76
##
## REML score: 248.5029or
birds_h <- pca(decostand(birds_spp, method = "hellinger"))
birds_h##
## Call: pca(X = decostand(birds_spp, method = "hellinger"))
##
## Inertia Rank
## Total 0.313
## Unconstrained 0.313 36
##
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.05735 0.03967 0.02744 0.02666 0.02018 0.01619 0.01463 0.01420
## (Showing 8 of 36 unconstrained eigenvalues)summary(eigenvals(birds_h))## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 0.05735 0.03967 0.02744 0.02666 0.02018 0.01619 0.01463
## Proportion Explained 0.18321 0.12672 0.08766 0.08516 0.06446 0.05171 0.04672
## Cumulative Proportion 0.18321 0.30993 0.39759 0.48275 0.54721 0.59892 0.64564
## PC8 PC9 PC10 PC11 PC12 PC13
## Eigenvalue 0.01420 0.01132 0.01078 0.009609 0.008787 0.007767
## Proportion Explained 0.04537 0.03617 0.03444 0.030695 0.028068 0.024810
## Cumulative Proportion 0.69101 0.72718 0.76162 0.792312 0.820380 0.845190
## PC14 PC15 PC16 PC17 PC18 PC19
## Eigenvalue 0.006586 0.005725 0.004716 0.004524 0.00399 0.003512
## Proportion Explained 0.021039 0.018288 0.015066 0.014451 0.01274 0.011218
## Cumulative Proportion 0.866229 0.884518 0.899584 0.914035 0.92678 0.937998
## PC20 PC21 PC22 PC23 PC24 PC25
## Eigenvalue 0.003175 0.002641 0.002359 0.002027 0.001833 0.001592
## Proportion Explained 0.010141 0.008437 0.007537 0.006475 0.005857 0.005086
## Cumulative Proportion 0.948139 0.956576 0.964112 0.970587 0.976444 0.981530
## PC26 PC27 PC28 PC29 PC30 PC31
## Eigenvalue 0.001363 0.001070 0.0008743 0.0007373 0.0006203 0.0004077
## Proportion Explained 0.004353 0.003418 0.0027928 0.0023552 0.0019816 0.0013022
## Cumulative Proportion 0.985883 0.989301 0.9920941 0.9944494 0.9964310 0.9977332
## PC32 PC33 PC34 PC35 PC36
## Eigenvalue 0.0002571 0.0002058 0.0001123 7.409e-05 6.044e-05
## Proportion Explained 0.0008212 0.0006573 0.0003586 2.367e-04 1.931e-04
## Cumulative Proportion 0.9985544 0.9992117 0.9995703 9.998e-01 1.000e+00plot(birds_h)
ev_h <- envfit(birds_h ~ Altit, data = birds_env)
ev_h##
## ***VECTORS
##
## PC1 PC2 r2 Pr(>r)
## Altit -0.86211 0.50673 0.4906 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(ev_h, add = TRUE)
ordisurf(birds_h ~ Altit, data = birds_env, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 5.54 total = 6.54
##
## REML score: 251.2566or
birds_nmds <- metaMDS(decostand(birds_spp, method = "hellinger"),
dist = "euclidean",
autotransform = FALSE, expand = FALSE)## Run 0 stress 0.2233509
## Run 1 stress 0.2231504
## ... New best solution
## ... Procrustes: rmse 0.03760199 max resid 0.2214757
## Run 2 stress 0.2237093
## Run 3 stress 0.2237343
## Run 4 stress 0.2237093
## Run 5 stress 0.2598602
## Run 6 stress 0.2273902
## Run 7 stress 0.2372302
## Run 8 stress 0.2400796
## Run 9 stress 0.2665373
## Run 10 stress 0.2603492
## Run 11 stress 0.2348269
## Run 12 stress 0.2341483
## Run 13 stress 0.271679
## Run 14 stress 0.2341093
## Run 15 stress 0.2275538
## Run 16 stress 0.2343994
## Run 17 stress 0.2701225
## Run 18 stress 0.2443406
## Run 19 stress 0.2231427
## ... New best solution
## ... Procrustes: rmse 0.0016735 max resid 0.008080747
## ... Similar to previous best
## Run 20 stress 0.226002
## *** Best solution repeated 1 timesbirds_nmds##
## Call:
## metaMDS(comm = decostand(birds_spp, method = "hellinger"), distance = "euclidean", autotransform = FALSE, expand = FALSE)
##
## global Multidimensional Scaling using monoMDS
##
## Data: decostand(birds_spp, method = "hellinger")
## Distance: euclidean
##
## Dimensions: 2
## Stress: 0.2231427
## Stress type 1, weak ties
## Best solution was repeated 1 time in 20 tries
## The best solution was from try 19 (random start)
## Scaling: centring, PC rotation
## Species: non-expanded scores based on 'decostand(birds_spp, method = "hellinger")'stressplot(birds_nmds)
ordipointlabel(birds_nmds)
ev_nmds <- envfit(birds_nmds ~ Altit, data = birds_env)
ev_nmds##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Altit 0.72538 0.68835 0.443 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(ev_nmds, add = TRUE)
ordisurf(birds_nmds ~ Altit, data = birds_env, add = TRUE)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 5.45 total = 6.45
##
## REML score: 255.986But does it make any difference?
layout(matrix(1:4, ncol = 2, byrow = TRUE))
plot(birds_pca, main = "PCA")
plot(ev_pca, add = TRUE)
ordisurf(birds_pca ~ Altit, data = birds_env, add = TRUE,
col = "forestgreen")##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 5.21 total = 6.21
##
## REML score: 257.1538plot(birds_ca, main = "CA")
plot(ev_ca, add = TRUE)
ordisurf(birds_ca ~ Altit, data = birds_env, add = TRUE,
col = "forestgreen")##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 6.76 total = 7.76
##
## REML score: 248.5029plot(birds_h, main = "PCA (Hellinger)")
plot(ev_h, add = TRUE)
ordisurf(birds_h ~ Altit, data = birds_env, add = TRUE,
col = "forestgreen")##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 5.54 total = 6.54
##
## REML score: 251.2566ordipointlabel(birds_nmds, main = "NMDS")
plot(ev_nmds, add = TRUE)
ordisurf(birds_nmds ~ Altit, data = birds_env, add = TRUE,
col = "forestgreen")
##
## Family: gaussian
## Link function: identity
##
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
##
## Estimated degrees of freedom:
## 5.45 total = 6.45
##
## REML score: 255.986layout(1)