RDA
Daniel Rodríguez García
2026-01-18
Análisis de redundancia
Un análisis de redundancia es una técnica estadística que permite ver las asociaciones que hay en una comunidad y el entorno. Utiliza las abundancias de individuos por especie y trata de relacionarlo con parámetros ambientales. Esto para explicar cual variable ambiental puede explicar la abundancia de especies de un espacio delimitado. Así se puede decidir sobre cómo puede afectar un cambio de estos parámetros en el ecosistema.
Utiliza una lógica en dónde se van a evaluando los parámetros uno por uno. Luego, se elige cual de ellos tiene mayor significancia. Eso se puede hacer ya sea agregando poco a poco variables o también teniendo todas las variables e ir quitando. Así, se observa si el modelo se vuelve más o menos significativo.
Pasos
Para poder correr el modelo se necesitan dos bases de datos combinadas. Una de ellas es la Matriz biológica que incluye la abundancia (cantidad de individuos) de cada especie.
library(readr)
Matriz_amb <- read_csv("Matriz_amb.csv")
Matriz_biol <- read_csv("Matriz_biol.csv")head(Matriz_amb)## # A tibble: 6 × 17
## Punto Temp Hum Elev DistSendero CobDosel AltDosel DensArboles DensArbustos
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 P1 22.5 85 1350 30 70 22 480 1.5
## 2 P2 23 82 1400 25 68 21 450 1.4
## 3 P3 21.8 88 1500 40 75 24 500 1.6
## 4 P4 24 79 1320 35 65 20 440 1.3
## 5 P5 23.5 80 1380 20 72 23 470 1.5
## 6 P6 22.8 87 1450 28 70 22 460 1.4
## # ℹ 8 more variables: RiqVeg <dbl>, BiomHerb <dbl>, Abertura <dbl>,
## # AreaBasal <dbl>, NDVI <dbl>, DistAgua <dbl>, Ruido <dbl>, Radiacion <dbl>head(Matriz_biol)## # A tibble: 6 × 13
## Punto SP1_Turdus SP2_Catharus SP3_Pitangus SP4_Cyanocorax SP5_Thraupis
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 P1 5 2 3 0 4
## 2 P2 3 1 4 1 5
## 3 P3 6 3 2 0 3
## 4 P4 2 0 5 1 6
## 5 P5 4 2 3 2 5
## 6 P6 5 3 4 1 4
## # ℹ 7 more variables: SP6_Cyanerpes <dbl>, SP7_Colibri <dbl>,
## # SP8_Ramphocelus <dbl>, SP9_Momotus <dbl>, SP10_Vireo <dbl>,
## # SP11_Saltator <dbl>, SP12_Cardinalis <dbl>Importante notar que:
Las matrices tienen múltiples columnas. En el caso de la biológica son las especies, en el caso de la ambiental son los parámetros. Cuando aparece 0 en un valor de la biológica, significa que en ese muestreo no se encontró esa especie.
Tanto la biológica como la ambiental tienen la misma cantidad de filas (Eso significa que son cuadradas).
No hay espacios entre los nombres
Transformación de los datos
Lo primero es estandarizar los datos. Esto se hace con varias metodologías estadísticas (sacando raíz cuadrada, logaritmo,…) y la importancia es para disminuir el peso que tienen números muy altos o muy bajos. Esto lo podemos hacer para poder comparar parámetros entre sí, por ejemplo el pH (valores entre 0 y 15) y la cobertura boscosa (valores entre 0 y 100). Lo más utilizado es la estandarización en ambientales tipo “standardize” y de “hellinger” para los biológicos
*Nota: consultar el help(decostand) donde vienen explicaciones de los métodos de estandarización
library(vegan)
library(dplyr)Matriz_biol <- Matriz_biol[,-1] #Elimino la primera columna
Matriz_amb <- Matriz_amb[,-1] #Elimino la primera columna
Matriz_biol_hell <- decostand(Matriz_biol, "total") #Hellinger
Matriz_amb_stand <- decostand(Matriz_amb, "standardize")Modelo
El modelo utiliza la matriz biológica estandarizada como variable dependiente (y) y la ambiental estandarizada como variable x (todas las columnas en conjunto).
modelo_rda<- rda(Matriz_biol_hell ~., data=Matriz_amb_stand) # 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) #significancia del modelo## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = Matriz_biol_hell ~ Temp + Hum + Elev + DistSendero + CobDosel + AltDosel + DensArboles + DensArbustos + RiqVeg + BiomHerb + Abertura + AreaBasal + NDVI + DistAgua + Ruido + Radiacion, data = Matriz_amb_stand)
## Df Variance F Pr(>F)
## Model 16 0.0095854 1.0212 0.566
## Residual 3 0.0017600El resultado muestra que no hay significancia estadística (Pr). También se encontró una baja varianza explicado por el ambiente (Variance). En este caso es menor que 1%. Eso quiere decir que los parámetros ambientales medidos no pueden explicar la composición de aves. Por ello, las aves no están respondiendo linealmente al gradiente ambiental o tal vez no se ven tan afectadas. También puede significar que el espacio medido (paisaje) es sumamente heterogéneo.
Afinación del modelo
Esto sucede cuando se toman un número grande de variables y pocos sitios de muestreo. Sin embargo, se puede reducir la cantidad de dichos parámetros escogiendo los que aportan más al modelo.
Puedo determinar la asociación (R2) de las dos matrices. Aquí puedo decidir cuáles son las variables que escojo para volver a generar el modelo, pero con menos variables. Este análisis siguiente me dice cuáles de los parámetros individualmente me aportan más al modelo.
mod.sel <- ordiR2step(rda(Matriz_biol_hell ~ 1, data=Matriz_amb_stand), # Corre el modelo sin utilizar columnas de la ambiental
rda(Matriz_biol_hell ~ ., data=Matriz_amb_stand), # 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_biol_hell ~ 1
##
## R2.adjusted
## + Temp 0.05137913
## + Radiacion 0.04835222
## + Abertura 0.04353432
## + Hum 0.02592560
## + CobDosel 0.02405079
## + DensArbustos 0.02106557
## + Elev 0.02085227
## + AreaBasal 0.02031515
## + BiomHerb 0.01905446
## + NDVI 0.01832043
## + DistSendero 0.01752951
## <All variables> 0.01752281
## + DensArboles 0.01627767
## + Ruido 0.01356588
## + AltDosel 0.01138445
## <none> 0.00000000
## + DistAgua -0.02532454
## + RiqVeg -0.02709746set.seed(1234)
rda1<- rda(Matriz_biol_hell ~1, data=Matriz_amb_stand)
rda2<- rda(Matriz_biol_hell ~., data=Matriz_amb_stand)
ordistep(rda1, scop= formula(rda2), direction= "forward", pstep=1000)##
## Start: Matriz_biol_hell ~ 1
##
## Df AIC F Pr(>F)
## + Abertura 1 -88.576 1.8648 0.040 *
## + Temp 1 -88.741 2.0291 0.120
## + Radiacion 1 -88.677 1.9654 0.120
## + Hum 1 -88.211 1.5057 0.175
## + AreaBasal 1 -88.097 1.3940 0.190
## + CobDosel 1 -88.173 1.4682 0.195
## + DensArbustos 1 -88.112 1.4089 0.230
## + NDVI 1 -88.056 1.3546 0.240
## + DistSendero 1 -88.040 1.3390 0.240
## + DensArboles 1 -88.014 1.3144 0.245
## + Elev 1 -88.107 1.4046 0.250
## + BiomHerb 1 -88.071 1.3691 0.265
## + Ruido 1 -87.959 1.2613 0.270
## + AltDosel 1 -87.915 1.2188 0.300
## + DistAgua 1 -87.186 0.5307 0.770
## + RiqVeg 1 -87.151 0.4987 0.825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: Matriz_biol_hell ~ Abertura
##
## Df AIC F Pr(>F)
## + BiomHerb 1 -88.663 1.8700 0.115
## + RiqVeg 1 -88.397 1.6206 0.120
## + AreaBasal 1 -88.653 1.8601 0.135
## + Elev 1 -88.436 1.6567 0.155
## + DistSendero 1 -88.347 1.5734 0.180
## + DistAgua 1 -88.004 1.2583 0.245
## + Temp 1 -87.840 1.1090 0.245
## + AltDosel 1 -87.887 1.1516 0.290
## + DensArbustos 1 -87.677 0.9620 0.420
## + NDVI 1 -87.653 0.9403 0.440
## + DensArboles 1 -87.570 0.8661 0.455
## + CobDosel 1 -87.614 0.9053 0.500
## + Radiacion 1 -87.498 0.8021 0.535
## + Ruido 1 -87.336 0.6587 0.690
## + Hum 1 -87.187 0.5274 0.775## Call: rda(formula = Matriz_biol_hell ~ Abertura, data =
## Matriz_amb_stand)
##
## Inertia Proportion Rank
## Total 0.011345 1.000000
## Constrained 0.001065 0.093875 1
## Unconstrained 0.010280 0.906125 11
## Inertia is variance
##
## Eigenvalues for constrained axes:
## RDA1
## 0.001065
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.004437 0.002108 0.001275 0.000716 0.000599 0.000383 0.000344 0.000205
## PC9 PC10 PC11
## 0.000113 0.000075 0.000024Este resultado me dice cómo va aumentando la significancia del modelo entre más vamos agregando parámetros. El modelo comienza sin parámetros, luego se agrega “Abertura” y tiene un valor de significancia (0.40). Luego, si le agrego otro “temp” tengo una significancia más baja (0.12). Esto signfica que el modelo va en dirección forward. Así puedo decir que la abertura es la única variable que puede explicar las diferencias entre la abundancia de especies.
Puede ser que la composición de especies sea explicada por varias variables. Pero en este caso, sólo fue una. Ya que si se agregan otras variables como temp, el valor de p en lugar de bajar, sube (p = 0.12). Recordar que entre más bajo es el p-value es porque la significancia es mayor.
Ya con esto decidido, puedo decir que el parámetro abertura es la que explica la composición de especies de aves del área estudiada. Es decir, que si yo quiero hacer algún cambio, es probable que si yo aumente o disminuya la cantidad de abertura, esto pueda afectar la comunidad de aves de la zona.
Visualización
Se puede realizar un plot que permita visualizar qué variables son más asociadas a ciertas especies.
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)