Carga de librerias.

library(sf)

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spatialreg)

## Loading required package: spData

## Loading required package: Matrix

library(spdep)

## 
## Attaching package: 'spdep'

## The following objects are masked from 'package:spatialreg':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(ggrepel)
library(gamlss)

## Loading required package: splines

## Loading required package: gamlss.data

## 
## Attaching package: 'gamlss.data'

## The following object is masked from 'package:datasets':
## 
##     sleep

## Loading required package: gamlss.dist

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: parallel

##  **********   GAMLSS Version 5.5-0  **********

## For more on GAMLSS look at https://www.gamlss.com/

## Type gamlssNews() to see new features/changes/bug fixes.

library(gamlss)
library(dplyr)

Vivienda Habilitada en Bogotá

El propósito de este estudio es analizar la estructura espacial de los desarrollos de vivienda nueva VIS en Bogotá y evaluar si las UPZ con mayor presencia de estos proyectos tienden a agruparse geográficamente, es decir, si existe dependencia o autocorrelación espacial en su distribución.

Carda de datos.

El conjunto de datos corresponde al registro de viviendas nuevas VIP, VIS y No VIS construidas en el Distrito Capital que cuentan con servicio definitivo de acueducto y alcantarillado. Estas viviendas incluyen soluciones habitacionales desarrolladas en predios de propiedad privada, sin estar asociadas a proyectos urbanísticos formales, y que representan la consolidación de nuevas unidades residenciales dentro del tejido urbano.

vivienda <- st_read("/Users/usermac/Documents/especializacion/estadistica espacial/vivienda-habilitada/Vivienda Habilitada.shp")

## Reading layer `Vivienda Habilitada' from data source 
##   `/Users/usermac/Documents/especializacion/estadistica espacial/vivienda-habilitada/Vivienda Habilitada.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 7976 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.2144 ymin: 4.4699 xmax: -74.0128 ymax: 4.8227
## Geodetic CRS:  MAGNA-SIRGAS

st_crs(vivienda)$wkt

## [1] "GEOGCRS[\"MAGNA-SIRGAS\",\n    DATUM[\"Marco Geocentrico Nacional de Referencia\",\n        ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank.\"],\n        BBOX[-4.23,-84.77,15.51,-66.87]],\n    ID[\"EPSG\",4686]]"

Vemos que el CRS del conjunto de datos de Vivienda Habilitada es de \(4686\)

unique(vivienda$año_habil)

##  [1] 2012 2014 2013 2015 2016 2018 2019 2017 2020 2021 2022 2023

Y que tenemos datos desde el 2012 al 2023.

ggplot(data = vivienda) +
  geom_sf(
    color =  "#FA5939",
    size = 0.2
  ) +
  labs(title = "Mapa de Vivienda Habilitada en Bogotá") +
  theme_light()

Como deseamos revisar por UPZ, desde la pagina de Datos Abiertos de Bogotá, podemos traer la geometrias de las UPZ.

url_upz <- "https://mapas.gobiernobogota.gov.co/waserver/rest/services/Mapa_Base/MapServer/8/query?where=1%3D1&outFields=%2A&f=geojson"
upz <- st_read(url_upz)

## Reading layer `OGRGeoJSON' from data source 
##   `https://mapas.gobiernobogota.gov.co/waserver/rest/services/Mapa_Base/MapServer/8/query?where=1%3D1&outFields=%2A&f=geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 116 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.44978 ymin: 3.73103 xmax: -74.00978 ymax: 4.836779
## Geodetic CRS:  WGS 84

ggplot(data = upz) +
  geom_sf(
    fill = NA,
    color =  "#FA5939",
    size = 0.2
  ) +
  labs(title = "UPZ - Bogotá") +
  theme_light()

Fijamos en ambos conjuntos de datos el mismo CRS

vivienda <- st_transform(vivienda, st_crs(upz))

Hacemos el join espacil.

vivienda_upz <- st_join(vivienda, upz[c("UPLCODIGO", "UPLNOMBRE")],
                        join = st_within, left = TRUE)

head(vivienda_upz)

## Simple feature collection with 6 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.1931 ymin: 4.5812 xmax: -74.0677 ymax: 4.7392
## Geodetic CRS:  WGS 84
##   id_viviend año_habil mes_habili     direccion_ estrato_vi numero_viv
## 1          4      2012          1  CL 5B S 3A 20          2          0
## 2          5      2012          1 CL 128D 86B 14          2          0
## 3          6      2012          1 KR 73B 146F 25          3          0
## 4         10      2012          1 KR 88C 62 09 S          2          0
## 5         15      2012          1  KR 36A 53A 96          5          0
## 6         16      2012          1   CL 152 72 35          4          0
##   numero_v_1 numero_v_2 total_vivi tipo_vivie iniciativa    lotcodigo UPLCODIGO
## 1         48          0         48          2          1 001101006045     UPZ33
## 2         80          0         80          2          1 009263079033     UPZ28
## 3        264          0        264          2          1 009128020002     UPZ23
## 4         84          0         84          2          1 004637079014     UPZ84
## 5          0         20         20          2          1 005106008030    UPZ106
## 6          0         50         50          2          1 009128015003     UPZ23
##          UPLNOMBRE                geometry
## 1          SOSIEGO POINT (-74.0815 4.5812)
## 2        EL RINCON POINT (-74.0856 4.7221)
## 3 CASA BLANCA SUBA POINT (-74.0703 4.7388)
## 4  BOSA OCCIDENTAL POINT (-74.1931 4.6255)
## 5     LA ESMERALDA POINT (-74.0813 4.6448)
## 6 CASA BLANCA SUBA POINT (-74.0677 4.7392)

resumen_upz_year <- vivienda_upz %>% 
  st_drop_geometry() %>%
  group_by( año_habil) %>% 
  summarise(
    sum_total_vis = sum(numero_v_1, na.rm = TRUE),
    sum_total_novis = sum(numero_v_2, na.rm = TRUE),
    estrato_moda = names(which.max(table(estrato_vi)))
  )

ggplot(resumen_upz_year, aes(x = año_habil, y = sum_total_vis)) +
  geom_line(color = "#FA5939", size = 1) + 
  geom_point(color = "#FA5939", size = 2) +
  scale_x_continuous(breaks = unique(resumen_upz_year$año_habil)) + 
  labs(
    title = "Número total de vivienda VIS habilitada por año",
    subtitle = "Serie de tiempo",
    x = "Año",
    y = "Total de Vivienda VIS"
  ) +
  theme_light() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

resumen_upz <- vivienda_upz %>% 
  st_drop_geometry() %>%
  group_by(UPLCODIGO, UPLNOMBRE) %>% 
  summarise(
    sum_total_vis = sum(numero_v_1, na.rm = TRUE),
    sum_total_novis = sum(numero_v_2, na.rm = TRUE),
    estrato_moda = names(which.max(table(estrato_vi)))
  )

## `summarise()` has grouped output by 'UPLCODIGO'. You can override using the
## `.groups` argument.

head(resumen_upz)

## # A tibble: 6 × 5
## # Groups:   UPLCODIGO [6]
##   UPLCODIGO UPLNOMBRE                 sum_total_vis sum_total_novis estrato_moda
##   <chr>     <chr>                             <dbl>           <dbl> <chr>       
## 1 UPZ1      PASEO DE LOS LIBERTADORES             1               0 2           
## 2 UPZ10     LA URIBE                            641            3859 3           
## 3 UPZ100    GALERIAS                            482            2267 4           
## 4 UPZ101    TEUSAQUILLO                         339             842 4           
## 5 UPZ102    LA SABANA                           534            1047 3           
## 6 UPZ103    PARQUE SALITRE                        0             325 3

Creamos el conjunto de datos final

upz_final <- upz %>%
  inner_join(resumen_upz, 
            by = c("UPLCODIGO", "UPLNOMBRE"))

head(upz_final)

## Simple feature collection with 6 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.18523 ymin: 4.59182 xmax: -74.04869 ymax: 4.740841
## Geodetic CRS:  WGS 84
##   OBJECTID UPLCODIGO UPLTIPO     UPLNOMBRE
## 1        1     UPZ37       1  SANTA ISABEL
## 2        2     UPZ49       1        APOGEO
## 3        3     UPZ19       1      EL PRADO
## 4        4     UPZ24       1          NIZA
## 5        5    UPZ100       1      GALERIAS
## 6        6     UPZ81       1 GRAN BRITALIA
##                                                               UPLAADMINI
## 1              Dcto. 349-15/08/2002 Mod.=Res 0071/2004 (Gaceta 290/2004)
## 2                                  Dcto. 180-8/06/2005 (Gaceta 366/2005)
## 3    Dcto. 299-10/07/2002 (Gaceta 251/2002) Mod.=Res 820/2007, 1673/2010
## 4 Dcto. 175-31/05/2006 Mod.=Dec 368/2008 (Gaceta 507/2008), Res 699/2009
## 5                                 Dcto. 621-29/12/2006 (Gaceta 456/2007)
## 6                                 Dcto. 217-13/07/2005 (Gaceta 373/2005)
##      LOCNOMBRE LOCCODIGO ENTCODIGO   SHAPE.AREA  SHAPE.LEN sum_total_vis
## 1 LOS MARTIRES        14        14 0.0001632389 0.05539622            52
## 2         BOSA        07         7 0.0001714811 0.07336394          1592
## 3         SUBA        11        11 0.0003529858 0.09591693             8
## 4         SUBA        11        11 0.0006162044 0.12220205           266
## 5  TEUSAQUILLO        13        13 0.0001934658 0.06056277           482
## 6      KENNEDY        08         8 0.0001465005 0.06921748           446
##   sum_total_novis estrato_moda                       geometry
## 1             599            3 MULTIPOLYGON (((-74.08478 4...
## 2               5            2 MULTIPOLYGON (((-74.15216 4...
## 3            3605            3 MULTIPOLYGON (((-74.04869 4...
## 4            2943            4 MULTIPOLYGON (((-74.07828 4...
## 5            2267            4 MULTIPOLYGON (((-74.06463 4...
## 6              35            2 MULTIPOLYGON (((-74.17064 4...

upz_final$proportion_vis <- upz_final$sum_total_vis/(upz_final$sum_total_vis+upz_final$sum_total_novis)

upz_final <- upz_final[!is.na(upz_final$proportion_vis), ]

Fijamos una proyección de CRS valida para encontrar distancias.

upz_final <- st_transform(upz_final, crs = 3116)

ggplot(data = upz_final) +
  geom_sf(aes(fill = proportion_vis), color = "white", size = 0.1) + 
  scale_fill_viridis_c(option = "magma", direction = -1, name = "Proporción Viviendas VIS") +
  labs(
    title = "Disponibilidad de Vivienda Nueva VIS",
    subtitle = "Distribución por UPZ en Bogotá",
    caption = "Fuente: Datos Abiertos de Bogotá"
  ) +
  theme_void() + 
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(color = "gray50")
  )

ggplot(data = upz_final) +
  geom_sf(aes(fill = sum_total_vis), color = "white", size = 0.1) + 
  scale_fill_viridis_c(option = "magma", direction = -1, name = "Total Viviendas VIS") +
  labs(
    title = "Disponibilidad de Vivienda Nueva VIS",
    subtitle = "Distribución por UPZ en Bogotá",
    caption = "Fuente: Datos Abiertos de Bogotá"
  ) +
  theme_void() + 
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(color = "gray50")
  )

ggplot(data = upz_final) +
  geom_sf(aes(fill = sum_total_novis), color = "white", size = 0.1) + 
  scale_fill_viridis_c(option = "magma", direction = -1, name = "Total Viviendas") +
  labs(
    title = "Disponibilidad de Vivienda Nueva No VIS",
    subtitle = "Distribución por UPZ en Bogotá",
    caption = "Fuente: Datos Abiertos de Bogotá"
  ) +
  theme_void() + 
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(color = "gray50")
  )

Calculamos los centroides de las UPZ.

centroids <- st_centroid(upz_final)

## Warning: st_centroid assumes attributes are constant over geometries

head(centroids)

## Simple feature collection with 6 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 989185.2 ymin: 1000244 xmax: 1002292 ymax: 1014152
## Projected CRS: MAGNA-SIRGAS / Colombia Bogota zone
##   OBJECTID UPLCODIGO UPLTIPO     UPLNOMBRE
## 1        1     UPZ37       1  SANTA ISABEL
## 2        2     UPZ49       1        APOGEO
## 3        3     UPZ19       1      EL PRADO
## 4        4     UPZ24       1          NIZA
## 5        5    UPZ100       1      GALERIAS
## 6        6     UPZ81       1 GRAN BRITALIA
##                                                               UPLAADMINI
## 1              Dcto. 349-15/08/2002 Mod.=Res 0071/2004 (Gaceta 290/2004)
## 2                                  Dcto. 180-8/06/2005 (Gaceta 366/2005)
## 3    Dcto. 299-10/07/2002 (Gaceta 251/2002) Mod.=Res 820/2007, 1673/2010
## 4 Dcto. 175-31/05/2006 Mod.=Dec 368/2008 (Gaceta 507/2008), Res 699/2009
## 5                                 Dcto. 621-29/12/2006 (Gaceta 456/2007)
## 6                                 Dcto. 217-13/07/2005 (Gaceta 373/2005)
##      LOCNOMBRE LOCCODIGO ENTCODIGO   SHAPE.AREA  SHAPE.LEN sum_total_vis
## 1 LOS MARTIRES        14        14 0.0001632389 0.05539622            52
## 2         BOSA        07         7 0.0001714811 0.07336394          1592
## 3         SUBA        11        11 0.0003529858 0.09591693             8
## 4         SUBA        11        11 0.0006162044 0.12220205           266
## 5  TEUSAQUILLO        13        13 0.0001934658 0.06056277           482
## 6      KENNEDY        08         8 0.0001465005 0.06921748           446
##   sum_total_novis estrato_moda                 geometry proportion_vis
## 1             599            3 POINT (997978.3 1000244)    0.079877112
## 2               5            2 POINT (990283.1 1000455)    0.996869130
## 3            3605            3  POINT (1002292 1014152)    0.002214226
## 4            2943            4  POINT (1000430 1013413)    0.082891867
## 5            2267            4  POINT (1000527 1005235)    0.175336486
## 6              35            2 POINT (989185.2 1002620)    0.927234927

Creamos la matriz de distancias.

wdist <- st_distance(centroids)
head(wdist)

## Units: [m]
##           1         2         3         4         5         6         7
## 1     0.000  7698.074 14561.489 13395.404  5604.496  9108.362  2274.606
## 2  7698.074     0.000 18216.068 16458.068 11304.774  2427.215  7809.227
## 3 14561.489 18216.068     0.000  2003.831  9089.550 17458.165 12469.953
## 4 13395.404 16458.068  2003.831     0.000  8178.593 15586.457 11209.257
## 5  5604.496 11304.774  9089.550  8178.593     0.000 11639.950  3842.096
## 6  9108.362  2427.215 17458.165 15586.457 11639.950     0.000  8631.774
##           8         9        10       11        12       13       14       15
## 1  9247.093 10207.800  8322.143 8089.698  7027.446 7168.595 11575.23 10871.84
## 2 15013.135  8801.915 12375.309 8715.304 11106.812 9857.396 13651.71 12023.08
## 3  6531.132 10979.460  6301.073 9675.951  7660.142 8414.114 25959.20 25415.58
## 4  6432.510  8992.458  5088.118 7817.465  6370.971 6797.331 24927.39 24262.14
## 5  3814.052  8722.902  3212.023 6066.236  2395.866 3876.963 16869.66 16363.77
## 6 15122.392  7196.500 12070.542 7782.245 10934.479 9406.411 16074.81 14450.04
##          16        17        18        19        20        21        22
## 1  3419.174  3639.052  4066.093  4070.145 11712.473  9515.321  5672.646
## 2  5744.101  9031.171 10566.443 11211.594  5463.410  4516.455 12001.821
## 3 13291.982 10952.239 10807.152 11577.156 17566.454 15606.930  9506.457
## 4 11778.598  9756.435  9901.262 10828.354 15594.834 13676.894  8803.810
## 5  5569.479  2340.643  1728.087  2723.747 13244.455 10832.914  1136.865
## 6  6393.907  9555.080 11247.917 12098.749  3063.225  2370.098 12503.036
##          23        24        25        26        27        28        29
## 1  5005.627  6400.546  2517.179  4140.594  1872.961  5165.492  2500.736
## 2  8477.453 12849.627  9163.846  7249.729  9020.586  2717.741 10101.192
## 3 10243.505  9207.861 16466.012 11562.022 12880.770 16008.087 13320.391
## 4  8783.447  8678.343 15473.476 10071.583 11833.180 14357.388 12442.692
## 5  3222.276  1819.103  7376.472  4143.721  3822.371  8634.449  4267.413
## 6  8521.546 13341.936 10949.768  7502.623 10087.662  3980.804 11306.825
##          30        31        32        33        34        35        36
## 1 16816.879 14195.268 13678.411 12137.692  9825.240  5612.414 10117.903
## 2 18671.623 14847.328 18178.082 16950.772 14146.027 10218.195  4327.762
## 3  4480.442  6550.745  2128.591  3343.231  4744.891  9008.902 16566.488
## 4  3947.590  4763.308  3441.782  3895.635  3917.459  7783.010 14629.284
## 5 11918.806 10216.676  8079.426  6533.217  4360.082  1599.344 11738.598
## 6 17365.842 13362.096 17699.167 16626.792 13797.841 10320.568  1941.297
##          37        38        39        40        41       42        43
## 1  7310.432  7326.658 18693.024 17052.475  9216.172 9591.678 16447.780
## 2 11169.410  9010.612 21663.961 20236.601 12201.309 9685.999 17061.349
## 3 21520.704 21887.800  4252.664  2589.355  6024.360 9251.441  6821.133
## 4 20541.425 20682.864  5370.308  3814.958  4399.108 7293.487  5507.135
## 5 12434.681 12890.795 13302.173 11640.855  4656.425 7210.811 12230.841
## 6 13464.308 11407.250 20595.082 19268.559 11574.386 8436.732 15479.262
##          44        45       46        47        48        49       50        51
## 1 12167.780 12864.611 24893.34  8126.789  8029.358  6126.225 19514.07  5765.852
## 2 14073.165 10463.437 27475.18  3325.761  5837.413  8786.664 23217.05 10303.986
## 3  4902.646 11960.543 10427.16 15604.232 21847.462 20678.827  5062.72 19951.392
## 4  2899.147  9969.387 11524.53 13740.085 20360.102 19515.826  6767.35 18966.658
## 5  7707.081 11276.874 19504.19  9985.442 13412.745 11649.567 13956.25 10863.683
## 6 12980.763  8439.059 26189.15  1857.735  8239.931 11114.798 22337.77 12499.636
##          52        53       54        55        56        57        58
## 1 11679.087  9233.475 13243.28 17076.158 18355.108  3686.609  7202.770
## 2 15799.089  5367.385 14620.05 20982.914 22352.495  8714.616  2565.659
## 3  2911.251 14351.424 27703.78  2770.106  4140.094 18090.684 15811.819
## 4  2524.050 12415.609 26628.99  4671.105  6021.546 17017.685 14006.643
## 5  6179.217  9931.571 18618.15 11511.395 12770.333  9015.515  9579.104
## 6 15292.963  3497.154 17044.85 20221.126 21586.799 10762.013  2080.835
##           59        60        61        62        63        64        65
## 1  7416.1320  5601.336  5072.309 15358.825 17240.941 12674.035  8525.508
## 2   987.6943  5708.989  4212.089 19773.905 21770.384 18034.857  3831.344
## 3 17288.3514 12525.373 14134.749  2294.532  3967.423  4441.266 21260.769
## 4 15512.7943 10813.768 12456.840  4214.638  5970.522  5401.318 19629.482
## 5 10607.2428  6248.975  7218.708  9757.647 11637.202  7148.253 13400.998
## 6  1922.6563  5467.400  4509.136 19214.432 21197.043 17868.027  6091.487
##          66        67        68       69        70        71        72       73
## 1  2890.633  6890.267  9283.037 10171.33 13179.095 13681.040 16362.606 6547.021
## 2  7352.519  4444.018 11103.870 13355.80 13365.405 12872.050 19104.010 8497.121
## 3 17450.546 20359.988 23805.174 24316.34  7575.528  9279.531  2605.938 9740.087
## 4 16258.275 18838.544 22678.438 23376.38  5642.612  7357.937  2968.255 8060.982
## 5  8477.692 12084.739 14741.629 15243.41  9676.209 10771.188 11116.664 4515.358
## 6  9352.760  6870.984 13520.695 15736.03 11831.441 11105.521 18037.154 8045.226
##         74        75        76       77        78        79        80        81
## 1 5477.602  8792.889  8375.531 11688.45  6552.281 22431.094  7147.704  6937.895
## 2 9719.567 10568.651  7310.633 11714.52  2334.191 25676.013  5262.819  7281.973
## 3 9277.683  7725.996 22600.104 26238.80 18879.751  7872.050 13153.185 10976.019
## 4 7959.360  5897.698 21196.823 24987.95 17255.920  9302.119 11317.270  9181.101
## 5 2113.114  5391.150 13922.524 17268.84 11110.206 16922.426  7787.294  6012.405
## 6 9778.670  9735.196  9727.017 14116.54  4747.271 24609.468  4340.428  6562.989
##          82        83        84        85        86        87        88
## 1  9109.443 11262.524  5025.802  3069.237  5141.847 14762.614 10701.092
## 2  6436.047  8394.215  7229.455  5531.175  3435.169 16526.953 12482.606
## 3 12997.438 12751.905 19541.525 17083.694 18108.682  4259.421  6118.397
## 4 11054.782 10751.321 18285.029 15724.882 16575.053  2805.980  4162.996
## 5  9049.855 10467.364 10630.006  8434.695  9982.813 10063.096  6626.992
## 6  4790.540  6400.968  9519.245  7505.041  5680.396 15272.247 11463.970
##          89        90        91        92        93        94        95
## 1  9930.334  1468.232  3409.547  2077.449 10817.510  5760.609  3843.634
## 2 10820.395  6586.141  4288.691  5896.384  3786.646  6118.352  6420.333
## 3  7944.341 15754.918 15848.321 14344.251 18463.427 20000.620 18255.478
## 4  5988.411 14488.083 14351.674 12949.496 16539.096 18631.415 16963.550
## 5  6779.638  6940.463  7728.609  5968.147 13175.317 11293.742  9424.844
## 6  9685.573  8254.974  5857.733  7091.437  1709.162  8494.312  8581.502
##          96        97        98        99       100       101       102
## 1  8469.775 11332.075  9302.936  7711.923  3109.422  3699.090  2632.593
## 2  1779.061 14490.442 14603.981 12758.532 10580.716 11200.643 10328.371
## 3 19899.084  3728.720  5832.138  6942.129 15426.500 12599.465 14207.845
## 4 18169.956  2092.955  5523.693  6118.584 14612.233 11873.108 13369.665
## 5 12653.878  6254.605  3713.136  2149.971  6434.373  3762.235  5191.083
## 6  3846.724 13749.944 14563.547 12749.031 12163.009 12285.596 11704.883
##         103       104       105       106       107
## 1  1658.928 10997.855  4305.937  1749.816 13167.463
## 2  8196.456  3357.981 10205.851  9266.090  5964.181
## 3 16053.193 19960.289 18108.243 15064.740 19986.073
## 4 14959.180 18071.176 17186.550 14093.767 18015.787
## 5  6998.660 14062.700  9037.259  5975.894 15338.540
## 6  9929.440  2533.324 12193.583 10804.259  4059.103

upz_final.nb <- poly2nb(upz_final, queen = TRUE)
head(upz_final.nb)

## [[1]]
## [1]   7  27  90  92 103 106
## 
## [[2]]
## [1]  28  59  78  96 104
## 
## [[3]]
## [1]  4 32 33 40 52 55 62 72
## 
## [[4]]
## [1]  3 30 52 72 87 97
## 
## [[5]]
## [1] 12 18 22 35 74 99
## 
## [[6]]
## [1]  21  36  47  58  59  93 104

Creamos una matriz de pesos binaria

UPLNOMBRE <- upz_final$UPLNOMBRE

n <- length(upz_final.nb)
MatW <- matrix(0, n, n, dimnames = list(UPLNOMBRE, UPLNOMBRE))

for (i in 1:n) {
  vecinos <- upz_final.nb[[i]]
  MatW[i, vecinos] <- 1
}

head(MatW)

##               SANTA ISABEL APOGEO EL PRADO NIZA GALERIAS GRAN BRITALIA
## SANTA ISABEL             0      0        0    0        0             0
## APOGEO                   0      0        0    0        0             0
## EL PRADO                 0      0        0    1        0             0
## NIZA                     0      0        1    0        0             0
## GALERIAS                 0      0        0    0        0             0
## GRAN BRITALIA            0      0        0    0        0             0
##               ZONA INDUSTRIAL EL REFUGIO FONTIBON DOCE DE OCTUBRE MODELIA
## SANTA ISABEL                1          0        0               0       0
## APOGEO                      0          0        0               0       0
## EL PRADO                    0          0        0               0       0
## NIZA                        0          0        0               0       0
## GALERIAS                    0          0        0               0       0
## GRAN BRITALIA               0          0        0               0       0
##               PARQUE SALITRE JARDIN BOTANICO ALFONSO LOPEZ COMUNEROS SAN RAFAEL
## SANTA ISABEL               0               0             0         0          0
## APOGEO                     0               0             0         0          0
## EL PRADO                   0               0             0         0          0
## NIZA                       0               0             0         0          0
## GALERIAS                   1               0             0         0          0
## GRAN BRITALIA              0               0             0         0          0
##               QUINTA PAREDES TEUSAQUILLO SAGRADO CORAZON EL PORVENIR
## SANTA ISABEL               0           0               0           0
## APOGEO                     0           0               0           0
## EL PRADO                   0           0               0           0
## NIZA                       0           0               0           0
## GALERIAS                   0           1               0           0
## GRAN BRITALIA              0           0               0           0
##               PATIO BONITO CHAPINERO CIUDAD SALITRE ORIENTAL PARDO RUBIO
## SANTA ISABEL             0         0                       0           0
## APOGEO                   0         0                       0           0
## EL PRADO                 0         0                       0           0
## NIZA                     0         0                       0           0
## GALERIAS                 0         1                       0           0
## GRAN BRITALIA            1         0                       0           0
##               SOSIEGO PUENTE ARANDA LA SABANA CARVAJAL LAS NIEVES SUBA BOLIVIA
## SANTA ISABEL        0             0         1        0          0    0       0
## APOGEO              0             0         0        1          0    0       0
## EL PRADO            0             0         0        0          0    0       0
## NIZA                0             0         0        0          0    1       0
## GALERIAS            0             0         0        0          0    0       0
## GRAN BRITALIA       0             0         0        0          0    0       0
##               COUNTRY CLUB SANTA BARBARA LOS ANDES LA ESMERALDA LAS MARGARITAS
## SANTA ISABEL             0             0         0            0              0
## APOGEO                   0             0         0            0              0
## EL PRADO                 1             1         0            0              0
## NIZA                     0             0         0            0              0
## GALERIAS                 0             0         0            1              0
## GRAN BRITALIA            0             0         0            0              1
##               LOS LIBERTADORES DANUBIO SAN JOSE DE BAVARIA BRITALIA LAS FERIAS
## SANTA ISABEL                 0       0                   0        0          0
## APOGEO                       0       0                   0        0          0
## EL PRADO                     0       0                   0        1          0
## NIZA                         0       0                   0        0          0
## GALERIAS                     0       0                   0        0          0
## GRAN BRITALIA                0       0                   0        0          0
##               CAPELLANIA TIBABUYES MINUTO DE DIOS FONTIBON SAN PABLO GUAYMARAL
## SANTA ISABEL           0         0              0                  0         0
## APOGEO                 0         0              0                  0         0
## EL PRADO               0         0              0                  0         0
## NIZA                   0         0              0                  0         0
## GALERIAS               0         0              0                  0         0
## GRAN BRITALIA          0         0              0                  0         0
##               CORABASTOS LUCERO DIANA TURBAY VERBENAL LA GLORIA LA ALHAMBRA
## SANTA ISABEL           0      0            0        0         0           0
## APOGEO                 0      0            0        0         0           0
## EL PRADO               0      0            0        0         0           1
## NIZA                   0      0            0        0         0           1
## GALERIAS               0      0            0        0         0           0
## GRAN BRITALIA          1      0            0        0         0           0
##               CALANDAIMA CIUDAD USME TOBERIN LA URIBE 20 DE JULIO
## SANTA ISABEL           0           0       0        0           0
## APOGEO                 0           0       0        0           0
## EL PRADO               0           0       1        0           0
## NIZA                   0           0       0        0           0
## GALERIAS               0           0       0        0           0
## GRAN BRITALIA          0           0       0        0           0
##               KENNEDY CENTRAL TIMIZA BAVARIA AMERICAS LOS CEDROS
## SANTA ISABEL                0      0       0        0          0
## APOGEO                      0      1       0        0          0
## EL PRADO                    0      0       0        0          1
## NIZA                        0      0       0        0          0
## GALERIAS                    0      0       0        0          0
## GRAN BRITALIA               1      1       0        0          0
##               SAN CRISTOBAL NORTE USAQUEN JERUSALEM SAN JOSE SAN FRANCISCO
## SANTA ISABEL                    0       0         0        0             0
## APOGEO                          0       0         0        0             0
## EL PRADO                        0       0         0        0             0
## NIZA                            0       0         0        0             0
## GALERIAS                        0       0         0        0             0
## GRAN BRITALIA                   0       0         0        0             0
##               GRAN YOMASA LA FLORA GARCES NAVAS ENGATIVA CASA BLANCA SUBA
## SANTA ISABEL            0        0            0        0                0
## APOGEO                  0        0            0        0                0
## EL PRADO                0        0            0        0                1
## NIZA                    0        0            0        0                1
## GALERIAS                0        0            0        0                0
## GRAN BRITALIA           0        0            0        0                0
##               CIUDAD SALITRE OCCIDENTAL PARQUE SIMON BOLIVAR - CAN
## SANTA ISABEL                          0                          0
## APOGEO                                0                          0
## EL PRADO                              0                          0
## NIZA                                  0                          0
## GALERIAS                              0                          1
## GRAN BRITALIA                         0                          0
##               SANTA CECILIA EL TESORO MONTE BLANCO ARBORIZADORA
## SANTA ISABEL              0         0            0            0
## APOGEO                    0         0            0            1
## EL PRADO                  0         0            0            0
## NIZA                      0         0            0            0
## GALERIAS                  0         0            0            0
## GRAN BRITALIA             0         0            0            0
##               PASEO DE LOS LIBERTADORES CASTILLA GRANJAS DE TECHO TINTAL NORTE
## SANTA ISABEL                          0        0                0            0
## APOGEO                                0        0                0            0
## EL PRADO                              0        0                0            0
## NIZA                                  0        0                0            0
## GALERIAS                              0        0                0            0
## GRAN BRITALIA                         0        0                0            0
##               ZONA FRANCA MARRUECOS QUIROGA VENECIA EL RINCON BOYACA REAL
## SANTA ISABEL            0         0       0       0         0           0
## APOGEO                  0         0       0       0         0           0
## EL PRADO                0         0       0       0         0           0
## NIZA                    0         0       0       0         1           0
## GALERIAS                0         0       0       0         0           0
## GRAN BRITALIA           0         0       0       0         0           0
##               ALAMOS RESTREPO MUZU CIUDAD MONTES BOSA OCCIDENTAL TUNJUELITO
## SANTA ISABEL       0        1    0             1               0          0
## APOGEO             0        0    0             0               0          0
## EL PRADO           0        0    0             0               0          0
## NIZA               0        0    0             0               0          0
## GALERIAS           0        0    0             0               0          0
## GRAN BRITALIA      0        0    0             0               1          0
##               MARCO FIDEL SUAREZ ISMAEL PERDOMO LA FLORESTA CHICO LAGO
## SANTA ISABEL                   0              0           0          0
## APOGEO                         0              1           0          0
## EL PRADO                       0              0           0          0
## NIZA                           0              0           1          0
## GALERIAS                       0              0           0          0
## GRAN BRITALIA                  0              0           0          0
##               LOS ALCAZARES LOURDES LA MACARENA LA CANDELARIA CIUDAD JARDIN
## SANTA ISABEL              0       0           0             0             1
## APOGEO                    0       0           0             0             0
## EL PRADO                  0       0           0             0             0
## NIZA                      0       0           0             0             0
## GALERIAS                  1       0           0             0             0
## GRAN BRITALIA             0       0           0             0             0
##               BOSA CENTRAL SAN BLAS LAS CRUCES TINTAL SUR
## SANTA ISABEL             0        0          1          0
## APOGEO                   1        0          0          0
## EL PRADO                 0        0          0          0
## NIZA                     0        0          0          0
## GALERIAS                 0        0          0          0
## GRAN BRITALIA            1        0          0          0

image(MatW, main = "Matriz de Adyacencia (Conexiones entre UPZ)", axes = FALSE)

Creamos diferentes versiones de la matriz de pesos.

coords <- st_coordinates(st_centroid(upz_final))

## Warning: st_centroid assumes attributes are constant over geometries

Criterio:

queen: en todas las direcciones.
rock: en lineas rectas

nb_queen <- poly2nb(upz_final, queen = TRUE)
nb_rook <- poly2nb(upz_final, queen = FALSE)

par(mfrow = c(1, 2), mar = c(1,1,3,1))
plot(st_geometry(upz_final), border = "grey80", main = "Criterio Queen (Reina)")
plot(nb_queen, coords, add = TRUE, col = "red", lwd = 1, pch = 19, cex = 0.5)
plot(st_geometry(upz_final), border = "grey80", main = "Criterio Rook (Torre)")
plot(nb_rook, coords, add = TRUE, col = "blue", lwd = 1, pch = 19, cex = 0.5)

par(mfrow = c(1, 1))

Criterio

-4 vecinos mas cercanos.

nb_knn <- knn2nb(knearneigh(coords, k = 4))
plot(st_geometry(upz_final), border = "grey80")
plot(nb_knn, coords, add = TRUE, col = "#FA5939", pch = 19, cex = 0.6)
title("Red de Vecinos Cercanos (K=4)")

Criterio

Distancia maxima.

nb_dist <- dnearneigh(coords, d1 = 0, d2 = 2000)

## Warning in dnearneigh(coords, d1 = 0, d2 = 2000): neighbour object has 17
## sub-graphs

plot(st_geometry(upz_final), border="grey80")
plot(nb_dist, coords, add=TRUE, col = "#FA5939", pch =19, cex = 0.6)
title("Red de distancia 2000 m")

### Criterio - Delunay: Intenta formar los triangulo mas equilateros posibles. - Gabriel Graph: Malla de Delunay cortar lineas debiles o lejanas. - No existe ningún punto C tal que \(Dist(A,C) < Dist(A,B)\) Y \(Dist(B,C) < Dist(A,B)\).

nb_tri <- tri2nb(coords)
nb_gab <- graph2nb(gabrielneigh(coords))
nb_rel <- graph2nb(relativeneigh(coords))

par(mfrow = c(1, 3), mar = c(1, 1, 3, 1))
plot(st_geometry(upz_final), border = "grey90", main = "1. Delaunay\n(Triangulación)")
plot(nb_tri, coords, add = TRUE, col = "blue", pch = 19, cex = 0.4, lwd = 0.8)

plot(st_geometry(upz_final), border = "grey90", main = "2. Gabriel Graph\n(Círculos Vacíos)")
plot(nb_gab, coords, add = TRUE, col = "#009E73", pch = 19, cex = 0.4, lwd = 0.8)


plot(st_geometry(upz_final), border = "grey90", main = "3. Relative Neighbor\n(Vecinos Relativos)")
plot(nb_rel, coords, add = TRUE, col = "#D55E00", pch = 19, cex = 0.4, lwd = 0.8)

par(mfrow = c(1, 1))

Criterio

Busca los 4 vecinos mas cercanos
Calcula su distancia
Usamos el inverso de la distancia.

k_nb <- knn2nb(knearneigh(coords, k = 4))
dist_list <- nbdists(k_nb, coords)
idw_weights <- lapply(dist_list, function(x) 1/x)
final_idw <- nb2listw(k_nb, glist = idw_weights, style = "B")
plot(st_geometry(upz_final), border = "grey85")
plot(final_idw$neighbours, coords, add = TRUE, col = "red", pch = 19, cex = 0.5)

title("Red IDW (Basada en K=6)")

upz_final.lww <- nb2listw(upz_final.nb, style = "W", zero.policy = TRUE)
upz_final.lwb <- nb2listw(upz_final.nb, style = "B", zero.policy = TRUE)
upz_final.lwc <- nb2listw(upz_final.nb, style = "C", zero.policy = TRUE)
upz_final.lwu <- nb2listw(upz_final.nb, style = "U", zero.policy = TRUE)
par(mfrow = c(2, 2), mar = c(2, 2, 2, 1))
image(listw2mat(upz_final.lww), main = "Estilo W (Promedio)\nFilas suman 1", axes = FALSE)
image(listw2mat(upz_final.lwb), main = "Estilo B (Binario)\nTodos valen 1", axes = FALSE)
image(listw2mat(upz_final.lwc), main = "Estilo C (Global)\nSuma total = N", axes = FALSE)
image(listw2mat(upz_final.lwu), main = "Estilo U (Universal)\nSuma total = 1", axes = FALSE)

list_of_w_1 <- list(
  "Queen"     = nb_queen,   
  "Rook"      = nb_rook,
  "KNN"       = nb_knn,         
  "Distancia" = nb_dist,       
  "Delaunay"  = nb_tri,        
  "Gabriel"   = nb_gab,         
  "Relativo"  = nb_rel        
)

df_resumen_1 <- data.frame()
for (nombre in names(list_of_w_1)) {
  print(nombre)
  objeto_actual <- list_of_w_1[[nombre]]
  w_para_test <- nb2listw(objeto_actual, style = "W", zero.policy = TRUE)
  test <- moran.test(upz_final$proportion_vis, listw = w_para_test, zero.policy = TRUE)
  fila_nueva <- data.frame(
    Matriz  = nombre,
    Moran_I = as.numeric(test$statistic),
    P_Value = test$p.value
  )
  df_resumen_1 <- rbind(df_resumen_1, fila_nueva)
}

## [1] "Queen"
## [1] "Rook"
## [1] "KNN"
## [1] "Distancia"
## [1] "Delaunay"
## [1] "Gabriel"
## [1] "Relativo"

df_resumen_1 <- df_resumen_1[order(-df_resumen_1$Moran_I), ]
print(df_resumen_1)

##      Matriz   Moran_I      P_Value
## 5  Delaunay 10.707670 4.684499e-27
## 3       KNN  9.439857 1.866436e-21
## 1     Queen  9.022941 9.155288e-20
## 2      Rook  8.298035 5.292368e-17
## 6   Gabriel  6.081651 5.947550e-10
## 4 Distancia  5.498352 1.916788e-08
## 7  Relativo  4.327664 7.534953e-06

list_of_w_2 <-list(
  "Estilo W"  = upz_final.lww,
  "Estilo B"  = upz_final.lwb,
  "Estilo C"  = upz_final.lwc,
  "Estilo U"  = upz_final.lwu
)

df_resumen_2 <- data.frame()
for (nombre in names(list_of_w_2)) {
  print(nombre)
  objeto_actual <- list_of_w_2[[nombre]]
  test <- moran.test(upz_final$proportion_vis, listw = objeto_actual, zero.policy = TRUE)
  fila_nueva <- data.frame(
    Matriz  = nombre,
    Moran_I = as.numeric(test$statistic),
    P_Value = test$p.value
  )
  df_resumen_2 <- rbind(df_resumen_2, fila_nueva)
}

## [1] "Estilo W"
## [1] "Estilo B"
## [1] "Estilo C"
## [1] "Estilo U"

df_resumen_2 <- df_resumen_2[order(-df_resumen_2$Moran_I), ]
print(df_resumen_2)

##     Matriz  Moran_I      P_Value
## 4 Estilo U 9.436802 1.921653e-21
## 2 Estilo B 9.436802 1.921653e-21
## 3 Estilo C 9.436802 1.921653e-21
## 1 Estilo W 9.022941 9.155288e-20

df_resumen <- rbind(df_resumen_1, df_resumen_2)
print(df_resumen)

##       Matriz   Moran_I      P_Value
## 5   Delaunay 10.707670 4.684499e-27
## 3        KNN  9.439857 1.866436e-21
## 1      Queen  9.022941 9.155288e-20
## 2       Rook  8.298035 5.292368e-17
## 6    Gabriel  6.081651 5.947550e-10
## 4  Distancia  5.498352 1.916788e-08
## 7   Relativo  4.327664 7.534953e-06
## 41  Estilo U  9.436802 1.921653e-21
## 21  Estilo B  9.436802 1.921653e-21
## 31  Estilo C  9.436802 1.921653e-21
## 11  Estilo W  9.022941 9.155288e-20

p_value_min <- df_resumen[which.min(df_resumen$P_Value), ]
print(p_value_min)

##     Matriz  Moran_I      P_Value
## 5 Delaunay 10.70767 4.684499e-27

w_to_use <- nb2listw(nb_tri, style = "W", zero.policy = TRUE)

var_stand <- scale(upz_final$proportion_vis) 
lag_stand <- lag.listw(w_to_use, var_stand)

df_moran <- data.frame(
  x = as.numeric(var_stand), 
  wx = as.numeric(lag_stand), 
  nombre = upz_final$UPLNOMBRE
)

df_moran$cuadrante <- NA
df_moran$cuadrante[df_moran$x >= 0 & df_moran$wx >= 0] <- "High-High"
df_moran$cuadrante[df_moran$x < 0  & df_moran$wx < 0]  <- "Low-Low"
df_moran$cuadrante[df_moran$x < 0  & df_moran$wx >= 0] <- "Low-High"
df_moran$cuadrante[df_moran$x >= 0 & df_moran$wx < 0]  <- "High-Low"

ggplot(df_moran, aes(x = x, y = wx)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  
  geom_point(aes(color = cuadrante), size = 2.5, alpha = 0.8) +
  
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "black", size = 0.8) +
  geom_text_repel(
    data = subset(df_moran, abs(x) > 1 | abs(wx) > 1),
    aes(label = nombre), 
    size = 3, 
    box.padding = 0.5,
    max.overlaps = 20
  ) +
  
  scale_color_manual(values = c("High-High" = "red", 
                                "Low-Low" = "blue", 
                                "High-Low" = "lightblue", 
                                "Low-High" = "pink")) +
  
  labs(
    title = "Dispersograma de Moran: Vivienda VIS (Bogotá 2012-2023)",
    subtitle = "Identificación de clusters espaciales (Estandarizado)",
    x = "Vivienda VIS (Desviaciones estándar)",
    y = "Rezago Espacial (Vecinos)",
    color = "Cuadrante"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

- High - High: UPZ con mucha vivienda VIS rodeadas de vecinos con mucha vivienda VIS - Low - Low: UPZ con poca vivienda nueva VIS rodeada de vecinos con poca vivienda nueva VIS

nearng = dnearneigh(st_coordinates(upz_final)[, 1:2], 0, 1000)
nearng

## Neighbour list object:
## Number of regions: 75024 
## Number of nonzero links: 368076518 
## Percentage nonzero weights: 6.539397 
## Average number of links: 4906.117 
## 1 region with no links:
## 39320

upz_final.lw.g = nb2listw(nearng, style="B", zero.policy = T)

localG = localG(upz_final$sum_total_vis, upz_final.lwb)
head(localG)

## [1] -1.0064902  4.1662335 -1.5051042  0.7068309 -1.2682862  1.7454371

upz_final$Gi_score <- as.numeric(localG)

ggplot(upz_final) +
  geom_sf(aes(fill = Gi_score), color = "white", size = 0.1) +
  scale_fill_distiller(palette = "RdBu", direction = -1) + 
  
  labs(title = "Análisis de Puntos Calientes (Getis-Ord Gi*)",
       subtitle = "Rojo = Aglomeración alta (Hotspot) | Azul = Aglomeración baja (Coldspot)",
       fill = "Z-Score") +
  theme_void()

sim.G = matrix(0,1000,length(upz_final$UPLCODIGO))
for(i in 1:1000) sim.G[i,] = localG(sample(upz_final$sum_total_vis), upz_final.lwb)
mc.pvalor.G = (colSums(sweep(sim.G,2,localG,">="))+1)/(nrow(sim.G)+1)
mc.pvalor.G

##   [1] 0.860139860 0.004995005 0.986013986 0.207792208 0.942057942 0.069930070
##   [7] 0.831168831 0.988011988 0.129870130 0.951048951 0.703296703 0.959040959
##  [13] 0.991008991 0.224775225 0.608391608 0.886113886 0.868131868 0.815184815
##  [19] 0.723276723 0.027972028 0.040959041 0.932067932 0.867132867 0.934065934
##  [25] 0.661338661 0.979020979 0.874125874 0.321678322 0.908091908 0.224775225
##  [31] 0.203796204 0.992007992 0.998001998 0.991008991 0.897102897 0.227772228
##  [37] 0.954045954 0.263736264 0.664335664 0.761238761 0.993006993 0.822177822
##  [43] 0.063936064 0.422577423 0.063936064 0.952047952 0.397602398 0.901098901
##  [49] 0.168831169 0.511488511 0.272727273 0.984015984 0.158841159 0.448551449
##  [55] 0.795204795 0.594405594 0.432567433 0.493506494 0.010989011 0.448551449
##  [61] 0.709290709 0.974025974 0.674325674 1.000000000 0.039960040 0.546453546
##  [67] 0.237762238 0.413586414 0.248751249 0.515484515 0.274725275 0.257742258
##  [73] 0.925074925 0.975024975 0.979020979 0.562437562 0.353646354 0.406593407
##  [79] 0.440559441 0.188811189 0.262737263 0.001998002 0.272727273 0.650349650
##  [85] 0.798201798 0.191808192 0.197802198 0.762237762 0.263736264 0.935064935
##  [91] 0.826173826 0.931068931 0.012987013 0.062937063 0.610389610 0.003996004
##  [97] 0.760239760 0.997002997 0.954045954 0.571428571 0.652347652 0.656343656
## [103] 0.656343656 0.013986014 0.593406593 0.464535465 0.007992008

ggplot(upz_final) + 
  geom_sf(aes(fill = mc.pvalor.G), color = "white", size = 0.1) +
  scale_fill_viridis_c(
    option = "magma", 
    direction = -1, 
    name = "P-Valor",
    trans = "sqrt", 
    breaks = c(0, 0.01, 0.05, 0.1, 0.5, 1)
  ) +
  
  labs(
    title = "Significancia Estadística (P-Values Monte Carlo)",
    subtitle = "Zonas oscuras indican alta probabilidad de aglomeración no aleatoria",
    caption = "Nota: Valores cercanos a 0 indican alta significancia."
  ) +
  theme_void() +
  theme(legend.position = "right")

## Regresión espacial

reg.eq1=sum_total_vis ~ sum_total_novis + estrato_moda
reg1=lm(reg.eq1,data=upz_final)   
summary(reg1)

## 
## Call:
## lm(formula = reg.eq1, data = upz_final)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1588.5  -660.8  -244.3   210.0  5217.5 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      6.843e+02  3.227e+02   2.121   0.0364 *
## sum_total_novis  1.094e-01  6.595e-02   1.659   0.1002  
## estrato_moda2    6.864e+02  4.025e+02   1.705   0.0912 .
## estrato_moda3   -6.200e+01  3.947e+02  -0.157   0.8755  
## estrato_moda4   -6.774e+02  4.770e+02  -1.420   0.1587  
## estrato_moda5   -8.861e+02  6.577e+02  -1.347   0.1809  
## estrato_moda6   -1.008e+03  6.735e+02  -1.497   0.1375  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1250 on 100 degrees of freedom
## Multiple R-squared:  0.1356, Adjusted R-squared:  0.08376 
## F-statistic: 2.615 on 6 and 100 DF,  p-value: 0.02137

plot(reg1)

SLX (Spatial Lag of X)

reg2=lmSLX(reg.eq1,data=upz_final, w_to_use,, zero.policy = TRUE)
summary(reg2)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                      Estimate    Std. Error  t value     Pr(>|t|)  
## (Intercept)           8.995e+02   4.352e+02   2.067e+00   4.149e-02
## sum_total_novis       1.274e-01   6.606e-02   1.929e+00   5.672e-02
## estrato_moda2         1.195e+03   4.479e+02   2.669e+00   8.973e-03
## estrato_moda3         1.357e+03   5.219e+02   2.600e+00   1.083e-02
## estrato_moda4         1.083e+03   6.244e+02   1.735e+00   8.610e-02
## estrato_moda5         8.600e+02   7.730e+02   1.113e+00   2.688e-01
## estrato_moda6         5.759e+02   8.761e+02   6.573e-01   5.126e-01
## lag.sum_total_novis   6.229e-02   1.347e-01   4.626e-01   6.447e-01
## lag.estrato_moda2     1.386e+02   7.831e+02   1.770e-01   8.599e-01
## lag.estrato_moda3    -2.395e+03   7.329e+02  -3.267e+00   1.516e-03
## lag.estrato_moda4    -2.030e+03   1.033e+03  -1.964e+00   5.243e-02
## lag.estrato_moda5    -2.941e+03   1.605e+03  -1.833e+00   7.000e-02
## lag.estrato_moda6    -1.979e+03   1.704e+03  -1.161e+00   2.485e-01

plot(reg2)

Modelo Lag (SAR - Spatial Autoregressive)

w_regresion <- nb2listw(upz_final.nb, style = "W", zero.policy = TRUE)
reg3 <- lagsarlm(reg.eq1, 
                 data = upz_final, 
                 listw = w_regresion,
                 zero.policy = TRUE)
summary(reg3)

## 
## Call:lagsarlm(formula = reg.eq1, data = upz_final, listw = w_regresion, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1381.58  -647.29  -329.40   201.83  4870.37 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                    Estimate  Std. Error z value Pr(>|z|)
## (Intercept)      235.372947  306.350450  0.7683  0.44230
## sum_total_novis    0.099759    0.060338  1.6533  0.09826
## estrato_moda2    585.363491  370.228753  1.5811  0.11386
## estrato_moda3    213.182871  361.310722  0.5900  0.55517
## estrato_moda4   -378.767691  436.435230 -0.8679  0.38547
## estrato_moda5   -469.768348  602.672998 -0.7795  0.43570
## estrato_moda6   -539.134461  616.291044 -0.8748  0.38168
## 
## Rho: 0.38429, LR test value: 8.3138, p-value: 0.0039344
## Asymptotic standard error: 0.11704
##     z-value: 3.2834, p-value: 0.0010258
## Wald statistic: 10.78, p-value: 0.0010258
## 
## Log likelihood: -907.0181 for lag model
## ML residual variance (sigma squared): 1306600, (sigma: 1143.1)
## Number of observations: 107 
## Number of parameters estimated: 9 
## AIC: 1832, (AIC for lm: 1838.3)
## LM test for residual autocorrelation
## test value: 2.6712, p-value: 0.10218

Modelo de Error Espacial (SEM - Spatial Error Model)

reg4 = errorsarlm(reg.eq1, data = upz_final, w_regresion, zero.policy = TRUE)
summary(reg4)

## 
## Call:errorsarlm(formula = reg.eq1, data = upz_final, listw = w_regresion, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1299.071  -677.093  -298.695    97.948  4765.600 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)     207.929807 419.569361  0.4956  0.62019
## sum_total_novis   0.083739   0.060980  1.3732  0.16968
## estrato_moda2   691.479247 429.236239  1.6110  0.10719
## estrato_moda3   792.930950 462.412203  1.7148  0.08639
## estrato_moda4   142.854006 531.958441  0.2685  0.78828
## estrato_moda5    88.991028 679.180946  0.1310  0.89575
## estrato_moda6   118.973646 699.248732  0.1701  0.86490
## 
## Lambda: 0.50625, LR test value: 7.0914, p-value: 0.0077457
## Asymptotic standard error: 0.1067
##     z-value: 4.7446, p-value: 2.0887e-06
## Wald statistic: 22.512, p-value: 2.0887e-06
## 
## Log likelihood: -907.6293 for error model
## ML residual variance (sigma squared): 1285700, (sigma: 1133.9)
## Number of observations: 107 
## Number of parameters estimated: 9 
## AIC: 1833.3, (AIC for lm: 1838.3)

Modelo de Durbin Espacial con Error (SDEM - Spatial Durbin Error Model)

reg5 = errorsarlm(reg.eq1, data = upz_final, w_regresion, etype = "emixed", zero.policy = TRUE)
summary(reg5)

## 
## Call:errorsarlm(formula = reg.eq1, data = upz_final, listw = w_regresion, 
##     etype = "emixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1725.24  -680.63  -131.15   322.34  3881.53 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                        Estimate  Std. Error z value Pr(>|z|)
## (Intercept)          7.1123e+02  4.7604e+02  1.4940 0.135166
## sum_total_novis      1.2715e-01  5.8301e-02  2.1810 0.029182
## estrato_moda2        1.0727e+03  4.0683e+02  2.6368 0.008370
## estrato_moda3        1.1535e+03  4.6044e+02  2.5051 0.012241
## estrato_moda4        6.8570e+02  5.4092e+02  1.2677 0.204919
## estrato_moda5        7.8873e+02  7.1835e+02  1.0980 0.272213
## estrato_moda6       -8.1603e+01  7.7283e+02 -0.1056 0.915907
## lag.sum_total_novis  1.7805e-01  1.2289e-01  1.4488 0.147396
## lag.estrato_moda2    5.0139e+02  6.7654e+02  0.7411 0.458629
## lag.estrato_moda3   -2.1120e+03  6.7632e+02 -3.1228 0.001791
## lag.estrato_moda4   -2.0747e+03  1.0267e+03 -2.0207 0.043306
## lag.estrato_moda5   -1.3530e+03  1.7052e+03 -0.7935 0.427494
## lag.estrato_moda6   -2.5488e+03  1.3313e+03 -1.9146 0.055548
## 
## Lambda: 0.30842, LR test value: 4.3977, p-value: 0.035988
## Asymptotic standard error: 0.12818
##     z-value: 2.4061, p-value: 0.016125
## Wald statistic: 5.7892, p-value: 0.016125
## 
## Log likelihood: -896.8688 for error model
## ML residual variance (sigma squared): 1094200, (sigma: 1046)
## Number of observations: 107 
## Number of parameters estimated: 15 
## AIC: 1823.7, (AIC for lm: 1826.1)

Modelo Durbin Espacial (SDM - Spatial Durbin Model)

reg6 = lagsarlm(reg.eq1, data = upz_final, w_regresion, type = "mixed", zero.policy = TRUE)
summary(reg6)

## 
## Call:lagsarlm(formula = reg.eq1, data = upz_final, listw = w_regresion, 
##     type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1737.36  -668.37  -101.07   347.60  3850.38 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                        Estimate  Std. Error z value Pr(>|z|)
## (Intercept)          4.1818e+02  3.6547e+02  1.1442 0.252530
## sum_total_novis      1.1612e-01  5.7878e-02  2.0062 0.044831
## estrato_moda2        1.0323e+03  4.2471e+02  2.4306 0.015074
## estrato_moda3        1.2686e+03  4.9186e+02  2.5793 0.009900
## estrato_moda4        7.9848e+02  5.6870e+02  1.4040 0.160309
## estrato_moda5        9.3495e+02  7.3165e+02  1.2779 0.201297
## estrato_moda6        2.2310e+02  7.9076e+02  0.2821 0.777840
## lag.sum_total_novis  1.2523e-01  1.0847e-01  1.1545 0.248312
## lag.estrato_moda2    3.4066e+02  6.3781e+02  0.5341 0.593268
## lag.estrato_moda3   -2.0051e+03  6.3995e+02 -3.1333 0.001729
## lag.estrato_moda4   -1.6509e+03  9.5435e+02 -1.7298 0.083660
## lag.estrato_moda5   -1.2212e+03  1.5587e+03 -0.7835 0.433358
## lag.estrato_moda6   -2.3867e+03  1.2344e+03 -1.9335 0.053180
## 
## Rho: 0.29345, LR test value: 4.5209, p-value: 0.033484
## Asymptotic standard error: 0.12403
##     z-value: 2.3659, p-value: 0.017986
## Wald statistic: 5.5975, p-value: 0.017986
## 
## Log likelihood: -896.8072 for mixed model
## ML residual variance (sigma squared): 1095200, (sigma: 1046.5)
## Number of observations: 107 
## Number of parameters estimated: 15 
## AIC: 1823.6, (AIC for lm: 1826.1)
## LM test for residual autocorrelation
## test value: 0.5426, p-value: 0.46136

Modelo Manski (Manski Model)

reg7 = sacsarlm(reg.eq1, data = upz_final, w_regresion, type = "sacmixed", zero.policy = TRUE)
summary(reg7)

## 
## Call:sacsarlm(formula = reg.eq1, data = upz_final, listw = w_regresion, 
##     type = "sacmixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1757.07  -658.55  -106.84   329.39  3872.49 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                        Estimate  Std. Error z value Pr(>|z|)
## (Intercept)          5.2738e+02  5.3636e+02  0.9832 0.325487
## sum_total_novis      1.2104e-01  6.0541e-02  1.9993 0.045580
## estrato_moda2        1.0633e+03  4.1505e+02  2.5618 0.010414
## estrato_moda3        1.2451e+03  5.4105e+02  2.3013 0.021374
## estrato_moda4        7.7842e+02  6.1030e+02  1.2755 0.202144
## estrato_moda5        8.9947e+02  7.6637e+02  1.1737 0.240525
## estrato_moda6        9.9309e+01  9.4277e+02  0.1053 0.916108
## lag.sum_total_novis  1.4781e-01  1.4292e-01  1.0342 0.301043
## lag.estrato_moda2    4.0703e+02  7.6492e+02  0.5321 0.594640
## lag.estrato_moda3   -2.0692e+03  6.5153e+02 -3.1760 0.001493
## lag.estrato_moda4   -1.8531e+03  1.0565e+03 -1.7541 0.079421
## lag.estrato_moda5   -1.2701e+03  1.6252e+03 -0.7815 0.434490
## lag.estrato_moda6   -2.4579e+03  1.3513e+03 -1.8190 0.068918
## 
## Rho: 0.17722
## Asymptotic standard error: 0.43979
##     z-value: 0.40297, p-value: 0.68697
## Lambda: 0.15187
## Asymptotic standard error: 0.46704
##     z-value: 0.32519, p-value: 0.74504
## 
## LR test value: 28.975, p-value: 0.00032035
## 
## Log likelihood: -896.6876 for sacmixed model
## ML residual variance (sigma squared): 1100700, (sigma: 1049.1)
## Number of observations: 107 
## Number of parameters estimated: 16 
## AIC: 1825.4, (AIC for lm: 1838.3)

Modelo SARAR (Kelejian-Prucha, Cliff-Ord, SAC)

reg8 = sacsarlm(reg.eq1, data = upz_final, w_regresion, type = "sac", zero.policy = TRUE)
summary(reg8)

## 
## Call:sacsarlm(formula = reg.eq1, data = upz_final, listw = w_regresion, 
##     type = "sac", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1439.48  -653.94  -246.98   279.98  4501.25 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                    Estimate  Std. Error z value Pr(>|z|)
## (Intercept)      143.191293  242.092868  0.5915  0.55420
## sum_total_novis    0.106809    0.053364  2.0015  0.04534
## estrato_moda2    523.400771  304.820782  1.7171  0.08597
## estrato_moda3      3.853743  272.299886  0.0142  0.98871
## estrato_moda4   -421.650447  359.522614 -1.1728  0.24087
## estrato_moda5   -446.760397  542.751334 -0.8231  0.41043
## estrato_moda6   -710.704681  532.707963 -1.3341  0.18216
## 
## Rho: 0.62399
## Asymptotic standard error: 0.15266
##     z-value: 4.0874, p-value: 4.3631e-05
## Lambda: -0.46284
## Asymptotic standard error: 0.27172
##     z-value: -1.7034, p-value: 0.088502
## 
## LR test value: 10.524, p-value: 0.005185
## 
## Log likelihood: -905.913 for sac model
## ML residual variance (sigma squared): 1151200, (sigma: 1072.9)
## Number of observations: 107 
## Number of parameters estimated: 10 
## AIC: 1831.8, (AIC for lm: 1838.3)

reg.eq2=sum_total_vis ~ sum_total_novis + estrato_moda
reg4=errorsarlm(reg.eq2,data=upz_final, w_regresion, zero.policy = TRUE)

s = summary
s(reg4)

## 
## Call:errorsarlm(formula = reg.eq2, data = upz_final, listw = w_regresion, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1299.071  -677.093  -298.695    97.948  4765.600 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)     207.929807 419.569361  0.4956  0.62019
## sum_total_novis   0.083739   0.060980  1.3732  0.16968
## estrato_moda2   691.479247 429.236239  1.6110  0.10719
## estrato_moda3   792.930950 462.412203  1.7148  0.08639
## estrato_moda4   142.854006 531.958441  0.2685  0.78828
## estrato_moda5    88.991028 679.180946  0.1310  0.89575
## estrato_moda6   118.973646 699.248732  0.1701  0.86490
## 
## Lambda: 0.50625, LR test value: 7.0914, p-value: 0.0077457
## Asymptotic standard error: 0.1067
##     z-value: 4.7446, p-value: 2.0887e-06
## Wald statistic: 22.512, p-value: 2.0887e-06
## 
## Log likelihood: -907.6293 for error model
## ML residual variance (sigma squared): 1285700, (sigma: 1133.9)
## Number of observations: 107 
## Number of parameters estimated: 9 
## AIC: 1833.3, (AIC for lm: 1838.3)

reg.eq4=sum_total_vis ~ sum_total_novis
reg4=errorsarlm(reg.eq4,data=upz_final, w_regresion, zero.policy = TRUE)

s(reg4)

## 
## Call:errorsarlm(formula = reg.eq4, data = upz_final, listw = w_regresion, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1391.230  -660.065  -372.106    54.759  4971.520 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)     702.796916 231.878539  3.0309 0.002438
## sum_total_novis   0.074748   0.059627  1.2536 0.209992
## 
## Lambda: 0.48749, LR test value: 16.028, p-value: 6.2411e-05
## Asymptotic standard error: 0.10903
##     z-value: 4.4711, p-value: 7.7819e-06
## Wald statistic: 19.991, p-value: 7.7819e-06
## 
## Log likelihood: -910.946 for error model
## ML residual variance (sigma squared): 1374700, (sigma: 1172.5)
## Number of observations: 107 
## Number of parameters estimated: 4 
## AIC: 1829.9, (AIC for lm: 1843.9)

Mapa estimado

upz_final$valores_ajustados <- as.numeric(reg2$fitted.values)
ggplot(upz_final) +
  geom_sf(aes(fill = valores_ajustados), color = "gray80", size = 0.05) +
  scale_fill_viridis_c(
    option = "magma", 
    direction = -1,
    name = "Vivienda Estimada",
  ) +
  labs(
    title = "Modelo SEM: Valores Ajustados de Vivienda VIS",
    subtitle = "Estimación espacial del total de vivienda nueva (Bogotá 2012-2023)",
    caption = "Fuente: Estimación propia mediante Modelo de Error Espacial (SEM)"
  ) +
  theme_void() +
  theme(
    plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
    plot.subtitle = element_text(size = 10, color = "gray40", hjust = 0.5),
    legend.position = "right"
  )

moran.test(reg4$residuals, w_regresion)

## 
##  Moran I test under randomisation
## 
## data:  reg4$residuals  
## weights: w_regresion    
## 
## Moran I statistic standard deviate = -0.061747, p-value = 0.5246
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.013089970      -0.009433962       0.003505728

n <- nrow(upz_final)
upz_final$prop_vis_trans <- (upz_final$proportion_vis * (n - 1) + 0.5) / n
summary(upz_final$prop_vis_trans)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.004673 0.056186 0.280070 0.408739 0.817688 0.995327

variables_modelo <- c("proportion_vis", "sum_total_novis", "estrato_moda")

upz_final_clean <- upz_final %>%
  st_drop_geometry() %>%
  dplyr::select(all_of(variables_modelo)) %>%
  na.omit()
N <- nrow(upz_final_clean)
upz_final_clean$proportion_vis_star <- (upz_final_clean$proportion_vis * (N - 1) + 0.5) / N
formula_beta_star <- proportion_vis_star ~ sum_total_novis + estrato_moda

reg_beta <- gamlss(
  formula = formula_beta_star,
  sigma.formula = ~1,
  family = BE,
  data = upz_final_clean
)

## GAMLSS-RS iteration 1: Global Deviance = -169.7498 
## GAMLSS-RS iteration 2: Global Deviance = -171.8445 
## GAMLSS-RS iteration 3: Global Deviance = -172.016 
## GAMLSS-RS iteration 4: Global Deviance = -172.0321 
## GAMLSS-RS iteration 5: Global Deviance = -172.0337 
## GAMLSS-RS iteration 6: Global Deviance = -172.0339

print(summary(reg_beta))

## ******************************************************************
## Family:  c("BE", "Beta") 
## 
## Call:  gamlss(formula = formula_beta_star, sigma.formula = ~1,  
##     family = BE, data = upz_final_clean) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.340e+00  2.972e-01   4.510 1.78e-05 ***
## sum_total_novis -2.784e-05  1.403e-04  -0.198    0.843    
## estrato_moda2   -2.706e-01  3.517e-01  -0.769    0.444    
## estrato_moda3   -2.368e+00  3.592e-01  -6.592 2.13e-09 ***
## estrato_moda4   -2.669e+00  4.212e-01  -6.337 6.98e-09 ***
## estrato_moda5   -3.332e+00  5.898e-01  -5.650 1.54e-07 ***
## estrato_moda6   -3.428e+00  5.857e-01  -5.852 6.31e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  logit
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.2823     0.1278   2.208   0.0295 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  107 
## Degrees of Freedom for the fit:  8
##       Residual Deg. of Freedom:  99 
##                       at cycle:  6 
##  
## Global Deviance:     -172.0339 
##             AIC:     -156.0339 
##             SBC:     -134.6512 
## ******************************************************************
##                      Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept)      1.340438e+00 0.2971930246  4.5103289 1.781573e-05
## sum_total_novis -2.783521e-05 0.0001402714 -0.1984383 8.431091e-01
## estrato_moda2   -2.705988e-01 0.3517347530 -0.7693262 4.435313e-01
## estrato_moda3   -2.367954e+00 0.3592126976 -6.5920653 2.130021e-09
## estrato_moda4   -2.668660e+00 0.4211527561 -6.3365615 6.983871e-09
## estrato_moda5   -3.332438e+00 0.5898191699 -5.6499312 1.544878e-07
## estrato_moda6   -3.427552e+00 0.5856988691 -5.8520713 6.307030e-08
## (Intercept)      2.822575e-01 0.1278122027  2.2083769 2.952525e-02

Métodos regresión espacial

Yeison Buitrago

2025-12-06

Carga de librerias.

Vivienda Habilitada en Bogotá

Carda de datos.

Criterio:

Criterio

Criterio

Criterio

SLX (Spatial Lag of X)

Modelo Lag (SAR - Spatial Autoregressive)

Modelo de Error Espacial (SEM - Spatial Error Model)

Modelo de Durbin Espacial con Error (SDEM - Spatial Durbin Error Model)

Modelo Durbin Espacial (SDM - Spatial Durbin Model)

Modelo Manski (Manski Model)

Modelo SARAR (Kelejian-Prucha, Cliff-Ord, SAC)

Mapa estimado