Cargar los siguientes paquetes:

library(raster)
library(maptools)
library(rgeos)
library(spgwr)
library(classInt)
library(RColorBrewer)

Vamos a usar los mismos datos y objetos generados que en el ejercicio anterior (GLS; Villalobos et al. 2013. Biological Conservation). Ya deben estar en su WD (i.e. cargar el “workspace” del ejercicio anterior).

load("~gls.RData")

OLS (ordinary least-squares)

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

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

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

Ajustar un modelo “local” con las mismas características que el anterior (“global”), considerando celdas de 0.5º

Primero, crear un raster con las características deseadas

r <- raster(mexico, res=0.5)
r <- rasterize(mexico, r)
#obtener las coordenadas de las celdas (i.e. centroide)
xy <- xyFromCell(r, cell=1:length(values(r)))
#obtener las coordenadas de las "comunidades"
crds <- aves.spt@coords

Crear una función para ajustar el modelo lineal y su aplicación para cada celda

#función general
modelo.local <- function(d)  {
  m <- lm(spp_richness ~ mean_annual_temperature+temperature_seazonality+total_precipitation+precipitation_seazonality, data=d)
  summary(m)$r.squared
}

#"Loop" para iterar sobre cada celda (que tenga datos disponibles)
local.r2s <- list()
for (i in 1:nrow(xy)) {
  d <- sqrt((xy[i,1]-crds[,1])^2 + (xy[i,2]-crds[,2])^2)
  j <- which(d < 1)
  if (length(j) > 5) {
    d <- aves.spt[j,]
    local.r2s[[i]] <- modelo.local(d)
  } else {
    local.r2s[[i]] <- NA
  }
}

Obtener y graficar los resultados de los modelos locales

local.r2 <- sapply(local.r2s, function(x) x)
r.lr2 <- setValues(r, local.r2)
plot(r.lr2)
plot(mexico, add=T)

Y ahora, ¿cómo se ve? ¿existe “no-estacionaridad”?

GWR (geographically weighted regression)

Obtener el “bandwidth” del kernel espacial mediante validación cruzada (“cross-validation”)

aves.bw <- gwr.sel(spp_richness ~ mean_annual_temperature+temperature_seazonality+total_precipitation+precipitation_seazonality, data=aves.spt, gweight = gwr.Gauss, verbose = FALSE)

NOTA: función de ponderación espacial es de tipo Gaussiana. Para este ejemplo, dada la dispersión espacial de los datos, no es posible ajustar un bandwidth adaptativo.

Ajustar el modelo GWR

aves.gwr <- gwr(spp_richness ~ mean_annual_temperature+temperature_seazonality+total_precipitation+precipitation_seazonality, data=aves.spt, bandwidth = aves.bw, gweight = gwr.Gauss,hatmatrix=T)

Checar los valores de R2 locales

#valor de R2 promedio
mean(aves.gwr$SDF$localR2)
#mínimo R2?
min(aves.gwr$SDF$localR2)
#máximo R2?
max(aves.gwr$SDF$localR2)
#desvio estándar
sd(aves.gwr$SDF$localR2)

Con esto, ¿qué podemos decir acerca de la “no-estacionaridad”?

Podemos hacer una ANOVA entre los modelos global (OLS) y local (GWR)

#ANOVA "común" (aunque la función es particular para un objeto tipo "gwr")
anova(aves.gwr)
#ANOVA test basado en el libro de GWR (Fotheringham et al. 2002)
BFC02.gwr.test(aves.gwr)

Mejor si lo graficamos!

#Primera forma, con spplot()
aves.gwr.dat<-aves.gwr$SDF
cols<-brewer.pal(n=9, name="RdBu")

spplot(aves.gwr.dat, "localR2", at=quantile(aves.gwr$SDF$localR2), col.regions=cols, main="local R^2")

#Segunda forma, rasterizando los resultados
aves.coms.gwr <- rasterize(aves.gwr$SDF,r,field="localR2")
plot(aves.coms.gwr)
plot(mexico,add=T)

Y ahora, ¿qué tal se ve? ¿En qué regiones el modelo global sería “malo”?

*NOTA: Guarden el “espacio de trabajo” (generar archivo .RData), pues será utilizado en el próximo ejercicio.