GWR

1 Libraries and raw data

library(tidyverse)
library(magrittr)
library(sf)
library(tmap)

library(spdep) # Spatial weights, Moran's I, and LISA
library(spatialreg) # Spatial Regression

library(spgwr) # GWR

load("../04_MexicoCity_HTS/Mexico_HTS.RData")
ls()

[1] "MexicoCity_HTS"

1.1 Remove districts

summary(MexicoCity_HTS$TotalPeople)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    634     968    1035    1027    1103    1333       1 

MexicoCity_HTS %>% filter(is.na(TotalPeople)) %>% select(Distrito)

Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -99.09278 ymin: 19.41697 xmax: -99.04731 ymax: 19.45018
Geodetic CRS:  WGS 84
  Distrito                       geometry
1      034 MULTIPOLYGON (((-99.08461 1...

MexicoCity_HTS %<>% filter(Distrito != "034")

summary(MexicoCity_HTS$Total_Homes)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  272.0   285.0   293.0   293.7   301.0   352.0 

1.2 Map preparation

MexicoCity_Fringe <- st_read("../04_MexicoCity_HTS/UrbanFringe_MexicoCity.shp")

Reading layer `UrbanFringe_MexicoCity' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\04_MexicoCity_HTS\UrbanFringe_MexicoCity.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -99.36492 ymin: 19.04824 xmax: -98.9403 ymax: 19.59276
Geodetic CRS:  WGS 84

BB_Districts <- st_bbox(MexicoCity_HTS)
BB_Districts[1] <- -99.8
BB_Districts

     xmin      ymin      xmax      ymax 
-99.80000  18.93534 -98.59687  20.06826 

2 GWR

MexicoCity_HTS %<>%
  mutate(CarTrips_PerPopulation = Trips_Automovil/TotalPeople,
         Prop_Male = Total_Male/TotalPeople,
         Prop_HomesWithCars = Homes_With_Cars/Total_Homes,
         Prop_HomesWithBicycles = Homes_With_Bicycles/Total_Homes,
         Prop_Education_Low = TotalEducation_Low/TotalPeople,
         Prop_Education_Medium = TotalEducation_Medium/TotalPeople,
         Prop_Education_High = TotalEducation_High/TotalPeople)

ModelEquation <- "CarTrips_PerPopulation ~ 
                          Prop_Male + 
                          Prop_HomesWithCars + 
                          Prop_HomesWithBicycles +
                          # Prop_Education_Low + 
                          Age"

st_is_longlat(MexicoCity_HTS)

[1] TRUE

Data_SP <- as(MexicoCity_HTS, "Spatial")
Data_Centroids <- st_centroid(MexicoCity_HTS)

Warning: st_centroid assumes attributes are constant over geometries of x

Data_Centroids_SP <- as(Data_Centroids, "Spatial")

Data_Centroids_SP <- as(Data_Centroids, "Spatial")

Bandwidth = gwr.sel(ModelEquation, 
                    data = Data_SP, adapt = T)

Adaptive q: 0.381966 CV score: 15.79799 
Adaptive q: 0.618034 CV score: 16.06607 
Adaptive q: 0.236068 CV score: 15.58982 
Adaptive q: 0.145898 CV score: 15.57813 
Adaptive q: 0.1791799 CV score: 15.52235 
Adaptive q: 0.1889397 CV score: 15.52933 
Adaptive q: 0.1776273 CV score: 15.52027 
Adaptive q: 0.1655078 CV score: 15.54149 
Adaptive q: 0.1729981 CV score: 15.51867 
Adaptive q: 0.1742408 CV score: 15.51841 
Adaptive q: 0.1742815 CV score: 15.5184 
Adaptive q: 0.1755595 CV score: 15.51829 
Adaptive q: 0.1759044 CV score: 15.51828 
Adaptive q: 0.1759451 CV score: 15.51828 
Adaptive q: 0.1758637 CV score: 15.51828 
Adaptive q: 0.1759044 CV score: 15.51828 

Model_GWR = gwr(ModelEquation, 
                data = Data_SP, 
                adapt = Bandwidth, 
                hatmatrix = T, 
                se.fit = T)

Warning: CRS object has comment, which is lost in output; in tests, see
https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html

Model_GWR

Call:
gwr(formula = ModelEquation, data = Data_SP, adapt = Bandwidth, 
    hatmatrix = T, se.fit = T)
Kernel function: gwr.Gauss 
Adaptive quantile: 0.1759044 (about 33 of 193 data points)
Summary of GWR coefficient estimates at data points:
                            Min.   1st Qu.    Median   3rd Qu.      Max.  Global
X.Intercept.           -4.596397 -3.308457 -2.509538 -1.814455 -1.476676 -1.9269
Prop_Male              -0.610742  0.462755  0.909108  1.463229  2.486808 -0.1161
Prop_HomesWithCars      1.456090  1.636658  1.906816  2.131098  2.370322  1.8743
Prop_HomesWithBicycles -0.711226 -0.435152 -0.222469 -0.098109 -0.011571 -0.0213
Age                     0.025141  0.042468  0.055391  0.069582  0.095840  0.0513
Number of data points: 193 
Effective number of parameters (residual: 2traceS - traceS'S): 19.61997 
Effective degrees of freedom (residual: 2traceS - traceS'S): 173.38 
Sigma (residual: 2traceS - traceS'S): 0.2714758 
Effective number of parameters (model: traceS): 14.23004 
Effective degrees of freedom (model: traceS): 178.77 
Sigma (model: traceS): 0.2673519 
Sigma (ML): 0.2573072 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 56.97799 
AIC (GWR p. 96, eq. 4.22): 37.95128 
Residual sum of squares: 12.77795 
Quasi-global R2: 0.7102361 

results <- as.data.frame(Model_GWR$SDF)
names(results)

 [1] "sum.w"                         "X.Intercept."                 
 [3] "Prop_Male"                     "Prop_HomesWithCars"           
 [5] "Prop_HomesWithBicycles"        "Age"                          
 [7] "X.Intercept._se"               "Prop_Male_se"                 
 [9] "Prop_HomesWithCars_se"         "Prop_HomesWithBicycles_se"    
[11] "Age_se"                        "gwr.e"                        
[13] "pred"                          "pred.se"                      
[15] "localR2"                       "X.Intercept._se_EDF"          
[17] "Prop_Male_se_EDF"              "Prop_HomesWithCars_se_EDF"    
[19] "Prop_HomesWithBicycles_se_EDF" "Age_se_EDF"                   
[21] "pred.se.1"                    

2.1 Homes with cars

Temp <- results %>% 
  select(Variable = Prop_HomesWithCars, SE = Prop_HomesWithCars_se) %>% 
  mutate(T_value = Variable/SE) %>% 
  mutate(Significance = cut(T_value,
                            breaks = c(-0.01 + min(T_value), 
                                       -1.96, 1.96, 
                                       0.01 + max(T_value)),
                            labels = c("sig","nonsig", "sig")))

Temp_GWR_Coefficient <- Data_SP

Temp_GWR_Coefficient$Variable <- Temp$Variable
Temp_GWR_Coefficient$SE <- Temp$SE
Temp_GWR_Coefficient$T_value <- Temp$T_value
Temp_GWR_Coefficient$Significance <- Temp$Significance

Temp_GWR_Coefficient <- as(Temp_GWR_Coefficient,"sf")

Temp_GWR_Coefficient %<>% select(Variable, SE, T_value, Significance) 

Temp_GWR_Coefficient_SIGNIGICANT <- Temp_GWR_Coefficient %>% 
  filter(Significance == "sig")
Temp_GWR_Coefficient_NON_SIGNIGICANT <- Temp_GWR_Coefficient %>% 
  filter(Significance != "sig")

summary(Temp_GWR_Coefficient_SIGNIGICANT$Variable)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.456   1.637   1.907   1.893   2.131   2.370 

Q <- quantile(Temp_GWR_Coefficient_SIGNIGICANT$Variable, 
         probs = c(0, 0.2, 0.4, 0.6, 0.8, 1))
Q

      0%      20%      40%      60%      80%     100% 
1.456090 1.596206 1.755150 2.015792 2.165692 2.370322 

Breaks <- c(
  as.numeric(Q[1]), as.numeric(Q[2]),
  as.numeric(Q[3]), as.numeric(Q[4]),
  as.numeric(Q[5]), as.numeric(Q[6]))

Labels <- c(
  paste(floor(Breaks[1]*1000)/1000, round(Breaks[2], 3), sep = " - "),
  paste(round(Breaks[2], 3), round(Breaks[3], 3), sep = " - "),
  paste(round(Breaks[3], 3), round(Breaks[4], 3), sep = " - "),
  paste(round(Breaks[4], 3), round(Breaks[5], 3), sep = " - "),
  paste(round(Breaks[5], 3), ceiling(Breaks[6]*1000)/1000, sep = " - "))

MyPalette <- c("#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494")

GWR_Map_Cars <- tm_shape(MexicoCity_Fringe, bbox = BB_Districts) + 
  tm_polygons(lwd = 2.5, border.col = "black") +
  tm_shape(Temp_GWR_Coefficient_SIGNIGICANT) + 
  tm_polygons("Variable",
          palette = MyPalette, lwd = 0, border.col = "black", border.alpha = 0,
          title = "Theft (quantiles)", 
          breaks = Breaks, labels = Labels) +
  # tm_shape(Temp_GWR_Coefficient_NON_SIGNIGICANT) +
  # tm_polygons(col = "white") +
  tm_compass(size = 2, type = "arrow", position = c(0.85,0.87)) +
  tm_scale_bar(position = c(0.2,0.02)) +
  tm_layout(title = "Car ownership", legend.position = c(0.01, 0.58)) +
  tm_add_legend(
    type = c("fill"),
    labels = c("Not significant"), 
    col = c("white"),
    title = "")


GWR_Map_Cars

GWR_Map_Cars <- tm_shape(Temp_GWR_Coefficient_SIGNIGICANT) + 
  tm_polygons("Variable",
          palette = MyPalette, lwd = 0, border.col = "black", border.alpha = 0,
          title = "Quantiles", 
          breaks = Breaks, labels = Labels) +
  tm_borders(alpha = 0.4) +
  tm_shape(MexicoCity_Fringe) + 
  tm_polygons(alpha = 0, border.col = "black", lwd = 1.4, lty = "solid") +
  tm_compass(size = 2, type = "arrow", position = c(0.85,0.87)) + 
  tm_scale_bar(position = c(0.4,0.02)) +
  tm_layout(title = "Cars", legend.position = c(0.01,0.15), scale = 1)


GWR_Map_Cars