First year results: Question 2
SP
2023-21-02
Explaining plant assemblies
Effect of soil and grazing-related variables on plant community in Castril, Santiago and Pontones (CSP) pasturelands.
Environmental variables
Grazing-related variables
Spatio-temporal livestock mobility:
Stocking Rate considering total days in the pastureland.
\[\text{Stocking Rate} = \frac{Heads * Days}{Hectares}\] We have these data for 24 pastoral units.
hist(soil_6$stock, xlab = "Stocking Rate")
Plant density i.e total plant touch by transect
hist(soil_6$pl_dens, xlab = "Plant density")
Soil variables
hist(soil_6$soil_dens, xlab = "Soil bulk density")
hist(soil_6$C, xlab = "C %")
Correlation matrix
res <- cor(soil_6%>% filter(stock>1), method = c("spearman"))
res <- round(res, 2)
upper<-res
upper[upper.tri(res,diag=FALSE)]<-""
upper<-as.data.frame(upper) ;kable(upper)| bare | litt | ston | rock | mos | ph | soil_dens | depth | hcl | chro | val | ur_t | ur_fh | dung_au | dung_sp | stock | t_day | e_day | pl_dens | om | N | C | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| bare | 1 | |||||||||||||||||||||
| litt | 0.21 | 1 | ||||||||||||||||||||
| ston | -0.47 | -0.6 | 1 | |||||||||||||||||||
| rock | -0.19 | 0.19 | -0.28 | 1 | ||||||||||||||||||
| mos | -0.21 | 0.28 | -0.35 | 0.24 | 1 | |||||||||||||||||
| ph | -0.1 | 0.1 | 0.23 | -0.23 | -0.14 | 1 | ||||||||||||||||
| soil_dens | 0.13 | 0.11 | 0 | -0.32 | -0.1 | 0.14 | 1 | |||||||||||||||
| depth | 0.51 | -0.22 | -0.07 | -0.09 | -0.43 | 0.09 | 0.04 | 1 | ||||||||||||||
| hcl | -0.11 | -0.15 | 0.48 | -0.53 | -0.12 | 0.48 | 0.11 | -0.1 | 1 | |||||||||||||
| chro | 0.21 | -0.07 | -0.02 | -0.09 | -0.03 | 0.1 | 0.31 | 0.45 | 0.09 | 1 | ||||||||||||
| val | 0.08 | 0.22 | -0.09 | -0.31 | 0.14 | 0.5 | 0.38 | 0.17 | 0.4 | 0.28 | 1 | |||||||||||
| ur_t | 0.31 | 0.1 | -0.42 | -0.06 | 0.27 | -0.04 | 0.23 | 0.12 | 0.15 | 0.37 | 0.39 | 1 | ||||||||||
| ur_fh | 0.01 | -0.06 | -0.18 | 0.01 | -0.04 | -0.11 | 0.23 | 0.12 | -0.07 | 0.03 | 0.13 | 0.57 | 1 | |||||||||
| dung_au | -0.04 | 0.22 | 0 | 0.14 | 0.31 | -0.15 | -0.01 | -0.3 | 0.15 | 0.23 | -0.01 | 0 | -0.22 | 1 | ||||||||
| dung_sp | -0.13 | 0.17 | -0.16 | 0.22 | 0.17 | -0.3 | -0.14 | -0.44 | -0.16 | -0.21 | -0.32 | 0.22 | 0.25 | 0.52 | 1 | |||||||
| stock | 0 | 0.26 | -0.23 | 0.26 | 0.16 | -0.2 | 0.01 | -0.3 | -0.32 | -0.38 | -0.15 | 0.01 | 0.19 | -0.01 | 0.48 | 1 | ||||||
| t_day | -0.19 | 0.16 | 0.27 | -0.02 | -0.03 | 0.06 | 0.26 | -0.35 | 0.23 | 0.09 | 0.15 | -0.32 | -0.55 | 0.31 | -0.07 | 0.06 | 1 | |||||
| e_day | 0.03 | -0.16 | -0.19 | 0.11 | 0.02 | -0.34 | -0.24 | 0.05 | -0.17 | 0 | -0.3 | 0.06 | 0.28 | 0.1 | 0.27 | -0.01 | -0.65 | 1 | ||||
| pl_dens | -0.46 | 0.15 | -0.07 | 0.13 | 0.42 | 0.11 | -0.29 | -0.53 | 0.13 | -0.54 | -0.11 | -0.11 | 0.18 | -0.23 | 0.08 | 0.27 | -0.1 | 0.03 | 1 | |||
| om | -0.33 | -0.25 | 0.4 | 0.44 | 0.05 | -0.2 | -0.65 | -0.31 | -0.17 | -0.47 | -0.6 | -0.56 | -0.35 | 0.18 | 0.2 | 0.06 | -0.03 | 0.19 | 0.19 | 1 | ||
| N | -0.33 | -0.23 | 0.38 | 0.42 | 0.11 | -0.25 | -0.61 | -0.34 | -0.24 | -0.48 | -0.59 | -0.52 | -0.3 | 0.21 | 0.26 | 0.08 | -0.04 | 0.18 | 0.13 | 0.98 | 1 | |
| C | -0.04 | 0.05 | 0.08 | 0.04 | 0.22 | 0.42 | -0.44 | -0.25 | 0.07 | -0.44 | 0.05 | -0.21 | -0.18 | 0.1 | 0.22 | 0.28 | -0.12 | 0.02 | 0.25 | 0.45 | 0.47 | 1 |
correlation between Soil organic matter and Total nitrogen
ggplot(soil_6,aes(x=N,y=om)) + geom_point() + xlab("Soil organic matter %") + ylab("Total nitrogen %") + annotate(geom="text", x=0.3, y=20, label="rho = 0.98",color="black")
Total nitrogen is highly related to soil organic matter. Also these two variables are slightly and negatively associated with Plant Utilisation Rate (PUR) and soil bulk density.
Due to multicollinearity, we just keep soil organic matter in our model.
Redundancy Analysis
Model with environmental variables; soil and grazing-related variables
Filter out rows S05s and S10l as we could not collect information related to stocking rate
csp_rda_1_bis <- capscale(community_data_presence_2_2[c(-11,-16),] ~ bare + litt + ston + rock + mos + ph + soil_dens + depth + hcl + chro + val + ur_fh + ur_t + dung_sp + stock + t_day + pl_dens + om + e_day + C + Condition(scores(coord[c(-11,-16),])),
data = soil_6%>% filter(stock>1), distance = "bray")
csp_rda_1_bis## Call: capscale(formula = community_data_presence_2_2[c(-11, -16), ] ~
## bare + litt + ston + rock + mos + ph + soil_dens + depth + hcl + chro +
## val + ur_fh + ur_t + dung_sp + stock + t_day + pl_dens + om + e_day + C
## + Condition(scores(coord[c(-11, -16), ])), data = soil_6 %>%
## filter(stock > 1), distance = "bray")
##
## Inertia Proportion Rank
## Total 4.646118 1.000000
## Conditional 0.570398 0.122769 2
## Constrained 3.947245 0.849579 20
## Unconstrained 0.167249 0.035998 1
## Imaginary -0.038775 -0.008346 3
## Inertia is squared Bray distance
## Species scores projected from '[' 'community_data_presence_2_2' 'c(-11, -16)' ''
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3 CAP4 CAP5 CAP6 CAP7 CAP8 CAP9 CAP10 CAP11
## 0.9645 0.5273 0.4305 0.3465 0.3239 0.2351 0.2321 0.1782 0.1481 0.1344 0.0835
## CAP12 CAP13 CAP14 CAP15 CAP16 CAP17 CAP18 CAP19 CAP20
## 0.0766 0.0688 0.0640 0.0529 0.0311 0.0240 0.0189 0.0065 0.0004
##
## Eigenvalues for unconstrained axes:
## MDS1
## 0.16725s_csp_rda_1_bis = summary(csp_rda_1_bis)
s_csp_rda_1_bis$cont$importance[,1:4]## CAP1 CAP2 CAP3 CAP4
## Eigenvalue 0.9645051 0.5272634 0.4305289 0.34652549
## Proportion Explained 0.2344164 0.1281478 0.1046371 0.08422067
## Cumulative Proportion 0.2344164 0.3625642 0.4672014 0.5514220412% of inertia explained by the geographical distance between transect. 85% of inertia is related to the explanatory variables.
But this R2 (85%) is biased by the number of variable included in the model. So R2 is commonly adjusted using the Ezekiel’s formula:
The adjustment of R2 using Ezekiel’s formula turns out:
RsquareAdj(csp_rda_1_bis)$adj.r.squared## [1] 0.3247823Axes one and txo of RDA
sp_sel = s_csp_rda_1_bis$species[,1:3] %>% as.data.frame %>% mutate(position = seq(1:209)) %>% filter(if_any(c(CAP1,CAP2,CAP3), ~ .x >= 0.1)|if_any(c(CAP1,CAP2,CAP3), ~ .x <= -0.1)) %>% pull(position)
ggtriplotRDA(csp_rda_1_bis,site.sc = "lc", plot.sites = FALSE,plot.spe = FALSE, scaling = 2,select.spe = sp_sel)##
## -----------------------------------------------------------------------
## Site constraints (lc) selected. To obtain site scores that are weighted
## sums of species scores (default in vegan), argument site.sc must be set
## to wa.
## -----------------------------------------------------------------------
##
## No factor, hence levels cannot be plotted with symbols;
## 'plot.centr' is set to FALSE
Figure 1: Scaling 2 i.e angles between variables reflect their correlation
##
## -----------------------------------------------------------------------
## Site constraints (lc) selected. To obtain site scores that are weighted
## sums of species scores (default in vegan), argument site.sc must be set
## to wa.
## -----------------------------------------------------------------------
##
## No factor, hence levels cannot be plotted with symbols;
## 'plot.centr' is set to FALSE
Figure 2: Type Scaling 1 which focuses on the distance relationships among species
Axes one and three of RDA
ggtriplotRDA(ax2 =3,csp_rda_1_bis,site.sc = "lc", plot.sites = FALSE,plot.spe = FALSE,scaling = 2,select.spe = sp_sel)##
## -----------------------------------------------------------------------
## Site constraints (lc) selected. To obtain site scores that are weighted
## sums of species scores (default in vegan), argument site.sc must be set
## to wa.
## -----------------------------------------------------------------------
##
## No factor, hence levels cannot be plotted with symbols;
## 'plot.centr' is set to FALSE
Figure 3: Type Scaling 1
Contribution of each variables in each canonical axes:
s_csp_rda_1_bis$biplot[,1:3]## CAP1 CAP2 CAP3
## bare -0.29862854 -0.072780923 -0.442263326
## litt -0.18427010 0.218293922 0.051400016
## ston 0.05835504 -0.223745400 -0.199632979
## rock 0.14667877 0.265936203 0.026797942
## mos 0.25687666 0.126568389 0.521405647
## ph -0.12259612 -0.135233144 -0.237293904
## soil_dens 0.12036683 -0.213944332 0.276170482
## depth 0.28709645 -0.247352206 -0.215339641
## hcl 0.13786338 0.073925579 -0.270868321
## chro 0.41546999 -0.239419127 -0.248944253
## val 0.01899064 -0.177363317 0.068978437
## ur_fh 0.08330649 -0.236039236 0.124186817
## ur_t 0.12339227 -0.186303358 -0.081231673
## dung_sp -0.07141911 0.053409322 -0.056408518
## stock -0.19123321 -0.078398860 0.028774691
## t_day 0.14526910 -0.032541712 -0.074857719
## pl_dens -0.08203157 0.329357571 0.378613842
## om -0.10724910 0.139984634 -0.055297345
## e_day 0.19545957 0.020099986 0.065018734
## C -0.51324264 0.007256343 0.005342516Variance Inflation Factors (VIF)
Linear dependencies can be explored by computing the X variables’ Variance Inflation Factors (VIF). VIF greater than 10 should be avoided.
vif.cca(csp_rda_1_bis)## scores(coord[c(-11, -16), ])x scores(coord[c(-11, -16), ])y
## 63.41917 114.19753
## bare litt
## 54.41976 11.95434
## ston rock
## 77.56150 21.03297
## mos ph
## 16.24591 34.81459
## soil_dens depth
## 34.33008 19.67644
## hcl chro
## 28.40969 25.82693
## val ur_fh
## 20.86687 16.75383
## ur_t dung_sp
## 65.54123 14.34461
## stock t_day
## 13.88128 145.97920
## pl_dens om
## 23.58503 36.85092
## e_day C
## 82.49753 17.38069So our first model features the problem of multicollinearity.Yet, we keep all the variables to have a general view of the relationship between variables.
Result interpretation
From this first model we can draw that Our results reveal that axis 1 is mainly related to soil variables. This axis discriminates plant communities characterized by Thymus serpylloides with a large surface occupied by bare soil, high total carbon from plant communities with Poa bulbosa having a deep soil depth, high chromaticity which corresponds to light soils, therefore not so organic. Interestingly this soil features a high presence of moss (Figure 1-2).
Axis 2 suggests that grazing pressure leads to the differentiation of high-presence communities of Festuca hystrix vs. Poa ligulata and Koeleria vallesiana plant communities. The positive size shows nutrient-rich soils (high organic matter and nitogen) and high density of plants. Whereas the negative size of axes 2 shows nutrient-poor soil, high Plant Utilization Rate (total and for Festuca hystrix), high chromaticity and soil bulk density (Figure 1-2).
The positive size of axis 1 and the negative size of axis 2 suggest plant communities that have been subjected to significant grazing disturbance (Figure 1). Moreover, the presence of Minuartia hybrida an annual plants also suggest this hypothesis. A significant grazing alteration could lead to a increase in soil bulk density and a depletion of soil nutrients (Pulido et al. 2016; Yand et al. 2023).
Plant acronyms
df_acry = tibble(colnames(community_data_presence_2)[1:28],acry)
kable(df_acry)| colnames(community_data_presence_2)[1:28] | acry |
|---|---|
| Thymus_serpylloides | TS |
| Festuca_hystrix | FH |
| Helianthemum_cinereum | HC |
| Poa_ligulata | PL |
| Koeleria_vallesiana | KV |
| Poa_bulbosa | PB |
| Minuartia_hybrida | MH |
| Pilosella_pseudopilosella | PP |
| Bromus_tectorum | BT |
| Helianthemum_appenninum | HA |
| Cerastium_brachypetalum-subsp-brachypetalum | CB |
| Aegilops_geniculata | AG |
| Arenaria_tetraquetra | AT |
| Festuca_larga/suave | FL |
| Vulpia_unilateralis | VU |
| Avenula_bromoides | AB |
| Sedum_amplexicaule | SAm |
| Arenaria_leptoclados | AL |
| Sedum_acre | SAc |
| Xeranthemum_inapertum | XI |
| Veronica_sp | Vsp |
| Asperula_aristata | AA |
| Eryngium_campestre | EC |
| Ononis_pusilla | OP |
| Ononis_spinosa | OS |
| Silene_legionensis | SL |
| Erinacea_anthyllis | EA |
| Helictotrichon_filifolium | HF |