dist_covid

datos <- sf::st_read("C:/Users/sixto/OneDrive/Documentos/Papers/Land distances/dist_covid/datos_agregados.gpkg")

## Reading layer `datos_agregados' from data source 
##   `C:\Users\sixto\OneDrive\Documentos\Papers\Land distances\dist_covid\datos_agregados.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 91122 features and 2 fields (with 1 geometry empty)
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13016370 ymin: -7331031 xmax: 1.0463e-306 ymax: 3850032
## Projected CRS: WGS 84 / Pseudo-Mercator

datos <- sf::st_read("C:/Users/sixto/OneDrive/Documentos/Papers/Land distances/dist_covid/datos_globales_por_anio.gpkg")

## Reading layer `datos_globales_por_anio' from data source 
##   `C:\Users\sixto\OneDrive\Documentos\Papers\Land distances\dist_covid\datos_globales_por_anio.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 117639 features and 4 fields (with 2 geometries empty)
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13016370 ymin: -7331031 xmax: 1.0463e-306 ymax: 3850032
## Projected CRS: WGS 84 / Pseudo-Mercator

datos$covid <- ifelse(datos$anio < 2022, 0, 1)

alcaldias = c("BOGOTA","MEDELLIN","CORDOBA","LA PLATA", "CURITIBA","BELO HORIZONTE", "BUENOS AIRES", "FORTALEZA", "LOJA", "EL SALVADOR")


lat_coord = c(-74.0797202,-75.5781105,-64.2007216,-57.9576927,-49.2859511, -43.9507902, -58.3720765,-38.5247386, -79.2021704,-89.1926163)
lon_coord = c(4.5954893,6.2439913,-31.4131208,-34.9200356,-25.437794, -19.9288533, -34.6084374,-3.7252513, -3.9960413, 13.6977517)

alcaldias <- data.frame(alcaldias,lat_coord, lon_coord)

colnames(alcaldias) <- c("alcaldia",'lon','lat')

calculadora = function(ciudad, ciudad_bb , nombre_ciudad, alcaldia) {


#Primero, decidimos una ciudad#
  
 osm = function(ciudad , llave, valor) {

  opq(bbox = ciudad) %>%
    add_osm_feature(key = llave, value = valor) %>%
    osmdata_sf ()}
 
 

filtrar_lista = function(dataset) {

 dataset$osm_polygons}

  #Sacamos a la alcaldia, La alcaldia se pone en formato lon, lat (yyy,xxx)#
  
alcaldia =  st_as_sf(x = alcaldia, coords = c("lon", "lat"),  crs = 4326)

#Sacamos todas las calles#

calles = osm(ciudad,"highway", c("motorway","primary","secondary","tertiary"))
calles = calles$osm_lines

calles_autopista = filter(calles, highway =="motorway")
calles_primarias  = filter(calles, highway =="primary")
calles_secundarias  = filter(calles, highway =="secondary")
calles_terciarias  = filter(calles, highway =="tertiary")

#Paradas de tren#
tren = osm(ciudad, "railway","station")

tren = tren$osm_polygons
#Acá preguntar que onda tema lines, points and polygons#


#Ahora los polígonos industriales#

usos_industriales =  osm(ciudad, "landuse","industrial")
usos_industriales = usos_industriales$osm_polygons


#Ahora los rios#

rio =  osm(ciudad, "water","river")
rio = filtrar_lista(rio)

#Para fortaleza ver cuales tienen playa#
#If not null#
#if null then nothing#
#IFNA() =NA sacarlo de la formula#

#Vamos con los parques

parques    = osm(ciudad, "leisure", "park")

parques = parques$osm_polygons

#Calcular distancias del dataset#

#Distancia de tamaño#
 
#El valor de la ciudad para hacer el filtro del dataset#


datos <- st_transform(datos, st_crs(ciudad_bb))
ciudad_suelo = st_intersection(datos, ciudad_bb)
ciudad_suelo = sample_n(ciudad_suelo,1500)


#Calculamos la distancia en metros lineales, en un futuro un ruteo puede ser mejor#


ciudad_suelo = ciudad_suelo  %>% mutate(dist_parque  = (0.00001 + apply(st_distance(ciudad_suelo, parques), 1 , function (x) min(x))))

ciudad_suelo = ciudad_suelo  %>% mutate(dist_autopista  = (0.00001 + apply(st_distance(ciudad_suelo, calles_autopista), 1 , function (x) min(x))))

ciudad_suelo = ciudad_suelo  %>% mutate(dist_calles_primarias  = (0.00001 + apply(st_distance(ciudad_suelo, calles_autopista), 1 , function (x) min(x))))


ciudad_suelo = ciudad_suelo  %>% mutate(dist_calles_secundarias  = (0.00001 + apply(st_distance(ciudad_suelo, calles_secundarias), 1 , function (x) min(x))))


ciudad_suelo = ciudad_suelo  %>% mutate(dist_calles_terciarias  = (0.00001 + apply(st_distance(ciudad_suelo, calles_terciarias), 1 , function (x) min(x))))


ciudad_suelo = ciudad_suelo  %>% mutate(dist_tren  = (0.00001 + apply(st_distance(ciudad_suelo, tren), 1 , function (x) min(x))))

ciudad_suelo = ciudad_suelo  %>% mutate(dist_usos_industriales  = (0.00001 + apply(st_distance(ciudad_suelo, usos_industriales), 1 , function (x) min(x))))


ciudad_suelo = ciudad_suelo  %>% mutate(dist_rio  = (0.00001 + apply(st_distance(ciudad_suelo, rio), 1 , function (x) min(x))))

ciudad_suelo = ciudad_suelo  %>% mutate(dist_alcaldia  = (0.00001 + apply(st_distance(ciudad_suelo, alcaldia), 1 , function (x) min(x))))


}

Empezamos por Bogotá

ciudad = "Bogota D.C, Colombia"

ciudad_bb =  getbb(ciudad, format_out = "sf_polygon") 

ciudad_bb = ciudad_bb$polygon

ciudad_bb_ <- ciudad_bb[1, ]

filtrar_lista = function(dataset) {

 dataset$osm_polygons}

acaldia_bogota = alcaldias %>% filter(alcaldia =="BOGOTA")

Intersectamos los datos para que estén en el polígono de Bogotá

# Transform CRS of datos to match ciudad_bb_
datos <- st_transform(datos, st_crs(ciudad_bb_))

# Perform intersection
datos_bogota <- st_intersection(datos, ciudad_bb_)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

mapview(datos_bogota , zcol = "v_suelo")

bogota_econometria  = calculadora(ciudad = "Bogotá, Colombia", ciudad_bb = ciudad_bb_ ,alcaldia = acaldia_bogota, nombre_ciudad = "BOGOTA"   ) 

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: There were 1500 warnings in `stopifnot()`.
## The first warning was:
## ℹ In argument: `dist_autopista = (...)`.
## Caused by warning in `min()`:
## ! ningún argumento finito para min; retornando Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1499 remaining warnings.

## Warning: There were 1500 warnings in `stopifnot()`.
## The first warning was:
## ℹ In argument: `dist_calles_primarias = (...)`.
## Caused by warning in `min()`:
## ! ningún argumento finito para min; retornando Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1499 remaining warnings.

colnames(bogota_econometria)

##  [1] "pais"                    "ciudad"                 
##  [3] "v_suelo"                 "anio"                   
##  [5] "covid"                   "geom"                   
##  [7] "dist_parque"             "dist_autopista"         
##  [9] "dist_calles_primarias"   "dist_calles_secundarias"
## [11] "dist_calles_terciarias"  "dist_tren"              
## [13] "dist_usos_industriales"  "dist_rio"               
## [15] "dist_alcaldia"

form = log(v_suelo) ~ log(dist_parque) + log(dist_calles_secundarias) + log(dist_calles_terciarias) + log(dist_tren) + log(dist_usos_industriales) + log(dist_rio) +
                  log(dist_alcaldia)    +  log(dist_alcaldia):covid 

lm = lm(form, data = st_drop_geometry(bogota_econometria))
options(scipen = 999 )
summary(lm)

## 
## Call:
## lm(formula = form, data = st_drop_geometry(bogota_econometria))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.69189 -0.24177  0.05146  0.29146  1.26995 
## 
## Coefficients:
##                               Estimate Std. Error t value             Pr(>|t|)
## (Intercept)                   8.620875   0.237420  36.311 < 0.0000000000000002
## log(dist_parque)             -0.018467   0.009717  -1.900               0.0576
## log(dist_calles_secundarias) -0.105092   0.010340 -10.163 < 0.0000000000000002
## log(dist_calles_terciarias)  -0.018118   0.011252  -1.610               0.1076
## log(dist_tren)               -0.168687   0.012407 -13.596 < 0.0000000000000002
## log(dist_usos_industriales)   0.021464   0.009602   2.235               0.0255
## log(dist_rio)                -0.203497   0.015348 -13.259 < 0.0000000000000002
## log(dist_alcaldia)            0.135878   0.022049   6.163     0.00000000091877
## log(dist_alcaldia):covid     -0.017396   0.002470  -7.044     0.00000000000285
##                                 
## (Intercept)                  ***
## log(dist_parque)             .  
## log(dist_calles_secundarias) ***
## log(dist_calles_terciarias)     
## log(dist_tren)               ***
## log(dist_usos_industriales)  *  
## log(dist_rio)                ***
## log(dist_alcaldia)           ***
## log(dist_alcaldia):covid     ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4263 on 1491 degrees of freedom
## Multiple R-squared:  0.2944, Adjusted R-squared:  0.2906 
## F-statistic: 77.74 on 8 and 1491 DF,  p-value: < 0.00000000000000022

mapview(bogota_econometria, zcol= "v_suelo")

mapview(bogota_econometria, zcol= "dist_parque")

mapview(bogota_econometria, zcol= "dist_calles_secundarias")

mapview(bogota_econometria, zcol= "dist_calles_terciarias")

mapview(bogota_econometria, zcol= "dist_tren")

mapview(bogota_econometria, zcol= "dist_usos_industriales")

mapview(bogota_econometria, zcol= "dist_alcaldia")

mapview(bogota_econometria, zcol= "dist_rio")

Ahora vamos con la Parte de análisis espacial

  datos_1 = sample_n(bogota_econometria,1500)
  datos_1 = st_transform(datos_1, crs = 3857)
  
  
cord <- st_coordinates(datos_1)
cord = cord[,1:2]
summary(cord)

##        X                  Y         
##  Min.   :-8261030   Min.   :500221  
##  1st Qu.:-8253817   1st Qu.:511018  
##  Median :-8250423   Median :514995  
##  Mean   :-8250825   Mean   :515625  
##  3rd Qu.:-8247815   3rd Qu.:521001  
##  Max.   :-8239212   Max.   :531188

d <- dnearneigh(cord, 0, 500)
dlist <- nbdists(d, coordinates(cord))
# idlist <- lapply(dlist, function(x) 1/x)
# lw <- nb2listw(d, glist=idlist ,style="W" , zero.policy = TRUE)
lw <- nb2listw(d, style="W" , zero.policy = TRUE)
print(lw, zero.policy = T)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1500 
## Number of nonzero links: 8870 
## Percentage nonzero weights: 0.3942222 
## Average number of links: 5.913333 
## 51 regions with no links:
## 57 113 156 164 194 199 208 214 222 225 254 265 277 322 351 364 384 388
## 414 435 441 484 487 617 629 647 653 682 702 766 831 912 945 985 1021
## 1036 1072 1095 1119 1122 1126 1224 1227 1228 1244 1276 1330 1361 1366
## 1436 1476
## 
## Weights style: W 
## Weights constants summary:
##      n      nn   S0       S1       S2
## W 1449 2099601 1449 681.8716 5936.388

lw2 = listw2mat(lw)
v_suelo = as.matrix(datos_1$v_suelo)
vecinos = lw2 %*% v_suelo
vecinos =  as.data.frame(vecinos)
summary(vecinos)

##        V1        
##  Min.   :   0.0  
##  1st Qu.: 400.7  
##  Median : 548.7  
##  Mean   : 518.8  
##  3rd Qu.: 657.8  
##  Max.   :1846.0

datos = cbind(datos_1, vecinos)
# datos$vecinos = NULL

matriz_vecindarios <- function(datos, distance_threshold, weight_style) {

  # Extract coordinates from datos_1

  datos = sample_n(datos,950)
  datos = st_transform(datos, crs = 3857)
  cord <- st_coordinates(datos)
  cord <- cord[, 1:2]
  
  # Calculate summary statistics of the coordinates
  summary_cord <- summary(cord)
  
  # Create neighbor object
  d <- dnearneigh(cord, 0, distance_threshold)
  
  # Calculate neighbor distances
  dlist <- nbdists(d, coordinates(cord))
  
  # Create spatial weights matrix
  lw <- nb2listw(d, style = weight_style, zero.policy = TRUE)
  
  # Print the spatial weights matrix
  print(lw, zero.policy = TRUE)
  
  # Convert the spatial weights matrix to a matrix
  lw2 <- listw2mat(lw)
  
  # Extract the desired variable from datos_1
  v_suelo <- as.matrix(datos$v_suelo)
  
  # Multiply the spatial weights matrix with the variable
  vecinos <- lw2 %*% v_suelo
  vecinos <- as.data.frame(vecinos)
  
  # Combine the original data with the calculated variable
  datos <- cbind(datos, vecinos)
  
  return(datos)
}

bogota_econometria_espacial = matriz_vecindarios(bogota_econometria, 750, "W")

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 950 
## Number of nonzero links: 6984 
## Percentage nonzero weights: 0.7738504 
## Average number of links: 7.351579 
## 31 regions with no links:
## 8 42 127 171 216 233 274 281 284 307 331 347 350 440 480 486 532 550
## 551 594 606 670 692 705 755 784 789 816 838 904 915
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0       S1       S2
## W 919 844561 919 360.8296 3748.257

summary(bogota_econometria_espacial$v_suelo)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    25.0   416.5   551.0   541.6   664.0  1874.0

summary(bogota_econometria_espacial$V1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   414.1   541.9   522.3   649.2  1773.0

mapview(bogota_econometria_espacial, zcol = "V1")

lm = lm(form, data = st_drop_geometry(datos_1))
options(scipen = 999 )
summary(lm)

## 
## Call:
## lm(formula = form, data = st_drop_geometry(datos_1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.69189 -0.24177  0.05146  0.29146  1.26995 
## 
## Coefficients:
##                               Estimate Std. Error t value             Pr(>|t|)
## (Intercept)                   8.620875   0.237420  36.311 < 0.0000000000000002
## log(dist_parque)             -0.018467   0.009717  -1.900               0.0576
## log(dist_calles_secundarias) -0.105092   0.010340 -10.163 < 0.0000000000000002
## log(dist_calles_terciarias)  -0.018118   0.011252  -1.610               0.1076
## log(dist_tren)               -0.168687   0.012407 -13.596 < 0.0000000000000002
## log(dist_usos_industriales)   0.021464   0.009602   2.235               0.0255
## log(dist_rio)                -0.203497   0.015348 -13.259 < 0.0000000000000002
## log(dist_alcaldia)            0.135878   0.022049   6.163     0.00000000091877
## log(dist_alcaldia):covid     -0.017396   0.002470  -7.044     0.00000000000285
##                                 
## (Intercept)                  ***
## log(dist_parque)             .  
## log(dist_calles_secundarias) ***
## log(dist_calles_terciarias)     
## log(dist_tren)               ***
## log(dist_usos_industriales)  *  
## log(dist_rio)                ***
## log(dist_alcaldia)           ***
## log(dist_alcaldia):covid     ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4263 on 1491 degrees of freedom
## Multiple R-squared:  0.2944, Adjusted R-squared:  0.2906 
## F-statistic: 77.74 on 8 and 1491 DF,  p-value: < 0.00000000000000022

lm.morantest(lm, lw, zero.policy = T) 

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## Moran I statistic standard deviate = 37, p-value < 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.6582952553    -0.0039777655     0.0003203816

No hay evidencia de dependencia espacial en los residuos. En principio, el modelo lineal ajustaría ok. Veamos un poquito más en detalle con multiplicadores de lagrange.

lm.LMtests(lm, lw, test=c("RLMerr","RLMlag"), zero.policy = T)

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## RLMerr = 1368.1, df = 1, p-value < 0.00000000000000022
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## RLMlag = 2.2781, df = 1, p-value = 0.1312

Acá si aparece dependencia espacial, tanto en el término de error como en la variable dependiente.

sac = spatialreg::sacsarlm(lm, listw = lw, datos_1, zero.policy = TRUE)
summary(sac, Nagelkerke = TRUE)

## 
## Call:spatialreg::sacsarlm(formula = lm, data = datos_1, listw = lw, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.907234 -0.105463  0.031033  0.125592  1.958760 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                                Estimate Std. Error z value
## (Intercept)                   8.1419284  0.4076996 19.9704
## log(dist_parque)             -0.0102900  0.0066298 -1.5521
## log(dist_calles_secundarias) -0.0705921  0.0104824 -6.7344
## log(dist_calles_terciarias)  -0.0204382  0.0082004 -2.4924
## log(dist_tren)               -0.1420919  0.0217244 -6.5407
## log(dist_usos_industriales)   0.0211629  0.0107720  1.9646
## log(dist_rio)                -0.1567550  0.0251888 -6.2232
## log(dist_alcaldia)            0.1154309  0.0381494  3.0258
## log(dist_alcaldia):covid     -0.0145522  0.0015963 -9.1165
##                                           Pr(>|z|)
## (Intercept)                  < 0.00000000000000022
## log(dist_parque)                           0.12065
## log(dist_calles_secundarias)      0.00000000001647
## log(dist_calles_terciarias)                0.01269
## log(dist_tren)                    0.00000000006124
## log(dist_usos_industriales)                0.04946
## log(dist_rio)                     0.00000000048708
## log(dist_alcaldia)                         0.00248
## log(dist_alcaldia):covid     < 0.00000000000000022
## 
## Rho: -0.028435
## Asymptotic standard error: 0.0078519
##     z-value: -3.6214, p-value: 0.00029298
## Lambda: 0.70841
## Asymptotic standard error: 0.017001
##     z-value: 41.67, p-value: < 0.000000000000000222
## 
## LR test value: 986.06, p-value: < 0.000000000000000222
## 
## Log likelihood: -352.0665 for sac model
## ML residual variance (sigma squared): 0.079715, (sigma: 0.28234)
## Nagelkerke pseudo-R-squared: 0.63433 
## Number of observations: 1500 
## Number of parameters estimated: 12 
## AIC: 728.13, (AIC for lm: 1710.2)

De hecho, y en línea con lo anterior, el modelo mejora marginalmente (ver Akaike = AIC), y sólo da significativo el parámetro del rezago espacial en el término de error. Ajustemos de vuelta el modelo teniendo esto en cuenta.

spatialreg::impacts(sac, listw = lw, zero.policy = TRUE)

## Impact measures (sac, exact):
##                                   Direct      Indirect       Total
## log(dist_parque)             -0.01029172  0.0002862648 -0.01000546
## log(dist_calles_secundarias) -0.07060416  0.0019638586 -0.06864030
## log(dist_calles_terciarias)  -0.02044171  0.0005685874 -0.01987313
## log(dist_tren)               -0.14211614  0.0039529680 -0.13816317
## log(dist_usos_industriales)   0.02116647 -0.0005887465  0.02057772
## log(dist_rio)                -0.15678176  0.0043608930 -0.15242086
## log(dist_alcaldia)            0.11545065 -0.0032112662  0.11223939
## log(dist_alcaldia):covid     -0.01455473  0.0004048407 -0.01414989

Repetimos todo pero con otro modelo lineal que analiza la variación de la pendiente

lm.morantest(lm, lw, zero.policy = T) 

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## Moran I statistic standard deviate = 37, p-value < 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.6582952553    -0.0039777655     0.0003203816

No hay evidencia de dependencia espacial en los residuos. En principio, el modelo lineal ajustaría ok. Veamos un poquito más en detalle con multiplicadores de lagrange.

lm.LMtests(lm, lw, test=c("RLMerr","RLMlag"), zero.policy = T)

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## RLMerr = 1368.1, df = 1, p-value < 0.00000000000000022
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## RLMlag = 2.2781, df = 1, p-value = 0.1312

Acá si aparece dependencia espacial, tanto en el término de error como en la variable dependiente.

sac = spatialreg::sacsarlm(lm, listw = lw, datos_1, zero.policy = TRUE)
summary(sac, Nagelkerke = TRUE)

## 
## Call:spatialreg::sacsarlm(formula = lm, data = datos_1, listw = lw, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.907234 -0.105463  0.031033  0.125592  1.958760 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                                Estimate Std. Error z value
## (Intercept)                   8.1419284  0.4076996 19.9704
## log(dist_parque)             -0.0102900  0.0066298 -1.5521
## log(dist_calles_secundarias) -0.0705921  0.0104824 -6.7344
## log(dist_calles_terciarias)  -0.0204382  0.0082004 -2.4924
## log(dist_tren)               -0.1420919  0.0217244 -6.5407
## log(dist_usos_industriales)   0.0211629  0.0107720  1.9646
## log(dist_rio)                -0.1567550  0.0251888 -6.2232
## log(dist_alcaldia)            0.1154309  0.0381494  3.0258
## log(dist_alcaldia):covid     -0.0145522  0.0015963 -9.1165
##                                           Pr(>|z|)
## (Intercept)                  < 0.00000000000000022
## log(dist_parque)                           0.12065
## log(dist_calles_secundarias)      0.00000000001647
## log(dist_calles_terciarias)                0.01269
## log(dist_tren)                    0.00000000006124
## log(dist_usos_industriales)                0.04946
## log(dist_rio)                     0.00000000048708
## log(dist_alcaldia)                         0.00248
## log(dist_alcaldia):covid     < 0.00000000000000022
## 
## Rho: -0.028435
## Asymptotic standard error: 0.0078519
##     z-value: -3.6214, p-value: 0.00029298
## Lambda: 0.70841
## Asymptotic standard error: 0.017001
##     z-value: 41.67, p-value: < 0.000000000000000222
## 
## LR test value: 986.06, p-value: < 0.000000000000000222
## 
## Log likelihood: -352.0665 for sac model
## ML residual variance (sigma squared): 0.079715, (sigma: 0.28234)
## Nagelkerke pseudo-R-squared: 0.63433 
## Number of observations: 1500 
## Number of parameters estimated: 12 
## AIC: 728.13, (AIC for lm: 1710.2)

De hecho, y en línea con lo anterior, el modelo mejora marginalmente (ver Akaike = AIC), y sólo da significativo el parámetro del rezago espacial en el término de error. Ajustemos de vuelta el modelo teniendo esto en cuenta.

spatialreg::impacts(sac, listw = lw, zero.policy = TRUE)

## Impact measures (sac, exact):
##                                   Direct      Indirect       Total
## log(dist_parque)             -0.01029172  0.0002862648 -0.01000546
## log(dist_calles_secundarias) -0.07060416  0.0019638586 -0.06864030
## log(dist_calles_terciarias)  -0.02044171  0.0005685874 -0.01987313
## log(dist_tren)               -0.14211614  0.0039529680 -0.13816317
## log(dist_usos_industriales)   0.02116647 -0.0005887465  0.02057772
## log(dist_rio)                -0.15678176  0.0043608930 -0.15242086
## log(dist_alcaldia)            0.11545065 -0.0032112662  0.11223939
## log(dist_alcaldia):covid     -0.01455473  0.0004048407 -0.01414989

lm = lm(form, data = st_drop_geometry(datos_1))
options(scipen = 999 )
summary(lm)

## 
## Call:
## lm(formula = form, data = st_drop_geometry(datos_1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.69189 -0.24177  0.05146  0.29146  1.26995 
## 
## Coefficients:
##                               Estimate Std. Error t value             Pr(>|t|)
## (Intercept)                   8.620875   0.237420  36.311 < 0.0000000000000002
## log(dist_parque)             -0.018467   0.009717  -1.900               0.0576
## log(dist_calles_secundarias) -0.105092   0.010340 -10.163 < 0.0000000000000002
## log(dist_calles_terciarias)  -0.018118   0.011252  -1.610               0.1076
## log(dist_tren)               -0.168687   0.012407 -13.596 < 0.0000000000000002
## log(dist_usos_industriales)   0.021464   0.009602   2.235               0.0255
## log(dist_rio)                -0.203497   0.015348 -13.259 < 0.0000000000000002
## log(dist_alcaldia)            0.135878   0.022049   6.163     0.00000000091877
## log(dist_alcaldia):covid     -0.017396   0.002470  -7.044     0.00000000000285
##                                 
## (Intercept)                  ***
## log(dist_parque)             .  
## log(dist_calles_secundarias) ***
## log(dist_calles_terciarias)     
## log(dist_tren)               ***
## log(dist_usos_industriales)  *  
## log(dist_rio)                ***
## log(dist_alcaldia)           ***
## log(dist_alcaldia):covid     ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4263 on 1491 degrees of freedom
## Multiple R-squared:  0.2944, Adjusted R-squared:  0.2906 
## F-statistic: 77.74 on 8 and 1491 DF,  p-value: < 0.00000000000000022

lm.morantest(lm, lw, zero.policy = T) 

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## Moran I statistic standard deviate = 37, p-value < 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     0.6582952553    -0.0039777655     0.0003203816

No hay evidencia de dependencia espacial en los residuos. En principio, el modelo lineal ajustaría ok. Veamos un poquito más en detalle con multiplicadores de lagrange.

lm.LMtests(lm, lw, test=c("RLMerr","RLMlag"), zero.policy = T)

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## RLMerr = 1368.1, df = 1, p-value < 0.00000000000000022
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = form, data = st_drop_geometry(datos_1))
## weights: lw
## 
## RLMlag = 2.2781, df = 1, p-value = 0.1312

Acá si aparece dependencia espacial, tanto en el término de error como en la variable dependiente.

sac = spatialreg::sacsarlm(lm, listw = lw, datos_1, zero.policy = TRUE)
summary(sac, Nagelkerke = TRUE)

## 
## Call:spatialreg::sacsarlm(formula = lm, data = datos_1, listw = lw, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.907234 -0.105463  0.031033  0.125592  1.958760 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                                Estimate Std. Error z value
## (Intercept)                   8.1419284  0.4076996 19.9704
## log(dist_parque)             -0.0102900  0.0066298 -1.5521
## log(dist_calles_secundarias) -0.0705921  0.0104824 -6.7344
## log(dist_calles_terciarias)  -0.0204382  0.0082004 -2.4924
## log(dist_tren)               -0.1420919  0.0217244 -6.5407
## log(dist_usos_industriales)   0.0211629  0.0107720  1.9646
## log(dist_rio)                -0.1567550  0.0251888 -6.2232
## log(dist_alcaldia)            0.1154309  0.0381494  3.0258
## log(dist_alcaldia):covid     -0.0145522  0.0015963 -9.1165
##                                           Pr(>|z|)
## (Intercept)                  < 0.00000000000000022
## log(dist_parque)                           0.12065
## log(dist_calles_secundarias)      0.00000000001647
## log(dist_calles_terciarias)                0.01269
## log(dist_tren)                    0.00000000006124
## log(dist_usos_industriales)                0.04946
## log(dist_rio)                     0.00000000048708
## log(dist_alcaldia)                         0.00248
## log(dist_alcaldia):covid     < 0.00000000000000022
## 
## Rho: -0.028435
## Asymptotic standard error: 0.0078519
##     z-value: -3.6214, p-value: 0.00029298
## Lambda: 0.70841
## Asymptotic standard error: 0.017001
##     z-value: 41.67, p-value: < 0.000000000000000222
## 
## LR test value: 986.06, p-value: < 0.000000000000000222
## 
## Log likelihood: -352.0665 for sac model
## ML residual variance (sigma squared): 0.079715, (sigma: 0.28234)
## Nagelkerke pseudo-R-squared: 0.63433 
## Number of observations: 1500 
## Number of parameters estimated: 12 
## AIC: 728.13, (AIC for lm: 1710.2)

De hecho, y en línea con lo anterior, el modelo mejora marginalmente (ver Akaike = AIC), y sólo da significativo el parámetro del rezago espacial en el término de error. Ajustemos de vuelta el modelo teniendo esto en cuenta.

spatialreg::impacts(sac, listw = lw, zero.policy = TRUE)

## Impact measures (sac, exact):
##                                   Direct      Indirect       Total
## log(dist_parque)             -0.01029172  0.0002862648 -0.01000546
## log(dist_calles_secundarias) -0.07060416  0.0019638586 -0.06864030
## log(dist_calles_terciarias)  -0.02044171  0.0005685874 -0.01987313
## log(dist_tren)               -0.14211614  0.0039529680 -0.13816317
## log(dist_usos_industriales)   0.02116647 -0.0005887465  0.02057772
## log(dist_rio)                -0.15678176  0.0043608930 -0.15242086
## log(dist_alcaldia)            0.11545065 -0.0032112662  0.11223939
## log(dist_alcaldia):covid     -0.01455473  0.0004048407 -0.01414989