En este capítulo vamos a dar una breve introducción a los principales conceptos y métodos básicos en análisis de series de tiempo y de datos espaciales.

0.1 Introducción al análisis de series de tiempo

Para el análisis de series de tiempo tenemos que tomar un enfoque distinto al que hemos venido observando. No podemos tomar varias observaciones en un momento dado, por lo que solo podemos tener una realización \(\{x_t\}\) de distribuciones conjuntas de una secuencia de variables aleatorias \(\{X_t \}\). De este modo, \(\{x_t\}\) viene a ser la serie de tiempo observada.

Una serie de tiempo \(( \{ x_t \})\) puede ser descompuesta en tres componentes:

  1. Tendencia (\(m_t\))
  2. Componente estacional (\(s_t\))
  3. Remanente (\(\epsilon_t\))

De tal forma que:

\(x_t=m_t+s_t+e_t\)

0.1.0.1 Tendencia (\(m_t\))

Se utiliza para extraer la señal de la serie. Podemos unas usar filtros lineales para hacer esto, por ejemplo, una media móvil.

\(m_t=\sum_{a}^{i=0}\frac{1}{2a+1}x_{t+1}\)

Si \(a=1\), entonces:

\(m_t=\frac{1}{3}(x_{t-1}+x_{t}+x_{t+1})\)

0.1.0.2 Componente estacional \((s_t)\)

Como su nombre lo menciona, nos da información de fluctuaciones repetitivas dentro de cada año.

0.1.0.3 Remanente \((e_t)\)

Para calcular el remanente, simplemente restamos los otros dos componentes de la serie de tiempo original \(x_t\):

\(e_t=x_t-m_t-s_t\)

0.1.1 White noise (ruido blanco)

Una serie de tiempo \(\{w_t\}\) es un ruido blanco discreto si sus componentes son independientes e identicamente distribuidos con media cero. Normalmene asumimos:

\(w_t \sim N(0,q)\)

Donde nos referimos al ruido blanco gaussiano. Esto es lo que esperamos obtener cuando removemos la autocorrelación de la serie (i.e. tendencias, efectos estacionales, etc). Entonces \(\epsilon_t\) será un ruido blanco con las propiedades:

\(\overline{x}=0\)

\[c_k=Cov(e_t,e_{t+k})= \begin{cases} q & \quad \text{si } k=0\\ 0 & \quad \text{si } k \neq 0 \end{cases}\]

\[r_k=Cor(e_t,e_{t+k})= \begin{cases} 1 & \quad \text{si } k=0\\ 0 & \quad \text{si } k \neq 0 \end{cases}\]

0.1.2 Random walk

Es una serie de tiempo $ { x_t }$:

\[x_t=x_{t-1}+w_t\] donde, \(w_t\) es ruido blanco, normalmente \(w_t \sim N(0,q)\)

Un random walk es una serie no estacionaria (i.e. sus propiedades dependen de donde nos ubiquemos en la serie de tiempo). Como también podemos notar, un ruido blanco \(w_t\) es estacionario.

0.1.3 Autocorrelación

0.1.3.1 Función de autocorrelación (ACF)

La autocorrelación es la correlación de una variable con si misma a diferentes lags o saltos \((k)\) de tiempos.

\[r_k=Cor(x_t,x_{t+k})\] El intervalo de confianza (95%) es calculado como:

\[-\frac{1}{n}\pm\frac{2}{\sqrt{n}}\] donde, \(n\) es el número de datos.

0.1.3.2 Función de autocorrelación parcial (PACF)

Mide la correlación lineal de una serie \(\{x_t\}\) y una versión “salteada” (lagged) de ella $ { x_{t+k} }$ con la dependencia lineal de $ { x_{t-1},x_{t-2},…,x_{t-(k-1)} }$ removida.

\[f_k=\begin{cases} Cor(x_1,x_0)=r_1 & \quad \text{si } k=1\\ Cor(x_k-x_k^{k-1},x_0-x_0^{k-1}) & \quad \text{si } k \geq 2 \end{cases}\]

0.1.3.3 Función de correlación cruzada

Empleada para dos o más series distintas, donde se quiere observar el grado de correlación entre ellas a diferentes ‘saltos’.

\[g_k^{xy}=\frac{1}{n}\sum_{t=1}^{n-k}(y_t - \overline{y})(x_{t+k} - \overline{x})\]

Luego:

\[r_k^{xy}=\frac{g_k^{xy}}{\sqrt{SD_xSD_y}}\]

0.1.4 Modelos autoregresivos (AR)

Describe el comportamiento de una serie de tiempo. Tiene un orden \(p\), abreviado como \(AR(p)\):

\[x_t=\phi_1 x_{t-1}+\phi_2 x_{t-2}+...+\phi_p x_{t-p}+w_t\] donde \(w_t\) es un ruido blanco.

0.1.5 Modelos de media móvil (MA)

Tiene un orden \(q\), abreviado como \(MA(q)\), es una suma ponderada del error aleatorio actual más los \(q\) errores más recientes:

\[x_t=w_t+\theta_1 w_{t-1}+\theta_2 w_{t-2}+...+\theta_q w_{t-q}\] donde \(w_t\) es un ruido blanco.

0.1.6 Modelos de media móvil autoregresivos (ARMA)

Es una mezcla de los dos modelos anteriores \(ARMA(p,q)\):

\[x_t=\phi_1 x_{t-1}+\phi_2 x_{t-2}+...+\phi_p x_{t-p}+w_t+ \theta_1 w_{t-1}+\theta_2 w_{t-2}+...+\theta_q w_{t-q}\]

donde \(w_t\) es un ruido blanco.

Ahora vamos a ver todos estos conceptos aprendidos en R:

#Leemos las librerias que vamos a utilizar para este capitulo
require(stats)
require(MARSS)
require(forecast)
require(datasets)
require(nlme)
require(pgirmess)
require(ncf)
require(spdep)
require(geoR)
require(gstat)
require(spatstat)         
require(raster)  

# Leer datos a utilizar
load('data/NHTemp.RData')
Temp = NHTemp
load('data/MLCO2.RData')
CO2 = MLCO2

# Veamos como se estructuran mis datos normalmente:
head(CO2)
##   year month    ppm
## 1 1958     3 315.71
## 2 1958     4 317.45
## 3 1958     5 317.50
## 4 1958     6 317.10
## 5 1958     7 315.86
## 6 1958     8 314.93
tail(CO2)
##     year month    ppm
## 704 2016    10 401.57
## 705 2016    11 403.53
## 706 2016    12 404.48
## 707 2017     1 406.13
## 708 2017     2 406.42
## 709 2017     3 407.18
# Vamos a transformalo en un formato ts
# Especificamos frecuencia (12 meses)
co2 = ts(data = CO2$ppm, frequency = 12, start = c(CO2[1, "year"], 
                                                    CO2[1, "month"]))

# Podemos hacer una figura de este objeto:
plot.ts(co2, ylab = expression(paste("CO"[2], " (ppm)")))

# Ahora para la variable temperatura
temp.ts = ts(data = Temp$Value, frequency = 12, start = c(1880, 1))


## Podemos juntar ambas series de tiempo
# Interseccion:
datI = ts.intersect(co2, temp.ts)

# Union:
datU = ts.union(co2, temp.ts)

# Plotear el objeto unido
plot(datI, main = "", yax.flip = TRUE)

# Usando la data de CO2, podemos descomponerla en los
# componentes: tendencia, estacional, remanente
co2.decomp = decompose(co2)

# Veamos esta descomposicion
plot(co2.decomp, yax.flip = TRUE)

# Creamos una funcion de autocorrelacion
acf(co2, lag.max = 36)

# Como interpretamos estos resultados?
# Serie altamente correlacionada en muchos 'lags' o saltos

# Creamos una funcion de autocorrelacion parcial
pacf(co2, lag.max = 36)

# Como interpretamos esto?
# Aproximacion de la correlacion entre el componente remanente.
# Importante para construir modelos AR MA ARMA (ver mas adelante)

# Veamos correlacion cruzada
# Leemos dos bases de datos: abun linces y ciclo solar
suns = ts.intersect(lynx, sunspot.year)[, "sunspot.year"]
lynx = ts.intersect(lynx, sunspot.year)[, "lynx"]

# Vamos a graficarlas:
plot(cbind(suns, lynx), yax.flip = TRUE)

# Veamos la correlacion cruzada
# x = ciclo solar
# y = lince
ccf(suns, log(lynx), ylab = "Cross-correlation")

# Como interpretamos esto?
# Mirar la correlacion cruzada a -3 y -5
# El numero de linces son relativamente mas bajos
# 3-5 anios despues de un ciclo solar alto 

# El numero de linces son relativamente mas altos 
# 11-13 anios antes de un ciclo solar alto

# Ruido blanco

# Vamos a simular ruido blanco:
set.seed(123)
GWN = rnorm(n = 100, mean = 5, sd = 0.2)

# Grafiquemos
plot.ts(GWN)
abline(h = 5, col = "blue", lty = "dashed")

# Veamos el ACF
acf(GWN, main = "", lag.max = 20)

# Vemos lo que esperabamos?
# No hay correlacion! Si tiene sentido!

# Random walk
# Vamos a simularlo
set.seed(123)
TT = 100 # n puntos a simular
# iniciamos ruido blanco y random walk
xx = ww = rnorm(n = TT, mean = 0, sd = 1)
# Loop:
for (t in 2:TT) {
  xx[t] = xx[t - 1] + ww[t]
}

# Ahora observemos lo simulado:
par(mfrow = c(1, 2))
plot.ts(xx, ylab = expression(italic(x[t])))
acf(xx)

# Tiene sentido el ACF?
# Si tiene sentido la ACF. Tiene una fuerte autocorrelacion

# AR models:

# Simulemos 2 AR time series:
set.seed(456)
# Especificamos el orden en 'order', los coeficientes phi en 'ar' 
# y la sd del ruido blanco:
AR.sm = list(order = c(1, 0, 0), ar = 0.1, sd = 0.1) 
AR.lg = list(order = c(1, 0, 0), ar = 0.9, sd = 0.1)
# Simulamos
AR1.sm = arima.sim(n = 50, model = AR.sm)
AR1.lg = arima.sim(n = 50, model = AR.lg)

# Veamos que obtenemos
par(mfrow = c(1, 2))
ylm = c(min(AR1.sm, AR1.lg), max(AR1.sm, AR1.lg))
plot.ts(AR1.sm, ylim = ylm, ylab = expression(italic(x)[italic(t)]), 
        main = expression(paste(phi, " = 0.1")))
plot.ts(AR1.lg, ylim = ylm, ylab = expression(italic(x)[italic(t)]), 
        main = expression(paste(phi, " = 0.9")))

# Ahora hagamos lo mismo para un orden mayor
set.seed(123)
ARp = c(0.7, 0.2, -0.1, -0.3)
AR.mods = list()
for (p in 1:4) {
  AR.mods[[p]] = arima.sim(n = 10000, list(ar = ARp[1:p]))
}

# Ploteamos
# par(mfrow = c(4, 3))
# for (p in 1:4) {
#   plot.ts(AR.mods[[p]][1:50], ylab = paste("AR(", p, ")", sep = ""))
#   acf(AR.mods[[p]], lag.max = 12)
#   pacf(AR.mods[[p]], lag.max = 12, ylab = "PACF")
# }

# Observemos el PACF
# El orden de proceso AR esta relacionado con
# el numero de autocorrelaciones significativas 
# a diferentes saltos.

# MA models:
# Hagamos lo mismo que hicimos para AR models:
set.seed(123)
MA.sm = list(order = c(0, 0, 1), ma = 0.2, sd = 0.1)
MA.lg = list(order = c(0, 0, 1), ma = 0.8, sd = 0.1)
MA.neg = list(order = c(0, 0, 1), ma = -0.5, sd = 0.1)
# Simulamos
MA1.sm = arima.sim(n = 50, model = MA.sm)
MA1.lg = arima.sim(n = 50, model = MA.lg)
MA1.neg = arima.sim(n = 50, model = MA.neg)


# Plot:
par(mfrow = c(1, 3))
plot.ts(MA1.sm, ylab = expression(italic(x)[italic(t)]), 
        main = expression(paste(theta, " = 0.2")))
plot.ts(MA1.lg, ylab = expression(italic(x)[italic(t)]), 
        main = expression(paste(theta,  " = 0.8")))
plot.ts(MA1.neg, ylab = expression(italic(x)[italic(t)]), 
        main = expression(paste(theta,  " = -0.5")))

# Analicemos el patron 
# No gran cambio: simples combinaciones de ruido blanco

# Ahora para diferentes ordenes:
set.seed(123)
MAq = c(0.7, 0.2, -0.1, -0.3)
MA.mods = list()
for (q in 1:4) {
  MA.mods[[q]] = arima.sim(n = 1000, list(ma = MAq[1:q]))
}

# Plot:
# par(mfrow = c(4, 3))
# for (q in 1:4) {
#   plot.ts(MA.mods[[q]][1:50], ylab = paste("MA(", q, ")", sep = ""))
#   acf(MA.mods[[q]], lag.max = 12)
#   pacf(MA.mods[[q]], lag.max = 12, ylab = "PACF")
# }

# ACF va a cero bastante rapido
# PACF: va a cero lentamente

# ARMA models:

# Vamos a estimar los parametros de un modelo ARMA

# Simulamos una serie de timepo:
set.seed(123)
ARMA22 = list(order = c(2, 0, 2), ar = c(-0.7, 0.2), ma = c(0.7, 0.2))
# la media del proceso
mu = 5
# simulamos usando la funcion arima.sim. Cambia 'n' para ver su efecto
ARMA.sim = arima.sim(n = 100, model = ARMA22) + mu
## estimate parameters
arima(x = ARMA.sim, order = c(2, 0, 2))
## 
## Call:
## arima(x = ARMA.sim, order = c(2, 0, 2))
## 
## Coefficients:
##           ar1     ar2     ma1     ma2  intercept
##       -0.3559  0.1967  0.3118  0.1743     4.9958
## s.e.   0.2809  0.2565  0.2826  0.2428     0.1196
## 
## sigma^2 estimated as 0.8785:  log likelihood = -135.6,  aic = 283.2
# Usamos esta funcion para estimar los parametros adecuados 
# para diferentes ordenes (lo que normalmente queremos hacer)
auto.arima(ARMA.sim, start.p = 0, max.p = 3, start.q = 0, max.q = 3)
## Series: ARMA.sim 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1     ar2    mean
##       -0.0804  0.3723  4.9978
## s.e.   0.0921  0.0938  0.1322
## 
## sigma^2 estimated as 0.919:  log likelihood=-136.3
## AIC=280.61   AICc=281.03   BIC=291.03
# Son buenos los estimados?

# Una aplicacion: --------------------------------
# Importante!: Usarlo cuando implementamos una regresion lineal
# con datos que tienen autocorrelacion temporal

# Leemos datos:
chum = read.csv('data/coastwide_chum_data.csv')

# Seleccionamos solo una parte de la base de datos:
stock1 = chum[chum$stock=='SE-SC AKPen',c('stock','entry.yr','spawners','recruits','ln.rs',
                                         'loc.sst','w.sst.1','w.sst.2','npgo2','npgo','pdo2','pdo','long','lat')]
head(stock1)
##           stock entry.yr spawners recruits       ln.rs     loc.sst    w.sst.1
## 764 SE-SC AKPen     1963   238600   576200  0.88166631  0.04182643 -0.2207373
## 765 SE-SC AKPen     1964   263000   251500 -0.04471104 -1.19076329  0.6549483
## 766 SE-SC AKPen     1965   160800   313100  0.66636127 -0.47777858  1.0118060
## 767 SE-SC AKPen     1966   203300   192000 -0.05718735 -1.02832160 -1.1020202
## 768 SE-SC AKPen     1967   354800   657883  0.61747285  0.59927106 -0.3414790
## 769 SE-SC AKPen     1968   132800   846192  1.85190205  0.20264870 -0.5955610
##         w.sst.2        npgo2       npgo   pdo2    pdo long  lat
## 764  0.21710551  0.164053783 -0.2557230 -0.659 -0.472  161 55.5
## 765  0.83337719 -0.007529911  0.5838306 -0.891 -0.846  161 55.5
## 766 -0.04510709 -0.714296310 -0.5988904 -0.735 -0.936  161 55.5
## 767 -0.72174962 -0.809721169 -0.8297022 -0.572 -0.534  161 55.5
## 768 -0.46852000 -0.679363840 -0.7897401 -0.516 -0.610  161 55.5
## 769 -0.32856791 -0.796891715 -0.5689875 -0.653 -0.422  161 55.5
# Implementamos una regresion lineal:
lm_loc_sst = lm(ln.rs~spawners+loc.sst,data=stock1)

# Observamos checking plots:
opar = par(mfrow=c(2,2)) # panel 2 por 2
plot(lm_loc_sst)

par(opar)

# Vemos si existe algun tipo de correlacion temporal
opar = par(mfrow=c(1,3)) # panel 2 por 2
plot(residuals(lm_loc_sst), type = 'b')
acf(residuals(lm_loc_sst))
pacf(residuals(lm_loc_sst))

par(opar)

# Vamos a tratar de tomar en cuenta esta correlacion temporal
gls_loc.sst = gls(ln.rs~spawners+loc.sst,data=stock1, method='ML')
summary(gls_loc.sst)
## Generalized least squares fit by maximum likelihood
##   Model: ln.rs ~ spawners + loc.sst 
##   Data: stock1 
##        AIC      BIC    logLik
##   50.96795 58.28252 -21.48398
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept)  1.5344406 0.15312740 10.020679   0e+00
## spawners    -0.0000023 0.00000049 -4.640995   0e+00
## loc.sst      0.2455959 0.06674652  3.679532   6e-04
## 
##  Correlation: 
##          (Intr) spwnrs
## spawners -0.923       
## loc.sst   0.262 -0.278
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2719752 -0.7088401  0.1435002  0.5538993  1.7688651 
## 
## Residual standard error: 0.3860086 
## Degrees of freedom: 46 total; 43 residual
gls_loc.sst_AR1 = gls(ln.rs~spawners+loc.sst,data=stock1,
                     correlation = corAR1(form=~entry.yr), method='ML')
summary(gls_loc.sst_AR1)
## Generalized least squares fit by maximum likelihood
##   Model: ln.rs ~ spawners + loc.sst 
##   Data: stock1 
##        AIC      BIC    logLik
##   43.24243 52.38564 -16.62122
## 
## Correlation Structure: AR(1)
##  Formula: ~entry.yr 
##  Parameter estimate(s):
##       Phi 
## 0.4590491 
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept)  1.5049862 0.17477859  8.610816  0.0000
## spawners    -0.0000022 0.00000051 -4.292580  0.0001
## loc.sst      0.1683614 0.05608123  3.002099  0.0045
## 
##  Correlation: 
##          (Intr) spwnrs
## spawners -0.836       
## loc.sst  -0.014  0.023
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4268110 -0.5899631  0.1828114  0.6460358  1.5441313 
## 
## Residual standard error: 0.3899023 
## Degrees of freedom: 46 total; 43 residual
# Cual de los dos modelos es mejor?
AIC(gls_loc.sst_AR1,gls_loc.sst)
##                 df      AIC
## gls_loc.sst_AR1  5 43.24243
## gls_loc.sst      4 50.96795
# Vemos si existe algun tipo de correlacion temporal
opar = par(mfrow=c(1,2)) # panel 2 por 2
pacf(residuals(gls_loc.sst,type='normalized'),main='Normalized residuals GLS_SST')
pacf(residuals(gls_loc.sst_AR1,type='normalized'),main='Normalized residuals GLS_SST_AR1')

par(opar)

0.2 Introducción al análisis espacial

Algunos conceptos importantes en este tema:

  • Ecologia espacial: se refiere al estudio y modelado del rol (o roles) del espacio sobre procesos ecológicos. Tiene aplicaciones en muchos campos como la conservación y manejo de recursos.

  • Escala espacial: importante cuando nos planteamos nuestra pregunta de investigación. Se refiere a la dimensión o dominio espacio-temporal de un proceso o patrón.

0.2.1 Definiciones:

  1. Grano: Unidad espacial más fina de medida para un proceso.

  2. Extensión: Área bajo investigación.

0.2.2 Datos o procesos puntuales

Diferentes tipos de datos se pueden colectar en campo:

  1. Presencia en un área determinada.

  2. Presencia/ausencia.

  3. Algún indicador de abundancia.

  4. Alguna covariable adicional.

Estudiar procesos puntuales nos ayuda a inferir comportamieno de la variable bajo estudio. Esto tiene una influencia importante sobre la ecología de la especie. Análisis del patrón de puntos examinan si hay regularidades en el proceso que representan (son agrupados? aleatorios? uniformes?).

  • Proceso puntual: son procesos aleatorios que resultan en realizaciones (observaciones) de puntos en tiempo y espacio.

  • Proceso puntual de Poisson: puntos son independientemente distribuidos (aleatorios) en el espacio.

0.2.2.1 Modelos nulos (null models)

Usados para evaluar la significancia del patrón observado de puntos. Normalmente son procesos aleatorios con los cuales vamos a comprar las observaciones.

0.2.2.2 Aleatoriedad espacial completa (CSR)

Asume que el número de puntos en una región sigue una distribución tipo Poisson con una media dada. Además, asume que los puntos son una muestra aleatoria independiente, y todos tienen la misma probabilidad de ocurrir. Nos ayudan a determinar si los puntos observados son aleatorios o no al compararlos con este tipo de proceso puntual.

Veamos ahora algunos métodos para estudiar el nivel de agregación de las observaciones.

0.2.2.3 K de Ripley

Calcula el grado de agregación espacial de puntos dentro de un círculo de radio r y lo contrasta con lo experado a partir de una distribución totalmente aleatoria (CSR).

0.2.2.4 Distancia entre vecinos: Función G

Estima la distribución acumulada de las distancias de los vecinos más cercanos a un punto dado. Útil para interpretar probabilidad de encontrar el vecino más cercano a una distancia r.

0.2.2.5 Marcas bivariadas o multivariadas

Usado para dos tipos de puntos (e.g. dos diferentes especies). Marcas son simplemente información acerca de los eventos (puntos), que puede ser categórica (e.g. presencia predador - presa) o continua (e.g. densidad de aves y árboles).

0.2.3 Efecto de borde

Que sucede cuando no tengo suficientes puntos cerca del borde? Algunos resultados serán sesgados. Los métodos anteriormente mencionados pueden incorporar correciones para lidiar con este efecto de borde.

Hasta el momento hemos asumido supuestos importantes en nuestros datos. Por ejemplo, isotropía, que quiere decir que los procesos espaciales no cambian en función de la dirección que se tome en cuenta. Además, estos métodos para estudiar la distribución de puntos asume que los datos provienen de un censo, lo cual dificilmente sucede.

0.2.4 Dependencia espacial y autocorrelación

Vamos a presentar definiciones importantes en este campo. Para esta parte vamos a trabajar con los valores de una variable en el espacio (ya no solamente procesos puntuales).

Primera ley de la geografía: Todo esta relacionado a todo lo demás, pero objetos cercanos lo están más que objetos distante.

  • Estadística espacial: tiene por objetivo cuantificar el patrón espacial y si significancia estadística.

  • Geoestadística: tiene por objetivo cuantificar la varianza espacial y usar esta información para interpolar datos.

  • Dependencia espacial: similaridad de una variable como una función del lugar o distancia geográfica.

Algunas otras definiciones:

  1. Mecanismos endógenos, provienen de la misma población.

  2. Mecanismos exógenos, provienen de factores externos a la población.

  3. Dependencia espacial, resultado de mecanismos endógenos y exógenos.

  4. Autocorrelación espacial, resultado solamente de mecanismos endógenos.

0.2.5 Cuantificación de correlación espacial

0.2.5.1 Correlogramas

Recordemos la correlación de Pearson:

\[r(z_1,z_2)=\frac{\sum_{i=1}^{n}(z_{1i}-\overline{z}_1)(z_{2i}-\overline{z}_2)}{\sqrt{\sum_{i=1}^n(z_{1i}-\overline{z}_1)^2\sum_{i=1}^n(z_{2i}-\overline{z}_2)^2}}=\frac{Cov(z_1,z_2)}{\sqrt{Var(z_1)Var(z_2)}}\]

Podemos extender esto al espacio (índice de Moran):

\[I(d)=\frac{n}{W(d)} \frac{\sum_{i=1}^n \sum_{j=1}^n w_{ij}(d)(z_i-\overline{z})(z_j-\overline{z})}{\sum_{i=1}^n(z_i-\overline{z})^2}\] Algunos puntos importantes:

  • La figura de distancia vs. \(I(d)\) es conocido como correlograma.

  • Valores cercanos a 0 indican ausencia de un patrón espacial.

  • Valores cercanos a 1 indican alta autocorrelación espacial.

  • Se necesitan normalmente más de 20 pares de puntos.

  • Asume isotropía, pero puede considerar anisotropía.

  • Sensible a outliers.

0.2.5.2 Variogramas

Otra forma de estudiar la autocorrelación espacial en una variable en el espacio. En geostadística, se centra en varianzas y covarianzas.

\[Var(z)=\frac{1}{n-1}\sum_{i=1}^n(z_i - \overline{z})^2\] \[Cov(z_1,z_2)=\frac{1}{n-1}\sum_{i=1}^n(z_{1,i} - \overline{z}_1)(z_{2,i} - \overline{z}_2)\] La semivarianza es estimada:

\[\gamma(d)=\frac{1}{2n(d)}\sum{i}^{n(d)}[z(x_i)-z(x_i+d)]^2\] Definiciones importantes:

  • Semivariograma (o variograma) es la representación gráfica de la semivarianza en función a la distancia.

  • Valores pequeños (cerca a 0) indican alta autocorrelación espacial.

  • Variogramas empíricos, provienen de los datos. Se representan simplemente por puntos.

  • Variogramas teóricos, vienen de un modelo que trata de ajustar a lo observado (variograma empírico). Es utilizado para la interpolación.

Existen tres parámetros comunes que son estimados en un variograma teórico:

  1. Nugget: Variabilidad en datos a distancias muy cortas.

  2. Rango: Distancia hasta donde se observa dependencia.

  3. Sill: Nivel de variabilidad que no se atribuye a autocorrelación espacial.

Existen varios tipos de modelos para implementar un variograma teórico, entre los más comunes son:

  • Esférico

  • Exponencial

  • Gausiano

  • Matérn

0.2.5.3 Kriging

Utilizado para inteprolar en una región dad basado en el grado de dependencia espacial. Usa información del semivariograma para realizar la interpolación espacial. El objetivo de esta interpolación es minimizar la varianza de la covarianza espacial en los datos. Normalmente modelado de la siguiente manera:

\[z=f(x_i)+v(x_i)+\varepsilon_i\]

donde, \(f(x_i)\) es una función de tendencia (en el espacio) determinística, \(x_i\) es un punto en el espacio, \(z\) son los valores predichos por el kriging, \(v(x_i)\) decribe la dependencia espacial, \(\varepsilon_i\) describe el error.

Para el kriging ordinario: \[f(x_i)=0\] Normalmente se usa kriging universal: \[f(x_i) \neq 0\]

Pasos importantes para realizar una interpolación espacial mediante kriging:

  1. Tener puntos de muestreo con variable a modelar.

  2. Obtener el variograma teórico (puede ser con un modelo esférico, exponencial, etc.)

  3. Definir el área y grilla de interpolación.

  4. Con esta información, puedo realizar la interpolación mediante el método de kriging.

A continuación vamos a ver algunas aplicaciones usando R:

# Procesos puntuales:
 
# Leemos datos:
cactus = read.csv("data/cactus.csv")
boundary = read.csv("data/cactus_boundaries.csv",header=T)

# Explorar datos
head(cactus)
##   PatchID     East   North       Area chelinidea
## 1       1 403377.3 3285723 0.19634954         26
## 2      10 403373.5 3285681 0.42552872          5
## 3      11 403382.0 3285680 0.27371126          9
## 4      12 403382.5 3285682 0.34714599          1
## 5      13 403385.3 3285681 0.67387162          5
## 6      14 403385.5 3285683 0.07696902          9
head(boundary)
##   Plot   Xmin     Xmax    Ymin    Ymax
## 1   10 403368 403424.6 3285673 3285728
# creamos spatstat objects
ppp.window = owin(xrange=c(boundary$Xmin, boundary$Xmax),
                 yrange=c(boundary$Ymin, boundary$Ymax))
ppp.cactus = ppp(cactus$East, cactus$North, window=ppp.window)

# Graficamos los puntos de muestreo:
par(mfrow = c(1,1))
plot(ppp.cactus)

# Un resumen de informacion
summary(ppp.cactus)
## Planar point pattern:  82 points
## Average intensity 0.0262668 points per square unit
## 
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
## 
## Window: rectangle = [403368, 403424.6] x [3285673, 3285728] units
##                     (56.57 x 55.18 units)
## Window area = 3121.81 square units
# Plots de densidad
plot(density(ppp.cactus))

# Elaborar cuadrantes:
Q = quadratcount(ppp.cactus, nx = 4, ny = 4) 

# Graficamos:
plot(ppp.cactus, cex = 2)
plot(Q, add = TRUE, cex = 1)

# Prueba chi-sq para verificar CSR:
quadrat.test(ppp.cactus, nx = 4, ny = 4, method="Chisq")
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp.cactus
## X2 = 35.463, df = 15, p-value = 0.004223
## alternative hypothesis: two.sided
## 
## Quadrats: 4 by 4 grid of tiles
# Mis datos no son aleatorios basados en esta prueba de hipotesis

#Ripley's K function:
#-----------------------#
Knone = Kest(ppp.cactus, correction="none")

# Graficamos el resultado
plot(Knone)

# Veamos en forma linearizada (mejor interpretacion):
Lnone = Lest(ppp.cactus, correction="none")
plot(Lnone, legend=F)

# Comportamiento de agregacion. r = radio del circulo

# Podemos ver la significancia del comportamiento de agregacion (MCMC)
Lcsr = envelope(ppp.cactus, Lest, nsim=99, correction="trans")
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
# Grafiquemos:
plot(Lcsr, . - r~r, shade=c("hi", "lo"), legend=F)

# Si la linea negra esta fuera de la zona gris, 
# el comportamiento agregativo es significativo.

#G function
#-----------------------#
Gtrans = Gest(ppp.cactus, correction="none")

# Grafiquemos
plot(Gtrans)

# Ver significancia de este indicador
Genv = envelope(ppp.cactus, Gest, nsim=99, rank=1, correction="rs", global=F)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
#plot G with envelope
plot(Genv, shade=c("hi", "lo"), legend=F)

# Interpretacion: este analisis sugiere que las
# distancias a los vecinos mas cercanos son aleatorias
# a pequenas escalas (<2m). Mientras que distancias
# son mas cercanas a las esperadas a escalas mayores.

# Marked point patterns
##############################################
# Ver la relacion entre cactus y un bicho

# Creamos una nueva columna
cactus$CheliPA = as.factor(ifelse(cactus$chelinidea>0, "presence", "absence"))
head(cactus)
##   PatchID     East   North       Area chelinidea  CheliPA
## 1       1 403377.3 3285723 0.19634954         26 presence
## 2      10 403373.5 3285681 0.42552872          5 presence
## 3      11 403382.0 3285680 0.27371126          9 presence
## 4      12 403382.5 3285682 0.34714599          1 presence
## 5      13 403385.3 3285681 0.67387162          5 presence
## 6      14 403385.5 3285683 0.07696902          9 presence
# Creamos un nuevo objeto
ppp.PA = ppp(cactus$East, cactus$North, window=ppp.window, marks=cactus$CheliPA)

# Resumen:
summary(ppp.PA)
## Marked planar point pattern:  82 points
## Average intensity 0.0262668 points per square unit
## 
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
## 
## Multitype:
##          frequency proportion   intensity
## absence         28  0.3414634 0.008969152
## presence        54  0.6585366 0.017297650
## 
## Window: rectangle = [403368, 403424.6] x [3285673, 3285728] units
##                     (56.57 x 55.18 units)
## Window area = 3121.81 square units
# Graficamos
plot(split(ppp.PA))

# Hacemos un analisis solamente considerando el bicho
cheli.data = subset(cactus, chelinidea>0)
ppp.bug = ppp(cheli.data$East, cheli.data$North, window=ppp.window)
Lbug = envelope(ppp.bug, Lest, nsim=99, rank=1, i="presence", global=F)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(Lbug, . - r~r, shade=c("hi", "lo"), legend=F)

# Ahora hacemos el analisis para las dos variables:
Lmulti = envelope(ppp.PA, Lcross, nsim=99, rank=1, i="presence", global=F, simulate=expression(rlabel(ppp.PA)))
## Generating 99 simulations by evaluating expression  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
# Graficamos
plot(Lmulti, . - r~r, shade=c("hi", "lo"), legend=F)

# Ambas especies tienen un comportamiento espacial
# similar (agregativo)

# -----------------------------
# Correlacion espacial y Geoestadistica:
# 

# Leemos los datos
matrix = read.csv('data/cactus_matrix.csv', header=T)
head(matrix)
##   x  y Height
## 1 0  0     35
## 2 0  2     65
## 3 0  4     75
## 4 0  6     60
## 5 0  8     80
## 6 0 10     20
# Hacemos una grafica:
plot(matrix[,"y"] ~ matrix[,"x"],
     pch=21, cex=1.2,
     bg=gray.colors(12)[cut(matrix[,"Height"], breaks = 12)])

# Veamos la distribucion de la variable a estudiar:
hist(matrix[,"Height"], xlab="Vegetation height (cm)",
     main="histogram of vegetation height")

# Distribucion normal es altamente recomendable

# Calculamos la matriz de distancias:
coords = cbind(matrix$x, matrix$y)
colnames(coords) = c("x", "y")
distmat = as.matrix(dist(coords))

#Maxima distancia para considerar en el correlograma
maxdist = 2/3*max(distmat)

# Construimos el correlograma
#--------------------------------------#
correlog.pgirmess = pgirmess::correlog(coords, matrix$Height, method="Moran",
                                        nbclass=14, alternative = "two.sided")

# Hacemos una grafica:
plot(correlog.pgirmess[,1], correlog.pgirmess[,2],
     xlab="Distance (m)", ylab="Moran's I")
abline(h=0)

# Podemos construit un correlograma con IC
spline.corr = spline.correlog(x = matrix$x, y = matrix$y, z = matrix$Height,
                               xmax = maxdist, resamp=99, type="boot")
## 10  of  99 
20  of  99 
30  of  99 
40  of  99 
50  of  99 
60  of  99 
70  of  99 
80  of  99 
90  of  99 
# Graficamos
plot(spline.corr)

# Variogramas
#############################################

# Crear un nuevo objeto:
geoR.veg = as.geodata(matrix)

# Grafiquemos
plot(geoR.veg)

# Semivariograma empirico:
emp = variog(geoR.veg, max.dist=maxdist)
## variog: computing omnidirectional variogram
# Grafiquemos
plot(emp)

# Podemos modificar el numero de lags tambien:
emp = variog(geoR.veg, max.dist=maxdist, breaks=c(seq(0,maxdist,by=3)))
## variog: computing omnidirectional variogram
# Grafiquemos
plot(emp)

# Verifiquemos anisotropia:

# Variograma direccional:
emp4 = variog4(geoR.veg, max.dist=maxdist)
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
# Grafiquemos
plot(emp4, legend = FALSE)

# Model-based variograms with likelihood fitting in geoR
#-------------------------------------------------------#

# Variograma teorico:
# Exponential variogram
mlexp = likfit(geoR.veg, cov.model="exp", ini=c(700,10))
## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
# Spherical variogram
mlsph = likfit(geoR.veg, cov.model="sph", ini=c(700,10))
## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
# Podemos comparar
summary(mlexp)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##    beta 
## 43.0713 
## 
## Parameters of the spatial component:
##    correlation function: exponential
##       (estimated) variance parameter sigmasq (partial sill) =  504.6
##       (estimated) cor. fct. parameter phi (range parameter)  =  5.883
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  732
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 17.62418
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
##  "-3298"      "4"   "6603"   "6621" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
##  "-3368"      "2"   "6739"   "6748" 
## 
## Call:
## likfit(geodata = geoR.veg, ini.cov.pars = c(700, 10), cov.model = "exp")
summary(mlsph)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##    beta 
## 43.3481 
## 
## Parameters of the spatial component:
##    correlation function: spherical
##       (estimated) variance parameter sigmasq (partial sill) =  451.5
##       (estimated) cor. fct. parameter phi (range parameter)  =  16.71
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  814.4
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 16.70779
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
##  "-3298"      "4"   "6604"   "6622" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
##  "-3368"      "2"   "6739"   "6748" 
## 
## Call:
## likfit(geodata = geoR.veg, ini.cov.pars = c(700, 10), cov.model = "sph")
AIC(mlexp,mlsph)
##       df      AIC
## mlexp  4 6603.375
## mlsph  4 6603.830
# Grafiquemos
plot(emp)
lines(mlexp, col="blue") 
lines(mlsph, col="red")

# Kriging e interpolation
############

# Creamos grilla de interpolacion:
new.grid.1m = expand.grid(0:max(matrix$x), 0:max(matrix$y))

# Ordinary kriging
krig.geoR.exp = krige.conv(geoR.veg,locations=new.grid.1m,
                          krige=krige.control(cov.pars=c(mlexp$cov.pars[1], mlexp$cov.pars[2]), nugget = mlexp$nugget,
                                              cov.model="exp", type.krige="OK"))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# Graficamos
image(krig.geoR.exp, main="krigged estimates")