Un ejemplo de los efectos de derrame
Yasna Cortés
2022-10-18
Objetivo
El objetivo de este taller es analizar los efectos de derrame (o efectos spillovers) que se originan a raíz del uso de modelos espaciales.
Ejemplo: Tiempos de viaje al CBD
Este es un ejemplo que aparece en Lesage and Pace (2009). Hay siete regiones que representan a una área metropolitana, donde la región 4 es el distrito central de negocios (CBD). Toda la región tiene una sola carretera, por tanto, todos los conmutantes usan la misma ruta para llegar al CBD.
Aquí, la variable dependiente son los tiempos de viaje al CBD (en minutos), mientras que las variables explicativas son la distancia (en millas) y la densidad poblacional.
A continuación, se presenta la estructura espacial de las observaciones/regiones:
Distribución Espacial de las Regiones
Variables del Ejemplo (Ver Lesage and Pace, 2009)
Ingresamos estas variables al código para trabajar con el ejemplo:
#Density population and distances
X<-cbind(c(10,20,30,50,30,20,10),
c(30,20,10,0,10,20,30))
colnames(X)<-c("Density", "Distance")
X## Density Distance
## [1,] 10 30
## [2,] 20 20
## [3,] 30 10
## [4,] 50 0
## [5,] 30 10
## [6,] 20 20
## [7,] 10 30#Travel times to CBD
Y<-cbind(c(42,37,30,26,30,37,42))
colnames(Y)<-c("Travel Times")
Y## Travel Times
## [1,] 42
## [2,] 37
## [3,] 30
## [4,] 26
## [5,] 30
## [6,] 37
## [7,] 42Construimos nuestra matriz de contiguidad, de acuerdo a la estructura espacial de las observaciones/regiones (matriz binaria):
#Contiguity matrix
C<-cbind(c(0,1,0,0,0,0,0),
c(1,0,1,0,0,0,0),
c(0,1,0,1,0,0,0),
c(0,0,1,0,1,0,0),
c(0,0,0,1,0,1,0),
c(0,0,0,0,1,0,1),
c(0,0,0,0,0,1,0))
C## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0 1 0 0 0 0 0
## [2,] 1 0 1 0 0 0 0
## [3,] 0 1 0 1 0 0 0
## [4,] 0 0 1 0 1 0 0
## [5,] 0 0 0 1 0 1 0
## [6,] 0 0 0 0 1 0 1
## [7,] 0 0 0 0 0 1 0También, la estandarizamos por filas
#Contiguity matrix
W<-C/rowSums(C)
W## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.0 1.0 0.0 0.0 0.0 0.0 0.0
## [2,] 0.5 0.0 0.5 0.0 0.0 0.0 0.0
## [3,] 0.0 0.5 0.0 0.5 0.0 0.0 0.0
## [4,] 0.0 0.0 0.5 0.0 0.5 0.0 0.0
## [5,] 0.0 0.0 0.0 0.5 0.0 0.5 0.0
## [6,] 0.0 0.0 0.0 0.0 0.5 0.0 0.5
## [7,] 0.0 0.0 0.0 0.0 0.0 1.0 0.0MODELO 1: Estimación vía OLS
Comenzaremos estimando el modelo vía OLS. Es decir,
\[ Y=X\beta+\epsilon \]
MODEL1.OLS<-lm(Y~X-1)
summary(MODEL1.OLS)##
## Call:
## lm(formula = Y ~ X - 1)
##
## Residuals:
## 1 2 3 4 5 6 7
## -0.9852 0.9926 0.9705 -1.5646 0.9705 0.9926 -0.9852
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## XDensity 0.55129 0.02064 26.71 1.37e-06 ***
## XDistance 1.24908 0.02839 43.99 1.15e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.284 on 5 degrees of freedom
## Multiple R-squared: 0.9991, Adjusted R-squared: 0.9987
## F-statistic: 2649 on 2 and 5 DF, p-value: 2.731e-08Por tanto, los parámetros de esta regresión son:
b_ols<-c(0.55129,1.24908)
b_ols## [1] 0.55129 1.24908Obtengamos el valor predicho de esta regresión (también, usar predict):
\[ Y = X\beta \]
y_hat_ols_S1<-crossprod(t(X), b_ols)
y_hat_ols_S1## [,1]
## [1,] 42.9853
## [2,] 36.0074
## [3,] 29.0295
## [4,] 27.5645
## [5,] 29.0295
## [6,] 36.0074
## [7,] 42.9853y_hat_ols_S1_2<-MODEL1.OLS$fitted.values
y_hat_ols_S1_2## 1 2 3 4 5 6 7
## 42.98524 36.00738 29.02952 27.56458 29.02952 36.00738 42.98524¿Qué ocurre si cambia la densidad poblacional en la región 2 (cambia en 20)?
Construimos una nueva matriz para incorporar este cambio:
X_T<-cbind(c(10,40,30,50,30,20,10),
c(30,20,10,0,10,20,30))
X_T## [,1] [,2]
## [1,] 10 30
## [2,] 40 20
## [3,] 30 10
## [4,] 50 0
## [5,] 30 10
## [6,] 20 20
## [7,] 10 30Obtengamos, nuevamente el valor predicho de esta regresión con la nueva variable de densidad:
y_hat_ols_S2<-crossprod(t(X_T), b_ols)
y_hat_ols_S2## [,1]
## [1,] 42.9853
## [2,] 47.0332
## [3,] 29.0295
## [4,] 27.5645
## [5,] 29.0295
## [6,] 36.0074
## [7,] 42.9853Resumen de resultados:
RES_OLS<-cbind(y_hat_ols_S1, y_hat_ols_S2, y_hat_ols_S2-y_hat_ols_S1)
RES_OLS## [,1] [,2] [,3]
## [1,] 42.9853 42.9853 0.0000
## [2,] 36.0074 47.0332 11.0258
## [3,] 29.0295 29.0295 0.0000
## [4,] 27.5645 27.5645 0.0000
## [5,] 29.0295 29.0295 0.0000
## [6,] 36.0074 36.0074 0.0000
## [7,] 42.9853 42.9853 0.0000MODELO 2: Estimación vía SAR
Introducimos nuestra matriz de pesos espaciales para incorporarla en nuestra estimación:
SW<-mat2listw(W)
SW## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 7
## Number of nonzero links: 12
## Percentage nonzero weights: 24.4898
## Average number of links: 1.714286
##
## Weights style: M
## Weights constants summary:
## n nn S0 S1 S2
## M 7 49 7 8.5 29Estimamos nuestro modelo usando un SAR:
\[ y=\rho W y +X\beta+\epsilon \]
MODEL1.SAR<-lagsarlm(Y~X-1, listw = SW)
summary(MODEL1.SAR)##
## Call:lagsarlm(formula = Y ~ X - 1, listw = SW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.073871 -0.059986 0.022529 0.051342 0.080155
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## XDensity 0.1351346 0.0090459 14.939 < 2.2e-16
## XDistance 0.5611967 0.0149665 37.497 < 2.2e-16
##
## Rho: 0.64298, LR test value: 38.047, p-value: 6.9063e-10
## Asymptotic standard error: 0.013907
## z-value: 46.233, p-value: < 2.22e-16
## Wald statistic: 2137.5, p-value: < 2.22e-16
##
## Log likelihood: 8.518579 for lag model
## ML residual variance (sigma squared): 0.0038434, (sigma: 0.061995)
## Number of observations: 7
## Number of parameters estimated: 4
## AIC: -9.0372, (AIC for lm: 27.01)
## LM test for residual autocorrelation
## test value: 4.6921, p-value: 0.030302Por tanto, los parámetros de esta regresión son:
b_sar<-c(0.1351346,0.5611967)
rho<-0.64298Obtengamos el valor predicho de esta regresión:
\[ Y=(I-\rho W)^{-1}X\beta \]
y_hat_sar_S1<-solve(diag(nrow(W))-rho*W)%*%crossprod(t(X),b_sar)
y_hat_sar_S1## [,1]
## [1,] 42.01651
## [2,] 37.06065
## [3,] 29.94227
## [4,] 26.00901
## [5,] 29.94227
## [6,] 37.06065
## [7,] 42.01651¿Qué ocurre si cambia la densidad poblacional en la región 2 (cambio en 20)?
Obtenemos los valores predichos considerando este cambio:
y_hat_sar_S2<-solve(diag(nrow(W))-rho*W)%*%crossprod(t(X_T),b_sar)
y_hat_sar_S2## [,1]
## [1,] 44.58643
## [2,] 41.05755
## [3,] 31.39798
## [4,] 26.54012
## [5,] 30.13860
## [6,] 37.14022
## [7,] 42.06766Resumen de resultados:
RES_SAR<-cbind(y_hat_sar_S1, y_hat_sar_S2, y_hat_sar_S2-y_hat_sar_S1)
RES_SAR## [,1] [,2] [,3]
## [1,] 42.01651 44.58643 2.56992136
## [2,] 37.06065 41.05755 3.99689160
## [3,] 29.94227 31.39798 1.45570804
## [4,] 26.00901 26.54012 0.53111252
## [5,] 29.94227 30.13860 0.19632629
## [6,] 37.06065 37.14022 0.07956368
## [7,] 42.01651 42.06766 0.05115785Matriz de derivadas parciales para la variable densidad:
\[ \begin{equation} \begin{bmatrix} \dfrac{\partial y_{1}}{\partial x_{1den}} & \dots & \dfrac{\partial y_{1}}{\partial x_{nden}} \\ \dots & \dots & \dots\\ \dfrac{\partial y_{n}}{\partial x_{1den}} & \dots & \dfrac{\partial y_{n}}{\partial x_{nden}} \end{bmatrix} \end{equation} =\beta_{den}(I-\rho W)^{-1} \]
b_den<-0.1351346
S_W<-b_den*(solve(diag(nrow(W))-rho*W))
S_W## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1764448009 0.128496068 0.046799558 0.01707474 0.006311694 0.002557893
## [2,] 0.0642480340 0.199844580 0.072785402 0.02655563 0.009816314 0.003978184
## [3,] 0.0233997789 0.072785402 0.179600648 0.06552698 0.024222116 0.009816314
## [4,] 0.0085373682 0.026555626 0.065526980 0.17726714 0.065526980 0.026555626
## [5,] 0.0031558469 0.009816314 0.024222116 0.06552698 0.179600648 0.072785402
## [6,] 0.0012789463 0.003978184 0.009816314 0.02655563 0.072785402 0.199844580
## [7,] 0.0008223369 0.002557893 0.006311694 0.01707474 0.046799558 0.128496068
## [,7]
## [1,] 0.0008223369
## [2,] 0.0012789463
## [3,] 0.0031558469
## [4,] 0.0085373682
## [5,] 0.0233997789
## [6,] 0.0642480340
## [7,] 0.1764448009Efectos Directos Promedio
ED<-sum(diag(S_W))/nrow(W)
ED## [1] 0.1841496Efectos Indirectos Promedio
d<-diag(S_W)
d1<-diag(d, 7,7)
PEI<-S_W-d1
EI<-sum(PEI)/nrow(W)
EI## [1] 0.1943575EFECTO TOTAL
ET<-ED+EI
ET## [1] 0.3785071Podemos estimar los impactos directamente:
impacts(MODEL1.SAR, listw = SW)## Impact measures (lag, exact):
## Direct Indirect Total
## XDensity 0.1841494 0.1943567 0.3785061
## XDistance 0.7647488 0.8071385 1.5718872Modelo 3: Estimación vía SLX
Vamos a realizar una modificación del modelo a estimar y usaremos solo la variable densidad:
DEN<-cbind(c(10,20,30,50,30,20,10))
DEN## [,1]
## [1,] 10
## [2,] 20
## [3,] 30
## [4,] 50
## [5,] 30
## [6,] 20
## [7,] 10DEN_T<-cbind(c(10,40,30,50,30,20,10))
DEN_T## [,1]
## [1,] 10
## [2,] 40
## [3,] 30
## [4,] 50
## [5,] 30
## [6,] 20
## [7,] 10MODEL1.SLX<-lmSLX(Y~DEN-1, listw = SW)
summary(MODEL1.SLX)##
## Call:
## lm(formula = formula(paste("y ~ 0 + ", paste(colnames(x), collapse = "+"))),
## data = as.data.frame(x), weights = weights)
##
## Residuals:
## 1 2 3 4 5 6 7
## 11.236 13.006 -15.374 3.549 -15.374 13.006 11.236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## DEN -0.6770 0.5775 -1.172 0.2939
## lag.DEN 1.8767 0.5975 3.141 0.0256 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.67 on 5 degrees of freedom
## Multiple R-squared: 0.8769, Adjusted R-squared: 0.8277
## F-statistic: 17.81 on 2 and 5 DF, p-value: 0.005317#Por tanto, los parámetros de esta regresión son:
b_slx<-c(-0.6770)
gam<-c(1.8767)y_hat_slx_S1<-crossprod(t(DEN),b_slx)+W%*%crossprod(t(DEN),gam)
y_hat_slx_S2<-crossprod(t(DEN_T),b_slx)+W%*%crossprod(t(DEN_T),gam)Resumen de resultados:
RES_SLX<-cbind(y_hat_slx_S1, y_hat_slx_S2, y_hat_slx_S2-y_hat_slx_S1)
RES_SLX## [,1] [,2] [,3]
## [1,] 30.7640 68.2980 37.534
## [2,] 23.9940 10.4540 -13.540
## [3,] 45.3745 64.1415 18.767
## [4,] 22.4510 22.4510 0.000
## [5,] 45.3745 45.3745 0.000
## [6,] 23.9940 23.9940 0.000
## [7,] 30.7640 30.7640 0.000impacts(MODEL1.SLX, listw = SW )## Impact measures (SlX, estimable):
## Direct Indirect Total
## DEN -0.6769953 1.876682 1.199687Matriz de derivadas parciales para la variable densidad:
\[ \begin{equation} \begin{bmatrix} \dfrac{\partial y_{1}}{\partial x_{1den}} & \dots & \dfrac{\partial y_{1}}{\partial x_{nden}} \\ \dots & \dots & \dots\\ \dfrac{\partial y_{n}}{\partial x_{1den}} & \dots & \dfrac{\partial y_{n}}{\partial x_{nden}} \end{bmatrix} \end{equation} =\beta_{den}I+W\gamma_{den} \]
b_den_slx<--0.6770
S_W_SLX<-b_den_slx*(diag(nrow(W)))+W*gam
S_W_SLX## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -0.67700 1.87670 0.00000 0.00000 0.00000 0.00000 0.00000
## [2,] 0.93835 -0.67700 0.93835 0.00000 0.00000 0.00000 0.00000
## [3,] 0.00000 0.93835 -0.67700 0.93835 0.00000 0.00000 0.00000
## [4,] 0.00000 0.00000 0.93835 -0.67700 0.93835 0.00000 0.00000
## [5,] 0.00000 0.00000 0.00000 0.93835 -0.67700 0.93835 0.00000
## [6,] 0.00000 0.00000 0.00000 0.00000 0.93835 -0.67700 0.93835
## [7,] 0.00000 0.00000 0.00000 0.00000 0.00000 1.87670 -0.67700Efectos Directos Promedio
ED_SLX<-sum(diag(S_W_SLX))/nrow(W)
ED_SLX## [1] -0.677Efectos Indirectos Promedio
d<-diag(S_W_SLX)
d1<-diag(d, 7,7)
PEI<-S_W_SLX-d1
EI_SLX<-sum(PEI)/nrow(W)
EI_SLX## [1] 1.8767EFECTO TOTAL
ET_SLX<-ED_SLX+EI_SLX
ET_SLX## [1] 1.1997