class: inverse, center, bottom background-image: url("img/initial.jpg") background-position: 50% 50% background-size: cover #Explorando la independencia espacial ## 🦌 para datos de fototrampeo 🦨 ### Gabriel Andrade Ponce; gpandradep@gmail.com --- # ¿Autocorrelación espacial (ACS)? En cualquier análisis de información espacial la ACS puede influir sobre nuestros resultados o inferencias. **¿Qué es la ACS?**: Correlación de una variable con ella misma, dada cierta distancia espacial. [(Fortin y Dale 2005)](https://doi.org/10.1017/CBO9780511542039). -- <img src="img/fig1.png" width="50%" style="display: block; margin: auto;" /> --- class: inverse, center ### Primera ley de la grografía <img src="img/law.png" width="50%" style="display: block; margin: auto;" /> *Tobler W., (1970) "A computer movie simulating urban growth in the Detroit region". Economic Geography, 46(Supplement): 234–240. * --- ## ¿ Qué puede causar ACS? ### Los factores más comúnes son [(Dorman *et al.* 2007)](https://onlinelibrary.wiley.com/doi/10.1111/j.2007.0906-7590.05171.x): .pull-left[ 1- **Procesos biológicos**: especiación, dispersión interacciones ecológicas entre otras, son fenómenos relacionados con el espacio. 2- **Especificación del muestreo**: Distancia de las unidades de muestreo respecto al movimiento de la especie (resolución o grano). 3- **Especificación del modelo**: Relaciones no lineales o modelos que no incluyen una variable ambiental determinante que causa la estructura espacial de la variable de interés.] .pull-right[ <img src="https://www.mdpi.com/entropy/entropy-22-00440/article_deploy/html/images/entropy-22-00440-g001.png" width="80%" style="display: block; margin: auto;" /> ] --- class: inverse, center ## La ACS es ¿mala o buena? <img src="https://c.tenor.com/EHVGbpNqP2oAAAAd/todo-depende-del-cristal-del-cual-se-mira-tadeo.gif" width="40%" /> --- class: inverse, center ## La ACS es ¿ mala o buena ? .pull-left[ ### Buena Si tu pregunta ecológica se relaciona con el espacio explicitamente. Entonces la autocorrelación espacial ayudará a informar sobre cómo ocurren los procesos ecológicos en el espacio. <img src="img/SCR.png" width="50%" style="display: block; margin: auto;" /> ] .pull-right[ ### No tan buena Es un gran problema para la prueba estadística de hipótesis y para las predicciones de las mismas, porque viola el supuesto (de casi todas las pruebas) de independencia. Lo que nos puede llevar a cometer errores **.white[Tipo I]** o incluso invertir la relación de la pendiente en algunos análisis. <img src="img/erroresTipoIyII.png" width="30%" /> ] --- ## Detectar y cuantificar la ACS El primer paso antes de empezar con pruebas o modelos más complejos para lidiar con la ACS, es identificar si es en efecto un problema. .pull-left[ ### Nota - 1- Si se realizan pruebas que asumen independencia, pero no son regresiones. La ACS debe verificarse en los datos "crudos". Ej. pruebas de t, patrones de actividad, entre otras. - 2- En regresiones donde se modela el efecto de posibles variables ambientales, la ACS se debe verificar en los residuales. ] .pull-right[ ### Procedimientos Existen diversos procedimientos, pero los más comúnes son la **I de Moran** o correlograma de Moran y los **semi-variogramas**. Para este ejercicio usaremos los correlogramas de Moran, pero recomiendo que exploren las ventajas que puede ofrecer un semi-variograma. ] --- ### Correlograma de I de Moran `$$I= \frac{n}{\sum_{i} \sum_{j}w_{ij}} \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{i,j}(X_{i}- \overline{X})(X_{j}- \overline{X})} {(X_{i}- \overline{X})^2}$$` donde `\(n\)` es el número de unidades `\(i\)` y `\(j\)`; `\(X\)` es la variable de interés; `\(\overline{X}\)` es la media de `\(X\)`; y `\(w_{íj}\)` es la matriz de pesos espaciales. El valor de `\(I\)` puede tomar valores de **1** (autocorrelación positiva), **-1** (autocorrelación positiva) o **0** (distribución aleatoria). .pull-left[ <img src="img/correl_ejem.png" width="700%" height="100%" style="display: block; margin: auto;" /> ] .pull-right[ <br> #### El correlograma es la inspección gráfica de la ACS a las diferentes distancias de pares de puntos de muestreo.] --- class: inverse, center ## Estudio de caso Tenemos un muestreo de cámaras trampa en agrupamientos, con una distancia mínima de ~ 500m. Cualquier director o fototrampero te diría que a esta distancia no hay independencia espacial. Particularmente, para especies grandes y que se mueven mucho.
--- ### Datos ```r # Tabla de funcionamiento de cámaras trampa con coordenadas CToperation <- read.csv("Data/CTtable.csv") # Usaremos el archivo que creamos en survey report, se llama events_by_station2.csv freq_reg <- read.csv("Data/surveyReport/events_by_station2.csv")%>% # Llamamos el csv filter(Species== "Odocoileus virginianus") %>% # Usaremos solo los datos de venados left_join(CToperation, by= "Station") # Unimos para agregar las coordenadas ``` <table> <thead> <tr> <th style="text-align:right;"> X.x </th> <th style="text-align:left;"> Station </th> <th style="text-align:left;"> Species </th> <th style="text-align:right;"> n_events </th> <th style="text-align:right;"> X.y </th> <th style="text-align:right;"> utm_x </th> <th style="text-align:right;"> utm_y </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:left;"> C1T2P1 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 693961 </td> <td style="text-align:right;"> 2009932 </td> </tr> <tr> <td style="text-align:right;"> 4 </td> <td style="text-align:left;"> C1T2P11 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 693553 </td> <td style="text-align:right;"> 2009643 </td> </tr> <tr> <td style="text-align:right;"> 7 </td> <td style="text-align:left;"> C1T2P21 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 10 </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 693144 </td> <td style="text-align:right;"> 2009355 </td> </tr> <tr> <td style="text-align:right;"> 10 </td> <td style="text-align:left;"> C1T4P1 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 4 </td> <td style="text-align:right;"> 693673 </td> <td style="text-align:right;"> 2010340 </td> </tr> <tr> <td style="text-align:right;"> 13 </td> <td style="text-align:left;"> C1T4P11 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 12 </td> <td style="text-align:right;"> 5 </td> <td style="text-align:right;"> 693264 </td> <td style="text-align:right;"> 2010052 </td> </tr> <tr> <td style="text-align:right;"> 16 </td> <td style="text-align:left;"> C1T4P21 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 6 </td> <td style="text-align:right;"> 692856 </td> <td style="text-align:right;"> 2009763 </td> </tr> </tbody> </table> --- ## Distribución de los registros de 🦌 <img src="index_files/figure-html/unnamed-chunk-11-1.png" width="80%" height="200%" style="display: block; margin: auto;" /> --- ## Modelando la frecuencia de captura del 🦌 ```r sp_glmdata <- read.csv("Data/covs.csv") %>% # Llamamos las covariables que vamos a usar right_join(freq_reg, by= "Station")# Las unimos con nuestra tabla de número de eventos ``` <table> <thead> <tr> <th style="text-align:right;"> X </th> <th style="text-align:left;"> Station </th> <th style="text-align:left;"> Cam </th> <th style="text-align:left;"> Habitat </th> <th style="text-align:right;"> Cluster </th> <th style="text-align:right;"> Vertcover_50 </th> <th style="text-align:right;"> Dcrops </th> <th style="text-align:right;"> MSAVI </th> <th style="text-align:right;"> Slope </th> <th style="text-align:right;"> Effort </th> <th style="text-align:right;"> Dpop_G </th> <th style="text-align:right;"> X.x </th> <th style="text-align:left;"> Species </th> <th style="text-align:right;"> n_events </th> <th style="text-align:right;"> X.y </th> <th style="text-align:right;"> utm_x </th> <th style="text-align:right;"> utm_y </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:left;"> C1T2P1 </td> <td style="text-align:left;"> MoultrieA30 </td> <td style="text-align:left;"> Scrub </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> -0.4744602 </td> <td style="text-align:right;"> 0.6765368 </td> <td style="text-align:right;"> 0.7316819 </td> <td style="text-align:right;"> -0.2195550 </td> <td style="text-align:right;"> -0.2080680 </td> <td style="text-align:right;"> -0.4260974 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 693961 </td> <td style="text-align:right;"> 2009932 </td> </tr> <tr> <td style="text-align:right;"> 2 </td> <td style="text-align:left;"> C1T2P11 </td> <td style="text-align:left;"> Moultrie </td> <td style="text-align:left;"> Scrub </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 0.6519157 </td> <td style="text-align:right;"> 1.4620448 </td> <td style="text-align:right;"> -0.8666152 </td> <td style="text-align:right;"> -0.9927301 </td> <td style="text-align:right;"> -1.4605398 </td> <td style="text-align:right;"> -0.3167450 </td> <td style="text-align:right;"> 4 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 693553 </td> <td style="text-align:right;"> 2009643 </td> </tr> <tr> <td style="text-align:right;"> 3 </td> <td style="text-align:left;"> C1T2P21 </td> <td style="text-align:left;"> MoultrieA30 </td> <td style="text-align:left;"> Scrub </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> -0.9243340 </td> <td style="text-align:right;"> 2.2193059 </td> <td style="text-align:right;"> -0.5838870 </td> <td style="text-align:right;"> -1.0258694 </td> <td style="text-align:right;"> -0.0719298 </td> <td style="text-align:right;"> -0.1805434 </td> <td style="text-align:right;"> 7 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 10 </td> <td style="text-align:right;"> 3 </td> <td style="text-align:right;"> 693144 </td> <td style="text-align:right;"> 2009355 </td> </tr> <tr> <td style="text-align:right;"> 4 </td> <td style="text-align:left;"> C1T4P1 </td> <td style="text-align:left;"> Primos </td> <td style="text-align:left;"> Columnar </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> -0.9243340 </td> <td style="text-align:right;"> 0.4603495 </td> <td style="text-align:right;"> -0.1430110 </td> <td style="text-align:right;"> 1.4183658 </td> <td style="text-align:right;"> -0.8887592 </td> <td style="text-align:right;"> -0.1922714 </td> <td style="text-align:right;"> 10 </td> <td style="text-align:left;"> Odocoileus virginianus </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 4 </td> <td style="text-align:right;"> 693673 </td> <td style="text-align:right;"> 2010340 </td> </tr> </tbody> </table> --- ## Modelando la frecuencia de captura del 🦌 Creamos algunos modelos que reflejan nuestras hipótesis sobre las variables que afectan la frecuencia de captura. ```r # Modelos lineales generalizados simples # sin variables m0 <- glm(n_events~ 1, data = sp_glmdata, family = "poisson") # la frecuencia de registro afectada por la distancia a cultivo m1 <- glm(n_events~ Dcrops, data = sp_glmdata, family = "poisson") # la frecuencia de registro afectada por el verdor de la vegetación m2 <- glm(n_events~ MSAVI, data = sp_glmdata, family = "poisson") # la frecuencia de registro afectada por la pendiente m3 <- glm(n_events~ Slope, data = sp_glmdata, family = "poisson") # la frecuencia de registro afectada por la distancia a poblados m4 <- glm(n_events~ Dpop_G, data = sp_glmdata, family = "poisson") # la frecuencia de registro afectada por el tipo de habitat m5 <- glm(n_events~ Habitat, data = sp_glmdata, family = "poisson" ) ``` --- ## Modelando la frecuencia de captura del 🦌 Según el criterio de información de AIC nuestro mejor modelo es aquel que incluye el hábitat. <table> <thead> <tr> <th style="text-align:left;"> </th> <th style="text-align:left;"> Modnames </th> <th style="text-align:right;"> K </th> <th style="text-align:right;"> AIC </th> <th style="text-align:right;"> Delta_AIC </th> <th style="text-align:right;"> ModelLik </th> <th style="text-align:right;"> AICWt </th> <th style="text-align:right;"> LL </th> <th style="text-align:right;"> Cum.Wt </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> 6 </td> <td style="text-align:left;"> freq~ Habitat </td> <td style="text-align:right;"> 4 </td> <td style="text-align:right;"> 440.026 </td> <td style="text-align:right;"> 0.000 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> -216.013 </td> <td style="text-align:right;"> 1 </td> </tr> <tr> <td style="text-align:left;"> 3 </td> <td style="text-align:left;"> freq~ MSAVI </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 497.635 </td> <td style="text-align:right;"> 57.610 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> -246.818 </td> <td style="text-align:right;"> 1 </td> </tr> <tr> <td style="text-align:left;"> 2 </td> <td style="text-align:left;"> freq~ D_cultivos </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 511.129 </td> <td style="text-align:right;"> 71.104 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> -253.565 </td> <td style="text-align:right;"> 1 </td> </tr> <tr> <td style="text-align:left;"> 5 </td> <td style="text-align:left;"> freq~ D_poblado </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 512.423 </td> <td style="text-align:right;"> 72.397 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> -254.211 </td> <td style="text-align:right;"> 1 </td> </tr> <tr> <td style="text-align:left;"> 4 </td> <td style="text-align:left;"> freq~ Slope </td> <td style="text-align:right;"> 2 </td> <td style="text-align:right;"> 512.494 </td> <td style="text-align:right;"> 72.469 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> -254.247 </td> <td style="text-align:right;"> 1 </td> </tr> <tr> <td style="text-align:left;"> 1 </td> <td style="text-align:left;"> freq~ 1 </td> <td style="text-align:right;"> 1 </td> <td style="text-align:right;"> 513.182 </td> <td style="text-align:right;"> 73.156 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> 0 </td> <td style="text-align:right;"> -255.591 </td> <td style="text-align:right;"> 1 </td> </tr> </tbody> </table> --- ## Procedemos a inspeccionar el modelo Debido al carácter relativo del AIC es necesario verificar que el mejor modelo es un buen modelo. Un mal ajuste puede ser causado por la existencia de autocorrelación en los residuales. ```r # Debido a que para glm poisson los residuales no se definen directamente, #usamos simulateResiduals del paquete DHARMa residuales <- simulateResiduals(fittedModel = m5, plot =F) ``` ```r # Verificamos visualmente que el modelo cumpla los requisitos de la distribución plotQQunif(residuales) ``` <img src="index_files/figure-html/unnamed-chunk-17-1.png" width="70%" height="100%" style="display: block; margin: auto;" /> --- class: inverse, center # ¿Será esta desviación del ajuste causada por la autocorrelación espacial? <img src="https://c.tenor.com/n7Q-cdm8ZLkAAAAd/suspicious-fry-futurama.gif" width="40%" height="40%" style="display: block; margin: auto;" /> --- ## Exploremos si existe autocorrelación espacial ```r # Creamos el data.frame para el análisis data_resm5 <- data.frame(res=residuals(residuales), # Residuales que creamos x= sp_glmdata$utm_x,# coordenadas en x y= sp_glmdata$utm_y) # coordenadas en y ``` ### Ahora usamos la función `correlog`, del paquete **nfc** ```r m5_cor <- correlog(x=data_resm5$x, # coordenadas en x y=data_resm5$y, # coordenadas en y z=data_resm5$res, # variable de interés na.rm=T, # en caso de NAs increment = 400, # Distancia mínima de unidades resamp=500) # número de iteraciones ``` ``` ## 50 of 500 100 of 500 150 of 500 200 of 500 250 of 500 300 of 500 350 of 500 400 of 500 450 of 500 500 of 500 ``` --- class: center ## Correlograma Ahora tenemos nuestro primer correlograma. Pueden hacerlo con la función `plot(m5_cor)` o con el código de ggplot del scrpit. <img src="index_files/figure-html/unnamed-chunk-21-1.png" width="100%" /> --- class: inverse ## La autocorrelación espacial no es el problema Ajustamos otro modelo que asumen una distribución de error *binomial negativa*. Con ello nos damos cuenta que el problema era el tipo de distribución. <img src="index_files/figure-html/unnamed-chunk-22-1.png" width="100%" /> --- class: inverse, center, midle # Esto apenas es el comienzo ###Esto fue un ejercicio sencillo, pero lidiar con la ACS merece realizar diversas lecturas en el tema, conocer los supuestos de las técnicas y aprender a interpretar los resultados. ###En caso de encontrar ACS y dependiendo de los objetivos es recomendable usar herramientas como modelos de mínimos cuadrados generalizados (GLS), modelos mixtos, considerar a las coordenadas como covariables o modelos autorregresivos. NOTA: *la mayoría de estas técnicas asumen error de probabilidad con distribución normal, por lo que en muchos casos los datos deben ser transformados* --- class: center ### Sería casi irresponsable decirte que estas listo para enfrentarte a la ACS en la vida real, así que te recomiendo leer: - 1- [Fox *et al.* 2015. Ecological Statistics: Contemporary theory and application](https://oxford.universitypressscholarship.com/view/10.1093/acprof:oso/9780199672547.001.0001/acprof-9780199672547)] - 2- [Plant 2019. Spatial Data Analysis in Ecology and Agriculture Using R](https://www.routledge.com/Spatial-Data-Analysis-in-Ecology-and-Agriculture-Using-R/Plant/p/book/9780367732325) - 3- [Dorman *et al.* 2007. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review](https://onlinelibrary.wiley.com/doi/10.1111/j.2007.0906-7590.05171.x) - 4- [Kuhn & Dorman 2012. Less than eight (and a half) misconceptions of spatial analysis](https://onlinelibrary.wiley.com/doi/full/10.1111/j.1365-2699.2012.02707.x) <img src="https://3.bp.blogspot.com/-feLSjKOc0nA/VspNAqbdzuI/AAAAAAAABAs/D1w5W-BtOSY/s1600/cat%2Bspeed%2Breader.gif" width="50%" height="50%" /> --- class: inverse background-image: url("img/end.JPG") background-position: 50% 50% background-size: cover ## Gracias