#indices de diversidad en el tiempo

Con la matriz de datos de aves obtenemos los indices de diversidad de Shannon-Weiner, Simpson y Equitatividad de Pielou

# Cargar datos
datos <- read_csv("Monitoreo Humedales_Aves.csv")
## Rows: 528 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): Sitio, Fecha
## dbl (11): Cisne cuello negro, Garza cuca, Garza grande, Garza chica, Tagua c...
## 
## ℹ 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.
# Transformar datos: crear matriz de especies por sitio
matriz <- datos %>%
  dplyr::select(-Fecha) %>% # Excluir la fecha para agrupar por sitio
  group_by(Sitio) %>% 
  summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>% # Sumar abundancias por sitio
  column_to_rownames("Sitio") %>%
  as.matrix()
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(where(is.numeric), sum, na.rm = TRUE)`.
## ℹ In group 1: `Sitio = "Humedal Angachilla"`.
## 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))
# Calcular índices de diversidad
shannon <- diversity(matriz, index = "shannon")  # Índice de Shannon
simpson <- diversity(matriz, index = "simpson")  # Índice de Simpson (1-D)
pielou <- shannon / log(specnumber(matriz))      # Índice de equitatividad de Pielou

# Crear un dataframe con los resultados
indices_diversidad <- tibble(
  Sitio = rownames(matriz),
  Shannon = shannon,
  Simpson = simpson,
  Pielou = pielou
)

# Mostrar resultados
kbl(indices_diversidad) %>%  kable_classic_2(full_width = F)
Sitio Shannon Simpson Pielou
Humedal Angachilla 1.3748043 0.6976914 0.7065097
Humedal Bueras 1.2518484 0.6135553 0.7778172
Humedal Krahmer 0.0000000 1.0000000 0.0000000
Humedal Miraflores 0.0000000 0.0000000 NaN
Humedal Prado Verde 0.0000000 1.0000000 0.0000000
Humedal Santa Inés 0.0000000 1.0000000 0.0000000
Humedal Santa Rosa 0.0000000 1.0000000 0.0000000
Laguna Llancahue 0.6730117 0.4800000 0.9709506
Laguna Los Lotos 0.0000000 0.0000000 NaN
Laguna Los Patos 0.1931592 0.0917405 0.2786698
Laguna Saval 0.0000000 0.0000000 NaN
Laguna de Alivio 0.0000000 0.0000000 NaN
Los Conquistadores 0.0000000 1.0000000 0.0000000
Los Fundadores 0.0000000 1.0000000 0.0000000
Piedra Blanca exterior 0.4634136 0.2265625 0.4218172
Piedra Blanca interior 0.0000000 0.0000000 NaN

Grafico de los indices de diversidad para cada humedal

indices_diversidad %>%
  pivot_longer(cols = -Sitio, names_to = "Indice", values_to = "Valor") %>%
  ggplot(aes(x = Sitio, y = Valor, fill = Indice)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  labs(title = "Índices de diversidad por humedal", x = "Humedal", y = "Valor del Índice") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_bar()`).

Lo malo es que varios sitios no presentan ninguna especie, por lo que los analisis relacionados con las variables ambientales no eran posibles de hacer. Sacando los sitios que no presentan especies y luego uniendo toda la informacion disponible para cada humedal, podemos seguir trabajando. Limpiamos las dos bases de datos, una con datos ambientales y otro solo con los datos de presencia de aves (solo las presentes y los muestreos que presentan valores positivos).

# Cargar los datos
df_env <- read_csv("Humedales_Valdivia_preparado_27ene.csv") %>% clean_names() %>% 
  dplyr::select(-muestra,-hora, -riqueza_sp_aves, -abundancia_sp_aves, -hum,-pp, -tipo)
## Rows: 528 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (2): sitio, tipo
## dbl  (14): muestra, profundidad_cm, temperatura_c, saturacion_o2_percent, co...
## dttm  (1): fecha
## time  (1): hora
## 
## ℹ 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.
df_aves <- read_csv("Monitoreo Humedales_Aves.csv") %>% clean_names() %>% 
  group_by(sitio) %>% 
  dplyr::select(-tagua_chica,- tagua_frente_roja, -taguita) %>% 
dplyr::filter(rowSums(across(where(is.numeric), ~.x!=0))>0) #quedan solo 10 humedales
## Rows: 528 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): Sitio, Fecha
## dbl (11): Cisne cuello negro, Garza cuca, Garza grande, Garza chica, Tagua c...
## 
## ℹ 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.
#para cortar los sitios que no tienen especies en los datos ambietnales
needed<-which(rownames(df_env) %in% rownames(df_aves))    
df_env2<-df_env[needed,]

#y limpar las bases dde datos para el analisis

df_aves2<-df_aves %>%  ungroup %>% dplyr::select(-fecha)
df_aves<- df_aves %>% ungroup %>%  dplyr::select(-sitio, -fecha)

df_env <- df_env2 %>% ungroup %>% dplyr::select(-fecha)

Con esta informacion vamos a hacer un grafico de barras para la abundancia de aves por especie presente por humedal durante todo el periodo

df_aves2 %>% pivot_longer(cols = -sitio, names_to = "Especies", values_to = "valor") %>%
ggplot(aes(x = sitio, y = valor, fill = Especies)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  labs(title = "Abundancia de especies de aves por humedal", x = "Humedal", y = "Numero de individuos") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Hagamos los analisis para las variables ambientales

Primero el PCA corremos el analisis: para explorar patrones en los datos rda es un analisis de redundancia que corre como pca en una matriz simple

#Principal Component Analysis PCA

df_env3 <- df_env %>% dplyr::select(-sitio)
pca1 <- rda(df_env3, scale=F)

#miramos la inercia total (varianza) y cuanto explica cada componente
#unconstrained (segun el numero de matrices hay presentes)
pca1
## Call: rda(X = df_env3, scale = F)
## 
## -- Model Summary --
## 
##               Inertia Rank
## Total           14143     
## Unconstrained   14143    9
## 
## Inertia is variance
## 
## -- Eigenvalues --
## 
## Eigenvalues for unconstrained axes:
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 
## 7740 4141  893  715  588   53   11    1    1
summary(pca1)
## 
## Call:
## rda(X = df_env3, scale = F) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total           14143          1
## Unconstrained   14143          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                             PC1       PC2       PC3       PC4      PC5
## Eigenvalue            7739.7151 4141.1518 892.93762 714.94161 588.2945
## Proportion Explained     0.5472    0.2928   0.06313   0.05055   0.0416
## Cumulative Proportion    0.5472    0.8400   0.90317   0.95372   0.9953
##                             PC6       PC7       PC8       PC9
## Eigenvalue            53.397969 10.905016 1.112e+00 0.8740895
## Proportion Explained   0.003775  0.000771 7.864e-05 0.0000618
## Cumulative Proportion  0.999089  0.999860 9.999e-01 1.0000000
#se puede usar la funcion "PCAsignificance" del paquete BiodiversityR
#esta funcion determina cuantos ejes considerar para interpretacion
#PCAsignificance(pca1) #no funciona BiodiversityR

#uso de grafico para ver
ordiplot(pca1)

ordiplot(pca1, type="t")

#usando graficos de ggvegan: funcion autoplot ver video https://www.youtube.com/watch?v=Tjxgd9FLeYc
#para modificar grafico

autoplot(pca1)

Análisis de Correspondencia Canónica (CCA)

El CCA nos ayuda a entender cómo las variables físico-químicas influyen en la comunidad de aves. Si una especie está cerca de una variable, indica una fuerte asociación. Sitios con condiciones ambientales similares deberían agruparse.

#Análisis de Correspondencia Canónica (CCA)


vare.cca <- cca(df_aves, df_env)
vare.cca
## Call: cca(X = df_aves, Y = df_env)
## 
## -- Model Summary --
## 
##               Inertia Proportion Rank
## Total          4.0356     1.0000     
## Constrained    1.4208     0.3521    7
## Unconstrained  2.6148     0.6479    7
## 
## Inertia is scaled Chi-square
## 
## -- Eigenvalues --
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7 
## 0.8631 0.3363 0.1311 0.0527 0.0202 0.0156 0.0017 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7 
## 0.7635 0.5625 0.5184 0.3169 0.2095 0.1895 0.0546
plot(vare.cca)

#usando graficos de ggvegan: funcion autoplot ver video https://www.youtube.com/watch?v=Tjxgd9FLeYc
#para modificar grafico

autoplot(vare.cca)

AICcPermanova

Este analisis (y el paquete en si) permite hacer una seleccion de modelos con el uso del paquete vegan. De esta forma se pueden analisis multiples modelos que podrian explicar las comunidades presentes (simplemente por la presencia de especies que existe por humedal) El valor de AICc para cada modelo poermite comparar la imporatancia que tienen, se hace una seleccion y se determina un umbral (en el caso de que varios modelos sean seleccionados).

###Calculating AICc
  
# armando todos los modelos posibles : full models
  #saco info incompleta como humedad y precipitacion
  

AllModels2 <- make_models(vars = c("profundidad_cm","temperatura_c","saturacion_o2_percent","concentracion_o2_mg_l",
                                  "conductividad_m_s_cm","tds_mg_l","p_h","orp_m_v","altura_marea_cm"))
## 1 of 9 ready 2025-03-07 15:29:51.640811
## 2 of 9 ready 2025-03-07 15:29:54.684511
## 3 of 9 ready 2025-03-07 15:29:57.769185
## 4 of 9 ready 2025-03-07 15:30:00.869666
## 5 of 9 ready 2025-03-07 15:30:03.99734
## 6 of 9 ready 2025-03-07 15:30:07.224955
## 7 of 9 ready 2025-03-07 15:30:10.476801
## 8 of 9 ready 2025-03-07 15:30:13.650863
## 9 of 9 ready 2025-03-07 15:30:15.292638

multicolinearidad y seleccion de modelos

La multicolinearidad algo que hay que evitar, ya que se presenta cuando dos o mas variables esan altamente relacionadas. Si hay variables altamente relacionadas, todas tienen un efecto sobrela respuesta, aunque no sea directo, por lo que sobreestima el modelo (le otorga una impoartancia mucho mayor porque suma estas importancias)

Chequeamos la multicolinearidad despues de armar los modelos y loego hacemos la seleccion. Para seleccionar lo modelos, usamos una funcion que elige los modelos dentro de ciertos criterios.

##########
# multicolinearidad  
#We can use the filter_vif function to filter out models that have a high degree of collinearity 
#(defined as having a maximum value of VIF of 5 or more):
  
NonColinear2 <- filter_vif(all_forms = AllModels2, env_data = df_env)


##Fittng the models

# After filtering out collinear models, we can fit all the remaining non-collinear models by using the fit_models function:
# sacando la variable "tipo"

  Fitted <- fit_models(
    all_forms = NonColinear2,
    veg_data = df_aves,
    env_data = df_env,
    ncores = 4,
    method = "bray"
  )
  

    #####################
#####    Model Selection
   
    
#    We can use the Fitted2 object generated in the previous section to select models:
      
      Selected <- select_models(Fitted)  
   kbl(Selected) %>%  kable_classic_2(full_width = F)
form max_vif AICc k N profundidad_cm temperatura_c saturacion_o2_percent concentracion_o2_mg_l conductividad_m_s_cm tds_mg_l p_h orp_m_v altura_marea_cm DeltaAICc AICWeight
Distance ~ profundidad_cm + temperatura_c + saturacion_o2_percent + tds_mg_l + p_h + orp_m_v + altura_marea_cm 1.393792 -179.0642 8 143 0.0184964 0.0376101 0.0279560 NA NA 0.1439834 0.0184943 0.0213685 0.0150578 0.0000000 0.2458672
Distance ~ profundidad_cm + temperatura_c + saturacion_o2_percent + tds_mg_l + p_h + orp_m_v 1.349981 -178.2813 7 143 0.0187292 0.0373907 0.0294074 NA NA 0.1423101 0.0205706 0.0228238 NA 0.7828973 0.1662251
Distance ~ profundidad_cm + temperatura_c + concentracion_o2_mg_l + tds_mg_l + p_h + orp_m_v + altura_marea_cm 1.413577 -178.0079 8 143 0.0186223 0.0352246 NA 0.0227389 NA 0.1406135 0.0174155 0.0206580 0.0155538 1.0563647 0.1449819
Distance ~ profundidad_cm + temperatura_c + saturacion_o2_percent + tds_mg_l + orp_m_v + altura_marea_cm 1.301462 -177.5992 7 143 0.0184695 0.0398309 0.0242666 NA NA 0.1394500 NA 0.0238115 0.0171341 1.4650325 0.1181878
Distance ~ temperatura_c + saturacion_o2_percent + tds_mg_l + p_h + orp_m_v + altura_marea_cm 1.393790 -177.5988 7 143 NA 0.0382933 0.0276772 NA NA 0.1426213 0.0184674 0.0252066 0.0152906 1.4654359 0.1181640
Distance ~ profundidad_cm + temperatura_c + saturacion_o2_percent + conductividad_m_s_cm + p_h + orp_m_v + altura_marea_cm 1.332925 -177.4953 8 143 0.0184965 0.0448333 0.0252403 NA 0.1362209 NA 0.0160240 0.0226643 0.0149949 1.5689365 0.1122045
Distance ~ profundidad_cm + temperatura_c + concentracion_o2_mg_l + tds_mg_l + p_h + orp_m_v 1.366479 -177.1491 7 143 0.0187506 0.0349070 NA 0.0236943 NA 0.1388052 0.0196018 0.0224105 NA 1.9151484 0.0943694

Resumen de pesos de las variables con un R2

###    Summary weighted by AICc
#    finally you can do a summarized r squared weighted by AICc using the akaike_adjusted_rsq 
#    function as seen bellow:
      
      Summary <- akaike_adjusted_rsq(Selected) 
kbl(Summary) %>%  kable_classic_2(full_width = F)
Variable Full_Akaike_Adjusted_RSq
profundidad_cm NA
temperatura_c 0.0381264
saturacion_o2_percent NA
concentracion_o2_mg_l NA
conductividad_m_s_cm NA
tds_mg_l 0.1258757
p_h 0.0163214
orp_m_v 0.0224934
altura_marea_cm 0.0114716