Spatial Autocorrelation

1 Libraries and raw data

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1── Conflicts ──────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(magrittr)

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract

library(sf)

Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

library(tmap)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

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

Loading required package: sp
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package
with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

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

[1] "MexicoCity_HTS"

The following would be the alternative:

# MexicoCity_HTS_Database <- read_delim("../03_MexicoCity_HTS/MexicoCity_HTS_Database.csv",
#                                       delim = ";")
# MexicoCity_Districts <- st_read("../03_MexicoCity_HTS/MexicoCity_Districts.shp")
# 
# 
# MexicoCity_HTS <- MexicoCity_Districts %>% 
#   left_join(MexicoCity_HTS_Database)

2 Neighbours

2.1 With polygons

from the “spdep” we will use the following functions:

• poly2nb(): to get a neighbours list based contiguous boundaries
• nb2listw(): to get the spatial weights

Neigh_ <- poly2nb(MexicoCity_HTS) # arugment queen = TRUE
WM_ <- nb2listw(Neigh_, style = "W")
# see https://r-spatial.github.io/spdep/reference/nb2listw.html for the style argument

class(Neigh_)

[1] "nb"

Neigh_

Neighbour list object:
Number of regions: 194 
Number of nonzero links: 1032 
Percentage nonzero weights: 2.742055 
Average number of links: 5.319588 

dim(MexicoCity_HTS)

[1] 194  58

length(Neigh_)

[1] 194

Neigh_[1]

[[1]]
[1] 2 3 4 7 8

Neigh_[5]

[[1]]
[1]  4 32 33

Important: Check this results with QGIS

2.2 With points

from the “spdep” we will use the following functions:

• knearneigh(): K nearest neighbours for spatial weights
• knn2nb(): Converts output from knearneigh() to a neighbours list of class nb
• nb2listw(): to get the spatial weights

Centroids_MexicoCity_HTS <- st_centroid(MexicoCity_HTS) %>% st_geometry()

Warning: st_centroid assumes attributes are constant over geometries of x

The knearneigh() function returns an intermediate form converted to an nb object by knn2nb(); knearneigh() can also take a longlat= argument to handle geographical coordinates

st_is_longlat(Centroids_MexicoCity_HTS)

[1] TRUE

KNN <- knearneigh(Centroids_MexicoCity_HTS,  k = 6) # longlat = NULL
Neigh_Centroids <- knn2nb(KNN)
WM_Centroids <- nb2listw(Neigh_Centroids, style = "W") 

class(Neigh_)

[1] "nb"

class(Neigh_Centroids)

[1] "nb"

class(Neigh_) == class(Neigh_Centroids)

[1] TRUE

class(WM_)

[1] "listw" "nb"   

class(WM_Centroids)

[1] "listw" "nb"   

class(WM_) == class(WM_Centroids)

[1] TRUE TRUE

Neigh_Centroids[1]

[[1]]
[1] 2 3 4 5 6 7

Neigh_Centroids[5]

[[1]]
[1]  4  6 32 33 34 36

KNN

$nn
       [,1] [,2] [,3] [,4] [,5] [,6]
  [1,]    2    3    4    5    6    7
  [2,]    1    3    7    8   16   21
  [3,]    1    2    4   21   28   32
  [4,]    1    3    5    6    7   32
  [5,]    4    6   32   33   34   36
  [6,]    1    4    5    7   11   12
  [7,]    1    2    6    8   10   11
  [8,]    1    2    7    9   10   16
  [9,]    7    8   10   14   15   55
 [10,]    7    8    9   11   13   14
 [11,]    6    7   10   12   13   14
 [12,]    6   11   13   36   37   38
 [13,]   10   11   12   14   40   46
 [14,]    9   10   11   13   15   47
 [15,]    8    9   10   14   47   55
 [16,]    2    8    9   17   21   57
 [17,]   16   18   20   21   91   96
 [18,]   17   19   20   91   95   96
 [19,]   17   18   20   96   99  100
 [20,]   17   18   19   21   25   26
 [21,]    1    2    3   16   17   20
 [22,]   23   24  101  115  126  135
 [23,]   22   24   25  100  101  126
 [24,]   22   23   25  100  125  126
 [25,]   23   24   26   27  100  125
 [26,]    3   20   25   27   28   32
 [27,]   26   28   29   30   32  125
 [28,]    3    4   26   27   29   32
 [29,]   27   30   31   32  125  138
 [30,]   27   29   31  137  138  141
 [31,]   29   30   32   33  138  139
 [32,]    3    5   27   28   29   31
 [33,]    5   31   32   34  138  139
 [34,]    5   33   35   36  162  163
 [35,]   34   36   37  162  163  164
 [36,]    5    6   12   34   35   37
 [37,]   12   35   36   38   85  163
 [38,]   12   36   37   39   40   85
 [39,]   38   40   41   42   84   85
 [40,]   12   13   38   39   41   46
 [41,]   38   39   40   42   43   46
 [42,]   39   40   41   43   77   78
 [43,]   41   42   44   45   69   73
 [44,]   43   45   49   50   68   69
 [45,]   41   43   44   46   49   50
 [46,]   13   14   40   41   45   49
 [47,]   14   15   46   48   49   51
 [48,]   45   47   49   50   51   64
 [49,]   44   45   46   47   48   50
 [50,]   44   45   48   49   51   64
 [51,]   47   48   50   52   53   64
 [52,]   48   51   53   54   56   63
 [53,]   51   52   54   55   56   57
 [54,]   52   53   55   56   57   58
 [55,]    9   15   53   56   57   58
 [56,]   53   54   55   57   58   87
 [57,]    9   53   54   55   56   58
 [58,]   54   55   56   57   87   88
 [59,]   52   54   56   60   61   87
 [60,]   52   59   61   62   63   65
 [61,]   52   59   60   62   63   65
 [62,]   52   60   61   63   65   66
 [63,]   51   52   61   62   64   65
 [64,]   48   49   50   51   65   66
 [65,]   51   52   62   63   64   66
 [66,]   44   50   63   64   65   68
 [67,]   62   64   65   66   68   70
 [68,]   44   50   64   66   69   70
 [69,]   43   44   45   50   68   73
 [70,]   66   67   68   69   71   72
 [71,]   68   69   70   72   73   75
 [72,]   69   70   73   74   75   76
 [73,]   43   69   72   74   76   77
 [74,]   72   73   76   77   78   79
 [75,]   72   74   76  184  185  186
 [76,]   74   77   78   79   80   81
 [77,]   41   42   73   76   78   79
 [78,]   39   42   76   77   79   84
 [79,]   76   77   78   80   81   84
 [80,]   76   78   79   81   83   84
 [81,]   74   76   79   80  179  184
 [82,]   83   85  164  166  167  168
 [83,]   80   82   84   85  167  168
 [84,]   39   78   79   82   83   85
 [85,]   39   82   83   84  164  167
 [86,]   88   89   90   93   94   97
 [87,]   56   57   58   88   89   90
 [88,]   58   87   89   90   91   93
 [89,]   87   88   90   91   92   93
 [90,]   87   88   89   91   92   93
 [91,]   17   88   89   92   95   96
 [92,]   89   90   91   93   95   96
 [93,]   88   89   90   91   92   94
 [94,]   90   91   92   93   95   97
 [95,]   18   91   92   93   94   96
 [96,]   17   18   19   91   92   95
 [97,]   94   95  103  104  105  106
 [98,]   99  100  101  102  103  104
 [99,]   18   19   96   98  100  101
[100,]   19   23   24   25   99  101
[101,]   22   23   98   99  100  115
[102,]   98   99  101  103  104  113
[103,]   97   98  102  104  105  113
[104,]  102  103  105  109  112  113
[105,]   97  103  104  106  109  110
[106,]   97  102  103  104  105  109
[107,]   97  103  105  106  108  109
[108,]  105  106  107  109  110  111
[109,]  103  104  105  106  110  112
[110,]  104  109  111  112  113  114
[111,]  109  110  112  114  116  117
[112,]  104  109  110  111  113  116
[113,]  102  104  110  112  116  124
[114,]  110  111  112  116  117  123
[115,]   22  101  113  116  124  127
[116,]  112  113  114  115  123  124
[117,]  114  116  122  123  130  131
[118,]  111  114  117  119  121  122
[119,]  110  111  114  117  118  121
[120,]  118  119  121  151  152  193
[121,]  117  118  122  131  151  152
[122,]  117  123  128  129  130  131
[123,]  116  122  124  127  128  130
[124,]   22  113  115  116  123  127
[125,]   24   27   29  126  135  136
[126,]   22   23   24  125  135  136
[127,]  123  124  128  129  132  133
[128,]  122  123  127  129  130  131
[129,]  127  128  130  131  132  150
[130,]  122  123  128  129  131  150
[131,]  122  128  129  130  132  150
[132,]  128  129  133  148  149  150
[133,]  127  129  132  134  135  144
[134,]  133  135  136  143  144  145
[135,]  125  126  133  134  136  144
[136,]   30  125  135  137  141  143
[137,]   30  136  138  141  142  143
[138,]   30   31  137  139  140  141
[139,]   30   31   33  138  140  141
[140,]  137  138  139  141  142  146
[141,]  136  137  138  140  142  143
[142,]  137  138  140  141  143  146
[143,]  136  137  141  142  145  146
[144,]  133  134  135  143  145  146
[145,]  134  143  144  146  147  148
[146,]  140  141  142  143  144  145
[147,]  144  145  148  149  154  155
[148,]  132  144  147  149  154  155
[149,]  132  147  148  150  154  155
[150,]  129  130  131  132  149  151
[151,]  130  131  149  150  153  154
[152,]  121  131  150  151  153  193
[153,]  148  149  150  151  154  155
[154,]  147  148  149  150  153  155
[155,]  145  147  148  149  154  156
[156,]  147  148  149  154  155  158
[157,]  152  153  154  156  158  193
[158,]  153  154  155  156  159  160
[159,]  140  141  142  146  172  173
[160,]  159  173  174  175  176  194
[161,]  158  160  174  175  182  194
[162,]   34   35  163  164  165  166
[163,]   35   37   85  162  164  167
[164,]   85  162  163  165  166  167
[165,]  162  164  166  167  169  172
[166,]   82  164  165  167  168  169
 [ reached getOption("max.print") -- omitted 28 rows ]

$np
[1] 194

$k
[1] 6

$dimension
[1] 2

$x
            X        Y
1   -99.13740 19.43337
2   -99.15832 19.44278
3   -99.13525 19.45262
4   -99.11631 19.43633
5   -99.10059 19.43750
6   -99.11947 19.41288
7   -99.14480 19.41449
8   -99.16832 19.41143
9   -99.17239 19.39336
10  -99.14686 19.39405
11  -99.13042 19.39160
12  -99.11070 19.39468
13  -99.12574 19.37219
14  -99.15352 19.36992
15  -99.17589 19.37559
16  -99.19272 19.43023
17  -99.19869 19.46129
18  -99.20752 19.48652
19  -99.19089 19.50085
20  -99.17018 19.48321
21  -99.16830 19.45841
22  -99.13043 19.56843
23  -99.14921 19.53886
24  -99.13142 19.53856
25  -99.14103 19.51295
26  -99.13901 19.49005
27  -99.10637 19.49802
28  -99.12413 19.47479
29  -99.09213 19.49323
30  -99.07574 19.49114
31  -99.07668 19.47299
32  -99.10698 19.46726
33  -99.07682 19.45479
34  -99.07180 19.43519
35  -99.06321 19.41941
36  -99.08725 19.41498
37  -99.07767 19.39733
38  -99.08538 19.38002
39  -99.06960 19.36130
40  -99.10608 19.36285
41  -99.09093 19.34480
42  -99.08073 19.33002
43  -99.10042 19.30897
44  -99.13094 19.29790
45  -99.12693 19.32541
46  -99.12462 19.34657
47  -99.17020 19.34874
48  -99.16899 19.32975
49  -99.15232 19.31927
50  -99.16085 19.30301
51  -99.19522 19.32222
52  -99.22690 19.33143
53  -99.21827 19.35622
54  -99.23954 19.36263
55  -99.20099 19.38189
56  -99.23713 19.37576
57  -99.21616 19.39384
58  -99.23262 19.40583
59  -99.29309 19.36052
60  -99.32819 19.29939
61  -99.28774 19.30552
62  -99.27522 19.26120
63  -99.22890 19.29703
64  -99.18227 19.29162
65  -99.22049 19.27396
66  -99.17999 19.25780
67  -99.21801 19.17601
68  -99.14105 19.25259
69  -99.10853 19.27055
70  -99.12394 19.20190
71  -99.05115 19.13951
72  -99.05479 19.23977
73  -99.07197 19.28122
74  -99.03723 19.28416
75  -98.98957 19.25017
76  -99.02825 19.30142
77  -99.05975 19.31900
78  -99.04703 19.33308
79  -99.02789 19.33118
80  -99.00581 19.33726
81  -98.99170 19.31691
82  -99.01841 19.37815
83  -99.02235 19.36845
84  -99.03926 19.35372
85  -99.05013 19.37985
86  -99.38270 19.46827
87  -99.26365 19.40619
88  -99.26397 19.42911
89  -99.26574 19.44086
90  -99.28033 19.44491
91  -99.24107 19.46047
92  -99.25378 19.46718
93  -99.28668 19.45974
94  -99.28332 19.49652
95  -99.24872 19.49980
96  -99.22339 19.49333
97  -99.29308 19.55098
98  -99.22059 19.55069
99  -99.19846 19.53216
100 -99.17253 19.53287
101 -99.17755 19.56628
102 -99.22670 19.57307
103 -99.25222 19.58164
104 -99.24630 19.59997
105 -99.29390 19.61652
106 -99.31230 19.59941
107 -99.40433 19.62720
108 -99.40821 19.72210
109 -99.26617 19.63697
110 -99.22985 19.65498
111 -99.21767 19.68483
112 -99.21634 19.63694
113 -99.20788 19.61320
114 -99.19143 19.68236
115 -99.15652 19.59268
116 -99.16838 19.64547
117 -99.14582 19.69994
118 -99.18829 19.76796
119 -99.23732 19.82269
120 -99.08422 19.93958
121 -99.08629 19.77902
122 -99.11797 19.68616
123 -99.13508 19.64643
124 -99.15306 19.61519
125 -99.10343 19.52537
126 -99.10877 19.54847
127 -99.10651 19.61263
128 -99.10614 19.64692
129 -99.08340 19.64297
130 -99.09299 19.66947
131 -99.08585 19.67980
132 -99.05911 19.62934
133 -99.07614 19.59222
134 -99.06258 19.56920
135 -99.08425 19.55313
136 -99.07231 19.52681
137 -99.06077 19.51401
138 -99.05715 19.48495
139 -99.04478 19.47077
140 -99.02975 19.49402
141 -99.04631 19.50873
142 -99.02566 19.50995
143 -99.04165 19.53551
144 -99.04408 19.56855
145 -99.02802 19.55814
146 -99.01995 19.53596
147 -99.00663 19.59737
148 -99.01949 19.60568
149 -99.01449 19.63166
150 -99.04451 19.66190
151 -99.02160 19.70280
152 -98.99054 19.77897
153 -98.96646 19.68240
154 -98.99249 19.64061
155 -98.98279 19.60416
156 -98.94136 19.58468
157 -98.78683 19.74491
158 -98.86314 19.58994
159 -98.96787 19.50332
160 -98.87044 19.48198
161 -98.75450 19.49393
162 -99.04181 19.42829
163 -99.04978 19.41043
164 -99.03594 19.40415
165 -99.01130 19.42050
166 -99.01065 19.40523
167 -99.02155 19.39383
168 -99.00148 19.38400
169 -98.98971 19.40083
170 -98.97977 19.37975
171 -98.96694 19.40083
172 -98.97537 19.42873
173 -98.94656 19.44317
174 -98.93969 19.42556
175 -98.90174 19.42825
176 -98.94295 19.40449
177 -98.96471 19.36905
178 -98.93927 19.35770
179 -98.93714 19.31382
180 -98.88981 19.34442
181 -98.89024 19.32954
182 -98.74829 19.34599
183 -98.89116 19.30383
184 -98.94940 19.30191
185 -98.94384 19.28344
186 -98.94887 19.25740
187 -98.89765 19.28179
188 -98.91053 19.25627
189 -98.83906 19.28392
190 -98.74895 19.23299
191 -98.88507 19.13936
192 -98.73891 19.05518
193 -98.96126 19.85208
194 -98.88028 19.39921

attr(,"class")
[1] "knn"
attr(,"call")
knearneigh(x = Centroids_MexicoCity_HTS, k = 6)

3 Spatial Autocorrelation

3.1 Preliminary exploration

We want to explore car-based trips. Should we analyze the total count? the share of total trips? or trips per inhabitants?

summary(MexicoCity_HTS$Trips_Automovil)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  110.0   300.8   438.0   540.9   642.8  2490.0 

tm_shape(MexicoCity_HTS) +
  tm_fill("Trips_Automovil") +
  tm_borders(alpha = 0.6)

MexicoCity_HTS %<>%
  mutate(CarTrips_Proportion = Trips_Automovil/TotalTrips,
         CarTrips_PerPopulation = Trips_Automovil/TotalPeople)

summary(MexicoCity_HTS$Trips_Automovil)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  110.0   300.8   438.0   540.9   642.8  2490.0 

summary(MexicoCity_HTS$CarTrips_Proportion)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.03783 0.08574 0.11236 0.12337 0.15241 0.41726 

summary(MexicoCity_HTS$CarTrips_PerPopulation)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.1008  0.2905  0.4219  0.5658  0.6500  3.8485       1 

Why do we have an NA?

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...

# View(MexicoCity_HTS %>% select(Distrito, TotalPeople))

MexicoCity_HTS %<>% replace_na(list(TotalPeople = 120)) # Any take on this?

MexicoCity_HTS %<>%
  mutate(CarTrips_Proportion = Trips_Automovil/TotalTrips,
         CarTrips_PerPopulation = Trips_Automovil/TotalPeople)

summary(MexicoCity_HTS$Trips_Automovil)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  110.0   300.8   438.0   540.9   642.8  2490.0 

summary(MexicoCity_HTS$CarTrips_Proportion)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.03783 0.08574 0.11236 0.12337 0.15241 0.41726 

summary(MexicoCity_HTS$CarTrips_PerPopulation)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1008  0.2916  0.4249  0.5680  0.6554  3.8485 

tm_shape(MexicoCity_HTS) +
  tm_fill("CarTrips_Proportion", style = "quantile") +
  tm_borders(alpha = 0.6)

tm_shape(MexicoCity_HTS) +
  tm_fill("CarTrips_PerPopulation", style = "quantile") +
  tm_borders(alpha = 0.6)

3.1.1 Ok, let’s fix that map

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 

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

st_crs(MexicoCity_HTS) == st_crs(MexicoCity_Fringe)

[1] TRUE

tm_shape(MexicoCity_HTS, bbox = BB_Districts) +
  tm_fill("CarTrips_PerPopulation", style = "quantile", title = "Quantiles") +
  tm_borders(alpha = 0.6) +
  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 = "Car trips / inhabitant", legend.position = c(0.01,0.15), scale = 1)

3.2 Moran’s I

We can use the two following functions from the spdep library:

• moran.test(): Analytical method. From the spdep
• moran.mc(): Montecarlo based method

moran.test(MexicoCity_HTS$Trips_Automovil, WM_)

    Moran I test under randomisation

data:  MexicoCity_HTS$Trips_Automovil  
weights: WM_    

Moran I statistic standard deviate = 10.616, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.459283602      -0.005181347       0.001914295 

moran.test(MexicoCity_HTS$CarTrips_Proportion, WM_)

    Moran I test under randomisation

data:  MexicoCity_HTS$CarTrips_Proportion  
weights: WM_    

Moran I statistic standard deviate = 11.238, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.487662076      -0.005181347       0.001923353 

moran.test(MexicoCity_HTS$CarTrips_Proportion, WM_)

    Moran I test under randomisation

data:  MexicoCity_HTS$CarTrips_Proportion  
weights: WM_    

Moran I statistic standard deviate = 11.238, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.487662076      -0.005181347       0.001923353 

moran.mc(MexicoCity_HTS$Trips_Automovil, WM_, nsim = 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$Trips_Automovil 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.45928, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

moran.mc(MexicoCity_HTS$CarTrips_Proportion, WM_, nsim = 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$CarTrips_Proportion 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.48766, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

moran.mc(MexicoCity_HTS$CarTrips_PerPopulation, WM_, nsim = 1000)

    Monte-Carlo simulation of Moran I

data:  MexicoCity_HTS$CarTrips_PerPopulation 
weights: WM_  
number of simulations + 1: 1001 

statistic = 0.4752, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

3.3 Moran’s I: the regression way

Variable.lag <- lag.listw(WM_, MexicoCity_HTS$CarTrips_PerPopulation) # What is this?

DataTemp <- data.frame(
  Lagged = Variable.lag, 
  Value = MexicoCity_HTS$CarTrips_PerPopulation)

MeanValues = mean(DataTemp$Value)
MeanLagged = mean(DataTemp$Lagged)

Model_Y_Ylagged <- lm(Lagged ~ Value, data = DataTemp)

DataTemp$Y_Est <- predict.lm(Model_Y_Ylagged, DataTemp)

summary(Model_Y_Ylagged)

Call:
lm(formula = Lagged ~ Value, data = DataTemp)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.01239 -0.15256 -0.05003  0.07945  0.99587 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.31404    0.02855   11.00   <2e-16 ***
Value        0.47520    0.03846   12.36   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2559 on 192 degrees of freedom
Multiple R-squared:  0.4429,    Adjusted R-squared:   0.44 
F-statistic: 152.7 on 1 and 192 DF,  p-value: < 2.2e-16

coef(Model_Y_Ylagged)[2] # How do you interpret this result?

    Value 
0.4751972 

ggplot() +
  geom_point(data = DataTemp, aes(x = Value, y = Lagged)) +
  geom_vline(xintercept = MeanValues, linetype = "dashed") +
  geom_hline(yintercept = MeanLagged, linetype = "dashed") +
  geom_smooth(data = DataTemp, aes(x = Value, y = Y_Est, method = "lm"), 
              color = "#f03b20") +
  # scale_x_continuous(limits = c(0, 4)) +
  # scale_y_continuous(limits = c(0, 2)) +
  xlab("Variable") + ylab("Lagged Variable") +
  theme_light()

Warning: Ignoring unknown aesthetics: method

3.4 Geary index

function geary.test() from the spdep library

geary.test(MexicoCity_HTS$CarTrips_PerPopulation, WM_)

    Geary C test under randomisation

data:  MexicoCity_HTS$CarTrips_PerPopulation 
weights: WM_ 

Geary C statistic standard deviate = 6.6255, p-value = 1.73e-11
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic       Expectation          Variance 
       0.59700354        1.00000000        0.00369967 

4 LISA

4.1 Moran’s I

Most of the work is handled by the localmoran() function from the spdep library/

MexicoCity_HTS$Variable <- MexicoCity_HTS$CarTrips_PerPopulation
ZonesAnalysis <- MexicoCity_HTS %>% select(Variable)

Variable.lag <- lag.listw(WM_, ZonesAnalysis$Variable) 

DataTemp <- data.frame(
  Lagged = Variable.lag, 
  Value = ZonesAnalysis$Variable)

MeanValues = mean(DataTemp$Value)
MeanLagged = mean(DataTemp$Lagged)

Model_Y_Ylagged <- lm(Lagged ~ Value, data = DataTemp)

DataTemp$Y_Est <- predict.lm(Model_Y_Ylagged, DataTemp)

LocalMoran_Temp <- localmoran(ZonesAnalysis$Variable, WM_)

LocalMoran_Temp <- as.data.frame(LocalMoran_Temp)
names(LocalMoran_Temp)[5] <- "PValue"

ZonesAnalysis %<>% bind_cols(LocalMoran_Temp)

ZonesAnalysis$Lagged <- lag.listw(WM_, ZonesAnalysis$Variable)
ZonesAnalysis$Value <- ZonesAnalysis$Variable

MeanValues = mean(ZonesAnalysis$Value)
MeanLagged = mean(ZonesAnalysis$Lagged)
MeanMoran = mean(ZonesAnalysis$Ii)

ZonesAnalysis %<>%
  mutate(
    Values_Centered = Value - MeanValues,
    Lagged_Centered = Lagged - MeanLagged,
    Moran_Centered = Ii - MeanMoran)

ZonesAnalysis %<>%
  mutate(Significance = if_else(PValue <= 0.05, 1, 0)) %>% 
  mutate(Cuadrants = 5) %>%
  mutate(Cuadrants = if_else(Values_Centered > 0 & Lagged_Centered > 0, 1, Cuadrants)) %>% 
  mutate(Cuadrants = if_else(Values_Centered > 0 & Lagged_Centered < 0, 2, Cuadrants)) %>%
  mutate(Cuadrants = if_else(Values_Centered < 0 & Lagged_Centered < 0, 3, Cuadrants)) %>%
  mutate(Cuadrants = if_else(Values_Centered < 0 & Lagged_Centered > 0, 4, Cuadrants)) %>% 
  mutate(Cuadrants = if_else(Significance == 0, 5, Cuadrants))

# 1: High-High - Spatial clusters
# 2: High-Low - Outliers
# 3: Low-low - Spatial clusters
# 4: Low-High - Outliers
# 5: Not significance
# c("High-High", "High-Low", "Low-low", "Low-high", "Not significant")
ZonesAnalysis$Cuadrants <- factor(as.character(ZonesAnalysis$Cuadrants), 
                                  levels = c("1", "2", "3", "4", "5"),
                                  ordered = TRUE)

ggplot() +
  geom_point(data = ZonesAnalysis, aes(x = Value, y = Lagged, color = Cuadrants)) +
  scale_color_manual(values = c("#1b9e77", "#7570b3", "#e7298a", "#d0d1e6"),
                    labels = c("High-High", "Low-low", "Low-high", "Not significant")) + 
  geom_vline(xintercept = MeanValues, linetype = "dashed") +
  geom_hline(yintercept = MeanLagged, linetype = "dashed") +
  geom_smooth(data = DataTemp, aes(x = Value, y = Y_Est, method = "lm"), 
              color = "#f03b20") +
  xlab("Variable") + ylab("Lagged Variable") +
  theme_light()

Warning: Ignoring unknown aesthetics: method

Breaks <- c(0, 1.5, 2.5, 3.5, 4.5, 5.5)
Labels <- c("High-high (clusters)", "High-low", 
            "Low-low (clusters)", "Low-high", "Not significant")
MyPalette <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#edf8fb")

Map_LISA <- tm_shape(ZonesAnalysis, bbox = BB_Districts) +
  tm_fill("Cuadrants",
           palette = MyPalette, breaks = Breaks, labels = Labels,
          title = "LISA cluster map") +
  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 = "Car trips / inhabitant", legend.position = c(0.01,0.15), scale = 1)

Map_LISA

```r
tmap_save(tm = Map_LISA, filename = \Maps_LISA.jpeg\, 
          width = 2500, height = 2500)

```