install.packages(c("tidyverse", # Manejo de datos y gráficas
"sp", # Manejo de datos espaciales
"sf", # Manejo de datos espaciales
"AICcmodavg", # selección de modelos
"DHARMa", # Diagnóstico de modelos
"MASS", # Extensión para glm nb
"pgirmess", #correlogramas
"ncf", # correlogramas
"gghighlight", # destacar figuras
"showtext", # Añadir texto
"ggtext", # markdown ggplot
"ggspatial", # objetos espaciales ggplot
"leaflet" #mapas interactivos
))Introducción a la autocorrelación espacial
@Gatorco_AP @gpandradep gpandradep@gmail.com
Bienvenidos 🤟🤟
En esta página encontrarán la presentación y los scripts para un pequeño ejercicio donde intentaremos entender de que se trata la autocorrelación espacial y como verificarla en nuestros datos.
Por favor sigue las instrucciones para descargar los paquetes y scripts.
Instrucciones 🧐
Instalación de paquetes
Usen el siguiente código parra instalar los paquetes necesarios. También pueden instalarlos uno por uno.
Descargar datos y script
Puedes descargar el proyecto en R que incluye los datos y los scripts. Para ello usa el siguiente boton y descomprime la carpeta donde lo desees
Presentación 📽️
Código 💻
Librerías
library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(sp) # Classes and Methods for Spatial Data
library(sf) # Simple Features for R
library(gghighlight) # Highlight Lines and Points in 'ggplot2'
library(showtext) # Using Fonts More Easily in R Graphs
library(ggtext) # Improved Text Rendering Support for 'ggplot2'
library(ggspatial) # Spatial Data Framework for ggplot2
library(leaflet) # Create Interactive Web Maps with the JavaScript 'Leaflet' Library
library(AICcmodavg) # Model Selection and Multimodel Inference Based on (Q)AIC(c)
library(DHARMa) # Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models
library(MASS) # Support Functions and Datasets for Venables and Ripley's MASS
# Paquetes para correlogramas
library(pgirmess) # Spatial Analysis and Data Mining for Field Ecologists
library(ncf) # Spatial Covariance Functions
font_add_google("Fira Sans", "firasans") # Familia de tipo de letra
showtext_auto()Cargar datos
CTtable <- read.csv("Data/CTtable.csv") %>% # Base de coordenadas de cámaras
dplyr::select(Station, utm_x, utm_y ) # Seleccionar columnas que usaremos| Station | utm_x | utm_y |
|---|---|---|
| C1T2P1 | 693961 | 2009932 |
| C1T2P11 | 693553 | 2009643 |
| C1T2P21 | 693144 | 2009355 |
| C1T4P1 | 693673 | 2010340 |
| C1T4P11 | 693264 | 2010052 |
| C1T4P21 | 692856 | 2009763 |
events_by_species <- read_csv("Data/surveyReport/events_by_species.csv") %>%
type.convert(n_events= col_double(), # Convertir tipo de columna
n_stations= col_double()) %>% # Convertir tipo de columna
mutate(Sp_n= paste("*",species,"*", " (n= ", n_events, " ) ", sep = "") ) # Columna para los nombres del eje| species | n_events | n_stations | Sp_n |
|---|---|---|---|
| Aves | 704 | 42 | Aves (n= 704 ) |
| Bassariscus astutus | 3 | 2 | Bassariscus astutus (n= 3 ) |
| Canis latrans | 6 | 4 | Canis latrans (n= 6 ) |
| Canis lupus familiaris | 4 | 4 | Canis lupus familiaris (n= 4 ) |
| Capra hircus | 13 | 6 | Capra hircus (n= 13 ) |
| Conepatus leuconotus | 79 | 24 | Conepatus leuconotus (n= 79 ) |
freq_reg <- read.csv("Data/surveyReport/events_by_station2.csv") %>% # Datos de registros por estación
filter(Species == "Odocoileus virginianus") %>% # Filtrar especie
left_join(CTtable, by= "Station") # Unir con la base da cámaras| Station | Species | n_events | utm_x | utm_y |
|---|---|---|---|---|
| C1T2P1 | Odocoileus virginianus | 0 | 693961 | 2009932 |
| C1T2P11 | Odocoileus virginianus | 3 | 693553 | 2009643 |
| C1T2P21 | Odocoileus virginianus | 10 | 693144 | 2009355 |
| C1T4P1 | Odocoileus virginianus | 1 | 693673 | 2010340 |
| C1T4P11 | Odocoileus virginianus | 12 | 693264 | 2010052 |
| C1T4P21 | Odocoileus virginianus | 0 | 692856 | 2009763 |
Gráfico de detecciones
(ndetectionplot <- ggplot(events_by_species, # Datos
aes(x= reorder(Sp_n,n_events), #Ordenar sp
y= n_events))+ # No. eventos
geom_bar(stat= "identity")+ # Geometria
gghighlight(species %in% c("Odocoileus virginianus"))+ # Señalar venados
coord_flip()+ # Girar ejes
labs(title = "Tabla de frecuencia de registro de especies", # Título
y= "Número de registros independientes", # Eje y
x= NULL)+ # Sin eje x
theme_minimal()+ # Tema
theme(text = element_text(family = 'firasans', size = 16), # tipo de letra
plot.title.position = 'plot', # Marco de posición
plot.title = element_text(face = 'bold', # Negrillas
margin = margin(t = 2, r = 0, b = 7, l = 0, unit = "mm"), # Margenes
hjust = 0.5), #posición central
axis.text.y= element_markdown())) # Itálica nombre de especies
Histograma de detecciónes
(hist <- ggplot(freq_reg, aes(x= n_events))+ # Datos
geom_histogram(binwidth = 3, color= "white")+ # Geometria
labs(title = " Histograma número de detecciones",
x= "Número de detecciones",
y= "Frecuencia") +
theme_bw()+ # Tema
theme(text = element_text(family = 'firasans', size = 16), # tipo de letra
plot.title.position = 'plot', # Marco de posición
plot.title = element_text(face = 'bold', # Negrillas
margin = margin(t = 2, r = 0, b = 7, l = 0, unit = "mm"), # Margenes
hjust = 0.5)))
Mapa de registros
UMA <- st_read("Shape/UMA.shp") # Leer el shape de la UMAReading layer `UMA' from data source
`D:\Inecol\Proyectos\Cursos\Independencia_espacial\Shape\UMA.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 688639 ymin: 2002353 xmax: 697655 ymax: 2014011
Projected CRS: WGS 84 / UTM zone 14N# Es importante verificar que este en la misma proyección (CRS)(map <- ggplot(UMA)+ # shape de uma
geom_sf()+ # Graficar objeto espacial
geom_point(data= freq_reg, aes(x=utm_x, y=utm_y, # Agregar puntos
size= n_events, color= n_events), # Tamaño y color
alpha=0.9)+ # Transparencia
scale_size(range = c(1,20))+ # Escala del radio para los puntos
labs(title= "Mapa de número detecciones venado cola blanca",
y= "Lat",
x= "Lon",
size= "Número de \ndetecciones",
color=NULL)+
annotation_scale(location = "bl", width_hint = 0.5, # Barra de escala
line_width=2, tick_height=10) +
annotation_north_arrow(location = "tr", # Flecha de norte
height = unit(2.5, "cm"), width = unit(2.5, "cm"),
pad_x = unit(1, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering)+
theme_bw()+
theme(text=element_text(size=15, family = "special"))+
theme(text = element_text(family = 'firasans', size = 16), # tipo de letra
plot.title.position = 'plot', # Marco de posición
plot.title = element_text(face = 'bold', # Negrillas
margin = margin(t = 2, r = 0, b = 7, l = 0, unit = "mm"), # Margenes
hjust = 0.5)))
Modelando la frecuencia de registro de venados
Vamos a modelar la frecuencia de registros del venado por medio de un GLM Poisson. Uno de los supuestos de los GLMs es que las unidades espaciales son independientes, de manera que vamos a verificarlo. Adicionalmente, vamos a probar algunas covariables que nos ayuden a explicar la distribución espacial de los registros. La no inclusión de una variable importante puede generar patrones en los residuales.
Vamos a cargar una base de datos de covariables
covs.data<-read.csv("Data/covars.csv", header=TRUE) %>% # covariables
dplyr::select(-X) # eliminar columna X
covs.data$Cluster<- as.character(covs.data$Cluster) # cluster es categórica
### Separar las variables númericas y categóricas
cov.num <- covs.data %>%
dplyr::select(where(is.numeric)) %>% # Seleccionar columnas numéricas
scale() %>% # estandarizar
as.data.frame()
cov.fac <- covs.data %>%
dplyr::select(where(is.character)) # Seleccionar variables de caracter
### Unir las bases
sp_glmdata <- data.frame(cov.fac, cov.num) %>% # Unir covariables
right_join(freq_reg, by= "Station")# Unir con la base de rgistros| Station | Cam | Habitat | Cluster | Vertcover_50 | Dcrops | MSAVI | Slope | Effort | Dpop_G | Species | n_events | utm_x | utm_y |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| C1T2P1 | MoultrieA30 | Scrub | 1 | -0.4744602 | 0.6765368 | 0.7316819 | -0.2195550 | -0.2080680 | -0.4260974 | Odocoileus virginianus | 0 | 693961 | 2009932 |
| C1T2P11 | Moultrie | Scrub | 1 | 0.6519157 | 1.4620448 | -0.8666152 | -0.9927301 | -1.4605398 | -0.3167450 | Odocoileus virginianus | 3 | 693553 | 2009643 |
| C1T2P21 | MoultrieA30 | Scrub | 1 | -0.9243340 | 2.2193059 | -0.5838870 | -1.0258694 | -0.0719298 | -0.1805434 | Odocoileus virginianus | 10 | 693144 | 2009355 |
| C1T4P1 | Primos | Columnar | 1 | -0.9243340 | 0.4603495 | -0.1430110 | 1.4183658 | -0.8887592 | -0.1922714 | Odocoileus virginianus | 1 | 693673 | 2010340 |
| C1T4P11 | MoultrieA30 | Scrub | 1 | -1.1509622 | 1.1723154 | 0.2129457 | -0.0597006 | -0.0719298 | -0.0948151 | Odocoileus virginianus | 12 | 693264 | 2010052 |
| C1T4P21 | Primos | Scrub | 1 | 0.6519157 | 1.9052628 | 0.5284530 | -1.0258694 | -0.0719298 | 0.0268169 | Odocoileus virginianus | 0 | 692856 | 2009763 |
Ahora ajustemos algunos modelos
## Modelos lineales generalizados
# sin variables
m0 <- glm(n_events~ 1, # Formula
data = sp_glmdata, # datos
family = "poisson") # tipo de distribución
# 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")
#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" )Seleccionamos el “mejor” con el criterio de información de Akaike
lista_mods <- list(m0, m1, m2, m3, m4, m5)
mod_names <- c("freq~ 1",
"freq~ D_cultivos",
"freq~ MSAVI",
"freq~ Slope",
"freq~ D_poblado",
"freq~ Habitat"
)
AIC <- aictab(lista_mods, # Lista de modelos
modnames = mod_names, # nombres
second.ord = F, # Seleccionar AIC
sort = T) # Ordenar por menor valorknitr::kable(AIC)| Modnames | K | AIC | Delta_AIC | ModelLik | AICWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 6 | freq~ Habitat | 4 | 440.0256 | 0.00000 | 1 | 1 | -216.0128 | 1 |
| 3 | freq~ MSAVI | 2 | 497.6351 | 57.60950 | 0 | 0 | -246.8175 | 1 |
| 2 | freq~ D_cultivos | 2 | 511.1291 | 71.10355 | 0 | 0 | -253.5646 | 1 |
| 5 | freq~ D_poblado | 2 | 512.4230 | 72.39737 | 0 | 0 | -254.2115 | 1 |
| 4 | freq~ Slope | 2 | 512.4941 | 72.46851 | 0 | 0 | -254.2471 | 1 |
| 1 | freq~ 1 | 1 | 513.1820 | 73.15641 | 0 | 0 | -255.5910 | 1 |
Podemos inspeccionar el mejor modelo
# para inspeccionar el modelo usar la función summary
summary(m5)
Call:
glm(formula = n_events ~ Habitat, family = "poisson", data = sp_glmdata)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5702 -2.4202 -0.7741 0.4265 7.9420
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.0745 0.1104 9.730 <2e-16 ***
Habitatcrop 0.0241 0.4229 0.057 0.955
Habitatsalt 2.3267 0.2134 10.904 <2e-16 ***
HabitatScrub 0.1203 0.1462 0.823 0.410
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 373.16 on 63 degrees of freedom
Residual deviance: 294.00 on 60 degrees of freedom
AIC: 440.03
Number of Fisher Scoring iterations: 6Dado 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.
Residuales del mejor modelo
Debido a que es un GLM de familia poisson los residuales no están definidos como en un modelo lineal simple. La función simulateResidual permite obtener una aproximación de residuales escalados para este tipo de modelos.
# Verificamos visualmente que el modelo cumpla los requisitos de la distribución
residuales <- simulateResiduals(fittedModel = m5, plot =F)
plotQQunif(residuales)
Creamos un data set con los residuales y las coordenadas
# Creamos el data.frame para el correlograma
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| res | x | y |
|---|---|---|
| 0.0005904 | 693961 | 2009932 |
| 0.4597318 | 693553 | 2009643 |
| 1.0000000 | 693144 | 2009355 |
| 0.2011827 | 693673 | 2010340 |
| 1.0000000 | 693264 | 2010052 |
| 0.0312027 | 692856 | 2009763 |
Mapa de los residuales
(map_res <- ggplot(UMA)+ # shape de uma
geom_sf()+ # Graficar objeto espacial
geom_point(data= data_resm5, aes(x=x, y=y, # Agregar puntos
size= res, color= res), # Tamaño y color
alpha=0.9)+ # Transparencia
scale_size(range = c(1,15))+ # Escala del radio para los puntos
labs(title= "Mapa de residuales modelo m5",
y= "Lat",
x= "Lon",
size= "Número de \ndetecciones",
color=NULL)+
annotation_scale(location = "bl", width_hint = 0.5, # Barra de escala
line_width=2, tick_height=10) +
annotation_north_arrow(location = "tr", # Flecha de norte
height = unit(2.5, "cm"), width = unit(2.5, "cm"),
pad_x = unit(1, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering)+
theme_bw()+
theme(text=element_text(size=15, family = "special"))+
theme(text = element_text(family = 'firasans', size = 16), # tipo de letra
plot.title.position = 'plot', # Marco de posición
plot.title = element_text(face = 'bold', # Negrillas
margin = margin(t = 2, r = 0, b = 7, l = 0, unit = "mm"), # Margenes
hjust = 0.5)))
Verificando autocorrelación espacial
El test de Moran I básicamente es la extensión de una prueba de correlación de Pearson para las diferentes distancias entre unidades de muestreo. El objetivo es estimar el grado de correlación espacial de una misma variable. Cuando ploteamos el valor de I para las diferentes distancias tenemos un correlograma espacial
Hay diferentes paquetes que permiten generar los correlogramas, cada una tiene sus virtudes y puede ser más conveniente o no, dependiendo de la capácidad de programación del usuario y los datos.
Primero vamos a generar la matriz de distancia. Para asegurar que cada categoría de distancia tenga suficientes pares se recomienda trabajar con 2/3 de la máxima distancia. En este caso la máxima distancia es ~ 8km
Wdist <- data.frame(x= sp_glmdata$utm_x,
y= sp_glmdata$utm_y) %>%
dist() %>%
as.matrix()
maxd <- 2/3* max(Wdist)
maxd[1] 5410.514Paquete pgirmess
pgirmess llama las funciones de spdep con una sintaxis mas amigable, lo que a su vez lo hace menos flexible. Esta aproximación asume que la variable de respuesta se distribuye de manera normal para calcular los p valores
pgirmess_correg <- pgirmess::correlog(coords = cbind(data_resm5$x, # coordenada x
data_resm5$y), # coordenada y
z= data_resm5$res, # residuales
method = "Moran", # Tipo demétodo
nbclass = NULL) # Automáticamente el número de pares| dist.class | coef | p.value | n |
|---|---|---|---|
| 721.7898 | 0.0294411 | 0.2755318 | 356 |
| 1364.7448 | -0.0651344 | 0.7410662 | 368 |
| 2007.6997 | -0.0740384 | 0.8288021 | 526 |
| 2650.6546 | -0.0101852 | 0.4619799 | 560 |
| 3293.6096 | 0.0479433 | 0.1357293 | 548 |
| 3936.5645 | 0.1084853 | 0.0224445 | 496 |
| 4579.5194 | 0.0807245 | 0.0852875 | 432 |
| 5222.4744 | -0.0904677 | 0.8024336 | 308 |
| 5865.4293 | -0.1115091 | 0.8297352 | 264 |
| 6508.3842 | 0.0223829 | 0.3466816 | 126 |
| 7151.3392 | -0.3394159 | 0.9428442 | 38 |
| 7794.2941 | -0.2052650 | 0.5430549 | 10 |
Paquete ncf
La ventaja del paquete ncf es que permite generar correlogramas usando un test de significancia no paramétrico, por medio simulaciones Monte Carlo.
ncf_correg <- ncf::correlog(x=data_resm5$x, # coordenadas en x
y=data_resm5$y, # coordenadas en y
z=data_resm5$res, # variable de interés,
latlon = FALSE, #formato lat lon
na.rm=T, # en caso de NAs
increment = 700, # Distancia mínima de unidades
resamp=500) # Número de remuestreos50 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 Modificamos data frame para la gráfica
corbase <- data.frame(coef=ncf_correg$correlation,
dist.class= ncf_correg$mean.of.class,
p.value= ncf_correg$p,
n= ncf_correg$n )%>%
bind_rows(as.data.frame(pgirmess_correg)) %>%
mutate(package= c(rep("ncf", 12), rep("pgirmess",12)),
p_valor= if_else(p.value<0.025, "significativo", "no-significativo"))| coef | dist.class | p.value | n | package | p_valor | |
|---|---|---|---|---|---|---|
| 1 | 0.0858466 | 527.4675 | 0.1696607 | 83 | ncf | no-significativo |
| 2 | 0.0267224 | 1044.7161 | 0.2574850 | 198 | ncf | no-significativo |
| 3 | -0.0938991 | 1756.1699 | 0.0898204 | 224 | ncf | no-significativo |
| 4 | -0.1348618 | 2438.5924 | 0.0159681 | 318 | ncf | significativo |
| 5 | 0.0710109 | 3154.0392 | 0.0558882 | 310 | ncf | no-significativo |
| 6 | 0.0799903 | 3842.9786 | 0.0499002 | 278 | ncf | no-significativo |
| 7 | 0.0749898 | 4526.0347 | 0.0718563 | 232 | ncf | no-significativo |
| 8 | -0.1415474 | 5243.6160 | 0.0598802 | 162 | ncf | no-significativo |
| 9 | -0.0848323 | 5908.4188 | 0.1916168 | 138 | ncf | no-significativo |
| 10 | 0.0546072 | 6567.7873 | 0.2774451 | 54 | ncf | no-significativo |
| 11 | -0.6117002 | 7212.9125 | 0.0079840 | 15 | ncf | significativo |
| 12 | -0.3202232 | 7906.3521 | 0.2634731 | 4 | ncf | no-significativo |
| …13 | 0.0294411 | 721.7898 | 0.2755318 | 356 | pgirmess | no-significativo |
| …14 | -0.0651344 | 1364.7448 | 0.7410662 | 368 | pgirmess | no-significativo |
| …15 | -0.0740384 | 2007.6997 | 0.8288021 | 526 | pgirmess | no-significativo |
| …16 | -0.0101852 | 2650.6546 | 0.4619799 | 560 | pgirmess | no-significativo |
| …17 | 0.0479433 | 3293.6096 | 0.1357293 | 548 | pgirmess | no-significativo |
| …18 | 0.1084853 | 3936.5645 | 0.0224445 | 496 | pgirmess | significativo |
| …19 | 0.0807245 | 4579.5194 | 0.0852875 | 432 | pgirmess | no-significativo |
| …20 | -0.0904677 | 5222.4744 | 0.8024336 | 308 | pgirmess | no-significativo |
| …21 | -0.1115091 | 5865.4293 | 0.8297352 | 264 | pgirmess | no-significativo |
| …22 | 0.0223829 | 6508.3842 | 0.3466816 | 126 | pgirmess | no-significativo |
| …23 | -0.3394159 | 7151.3392 | 0.9428442 | 38 | pgirmess | no-significativo |
| …24 | -0.2052650 | 7794.2941 | 0.5430549 | 10 | pgirmess | no-significativo |
Generamos nuestro gráfico
(ggplot(corbase ,aes(x=dist.class, y=coef, group= package))+
geom_hline(yintercept = 0, linetype= "dashed")+
geom_line(size=0.9, colour="black")+
geom_point(aes(fill= package, shape=p_valor),
shape= 21, size=5, stroke = 1.5)+
ylim(-1,1)+
scale_x_continuous(breaks = seq(0,8000, by=500),limits = c(500,6000))+
labs(x= "Unidades de distancia (m)", y= " Moran I",
title = " Correlograma de residuales",
colour= "Librería")+
theme_classic()+
theme(text = element_text(family = 'firasans', size = 20), # tipo de letra
plot.title.position = 'plot', # Marco de posición
plot.title = element_text(face = 'bold', # Negrillas
margin = margin(t = 2, r = 0,
b = 7, l = 0,
unit = "mm"), # Margenes
hjust = 0.5),
legend.position = c(0.8, 0.8))) # posición de la legendaWarning: Removed 6 row(s) containing missing values (geom_path).Warning: Removed 6 rows containing missing values (geom_point).
SP lines con ncf
ncf también permite ajustar splines cubicas por medio de bootstrap, en este caso no es necesario definir cortes de pares de distancias y permite obtener intervalos de confianza
spline <- ncf::spline.correlog(x=data_resm5$x, # coordenadas en x
y=data_resm5$y, # coordenadas en y
z=data_resm5$res,
xmax = maxd, # distancia máxima
type = "boot",
resamp = 500)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 Podemos hacer un plot de manera muy sencilla con plot(spline), pero como yo soy terco la quiero hacer con ggplot.
Primero extraigo los valores de la gráfica
spl_dist <- as.data.frame(spline[["boot"]][["boot.summary"]][["predicted"]][["x"]]) %>%
t()
sp_base <- spline[["boot"]][["boot.summary"]][["predicted"]][["y"]][c(2,6,10),] %>% t() %>%
as_tibble() %>%
mutate(dist= spl_dist[,1])| 0.025 | 0.5 | 0.975 | dist |
|---|---|---|---|
| -0.3358756 | 0.2151855 | 0.8007546 | 0.00000 |
| -0.3311834 | 0.2119576 | 0.7849898 | 18.09537 |
| -0.3264913 | 0.2081718 | 0.7691167 | 36.19073 |
| -0.3217992 | 0.2043483 | 0.7532436 | 54.28610 |
| -0.3152821 | 0.2007483 | 0.7373705 | 72.38146 |
| -0.3086233 | 0.1961292 | 0.7215707 | 90.47683 |
Ahora si el gráfico
(splineplot <- ggplot(sp_base, aes(x= dist, y= `0.5`))+
geom_ribbon(aes(ymin= `0.025`, ymax= `0.975`),
fill= "lightgray")+
geom_line(size=1)+
ylim(-1,1)+
scale_x_continuous(breaks = seq(0,8000, by=500),
limits = c(min(sp_base$dist), max(sp_base$dist)),
expand = c(0,0))+
geom_hline(yintercept = 0, linetype=2)+
geom_vline(xintercept = 500, linetype=2, colour= "blue", size=0.8)+
annotate(geom = "curve", x = 900, y = 0.7, xend = 530, yend = 0.5,
curvature = .3, arrow = arrow(length = unit(4, "mm")), size=0.8) +
annotate(geom = "text", x = 930, y = 0.7, size= 5,
label = "Distancia mínima", hjust = "left")+
labs(x= "Unidades de distancia (m)", y= " Moran I",
title = " Spatial line de residuales")+
theme_classic()+
theme(text = element_text(family = 'firasans', size = 20), # tipo de letra
plot.title.position = 'plot', # Marco de posición
plot.title = element_text(face = 'bold', # Negrillas
margin = margin(t = 2, r = 0,
b = 7, l = 0,
unit = "mm"), # Margenes
hjust = 0.5),
legend.position = c(0.8, 0.8))) # posición de la legenda
Lo que podemos concluir con todas las aproximaciones vistas es que parece ser que no hay una autocorrelación espacial fuerte de los residuales de nuestro modelo m5. Esto quiere decir que las desviaciones de la normalidad se deben a otro factor.
Cambiando especificación del modelo
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.
m5bn <- glm.nb(n_events~ Habitat, data = sp_glmdata)testDispersion(m5bn)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.52823, p-value = 0.624
alternative hypothesis: two.sidedresiduales.bn <- simulateResiduals(fittedModel = m5bn, plot =F)
plotQQunif(residuales.bn)
Sorpresa 👽
Esto es el inicio de la aventura
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.
Lecturas recomendadas
