9 Regresión lineal

Librerías usadas

• rgdal
• spdep
• ggplot2
• lmtest
• RColorBrewer
• classInt
• spatialreg

9.1 Una tipología de modelos econométricos

Podemos agrupar los modelos en función de sus características. Aquí se propone, restringido a modelos de regresión lineal, una clasificación de modelos econométricos espaciales.

9.1.1 Una clasificación general

Recordemos que la característica principal de la econometría espacial es cómo los efectos espaciales son tomados en cuenta.

Esta característica presupone que el espacio ha sido formalizado de una forma u otra.

En general la matriz de pesos espaciales logra ser aplicable en varios contextos empíricos siempre que:

• los pesos recojan adecuadamente la dependencia espacial
• la heterogeneidad espacial sea tomada en cuenta en la especificación del modelo.

Una distinción importante en términos de estrategias de estimación y pruebas de hipótesis es de procesos espaciales:

• simultáneos: expresado como probabilidad conjunta
• condicionales: expresado como probabilidad condicional.

Ambos se suelen presentar de manera autor regresiva: donde el valor de una variable en un punto se relaciona con sus valores en el resto del sistema espacial.

Un modelo simultáneo se remonta al trabajo de Whittle (1954) que en notación matricial es equivalente a un modelo espacial autorregresivo de primer orden:

\[ y = \rho W y+\epsilon \] con \(y\) como el vector de la variable y \(\epsilon\) el vector del término de error, y \(W\) como la matriz de pesos espaciales correspondiente. La estimación de este modelo implica la especificación de una distribución de probabilidad conjunta y optimización no lineal.

El modelo condicional fue propuesto originalmente por Besag (1974) como alternativa a la estimación no lineal de Whittle.

En esencia, consiste en una relación lineal entre la esperanza condicional de la variable dependiente y sus valores en el resto del sistema (\(S\)):

\[ E[y_i|y_j \forall j\in S, j \neq i] = \rho W y \] La estructura condicional del modelo hace posible míninos cuadrados como técnica de estimación siempre que los datos se codifiquen de un modo particular.

La codificación consiste en eliminar observaciones tal que las variables dependientes que quedan se puedan considerar como independientes (no contiguas). Si hay pocas observaciones no se recomienda hacerlo, también debe considerarse que este proceso no es único por lo que pueden obtenerse resultados diferentes según la codificación.

Besag (1975), y Lele and Ord (1986) han desarrollado -a pseudo-máxima verosimilitud para el modelo condicional que evita la codificación; pero es más compleja matemáticamente.

La mayoría de situaciones en ciencias regionales empíricas se expresan más fácilmente en forma simultánea. Este enfoque también es más cercano al enfoque econométrico tradicional.

De forma análoga al enfoque Box-Jenkins en series de tiempo, los especificaciones de modelos espaciales combinan procesos auto regresivos y de media móvil. Un proceso espacial de media móvil tiene la forma:

\[ y =\epsilon+ \rho W \epsilon \]

donde \(y\) es un vector de variables, o un vector de términos aleatorios independientes y \(W\) es la matriz de peros espaciales.

Otro punto importante pone la diferenciación de modelos son los datos:

• En observaciones transversales (cross-sectional) hay información insuficiente en los datos para extraer un patrón simultáneo total de la interacción, y a reparametrización del modelo en términos auto regresivos es una necesidad.

• Si se tienen datos cross-sectional y de series de tiempo, la información temporal puede relajar esa necesidad.

9.1.2 Una taxonomía de Modelos de Regresión Lineal Espacial

Esta clasificación asume observaciones de unidades disponibles en corte transversal en un punto del tiempo.

Se deriven varios modelos mediante imposición de restricciones en el modelo general.

El punto de partida es:

\[ y = \rho W_1 y+X\beta+\epsilon\tag{9.1} \]

\[ \epsilon = \lambda W_2 \epsilon+\mu \]

con \(\mu\sim N(0,\Omega)\), y los elementos de la diagonal de la matriz de covarianza del error \(\Omega\) como

\[ \Omega_{jj} = h_{i}(z\alpha) \quad h_i>0 \] En esta especificación:

• \(\beta_{K\times 1}\) es el vector de parámetros de las variables exógenas (no rezagadas) \(X_{N\times K}\).
• \(\rho\) es el coeficiente de la variable dependiente espacialmente rezagada.
• \(\lambda\) es el coeficiente de la estructura autorregresiva espacial para el error \(\epsilon\).

La perturbación \(\mu\) se asume normalmente distribuida con matriz de covarianza \(\Omega\). Los elementos de la diagonal de \(\Omega\) permiten heterocedasticidad como una función de \(P+1\) variables exógenas \(z\) que incluyen el término constante.

Los \(P\) parámetros \(\alpha\) están asociados con los términos no constantes, tal que, para \(\alpha=0\), se tiene que

\[ h = \sigma^2, \]

(la situación clásica de homocedasticidad).

Las matrices \(N\times N\) (\(W_1\) y \(W_2\)) son matrices de pesos espaciales estandarizadas o no estandarizados.

\(W_1\) se asocia con el proceso autorregresivo de la variable dependiente. \(W_2\) se asocia con el proceso autorregresivo del error.

Esto sugiere que los dos procesos pueden tener diferentes estructuras espaciales. Así propuesto, los modelos pueden tener \(3+K+P\) parámetros desconocidos, en forma vectorial:

\[ \theta = [\rho,\beta',\lambda,\sigma^2,\alpha']'. \tag{9.2} \] Varias estructuras espaciales se derivan cuando sub vectores de \(\theta\) se hacen cero.

9.1.2.1 Modelo clásico de regresión lineal, sin efectos espaciales

para \(\rho=0\), \(\lambda=0\), \(\alpha=0\), (\(P+2\) restricciones):

\[ y = X\beta + \epsilon \tag{9.3} \]

9.1.2.2 Modelo autorregresivo espacial (SLM)

incluye especificaciones de factor común (con \(WX\) incluido en las variables explicativas), como caso especial (si se incluye \(WX\) se tiene el modelo de Durbin)

para \(\rho=0\), \(\alpha=0\), (\(P+1\) restricciones):

\[ y = \rho W_1 y + X\beta+\epsilon \tag{9.4} \]

9.1.2.3 Modelo de regresión lineal con error autorregresivo espacial (SEM)

para \(\rho=0\), \(\alpha=0\), (\(P+1\) restricciones):

\[ y = X\beta+(I-\lambda W_2)^{-1}\mu \tag{9.5} \] A partir del modelo (9.1), nota que:

\[ \epsilon -\lambda W_2 \epsilon = \mu \\ (I-\lambda W_2)\epsilon = \mu \\ (I-\lambda W_2)^{-1}(I-\lambda W_2)\epsilon = (I-\lambda W_2)^{-1}\mu \\ \epsilon=(I-\lambda W_2)^{-1}\mu \]

9.1.2.4 Modelo autorregresivo espacial mixto regresivo con error espacial autorregresivo (SARMA)

para \(\alpha=0\), (\(P\) restricciones):

\[ y = \rho W_1 y+X\beta+(I-\lambda W_2)^{-1}\mu \tag{9.6} \] Se puede obtener cuatro especificaciones más si se permite heterocedasticidad (una forma para \(h(z\alpha)\)) en (9.3)-(9.6).

Los modelos que estimaremos en las siguientes secciones son SAR, SEM y SDM. La figura ?? muestra la relación entre ellos.

9.2 Modelo de rezagos espacial (SLM)

El modelo de rezagos espacial (Spatial lag model o SLM) es:

\[ y = \rho Wy+X\beta+\epsilon \tag{9.7} \] donde \(y\) es un vector \(n\times 1\) de observaciones de la variable dependiente, \(X\) es una matriz \(n\times K\) de observaciones de las variables explicativas, \(\beta\) es el vector de parámetros \(K\times 1\).

La siguiente figura ?? representa el modelo autorregresivo espacial en (9.7) para dos regiones. Las variables (\(x_1\), \(x_2\)) y los términos no observados (\(\epsilon_1\), \(\epsilon_2\)) tienen un efecto directo sobre \(y\) para ambas regiones. Tenga en cuenta que el modelo incorpora efectos de desbordamiento espacial por el efecto de \(y_1\) sobre \(y_2\) y viceversa. Es decir, el modelo refleja la simultaneidad inherente a la autocorrelación espacial.

La ecuación (9.7) expresa el modelo SLM en su forma estructural o de comportamiento donde se relaciona tanto las variables exógenas como endógenas con la dependiente \(y\).

Otra forma de expresar el modelo es en su forma reducida:

\[ y = (I-\rho W)^{-1}X\beta +(I-\rho W)^{-1}\epsilon \tag{9.8} \] que ya no contiene ninguna variable dependiente espacialmente rezagada en el lado de la derecha. La ecuación (9.8) expresa la naturaleza simultánea del proceso autorregresivo espacial.

La forma estructural es un modelo de la forma \(f(y,X,\epsilon)\), mientras la forma reducida es de la forma \(y = g(X,\epsilon)\).

Notemos que \((I-\rho W)\) debe ser invertible, es decir \(det(I-\rho W)\neq 0\). Entonces ¿qué valores de \(\rho\) hacen que \((I-\rho W)\) sea invertible? Se puede demostrar que para matrices simétricas, \((I-\rho W)\) es invertible si

\[ \omega_{min}^{-1}<\rho <\omega_{max}^{-1} \] donde \(\omega_{min}^{-1}\) y \(\omega_{max}^{-1}\) son los valores propios mínimo y máximo de \(W\).

Recuerda que es común estandarizar \(W\) por fila tal que sus filas sumen \(1\). Dado que \(W\) es no negativa, esto asegura que los pesos estén entre \(0\) y \(1\) y la estandarización se puede interpretar como promediar los valores de los vecinos.

Si \(W\) es estandarizada por filas, esto implica que \(-1<\rho <1\) pero recuerda que esto no significa que \(\rho\) sea la correlación convencional, esto solo sucede porque \(W\) es estandarizada por filas. Un problema de estandarizar por filas es que la comparación entre filas puede resultar complicado.

Otra forma de estandarizar \(W\) es multiplicarla por un número, digamos \(a\) tal que se remueva el efecto de las unidades de medida. En particular:

\[ a = min\{r,c\}\\ r = \underset{i}{max}\sum_j|w_{ij}| \quad \text{máximo de la suma for filas en valor absoluto}\\ c = \underset{j}{max}\sum_i|w_{ij}| \quad \text{máximo de la suma for columnas en valor absoluto} \] Como resultado, \((I-\rho W)\) será invertible para \(|\rho|=\frac{1}{a}\). Ten en cuenta que esta normalización tiene la ventaja de garantizar que las ponderaciones espaciales resultantes, \(w_{ij}\), estén todas entre \(0\) y \(1\) y, por lo tanto, aún puedan interpretarse como intensidades de influencia relativas.

Los momentos del modelo en forma reducida son:

Esperanza

\[ \mathbb{E}(y|X,W) = \mathbb{E}\left[(I-\rho W)^{-1}X\beta+ (I-\rho W)^{-1}\epsilon|X,W \right]\\ = (I-\rho W)^{-1}X\beta \]

[Expansión de Leontief]. Si \(|\rho|<1\), entonces \[(I-\rho W)^{-1}=\sum_{i=0}^{\infty}(\rho W)^i\]

• En promedio, el valor de \(y\) en una ubicación \(i\) no solo se explica por los valores de las variables explicativas asociadas a esta ubicación, sino también por las asociadas a todas las demás ubicaciones (vecinas o no) a través de la transformación espacial inversa \((I-\rho W)^{-1}\).

• Este efecto multiplicador espacial disminuye con la distancia, es decir, las potencias de \(W\) en la expansión en serie de \((I-\rho W)^{-1}\)

• Un choque aleatorio en una ubicación \(i\) no solo afecta el valor de \(y\) en esta ubicación, sino que también tiene un impacto en los valores de \(y\) en todas las demás ubicaciones a través de la misma transformación inversa espacial \((I-\rho W)^{-1}\).

Varianza

\[ Var(y|X,W) = \mathbb{E}(yy^T|W,X)\\ =(I-\rho W)^{-1}\mathbb{E}(\epsilon\epsilon^T|W,X)(I-\rho W^T) \]

• Esta matriz de varianza-covarianza es completa, lo que implica que cada ubicación está correlacionada con todas las demás ubicaciones del sistema.

• Sin embargo, esta correlación disminuye con la distancia. Esta ecuación también muestra que la covarianza entre cada par de términos de error no es nula y disminuye con el orden de proximidad.

• Además, los elementos de la diagonal de \(\mathbb{E}(\epsilon\epsilon^T|W,X)\) no son constantes. Esto implica heterocedasticidad de error de \(\epsilon\).

• Dado que no hemos asumido nada sobre la varianza del error, podemos decir que \(\mathbb{E}(\epsilon\epsilon^T|W,X)\) es una matriz completa, digamos \(\Omega_\epsilon\). Esto cubre la posibilidad de heterocedasticidad, autocorrelación espacial o ambas. En ausencia de cualquiera de estas complicaciones, la matriz de varianza se simplifica a la \(\sigma^2I\) habitual.

9.2.0.1 Ejemplo

Realizaremos una regresión donde la variable dependiente es la pobreza por NBI (personas) y la independiente es el valor agregado bruto no petrolero per cápita.

1. Preparamos los datos: esto incluye el cálculo de la matriz de ponderaciones espaciales.

# **** Importo datos: VAB no petrolero
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/VABNoPetroleroCantones2007-2019.csv"
datos <- read.csv(uu,sep = ",")
datos$COD_CANT[nchar(datos$COD_CANT)==3] = paste("0",datos$COD_CANT[nchar(datos$COD_CANT)==3],sep = "")

datos <- subset(datos,datos$YEAR == 2019)

# **** Importo datos: Población por cantones
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/proyeccion_cantonal_total_2010-2020.csv"

poblacion <- read.csv(uu,sep = ";")

poblacion$CODIGO[nchar(poblacion$CODIGO)==3] = paste("0",poblacion$CODIGO[nchar(poblacion$CODIGO)==3],sep = "")

# **** Importo datos: NBI
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/NBI_PER_CANT.csv"
NBI <- read.csv(uu, sep = ";")

NBI$CODIGO[nchar(NBI$CODIGO)==3] = paste("0",NBI$CODIGO[nchar(NBI$CODIGO)==3],sep = "")

datos <- merge(datos,poblacion[,c("CODIGO","A_2020","SuperficieKM2")],by.x = "COD_CANT",by.y = "CODIGO",all.x = TRUE)

datos <- merge(datos,NBI[,c("CODIGO","POBRES_P")],by.x = "COD_CANT",by.y = "CODIGO",all.x = TRUE)

datos$VAB_PC <- datos$VAB/datos$A_2020 #miles de USD / pob
datos$DensPob <- datos$A_2020/datos$SuperficieKM2 

# **** Importo datos: NBICOVID (https://www.covid19ecuador.org/cantones)

uu <- "https://raw.githubusercontent.com/andrab/ecuacovid/master/datos_crudos/positivas/por_fecha/cantones_por_dia_acumuladas.inec.csv"

c19 <- read.csv(uu, sep = ",")
c19 <- c19[,c("inec_canton_id","X24.03.2022")]
c19$inec_canton_id[nchar(c19$inec_canton_id)==3] = paste("0",c19$inec_canton_id[nchar(c19$inec_canton_id)==3],sep = "")
names(c19) <- c("COD_CANT","C19")

# names(datos)
datos <- merge(datos,c19, by.x = "COD_CANT",all.x = TRUE)
datos$C19_pc <- (datos$C19/datos$A_2020)*1000

# cor.test(datos$C19_pc,datos$DensPob,use="na.or.complete")
# cor.test(datos$C19_pc,datos$DensPob,use="na.or.complete",method = "spearman")
# cor.test(datos$C19_pc,datos$POBRES_P,use="na.or.complete",method = "spearman")
# cor.test(datos$DensPob,datos$POBRES_P,use="na.or.complete",method = "spearman")
# plot(datos$DensPob,datos$POBRES_P)


library(rgdal)
uu <- "https://github.com/vmoprojs/DataLectures/raw/master/SpatialData/2012_nxcantones.zip"
poligonos <- read_git_shp(uu)

## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/0p/n_r_hl095sv7nktfp_8n9n_80000gn/T/Rtmp4Gvwcm/file28d2799ba36/nxcantones.shp", layer: "nxcantones"
## with 224 features
## It has 6 fields

poligonos <- poligonos[poligonos$DPA_PROVIN!="90",]
poligonos <- poligonos[poligonos$DPA_PROVIN!="20",]
poligonos@data = merge(poligonos@data,datos,by.x  ="DPA_CANTON",by.y = "COD_CANT",all.x=TRUE)



# Me quedo con los datos que no son perdidos:
poligonos <- poligonos[!is.na(poligonos@data$POBRES_P),]

library(spdep)
set.ZeroPolicyOption(TRUE)

## [1] TRUE

# Se construye la lista de vecinos
cont.nb <- poly2nb(poligonos)
# Se construye la matriz de pesos espaciales
cont.listw <- nb2listw(cont.nb, style="W",zero.policy = TRUE)

2. Revisamos si el NBI tiene autocorrelación espacial.

y <- poligonos$POBRES_P 
x <- (poligonos$VAB_PC)
# x <- poligonos$CONSUMO_INTERMEDIO/poligonos$A_2019
# x <- poligonos$PRODUCCION/poligonos$A_2019
# cor(x,y)
# plot(y~x)
MImc <- moran.mc(y, cont.listw, 999)
MImc

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  y 
## weights: cont.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.084924, observed rank = 969, p-value = 0.031
## alternative hypothesis: greater

# (Zmi <- (MImc$statistic-mean(MImc$res))/sd(MImc$res))

Se rechaza la hipótesis nula de independencia espacial.

3. Ajustamos un modelo de regresión lineal.

mod.lm <- lm(y~x)
summary(mod.lm)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37428 -0.08681  0.00911  0.10271  0.40619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.820727   0.013324  61.597  < 2e-16 ***
## x           -0.018275   0.003332  -5.484 1.17e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1284 on 213 degrees of freedom
## Multiple R-squared:  0.1237, Adjusted R-squared:  0.1196 
## F-statistic: 30.07 on 1 and 213 DF,  p-value: 1.173e-07

En promedio, por cada unidad que aumenta el VAB per cápita, la pobreza disminuye 0.018275.

4. Revisamos si los residuos de la regresión tienen autocorrelación espacial.

library(RColorBrewer)
library(classInt)

res <- mod.lm$residuals

res.palette <- colorRampPalette(c("red","orange","white", "lightgreen","green"), space = "rgb")
pal <- res.palette(5)

classes_fx <- classIntervals(round(res,2), n=5, style="fixed", fixedBreaks=as.numeric(quantile(res)), rtimes = 1)

## Warning in classIntervals(round(res, 2), n = 5, style = "fixed", fixedBreaks =
## as.numeric(quantile(res)), : variable range greater than fixedBreaks

cols <- findColours(classes_fx,pal)


sp::plot(poligonos,col=cols, border="grey",main = "Residuos del modelo OLS")
legend(x="bottomleft",cex=0.8,fill=attr(cols,"palette"),bty="n",legend=names(attr(cols, "table")),ncol=5)

## Residual Autocorrelation

moran.mc(res, listw=cont.listw, zero.policy=T,nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res 
## weights: cont.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.10837, observed rank = 991, p-value = 0.009
## alternative hypothesis: greater

Se rechaza la hipótesis nula de independencia.

5. Caculamos el modelo SLM.

La función lagsarlm del paquete spatialreg (R. Bivand and Piras (2015)) nos permite hacer la estimación:

mod.slm <- spatialreg::lagsarlm(y~x,  listw=cont.listw, zero.policy=T, na.action=na.exclude)
summary(mod.slm)

## 
## Call:spatialreg::lagsarlm(formula = y ~ x, listw = cont.listw, na.action = na.exclude, 
##     zero.policy = T)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.376239 -0.083052  0.010570  0.100231  0.393399 
## 
## Type: lag 
## Regions with no neighbours included:
##  61 96 102 104 108 118 119 123 135 189 190 220 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  0.7699413  0.0350445 21.9704 < 2.2e-16
## x           -0.0178382  0.0033108 -5.3879  7.13e-08
## 
## Rho: 0.06787, LR test value: 2.4927, p-value: 0.11437
## Asymptotic standard error: 0.043284
##     z-value: 1.568, p-value: 0.11688
## Wald statistic: 2.4587, p-value: 0.11688
## 
## Log likelihood: 138.5559 for lag model
## ML residual variance (sigma squared): 0.016117, (sigma: 0.12695)
## Number of observations: 215 
## Number of parameters estimated: 4 
## AIC: -269.11, (AIC for lm: -268.62)
## LM test for residual autocorrelation
## test value: 2.478, p-value: 0.11545

6. Revisamos si los residuos siguen presentando autocorrelación espacial:

res.sar <- resid(mod.slm)
moran.mc(res.sar, listw=cont.listw, zero.policy=T,nsim = 9999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res.sar 
## weights: cont.listw  
## number of simulations + 1: 10000 
## 
## statistic = 0.068521, observed rank = 9261, p-value = 0.0739
## alternative hypothesis: greater

No se rechaza la hipótesis nula de independencia.

9.3 Modelo de Durbin espacial (SDM)

El modelo de Durbin así como su forma estructural (9.9) y su forma reducida (9.10):

\[ y = \rho Wy+X\beta+WX\gamma+\epsilon \tag{9.9} \]

\[ y = (I-\rho W)^{-1}(X\beta+WX\gamma+\epsilon) \tag{9.10} \]

El SDM da como resultado un modelo espacial autorregresivo de una forma especial, que incluye no solo la variable dependiente espacialmente rezagada y las variables explicativas, sino también las variables explicativas espacialmente rezagadas, \(WX\). \(y\) depende de factores de si misma de la matriz \(X\), más los mismos factores promediados en las \(n\) regiones vecinas.

Ten en cuenta en la siguiente figura ?? que la Región 1 no solo ejerce un impacto en la Región 2 (y viceversa) a través de \(y\), sino también a través de la variable independiente \(x\).

Este modelo también tiene muy buenas propiedades en términos de cálculo de efectos marginales que exploraremos más adelante.

Siguiendo con nuestro ejemplo, la estimación de Durbin en R se puede usar el argumento type = "mixed":

mod.sdm <- spatialreg::lagsarlm(y~x,  listw=cont.listw, zero.policy=T, na.action=na.exclude,type = "mixed")
summary(mod.sdm)

## 
## Call:spatialreg::lagsarlm(formula = y ~ x, listw = cont.listw, na.action = na.exclude, 
##     type = "mixed", zero.policy = T)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.375220 -0.083068  0.011459  0.097578  0.391961 
## 
## Type: mixed 
## Regions with no neighbours included:
##  61 96 102 104 108 118 119 123 135 189 190 220 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  0.7653862  0.0358089 21.3742 < 2.2e-16
## x           -0.0177758  0.0033091 -5.3718 7.793e-08
## lag.x        0.0038194  0.0083162  0.4593     0.646
## 
## Rho: 0.059768, LR test value: 1.6689, p-value: 0.19641
## Asymptotic standard error: 0.04792
##     z-value: 1.2472, p-value: 0.21231
## Wald statistic: 1.5556, p-value: 0.21231
## 
## Log likelihood: 138.6673 for mixed model
## ML residual variance (sigma squared): 0.016105, (sigma: 0.1269)
## Number of observations: 215 
## Number of parameters estimated: 5 
## AIC: -267.33, (AIC for lm: -267.67)
## LM test for residual autocorrelation
## test value: 3.325, p-value: 0.068234

9.4 Modelo con error espacial (SEM)

Otra forma de dependencia espacial ocurre cuando la dependencia opera a través del proceso de error, ya que los errores de diferentes áreas pueden mostrar autocorrelación espacial. El modelo SEM se formula como:

\[ y = X\beta+u \tag{9.11}\\ u = \lambda Wu+\epsilon \] donde \(\lambda\) es el parámetro autorregresivo del rezago espacial del error \(Wu\) y \(\epsilon\) es \(\text{i.i.d.}\).

En su forma estructural tenemos:

\[ y = X\beta+(I-\rho W)^{-1}\epsilon \] Esta expresión conduce a un efecto de difusión espacial global, pero no hay un efecto multiplicador espacial.

La siguiente figura ?? visualiza el \(SEM\) para dos regiones. Ten en cuenta que el término de error de ambas regiones está relacionado y el único efecto espacial va de \(\epsilon_1\) a \(\epsilon_2\) y viceversa.

Como afirman Luc Anselin and Bera (1998), la dependencia del error espacial puede interpretarse como una nuisance (\(\lambda\) como un parámetro de nuisance) en el sentido de que refleja la autocorrelación espacial en los errores de medición o en variables que de otro modo no son cruciales para el modelo (es decir, las variables ignoradas se propagan a través de las unidades espaciales de observación).

A diferencia de los modelos anteriores, los efectos de las interacciones entre los términos de error no requieren un modelo teórico para un proceso de interacción espacial o social, sino que son consistentes con una situación en la que los determinantes de la variable dependiente omitidos del modelo están espacialmente autocorrelacionados, o con una situación donde los choques no observados siguen un patrón espacial.

Los efectos de interacción entre los términos no observados también pueden interpretarse como un mecanismo para corregir a los políticos que buscan rentas por cambios imprevistos en la política fiscal. Ver por ejemplo, Allers and Elhorst (2005).

Siguiendo con nuestro ejemplo, la estimación de \(SEM\) en R se puede usar la función errorsarlm:

mod.sem <- spatialreg::errorsarlm(y~x,  listw=cont.listw, zero.policy=T)
summary(mod.sem)

## 
## Call:spatialreg::errorsarlm(formula = y ~ x, listw = cont.listw, zero.policy = T)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.362896 -0.084167  0.011036  0.093360  0.418968 
## 
## Type: error 
## Regions with no neighbours included:
##  61 96 102 104 108 118 119 123 135 189 190 220 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.820608   0.014498 56.6010  < 2e-16
## x           -0.018612   0.003257 -5.7145  1.1e-08
## 
## Lambda: 0.19491, LR test value: 4.6304, p-value: 0.031411
## Asymptotic standard error: 0.088053
##     z-value: 2.2135, p-value: 0.026862
## Wald statistic: 4.8997, p-value: 0.026862
## 
## Log likelihood: 139.6248 for error model
## ML residual variance (sigma squared): 0.015829, (sigma: 0.12581)
## Number of observations: 215 
## Number of parameters estimated: 4 
## AIC: -271.25, (AIC for lm: -268.62)

9.5 El mejor modelo espacial

¿Cómo decidir? Para el caso simple de elegir entre una alternativa SLM o SEM, existe evidencia de que el modelo adecuado es probablemente el que tiene el valor de prueba LM significativo más grande (Luc Anselin and Rey (1991)).

• Cuando el test \(\text{LM}_{LAG}\) es significativo y el \(\text{LM}_{ERR}\) es no es significativo, el modelo más apropiado es el modelo \(\text{SLM}\);

• en la misma línea, cuando la prueba \(\text{LM}_{ERR}\) es significativa y el valor \(\text{LM}_{LAG}\) no es significativo, el modelo más apropiado es el modelo \(\text{SEM}\).

Como puede adivinar, a veces es posible encontrar que ambas pruebas estadísticas son significativas. En este caso, una regla de decisión puede ser la siguiente:

• cuando el valor \(\text{LM}_{LAG}\) de prueba es mayor que el valor \(\text{LM}_{ERR}\) de prueba, sería mejor considerar el modelo \(\text{SLM}\);

• cuando el valor \(\text{LM}_{ERR}\) de prueba es mayor que el valor \(\text{LM}_{LAG}\) de prueba, sería mejor considerar el modelo \(\text{SEM}\)

Por supuesto, si ambas estadísticas son significativas, también podría ser apropiado estimar un modelo autorregresivo general (\(\text{SAC}\)).

Siguiendo con nuestro ejemplo, la prueba del multiplicador de Lagrange en R se puede hacer con la función lm.LMtests:

LM <- lm.LMtests(mod.lm, cont.listw, test="all")

## Please update scripts to use lm.RStests in place of lm.LMtests

print(LM)

## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = y ~ x)
## test weights: listw
## 
## RSerr = 4.8505, df = 1, p-value = 0.02764
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = y ~ x)
## test weights: listw
## 
## RSlag = 2.5537, df = 1, p-value = 0.11
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = y ~ x)
## test weights: listw
## 
## adjRSerr = 2.7219, df = 1, p-value = 0.09898
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = y ~ x)
## test weights: listw
## 
## adjRSlag = 0.42507, df = 1, p-value = 0.5144
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = y ~ x)
## test weights: listw
## 
## SARMA = 5.2756, df = 2, p-value = 0.07152

Las pruebas \(LM\) robustas tienen en cuenta la posibilidad alternativa, es decir, la prueba LMerr responderá tanto a una variable dependiente rezagada espacialmente omitida como a los residuos autocorrelacionados espacialmente, mientras que la prueba RLMerr robusta está diseñada para probar residuos autocorrelacionados espacialmente en la posible presencia de una variable dependiente espacialmente rezagada omitida.

9.6 Interpretación de los modelos espaciales

9.6.1 Propagación (spillovers) espaciales

Un enfoque importante de la ciencia regional es el desbordamiento espacial. Una definición básica de derrames en un contexto espacial sería que los cambios que ocurren en una región ejercen impactos en otras regiones (J. P. LeSage and Pace (2014)). Algunos ejemplos son:

• Los cambios en la tasa impositiva por una unidad espacial podrían tener un impacto en las decisiones de fijación de la tasa impositiva de las regiones cercanas, un fenómeno que ha sido etiquetado como imitación de impuestos y competencia de criterios entre gobiernos locales.

• Situaciones en las que las mejoras a la vivienda realizadas por un propietario ejercen un impacto beneficioso en los precios de venta de las viviendas vecinas.

• La innovación de los investigadores universitarios se difunde a las empresas cercanas.

• La contaminación del aire o del agua generada en una región se extiende a las regiones cercanas.

9.6.1.1 Spillovers globales

Los efectos de propagación globales surgen cuando los cambios en una característica de una región afectan los resultados de todas las regiones. Esto se aplica incluso a la propia región, ya que los impactos pueden pasar a los vecinos y regresar a la propia región (retroalimentación). Específicamente, los spillovers globales impactan a los vecinos, los vecinos a los vecinos, los vecinos a los vecinos de los vecinos, etc.

9.6.1.2 Spillovers locales

Los efectos de propagación locales representan una situación en la que el impacto recae solo en los vecinos cercanos o inmediatos, desapareciendo antes de que afecten a las regiones vecinas de los vecinos.

Como se puede observar en las definiciones anteriores, la principal diferencia es que la retroalimentación o la interacción endógena solo es posible para los efectos de contagio globales.

9.6.2 Efectos marginales

Matemáticamente, la noción de spillover puede pensarse como la derivada \(\partial y_i/\partial x_j\). Esto significa que los cambios en las variables explicativas en la región \(i\) impactan en la variable dependiente en la región \(j \neq i\).

En general, se puede notar que el cambio de cada variable en cada región implica \(n^2\) efectos marginales potenciales. Si tenemos \(K\) variables en nuestro modelo, esto implica \(K\times n^2\) medidas potenciales. Incluso para valores pequeños de \(n\) y \(K\), puede resultar bastante difícil informar estos resultados de forma compacta.

Para superar este problema, J. LeSage and Pace (2009) (p. 36-37) proponen las siguientes medidas de resumen escalares:

9.6.2.1 Impacto directo promedio (ADI)

Sea \(S_r = A(W)^{-1}(I_n\beta_r+W\theta_r)\) para la variable \(r\) (\(A(W)=(I_n-\rho W)\)). El impacto de los cambios en la i-ésima observación de \(x_r\), denotada por \(x_{ir}\), en \(y_i\) puede ser resumida midiendo el promedio \(S_r(W)_{ii}\):

\[ \text{ADI}=\frac{1}{n}\text{tr}(S_r(W)) \tag{9.12} \]

El promedio sobre el impacto directo asociado con todas las observaciones \(i\) es similar a las interpretaciones típicas de coeficientes de regresión que representan la respuesta promedio de las variables dependientes a independientes.

Por ejemplo, si la región \(i\) aumenta el número de viajeros que usan el transporte público, ¿cuál será el impacto promedio en los tiempos de viaje en la región \(i\)? Esta medida tendrá en cuenta los efectos de retroalimentación que surgen del cambio en el uso del transporte público de la \(i\)-ésima región en los tiempos de desplazamiento de las regiones vecinas en el sistema de regiones espacialmente dependientes.

9.6.2.2 Impacto total promedio hacia una observación (ATIT)

Sea \(S_r = A(W)^{-1}(I_n\beta_r+W\theta_r)\) para la variable \(r\) (\(A(W)=(I_n-\rho W)\)). La suma en la \(i\)-ésima fila de \(S_r (W)\) representaría el impacto total en la observación individual \(y_i\) resultante de cambiar la \(r\)-ésima variable explicativa en la misma cantidad en todas las \(n\) observaciones. Hay \(n\) de estas sumas dadas por el vector de columna \(c_r = S_r (W)\iota_n\), por lo que un promedio de estos impactos totales es:

\[ \text{ATIT} = \frac{1}{n}\iota'_nc_r \tag{9.13} \]

donde \(\iota_n\) denota un vector columna compuesto de unos, esto es: \(\iota_n = (1,1,\ldots,1)'\).

9.6.2.3 Impacto total promedio desde una observación (ATIF)

Sea \(S_r = A(W)^{-1}(I_n\beta_r+W\theta_r)\) para la variable \(r\) (\(A(W)=(I_n-\rho W)\)). La suma de la \(j\)-ésima columna de \(S_r (W)\) produciría el impacto total sobre todo \(y_i\) al cambiar la \(r\)-ésima variable explicativa en una cantidad en la \(j\)-ésima observación. Hay \(n\) de estas sumas dadas por el vector fila \(r_r\iota'_nS_r(W)\), por lo que un promedio de estos impactos totales es:

\[ \text{ATIF} = \frac{1}{n}r_r\iota_n \tag{9.14} \]

La ecuación (9.14) relaciona cómo los cambios en una sola observación \(j\) influyen en todas las observaciones. Por el contrario, la ecuación (9.13) considera cómo los cambios en todas las observaciones influyen en una sola observación \(i\). En ambos casos, promediar todas las \(n\) observaciones conduce al mismo resultado numérico. La implicación de este interesante resultado es que el impacto total promedio es el promedio de todas las derivadas de \(y_i\) con respecto a \(x_{jr}\) para cualquier \(i\), \(j\).

Entonces

\[ \bar{M}(r)_{\text{directo}} = n^{-1}\text{tr}(S_r(W)) \] \[ \bar{M}(r)_{\text{total}} = n^{-1}\iota'_n\text{tr}(S_r(W))\iota_n \]

\[ \bar{M}(r)_{\text{indirecto}} = \bar{M}(r)_{\text{total}}-\bar{M}(r)_{\text{directo}} \]

Para el modelo SEM, las propagaciones espaciales (spillovers) son cero por diseño, como en el caso de los modelos de regresión no espacial. Es decir, los coeficientes se interpretan como en regresión lineal clásica: la estimación del parámetro \(\beta_r\) es igual al efecto marginal promedio de la variable sobre la variable dependiente \(y\).

Siguiendo con nuestro ejemplo, el modelo elegido es el \(\text{SEM}\), en R estimamos sus impactos con la función impacts:

spatialreg::impacts(mod.slm, listw=cont.listw)

## Impact measures (lag, exact):
##        Direct     Indirect       Total
## x -0.01785774 -0.001279337 -0.01913707

Tras bastidores:

listw = cont.listw
slm <- mod.slm

rho <- slm$rho # rho estimado del modelo SLM
beta_hat <- coef(slm)[-1] # Parámtros estimados
A <- spatialreg::invIrW(listw, rho = rho) # (I - rho*W)^{-1}
X <- cbind(1, x) # Matriz de covariables
y_hat_pre <- A %*% crossprod(t(X), beta_hat) # y estimado

## Construimos: S_r(W) = A(W)^-1 (I * beta_r + W * theta_r)
Ibeta <- diag(length(listw$neighbours)) * coef(slm)["x"]
S <- A %*% Ibeta
ADI <- sum(diag(S)) / nrow(A) # Directo

n <- length(listw$neighbours)
Total <- crossprod(rep(1, n), S) %*% rep(1, n) / n

Indirect <- Total - ADI

c(Direct = ADI, Indirect=Indirect,Total = Total)

##       Direct     Indirect        Total 
## -0.017857738 -0.001206843 -0.019064581

listw = cont.listw
slm <- mod.sem

rho <- slm$lambda # rho estimado del modelo SLM
beta_hat <- coef(slm)[-1] # Parámtros estimados
A <- spatialreg::invIrW(listw, rho = rho) # (I - rho*W)^{-1}
X <- cbind(1, x) # Matriz de covariables
y_hat_pre <- A %*% crossprod(t(X), beta_hat) # y estimado

## Construimos: S_r(W) = A(W)^-1 (I * beta_r + W * theta_r)
Ibeta <- diag(length(listw$neighbours)) * coef(slm)["x"]
S <- A %*% Ibeta
ADI <- sum(diag(S)) / nrow(A) # Directo

n <- length(listw$neighbours)
Total <- crossprod(rep(1, n), S) %*% rep(1, n) / n

Indirect <- Total - ADI

c(Direct = ADI, Indirect=Indirect,Total = Total)

##       Direct     Indirect        Total 
## -0.018788673 -0.004078085 -0.022866758

Pruebe los modelos \(SLM\), \(SDM\) y \(SEM\), recuerde agregar covariables al modelo.

9.6.3 Una interpretación más práctica

• Difusión y efectos externos (desbordamiento): la mayoría de las veces se derivan de efectos directos e indirectos, que miden el impacto de variables exógenas en una ubicación y vecindario determinados. Muestran la interdependencia de las unidades espaciales: cuando es positivo, la mayoría de las veces se interpreta como cooperación y cuando es negativo, como separación/rivalidad/competencia.

	Efecto directo	Efecto Indirecto (spillover)
OLS/SEM	\(\beta_k\)	0
SAR/SAC	Elementos diagonales de la matriz \((I-\rho W)^{-1}\beta_k\)	Elementos NO diagonales de la matriz \((I-\rho W)^{-1}\beta_k\)
SXL/SDEM	\(\beta_k\)	\(\theta_k\)
SDM/GNS	Elementos diagonales de la matriz \((I-\rho W)^{-1}[\beta_k+W\theta_k]\)	Elementos NO diagonales de la matriz \((I-\rho W)^{-1}[\beta_k+W\theta_k]\)

donde,

• Modelos con dos componentes espaciales:

  • \(SAC\) (spatial autocorrelation model): \(Y = \beta_0+\rho W Y+X\beta+u \text{ donde } u = \lambda Wu+e\)
  • \(SDM\) (spatial Durbin model): \(Y = \beta_0+\rho W Y+X\beta+WX\theta+e\)
  • \(SDEM\) (spatial Durbin error model): \(Y = \beta_0+X\beta+WX\theta+u \text{ donde } u = \lambda Wu+e\)

• Modelos con un componente espacial

  • \(SAR\) (spatial autoregressive model): \(Y = \beta_0+\rho W Y+X\beta+e\) - \(SEM\) (spatial error model): \(Y = \beta_0+X\beta+u \text{ donde } u = \lambda Wu+e\)
    
  • \(SLX\) (explanatory variables spatially lagged): \(Y = \beta_0+X\beta++WX\theta+e\)

9.7 Ejercicios

1. Realizar un modelo de regresión entre NBI personas (2010) y el número de robos totales 2019 en Ecuador a nivel cantonal.

2. Usando un análisis de regresión espacial completo usando datos a nivel provincial de la Encuesta Nacional de Actividades de Innovación (AI) 2012-2014:

https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/ACTIprov.csv

Contiene las siguientes variables:

• DPAprov: DPA de la provincia
• exitoso: Porcentaje de innovaciones exitosas. Igual a 1 cuando la innovación se realiza en procesos y productos.
• actipor: Innovación interna. Valor destinado hacia actividades de investigación y desarrollo al interior de la empresa
• ImasD: Porcentaje de empresas que cooperaron en actividades de Investigación y Desarrollo.
• IngDis: Porcentaje de empresas que cooperaron en actividades de Ingeniería y Diseño.
• Capa: Porcentaje de empresas que cooperaron en actividades de Capacitación.
• Asistec: Porcentaje de empresas que cooperaron en actividades de Asistencia Técnica.
• Info: Porcentaje de empresas que cooperaron en actividades de Información.
• Pruebas: Porcentaje de empresas que cooperaron en actividades de Pruebas de Productos.
• Público: Porcentaje de gasto en Investigación y Desarrollo financiado por el Gobierno sobre el total de gasto en I+D.
• Apoyo: Monto de apoyo no reembolsable en Investigación y Desarrollo.
• Financiero: Promedio de financiamiento realizado por la Banca Privada.
• GDPP: Logaritmo del valor agregado bruto

Estime el modelo base:

\[ exitoso = \beta_0+\beta_1ImasD+\beta_2IngDis+\beta_3Capa+\beta_4AsisTec+\beta_5Info+\beta_6Pruebas+\beta_7Publico+\beta_8Apoyo+\beta_9Financiero+\beta_10log(GDPP) \] Pruebe los modelos \(SLM\), \(SDM\) y \(SEM\).

uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/ACTIprov.csv"