library(readr)
## Warning: package 'readr' was built under R version 4.5.2
library(vegan)
## Warning: package 'vegan' was built under R version 4.5.2
## Warning: package 'permute' was built under R version 4.5.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.2
# Matriz biológica: abundancia porcentual de líquenes
Matriz_liquenes <- read_csv("~/Ecología del Paisaje/Semana 2/Matriz_liquenes.csv")
# Matriz ambiental
Matriz_amb_liq <- read_csv("~/Ecología del Paisaje/Semana 2/Matriz_amb_liq.csv")
head(Matriz_liquenes)
## # A tibble: 6 × 21
## Sitio Sp1 Sp2 Sp3 Sp4 Sp5 Sp6 Sp7 Sp8 Sp9 Sp10 Sp11 Sp12
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 S1 5 0 8 2 6 1 9 4 7 3 10 2
## 2 S2 3 1 7 3 5 0 8 3 6 2 9 3
## 3 S3 6 2 9 1 7 2 7 5 8 4 11 1
## 4 S4 2 1 6 2 4 1 8 4 7 3 8 2
## 5 S5 7 3 8 3 8 3 9 6 9 5 10 3
## 6 S6 8 2 9 4 7 2 7 5 8 4 12 2
## # ℹ 8 more variables: Sp13 <dbl>, Sp14 <dbl>, Sp15 <dbl>, Sp16 <dbl>,
## # Sp17 <dbl>, Sp18 <dbl>, Sp19 <dbl>, Sp20 <dbl>
head(Matriz_amb_liq)
## # A tibble: 6 × 16
## Sitio Temp HumRel Radiacion Elevacion CobDosel Abertura NDVI BiomHerb
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 S1 18.3 72 105 1870 85 12 0.71 230
## 2 S2 17.8 74 110 1890 80 15 0.74 220
## 3 S3 18.5 70 120 1905 78 18 0.78 240
## 4 S4 19.1 68 125 1920 82 10 0.7 210
## 5 S5 17.6 75 130 1950 76 20 0.82 260
## 6 S6 18 71 118 1885 79 14 0.72 240
## # ℹ 7 more variables: DensArboles <dbl>, DensArbustos <dbl>, AreaBasal <dbl>,
## # Rugosidad <dbl>, DistAgua <dbl>, Ruido <dbl>, TempCorteza <dbl>
library(vegan)
library(dplyr)
Matriz_liquenes <- Matriz_liquenes[,-1] #Elimino la primera columna
Matriz_amb_liq <- Matriz_amb_liq[,-1] #Elimino la primera columna
Matriz_liquenes_hell <- decostand(Matriz_liquenes, "total") #Hellinger
Matriz_amb_liq <- decostand(Matriz_amb_liq, "standardize")
modelo_rda<- rda(Matriz_liquenes_hell ~., data=Matriz_amb_liq) # El punto luego de la virgulilla significa que está usando todas las columnas de la Matriz ambiental
r2adj<-RsquareAdj(modelo_rda)$adj.r.squared
anova.cca(modelo_rda, step= 1000)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Matriz_liquenes_hell ~ Temp + HumRel + Radiacion + Elevacion + CobDosel + Abertura + NDVI + BiomHerb + DensArboles + DensArbustos + AreaBasal + Rugosidad + DistAgua + Ruido + TempCorteza, data = Matriz_amb_liq)
## Df Variance F Pr(>F)
## Model 15 0.00067109 1.0255 0.478
## Residual 14 0.00061077
mod.sel <- ordiR2step (rda(Matriz_liquenes_hell ~ 1, data=Matriz_amb_liq), # Corre el modelo sin utilizar columnas de la ambiental
rda(Matriz_liquenes_hell ~ ., data=Matriz_amb_liq), # Corre el modelo utilizando todas las columnas
direction = "both") # Va añadiendo parámetros uno por uno y mide la significancia del modelo en cada paso (forward regression). Luego, agarra todos los parámetros y los va quitando uno por uno (backward regression).
## Step: R2.adj= 0
## Call: Matriz_liquenes_hell ~ 1
##
## R2.adjusted
## + HumRel 0.0372131522
## <All variables> 0.0130231810
## + AreaBasal 0.0051961406
## + TempCorteza 0.0016879330
## + Temp 0.0010094345
## <model> 0.0000000000
## + Rugosidad -0.0003531162
## + DensArbustos -0.0056694274
## + BiomHerb -0.0076732947
## + DensArboles -0.0085270793
## + Abertura -0.0092063597
## + CobDosel -0.0113869066
## + Elevacion -0.0113945506
## + NDVI -0.0140102544
## + Radiacion -0.0179338310
## + Ruido -0.0182783116
## + DistAgua -0.0183333021
Ver la significancia de cada variable (p-value) para el modelo.
Valores positivos significa que aporta al modelo
Valores negativos significan que deben desecharse porque no empeoran el modelo.
set.seed(1234)
rda1<- rda(Matriz_liquenes_hell ~1, data=Matriz_amb_liq)
rda2<- rda(Matriz_liquenes_hell ~., data=Matriz_amb_liq)
ordistep(rda1, scop= formula(rda2), direction= "forward", pstep=1000)
##
## Start: Matriz_liquenes_hell ~ 1
##
## Df AIC F Pr(>F)
## + HumRel 1 -198.99 2.1209 0.020 *
## + AreaBasal 1 -198.01 1.1515 0.245
## + TempCorteza 1 -197.90 1.0490 0.355
## + Temp 1 -197.88 1.0293 0.385
## + Rugosidad 1 -197.84 0.9898 0.420
## + DensArbustos 1 -197.68 0.8365 0.560
## + BiomHerb 1 -197.62 0.7792 0.640
## + DensArboles 1 -197.60 0.7548 0.660
## + Abertura 1 -197.58 0.7355 0.695
## + CobDosel 1 -197.51 0.6735 0.700
## + Elevacion 1 -197.51 0.6733 0.750
## + NDVI 1 -197.44 0.5993 0.810
## + Radiacion 1 -197.32 0.4891 0.895
## + Ruido 1 -197.31 0.4794 0.940
## + DistAgua 1 -197.31 0.4779 0.950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: Matriz_liquenes_hell ~ HumRel
##
## Df AIC F Pr(>F)
## + TempCorteza 1 -198.63 1.5130 0.120
## + Temp 1 -198.65 1.5329 0.155
## + AreaBasal 1 -198.37 1.2736 0.230
## + Rugosidad 1 -198.11 1.0232 0.395
## + DensArboles 1 -197.91 0.8446 0.560
## + DensArbustos 1 -197.89 0.8241 0.580
## + BiomHerb 1 -197.91 0.8373 0.585
## + Abertura 1 -197.89 0.8223 0.625
## + NDVI 1 -197.75 0.6957 0.715
## + CobDosel 1 -197.74 0.6822 0.755
## + Elevacion 1 -197.61 0.5669 0.840
## + Ruido 1 -197.60 0.5533 0.880
## + Radiacion 1 -197.55 0.5079 0.890
## + DistAgua 1 -197.51 0.4706 0.895
##
## Call: rda(formula = Matriz_liquenes_hell ~ HumRel, data = Matriz_amb_liq)
##
## Inertia Proportion Rank
## Total 1.282e-03 1.000e+00
## Constrained 9.026e-05 7.041e-02 1
## Unconstrained 1.192e-03 9.296e-01 19
##
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1
## 9.026e-05
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## 0.00027355 0.00025831 0.00015498 0.00013120 0.00007835 0.00007145 0.00005446
## PC8
## 0.00003879
## (Showing 8 of 19 unconstrained eigenvalues)
plot(modelo_rda, type = "n", scaling = 2)
# Sitios
points(modelo_rda, display = "sites", pch = 21, bg = "gray90", col = "black")
# Especies
text(modelo_rda, display = "species", col = "blue", cex = 0.8)
# Variables ambientales
bp <- scores(modelo_rda, display = "bp", scaling = 2)
arrows(0, 0, bp[,1], bp[,2], col = "red", length = 0.1)
text(bp[,1], bp[,2], labels = rownames(bp), col = "red", pos = 4)