Parte del código para el ajuste de modelos GLS está basado en código de Jérôme Guélat (https://rpubs.com/jguelat) y Marcus W. Beck (https://beckmw.wordpress.com/)

Cargar los paquetes necesarios:

library(gstat)
library(sp)
library(nlme) 
library(letsR)
library(classInt)
library(RColorBrewer)

En los ejemplos siguientes vamos a usar datos de las aves de México. Datos usados en Villalobos et al. 2013. Biological Conservation. (ya debe estar en su WD)

#para ver los datos "originales", cargar el shapefile
aves.endemic <- readShapePoly("spatialData/aves_endemics7887.shp")
#definir la proyección (ya debemos conocerla)
projection(aves.endemic) <- "+proj=longlat +ellps=clrk66 +no_defs"

#cargar un shapefile de México
mexico <- readShapePoly("spatialData/destdv1gw/destdv1gw.shp")
projection(mexico) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#"guardar" la referencia del objeto "mexico"geográfica en otro objeto, que será usada más adelante
crs_mexico <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#proyectar los dos shapefiles (aves endémicas y México) en la misma proyección
aves.geo <- spTransform(aves.endemic, proj4string(mexico))

Graficar el mapa de riqueza de las aves endémicas de México

#definir los colores y la escala
colours <- brewer.pal(9, "Blues")
brks<-classIntervals(aves.geo$spp_rich, n=9, style="quantile")
brks<- brks$brks
plot(aves.geo, col=colours[findInterval(aves.geo$spp_rich, brks, all.inside=TRUE)], axes=T, border=NA)

legend(x=-95, y=35, legend=leglabs(brks), fill=colours, bty="n")

Cargar datos de “comunidades” (“” porque fueron muestreadas aleatoriamente a partir del shapefile de aves endémicas mencionado arriba).

#Cargar el objeto de R proporcionado
load("spatialData/aves_comunidades.Rdata")

Graficar las “comunidades” en el mapa de México

#Checar que solo queden aquellas "comunidades" dentro del país
plot(mexico)
plot(aves_comunidades,add=T,col="blue",pch=19)

Cargar datos ambientales para México (datos de Cuervo-Robayo A. P., O. Téllez-Valdés, M. Gómez, C. Venegas-Barrera, J. Manjarrez y E. Martínez-Meyer. (2013). An update of high-resolution monthly climate surfaces for Mexico. International Journal of Climatology. Doi: 10.1002/joc.3848. http://idrisi.uaemex.mx/ligas/geogadatos?id=8)

#temperatura promedio anual
mex.temp <- raster("spatialData/bio1.asc")
#estacionalidad de temperatura
mex.temp.sz <- raster("spatialData/bio4.asc")
#precipitación total anual
mex.pp <- raster("spatialData/bio12.asc")
#estacionalidad de precipitación
mex.pp.sz <- raster("spatialData/bio15.asc")

#asignar un sistema de coordenadas (el mismo que el mapa de México)
crs(mex.temp) <- crs_mexico
crs(mex.temp.sz) <- crs_mexico
crs(mex.pp) <- crs_mexico
crs(mex.pp.sz) <- crs_mexico

Obtener los valores ambientales correspondientes a las “comunidades” de aves

aves_coms_temp <- as.data.frame(extract(mex.temp, aves_comunidades))
aves_coms_temp_sz <- as.data.frame(extract(mex.temp.sz, aves_comunidades))

aves_coms_pp <- as.data.frame(extract(mex.pp, aves_comunidades))
aves_coms_pp_sz <- as.data.frame(extract(mex.pp.sz, aves_comunidades))

Convertir los datos en objetos espaciales, preparando un único conjunto de datos que incluya la riqueza de especies y las variables ambientales

#crear objetos faltantes
long <- aves_comunidades$longitude
lat <-   aves_comunidades$latitude
aves_com_rich_coords <- cbind(long,lat)
aves_coms_rich <- aves_comunidades$aves_coms_rich

#Juntar los datos en un data.frame espacial
aves.data <- as.data.frame(cbind(aves_com_rich_coords,aves_coms_rich,aves_coms_pp,aves_coms_pp_sz,aves_coms_temp,aves_coms_temp_sz))
colnames(aves.data) <- c("longitude","latitude","spp_richness","total_precipitation","precipitation_seazonality","mean_annual_temperature","temperature_seazonality")
#Checar que no haya NAs, si los hay, eliminarlos
aves.data.comp <- na.omit(aves.data)
#Convertir el data.frame en un objeto "SpatialPointsDataFrame"
aves.sp <- aves.data.comp
coordinates(aves.sp) = ~ longitude + latitude
crs(aves.sp) <- crs_mexico
aves.spt <- spTransform(aves.sp, crs_mexico)

OLS

Ajustar un modelo lineal tradicional (OLS) entre riqueza de aves endémicas y ambiente (temperatura, precipitación y sus estacionalidades)

aves.ols <- lm(spp_richness ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, data=aves.spt)

¿Cómo se ve? ¿Es un buen modelo?

Vamos a ver los residuos y checar si existe autocorrelación espacial

#Extraer los residuos del modelo
rs <- aves.ols$residuals
#Crear un objeto espacial (SpatialPointsDataFrame) para fines practicos
spdata <- data.frame(resid = rs, x = aves.spt@coords[,1], y = aves.spt@coords[,2])
coordinates(spdata) <- c("x", "y")

#Graficar los residuos usando el objeto espacial
bubble(spdata, "resid", col = c("blue", "orange"), main = "Residuals", xlab = "X-coordinates", 
       ylab = "Y-coordinates")

#Generar un correlograma espacial
distmat <- lets.distmat(aves.spt@coords)
ols.correl <- lets.correl(spdata$resid,distmat,10)

#Generar un semi-variograma
ols.var <- variogram(resid~1,data=spdata)
plot(ols.var)

¿Existe autocorrelación espacial en los residuos del modelo lineal?

GLS (no espacial)

Cuando existe autocorrelación espacial en los residuos de un modelos, la mejor solución sería encontrar la (co)variable faltante (la que causa que no capturemos esa estructura y por tanto se genere autocorrelación espacial) e incluirla en el modelo. Sin embargo, esto no siempre es posible y tenemos que buscar otras soluciones. Una de estas soluciones es modelar la autocorrelación espacial de los residuos e incluirla en el modelo linear.

Para esto, usaremos un modelo GLS (generalised least-squares) ajustado con el paquete nlme. Un modelo GLS permite modelar directamente la estructura de covarianza espacial contenida en la matriz de varianza-covarianza. Esto se hace incorporando la estructura de autocorrelación directamente en los residuos del modelo (i.e. los residuos van a quedar espacialmente autocorrelacionados, pero no así el error de nuestro modelo, que finalmente es lo queremos para estimar de manera correcta los paraámetros).

Primero, para comparar los modelos que vamos a generar, necesitamos ajustar nuevamente el modelo linear inicial pero con una estructura GLS. Esto va a permitir que usemos el criterio de información (AIC) para compararlos (no podemos comparar los AICs de modelos ols(lm) contra gls).

Ajustar el modelo original como GLS

aves.gls <- gls(spp_richness ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, data=aves.spt)

Checar los residuos del modelo original pero ajustado GLS

#Extraer los residuos del modelo
rs.gls <- aves.gls$residuals
#Crear un objeto espacial (SpatialPointsDataFrame) para fines practicos
spdata.gls <- data.frame(resid = rs.gls, x = aves.spt@coords[,1], y = aves.spt@coords[,2])
coordinates(spdata.gls) <- c("x", "y")

#Graficar los residuos usando el objeto espacial
bubble(spdata.gls, "resid", col = c("blue", "orange"), main = "Residuals", xlab = "X-coordinates", 
       ylab = "Y-coordinates")

#Generar un correlograma espacial
gls.resid <- as.vector(spdata.gls$resid)
gls.correl <- lets.correl(gls.resid,distmat,10)

#Generar semi-variograma
gls.var <- variogram(resid~1,data=spdata.gls)
plot(gls.var)

gls.variocloud <- variogram(resid ~ 1, spdata.gls, cloud = TRUE)
plot(gls.variocloud)

GLS espacial

Ahora, ajustemos modelos GLS que incorporan la estructura de autocorrelación en los residuos. Dicha estructura de correlación se introduce como un argumento de la función de ajuste del modelo (la cual será usada para generar la matriz de varianza-covarianza). Recordemos que existen varias estructuras de correlación espacial y tenemos que escoger la que mejor se ajuste/modele los residuos del modelo. Para esto, vamos a ajustar diferentes modelos GLS para cada estructura y los vamos a comparar a través del criterio de información de Akaike (AIC). Noten que necesitamos las coordenadas de los sitios, pero esto ya lo tenemos en nuestro objeto espacial.

#GLS con correlación exponencial
gls.exp <- gls(spp_richness ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, correlation = corExp(form = ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, nugget = TRUE), data=aves.spt)

#GLS con correlación Gaussiana
gls.gaus <- gls(spp_richness ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, correlation = corGaus(form = ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, nugget = TRUE), data=aves.spt)

#GLS con correlacón esférica
gls.spher <- gls(spp_richness ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, correlation = corSpher(form = ~ mean_annual_temperature + temperature_seazonality + total_precipitation + precipitation_seazonality, nugget = TRUE), data=aves.spt)

Comparar los modelos por medio de sus AICs

AIC(aves.gls, gls.exp, gls.gaus, gls.spher)

¿Qué modelo es el “mejor”? Para corroborar:

anova(aves.gls,gls.gaus)

Ahora, vamos a checar los coeficientes del “mejor” modelo y compararlos con el modelo inicial

summary(aves.gls)

summary(gls.gaus)

¿Y cómo se ven los residuos del nuevo modelo?

#Extraer los residuos del modelo
rs.gls.gaus <- gls.gaus$residuals
#Crear un objeto espacial (SpatialPointsDataFrame) para fines practicos
spdata.gaus <- data.frame(resid = rs.gls.gaus, x = aves.spt@coords[,1], y = aves.spt@coords[,2])
coordinates(spdata.gaus) <- c("x", "y")

#Graficar los residuos usando el objeto espacial
bubble(spdata.gaus, "resid", col = c("blue", "orange"), main = "Residuals", xlab = "X-coordinates",
       ylab = "Y-coordinates")

#Generar un correlograma espacial
gaus.resid <- as.vector(spdata.gaus$resid)
gaus.correl <- lets.correl(gaus.resid,distmat,10)

#Generar semi-variograma
gaus.var <- variogram(resid~1,data=spdata.gaus)
plot(gaus.var)

gaus.variocloud <- variogram(resid ~ 1, spdata.gaus, cloud = TRUE)
plot(gaus.variocloud)

Listo! Y ahora, ¿qué podemos decir de los modelos GLS? ¿y de estos modelos para nuestros datos particulares?