Take-home_Exercise 02
Erika Aldisa
9/28/2020 (updated: 2020-10-14)
IS415 Geospatial Analytics and Application
Instructor: Dr. Kam Tin Seong.
Assoc. Professor of Information Systems (Practice)
1. Overview
Task - Analyze the spatio-temporal patterns of COVID-19 cases at Central Mexico (Mexico City (9), Mexico State (15), Morelos State (17)) by using localized spatial statistics methods.
- Geospatial data wrangling
Extracting the study area of Central Mexico & calculate the COVID-19 rate (i.e. cases per 10000 population) from e-week 13 until e-week 32
Thematic mapping technique: show spatio-temporal distribution of COVID-19 rates & describe the spatio-temporal patterns
Perform Local Moran’s I analysis
Display results by appropriate thematic mapping techniques & describe the spatio-temporal patterns revealed
- Perform Local Getis-Ord Gi* analysis
Display results by appropriate thematic mapping techniques & describe the spatio-temporal patterns revealed
1.1. Data sets used
With a special focus on Mexico City (9), Mexico State (15), Morelos State (17) and epidemiological week (e-week) 13 to e-week 32.
1.2. Install & Load R packages
R packages used:
tidyverse - for reading, handling and visualizing attribute data (readr, ggplot2 & dplyr)
tmap - for choropleth mapping
rgdal - for handling spatial data
spdep - for spatial statistical analysis
spdplyr - ‘dplyr’ methods for Spatial classes
raster - for handling spatial data
packages = c('tidyverse', 'tmap','rgdal', 'spdep', 'spdplyr','raster')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p, character.only = T)
}2. Geospatial data wrangling
2.1 Importing geospatial data
We will use readOGR() function to import the data as Spatial object.
mexico <- readOGR(dsn = "data/geospatial", layer = "municipalities_COVID")## OGR data source with driver: ESRI Shapefile
## Source: "D:\SMU Year 3 Sem 1\IS415 Geospatial Analytics and Application\Take-home_Ex02\data\geospatial", layer: "municipalities_COVID"
## with 2465 features
## It has 198 fieldshead(mexico@data, n = 10)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 new1 new2 new3
## 0 01001 01 001 Aguascalientes 797010 961977 0 0 0
## 1 01002 01 002 Asientos 45492 50864 0 0 0
## 2 01003 01 003 Calvillo 54136 60760 0 0 0
## 3 01004 01 004 CosÃo 15042 16918 0 0 0
## 4 01005 01 005 Jesús MarÃa 99590 130184 0 0 0
## 5 01006 01 006 Pabellón de Arteaga 41862 50032 0 0 0
## 6 01007 01 007 Rincón de Romos 49156 57981 0 0 0
## 7 01008 01 008 San José de Gracia 8443 9661 0 0 0
## 8 01009 01 009 Tepezalá 19668 22743 0 0 0
## 9 01010 01 010 El Llano 18828 21947 0 0 0
## new4 new5 new6 new7 new8 new9 new10 new11 new12 new13 new14 new15 new16 new17
## 0 0 0 0 0 0 0 0 1 7 30 8 10 51 60
## 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0 0 8
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 1 0 0 6 3
## 5 0 0 0 0 0 0 0 0 0 4 1 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 1 1 1 2
## 7 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## 9 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## new18 new19 new20 new21 new22 new23 new24 new25 new26 new27 new28 new29 new30
## 0 87 101 101 99 185 257 312 280 258 234 295 307 267
## 1 1 0 0 5 2 1 8 6 4 9 11 6 7
## 2 1 1 1 0 0 3 3 8 11 8 20 35 13
## 3 0 0 6 7 3 7 6 12 1 0 2 10 9
## 4 3 4 5 5 6 5 11 10 15 14 15 13 10
## 5 6 3 2 14 6 17 8 19 13 7 8 21 14
## 6 4 6 8 14 12 18 17 25 20 29 42 43 26
## 7 4 0 0 0 0 2 3 7 6 1 5 12 4
## 8 1 1 1 13 6 15 9 10 7 6 11 7 17
## 9 0 1 0 0 2 1 0 0 1 1 2 2 0
## new31 new32 cumul1 cumul2 cumul3 cumul4 cumul5 cumul6 cumul7 cumul8 cumul9
## 0 219 23 0 0 0 0 0 0 0 0 0
## 1 11 0 0 0 0 0 0 0 0 0 0
## 2 13 3 0 0 0 0 0 0 0 0 0
## 3 4 0 0 0 0 0 0 0 0 0 0
## 4 15 2 0 0 0 0 0 0 0 0 0
## 5 19 0 0 0 0 0 0 0 0 0 0
## 6 23 7 0 0 0 0 0 0 0 0 0
## 7 4 0 0 0 0 0 0 0 0 0 0
## 8 11 0 0 0 0 0 0 0 0 0 0
## 9 2 1 0 0 0 0 0 0 0 0 0
## cumul10 cumul11 cumul12 cumul13 cumul14 cumul15 cumul16 cumul17 cumul18
## 0 0 1 8 38 46 56 107 167 254
## 1 0 0 1 1 1 1 1 1 2
## 2 0 0 0 0 0 0 0 8 9
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 1 1 1 7 10 13
## 5 0 0 0 4 5 5 5 5 11
## 6 0 0 0 0 1 2 3 5 9
## 7 0 0 0 0 0 0 0 1 5
## 8 0 0 0 1 1 1 1 1 2
## 9 0 0 0 0 0 1 1 1 1
## cumul19 cumul20 cumul21 cumul22 cumul23 cumul24 cumul25 cumul26 cumul27
## 0 355 456 555 740 997 1309 1589 1847 2081
## 1 2 2 7 9 10 18 24 28 37
## 2 10 11 11 11 14 17 25 36 44
## 3 0 6 13 16 23 29 41 42 42
## 4 17 22 27 33 38 49 59 74 88
## 5 14 16 30 36 53 61 80 93 100
## 6 15 23 37 49 67 84 109 129 158
## 7 5 5 5 5 7 10 17 23 24
## 8 3 4 17 23 38 47 57 64 70
## 9 2 2 2 4 5 5 5 6 7
## cumul28 cumul29 cumul30 cumul31 cumul32 activ1 activ2 activ3 activ4 activ5
## 0 2376 2683 2950 3169 3192 0 0 0 0 0
## 1 48 54 61 72 72 0 0 0 0 0
## 2 64 99 112 125 128 0 0 0 0 0
## 3 44 54 63 67 67 0 0 0 0 0
## 4 103 116 126 141 143 0 0 0 0 0
## 5 108 129 143 162 162 0 0 0 0 0
## 6 200 243 269 292 299 0 0 0 0 0
## 7 29 41 45 49 49 0 0 0 0 0
## 8 81 88 105 116 116 0 0 0 0 0
## 9 9 11 11 13 14 0 0 0 0 0
## activ6 activ7 activ8 activ9 activ10 activ11 activ12 activ13 activ14 activ15
## 0 0 0 0 0 1 4 16 45 53 75
## 1 0 0 0 0 0 0 1 1 1 1
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 1 1 3
## 5 0 0 0 0 0 0 0 5 5 5
## 6 0 0 0 0 0 0 0 0 1 2
## 7 0 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 1 1 1
## 9 0 0 0 0 0 0 0 0 1 1
## activ16 activ17 activ18 activ19 activ20 activ21 activ22 activ23 activ24
## 0 139 203 309 393 500 622 848 1133 1424
## 1 1 1 2 2 4 9 10 14 19
## 2 0 9 9 11 11 11 12 15 18
## 3 0 0 0 1 7 14 18 27 34
## 4 7 12 14 19 23 28 34 44 53
## 5 5 5 13 14 22 34 41 60 70
## 6 3 5 9 15 26 41 52 75 86
## 7 1 1 5 5 5 5 5 7 12
## 8 1 2 2 3 5 18 26 39 48
## 9 1 1 2 2 2 2 5 5 5
## activ25 activ26 activ27 activ28 activ29 activ30 activ31 activ32 death1 death2
## 0 1713 1972 2248 2544 2845 3079 3189 3192 0 0
## 1 27 31 41 51 56 70 72 72 0 0
## 2 28 39 53 90 108 124 127 128 0 0
## 3 42 42 43 50 56 66 67 67 0 0
## 4 69 83 96 109 122 135 143 143 0 0
## 5 88 95 104 119 135 148 162 162 0 0
## 6 116 141 171 220 255 282 297 299 0 0
## 7 22 23 26 30 43 46 49 49 0 0
## 8 60 68 79 83 97 110 116 116 0 0
## 9 6 6 8 10 11 13 14 14 0 0
## death3 death4 death5 death6 death7 death8 death9 death10 death11 death12
## 0 0 0 0 0 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0 0 0
## death13 death14 death15 death16 death17 death18 death19 death20 death21
## 0 0 0 1 2 2 3 7 6 4
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 2
## 5 0 0 0 0 0 0 0 0 1
## 6 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0 0
## death22 death23 death24 death25 death26 death27 death28 death29 death30
## 0 17 18 27 25 24 19 19 19 20
## 1 0 0 0 0 2 0 2 2 0
## 2 0 0 0 0 1 0 0 0 1
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 1 3
## 5 1 0 1 1 1 0 1 0 1
## 6 0 1 0 1 0 3 0 0 1
## 7 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 1 0 0 0
## 9 0 0 0 0 0 0 1 1 0
## death31 death32 actvrt1 actvrt2 actvrt3 actvrt4 actvrt5 actvrt6 actvrt7
## 0 13 7 0 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0 0
## actvrt8 actvrt9 actvr10 actvr11 actvr12 actvr13 actvr14 actvr15
## 0 0 0 0.1039526 0.4158104 1.663241 4.5739139 5.0936769 6.133203
## 1 0 0 0.0000000 0.0000000 1.966027 1.9660271 1.9660271 0.000000
## 2 0 0 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.000000
## 3 0 0 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.000000
## 4 0 0 0.0000000 0.0000000 0.000000 0.7681436 0.7681436 2.304431
## 5 0 0 0.0000000 0.0000000 0.000000 9.9936041 9.9936041 9.993604
## 6 0 0 0.0000000 0.0000000 0.000000 0.0000000 1.7247029 3.449406
## 7 0 0 0.0000000 0.0000000 0.000000 0.0000000 0.0000000 0.000000
## 8 0 0 0.0000000 0.0000000 0.000000 4.3969573 4.3969573 4.396957
## 9 0 0 0.0000000 0.0000000 0.000000 0.0000000 4.5564314 4.556431
## actvr16 actvr17 actvr18 actvr19 actvr20 actvr21 actvr22
## 0 9.771543 15.592888 24.324906 26.403958 30.873919 32.537160 47.29843
## 1 0.000000 0.000000 1.966027 1.966027 5.898081 13.762189 15.72822
## 2 0.000000 14.812377 14.812377 18.104016 3.291639 3.291639 1.64582
## 3 0.000000 0.000000 0.000000 5.910864 41.376049 82.752098 100.48469
## 4 4.608861 8.449579 8.449579 9.217723 8.449579 10.754010 11.52215
## 5 0.000000 0.000000 15.989767 17.988487 33.978254 41.973137 53.96546
## 6 5.174109 6.898812 12.072920 20.696435 36.218761 55.190493 63.81401
## 7 10.350895 10.350895 51.754477 41.403581 41.403581 0.000000 0.00000
## 8 0.000000 4.396957 4.396957 8.793915 13.190872 70.351317 101.13002
## 9 4.556431 0.000000 4.556431 4.556431 4.556431 0.000000 13.66929
## actvr23 actvr24 actvr25 actvr26 actvr27 actvr28 actvr29
## 0 65.801989 83.36998 89.918990 87.216222 85.65693 86.38460 90.75061
## 1 19.660271 19.66027 33.422460 33.422460 43.25260 47.18465 49.15068
## 2 6.583278 11.52074 26.333114 39.499671 57.60369 102.04082 113.56155
## 3 118.217283 118.21728 141.860740 88.662963 53.19778 47.28691 82.75210
## 4 16.131015 19.20359 26.885024 29.957598 33.03017 30.72574 29.95760
## 5 75.951391 71.95395 93.939878 69.955229 67.95651 61.96035 79.94883
## 6 84.510443 77.61163 110.380987 113.830393 146.59975 179.36910 196.61613
## 7 20.701791 72.45627 175.965221 165.614326 144.91253 82.80716 207.01791
## 8 149.496548 131.90872 149.496548 127.511762 136.30568 101.13002 127.51176
## 9 13.669294 13.66929 4.556431 4.556431 13.66929 18.22573 22.78216
## actvr30 actvr31 actvr32 dethrt1 dethrt2 dethrt3 dethrt4 dethrt5 dethrt6
## 0 86.38460 67.04942 36.07155 0 0 0 0 0 0
## 1 57.01478 41.28657 31.45643 0 0 0 0 0 0
## 2 116.85319 60.89533 32.91639 0 0 0 0 0 0
## 3 135.94988 100.48469 65.01951 0 0 0 0 0 0
## 4 29.95760 26.11688 16.13101 0 0 0 0 0 0
## 5 87.94372 85.94500 53.96546 0 0 0 0 0 0
## 6 191.44202 132.80212 75.88693 0 0 0 0 0 0
## 7 207.01791 196.66701 62.10537 0 0 0 0 0 0
## 8 136.30568 145.09959 83.54219 0 0 0 0 0 0
## 9 22.78216 18.22573 13.66929 0 0 0 0 0 0
## dethrt7 dethrt8 dethrt9 dthrt10 dthrt11 dthrt12 dthrt13 dthrt14 dthrt15
## 0 0 0 0 0 0 0 0 0 0.1039526
## 1 0 0 0 0 0 0 0 0 0.0000000
## 2 0 0 0 0 0 0 0 0 0.0000000
## 3 0 0 0 0 0 0 0 0 0.0000000
## 4 0 0 0 0 0 0 0 0 0.0000000
## 5 0 0 0 0 0 0 0 0 0.0000000
## 6 0 0 0 0 0 0 0 0 0.0000000
## 7 0 0 0 0 0 0 0 0 0.0000000
## 8 0 0 0 0 0 0 0 0 0.0000000
## 9 0 0 0 0 0 0 0 0 0.0000000
## dthrt16 dthrt17 dthrt18 dthrt19 dthrt20 dthrt21 dthrt22 dthrt23
## 0 0.2079052 0.2079052 0.3118578 0.7276681 0.6237155 0.4158104 1.767194 1.871147
## 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000
## 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000
## 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000
## 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.5362871 0.000000 0.000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.9987208 1.998721 0.000000
## 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 1.724703
## 7 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000
## 8 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000
## 9 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000
## dthrt24 dthrt25 dthrt26 dthrt27 dthrt28 dthrt29 dthrt30 dthrt31
## 0 2.806720 2.598815 2.494862 1.975099 1.975099 1.9750992 2.079052 1.351384
## 1 0.000000 0.000000 3.932054 0.000000 3.932054 3.9320541 0.000000 0.000000
## 2 0.000000 0.000000 1.645820 0.000000 0.000000 0.0000000 1.645820 0.000000
## 3 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000 0.000000
## 4 0.000000 0.000000 0.000000 0.000000 0.000000 0.7681436 2.304431 0.000000
## 5 1.998721 1.998721 1.998721 0.000000 1.998721 0.0000000 1.998721 0.000000
## 6 0.000000 1.724703 0.000000 5.174109 0.000000 0.0000000 1.724703 0.000000
## 7 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000 0.000000
## 8 0.000000 0.000000 0.000000 4.396957 0.000000 0.0000000 0.000000 0.000000
## 9 0.000000 0.000000 0.000000 0.000000 4.556431 4.5564314 0.000000 0.000000
## dthrt32
## 0 0.7276681
## 1 0.0000000
## 2 0.0000000
## 3 0.0000000
## 4 0.0000000
## 5 0.0000000
## 6 0.0000000
## 7 0.0000000
## 8 0.0000000
## 9 0.0000000We can see that CVGEO variable is the unique ID of NOMGEO which is the concatenation of CV_ENT and CV_MUN.
Next, we should check if there are any missing/NA values in the mexico SpatialPolygonsDataFrame.
any(is.na(as.data.frame(mexico)))## [1] FALSEWe can see that there are no NA values, although there are 0 values.
Let’s check the Coordinate Reference System (CRS) of the mexico SpatialPolygonsDataFrame.
crs(mexico)## CRS arguments:
## +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defsLet’s transform mexico’s CRS to WGS 84 for standard use and so that our coordinates can be translated better.
mexico_wgs84 <- spTransform(mexico, CRS("+init=epsg:4326"))crs(mexico_wgs84)## CRS arguments: +proj=longlat +datum=WGS84 +no_defs2.2 Plotting the Base map
tm_shape(mexico_wgs84)+
tm_polygons("CVE_ENT")+
tm_borders(lwd=0.5)+
tm_layout(legend.height = 0.45,
legend.width = 0.35,
frame = FALSE) +
tm_compass(type="8star", size = 2, position = c("left", "bottom")) +
tm_scale_bar(width = 0.15, position = c("left", "bottom")) +
tm_credits("Secretaria de Salud (Secretary of Health, Mexico)",
position = c("left", "bottom"))
2.3 Extracting & Plotting the Study Area
We will be extracting Mexico City(9), Mexico State (15), Morelos State (17) by filtering the CVE_ENT value.
MEXICO = mexico_wgs84[(mexico_wgs84$CVE_ENT == "09" | mexico_wgs84$CVE_ENT == "15" | mexico_wgs84$CVE_ENT == "17"),]
MEX = mexico_wgs84[mexico_wgs84$CVE_ENT == "09",]
MES = mexico_wgs84[mexico_wgs84$CVE_ENT == "15",]
MOS = mexico_wgs84[mexico_wgs84$CVE_ENT == "17",]Let’s plot each municipality first and then plot the study area.
par(mfrow= c(1,3))
plot(MEX, main = "Mexico City")
plot(MES, main = "Mexico State")
plot(MOS, main = "Morelos State")
We can see that we cannot identify how big the area is if we plot it one by one.
qtm(MEXICO, fill = "CVE_ENT", frame = FALSE)
We can see that the 3 state actually make a full study area without any gaps. We can also figure out that Mexico City (CVE_ENT 9) is a small area compared to Mexico State (CVE_ENT 15).
2.4 Compute COVID-19 rate
We will calculate the COVID-19 rate (cases per population) from e-week 13 until e-week 32. So, we will be deselecting the variables for e-week 1 to e-week 12.
MXCO <- MEXICO[,-(7:18)]
MXCO <- MXCO[,-(27:38)]
MXCO <- MXCO[,-(47:58)]
MXCO <- MXCO[,-(67:78)]
MXCO <- MXCO[,-(87:98)]
MXCO <- MXCO[,-(107:118)]
head(MXCO@data)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 new13 new14
## 270 09002 09 002 Azcapotzalco 414711 408441 18 27
## 271 09003 09 003 Coyoacán 620416 621952 15 19
## 272 09004 09 004 Cuajimalpa de Morelos 186391 199809 18 9
## 273 09005 09 005 Gustavo A. Madero 1185772 1176967 21 68
## 274 09006 09 006 Iztacalco 384326 393821 9 16
## 275 09007 09 007 Iztapalapa 1815786 1815551 32 59
## new15 new16 new17 new18 new19 new20 new21 new22 new23 new24 new25 new26
## 270 41 69 98 137 224 316 336 323 349 310 345 328
## 271 63 106 162 179 277 319 342 374 389 299 382 272
## 272 18 29 35 45 58 63 91 159 141 137 140 126
## 273 97 204 337 373 546 677 778 735 831 692 611 525
## 274 38 123 155 209 254 319 305 276 283 220 245 206
## 275 125 301 535 693 863 966 1076 873 884 818 747 665
## new27 new28 new29 new30 new31 new32 cumul13 cumul14 cumul15 cumul16 cumul17
## 270 316 301 300 332 271 24 19 46 87 156 254
## 271 382 313 490 427 396 30 25 44 107 213 375
## 272 168 159 171 220 114 7 66 75 93 122 157
## 273 573 512 625 663 505 59 28 96 193 397 734
## 274 165 172 186 233 196 62 11 27 65 188 343
## 275 762 670 744 728 620 135 41 100 225 526 1061
## cumul18 cumul19 cumul20 cumul21 cumul22 cumul23 cumul24 cumul25 cumul26
## 270 391 615 931 1267 1590 1939 2249 2594 2922
## 271 554 831 1150 1492 1866 2255 2554 2936 3208
## 272 202 260 323 414 573 714 851 991 1117
## 273 1107 1653 2330 3108 3843 4674 5366 5977 6502
## 274 552 806 1125 1430 1706 1989 2209 2454 2660
## 275 1754 2617 3583 4659 5532 6416 7234 7981 8646
## cumul27 cumul28 cumul29 cumul30 cumul31 cumul32 activ13 activ14 activ15
## 270 3238 3539 3839 4171 4442 4466 34 70 130
## 271 3590 3903 4393 4820 5216 5246 41 82 166
## 272 1285 1444 1615 1835 1949 1956 70 86 112
## 273 7075 7587 8212 8875 9380 9439 68 173 321
## 274 2825 2997 3183 3416 3612 3674 20 58 129
## 275 9408 10078 10822 11550 12170 12305 75 174 427
## activ16 activ17 activ18 activ19 activ20 activ21 activ22 activ23 activ24
## 270 199 332 516 776 1126 1456 1795 2143 2451
## 271 321 513 731 1054 1398 1751 2117 2496 2841
## 272 140 182 237 298 385 498 660 801 917
## 273 571 940 1438 2015 2755 3557 4275 5115 5756
## 274 286 482 722 1008 1328 1618 1893 2147 2394
## 275 854 1492 2315 3207 4207 5208 6095 6979 7723
## activ25 activ26 activ27 activ28 activ29 activ30 activ31 activ32 death13
## 270 2810 3117 3444 3746 4034 4340 4460 4466 0
## 271 3169 3496 3862 4227 4681 5075 5245 5246 0
## 272 1080 1181 1384 1540 1677 1885 1954 1956 0
## 273 6329 6842 7459 8022 8636 9182 9427 9439 2
## 274 2591 2766 2953 3137 3369 3553 3666 3674 1
## 275 8434 9100 9874 10528 11287 11950 12283 12305 2
## death14 death15 death16 death17 death18 death19 death20 death21 death22
## 270 2 4 6 15 14 38 33 44 47
## 271 1 7 8 19 11 26 30 36 34
## 272 2 1 3 5 1 4 11 8 9
## 273 5 26 19 37 74 82 98 117 114
## 274 1 3 7 26 20 30 46 49 37
## 275 6 15 33 76 115 127 174 176 127
## death23 death24 death25 death26 death27 death28 death29 death30 death31
## 270 44 43 47 38 27 19 31 20 19
## 271 41 34 37 33 31 25 23 20 16
## 272 11 18 6 5 5 9 10 7 8
## 273 121 148 112 77 72 56 57 52 33
## 274 29 35 26 22 19 25 23 17 12
## 275 120 110 103 63 67 65 52 29 29
## death32 actvr13 actvr14 actvr15 actvr16 actvr17 actvr18 actvr19
## 270 8 8.324336 16.893505 29.86968 40.39751 64.14635 94.50569 141.26887
## 271 6 6.592149 12.058808 23.63526 45.01955 69.29795 90.84302 117.85475
## 272 1 32.531067 26.525332 26.02485 35.03346 48.04588 62.55974 79.07552
## 273 18 5.607634 14.358941 25.57421 42.73697 65.16750 94.90495 122.68823
## 274 3 5.078449 14.219658 31.23246 67.54338 107.66313 150.57602 183.33202
## 275 17 4.020818 9.308469 22.08696 42.90708 72.59504 103.99047 129.60253
## actvr20 actvr21 actvr22 actvr23 actvr24 actvr25 actvr26 actvr27
## 270 194.3977 230.1434 249.4852 248.9956 243.6092 248.5059 238.4677 243.1196
## 271 142.2939 163.9998 170.9135 176.5410 175.2547 169.1449 160.7841 164.1606
## 272 101.5970 130.6247 181.1730 208.1988 209.7003 210.2007 190.1816 233.7232
## 273 154.2099 180.0390 192.0190 200.5154 186.8362 174.5164 146.7331 144.6939
## 274 214.8184 227.5145 224.7214 207.9625 197.0438 177.2379 157.1780 141.9427
## 275 149.5414 159.3456 159.0702 152.6809 138.5254 128.8314 116.8240 118.4764
## actvr28 actvr29 actvr30 actvr31 actvr32 dthrt13 dthrt14
## 270 229.1641 224.5122 219.3707 174.81105 105.76803 0.0000000 0.4896668
## 271 170.1096 190.5292 195.0311 163.67823 90.84302 0.0000000 0.1607841
## 272 230.2199 248.2371 250.7395 207.19787 139.63335 0.0000000 1.0009559
## 273 143.8443 152.4257 146.3932 119.37463 68.22621 0.1699283 0.4248207
## 274 138.6417 153.1152 152.3535 134.32499 77.44635 0.2539225 0.2539225
## 275 115.3369 120.4593 114.3455 96.66487 56.07113 0.1101594 0.3304782
## dthrt15 dthrt16 dthrt17 dthrt18 dthrt19 dthrt20 dthrt21 dthrt22
## 270 0.9793336 1.469000 3.672501 3.427668 9.303669 8.079502 10.772670 11.507170
## 271 1.1254888 1.286273 3.054898 1.768625 4.180387 4.823523 5.788228 5.466660
## 272 0.5004780 1.501434 2.502390 0.500478 2.001912 5.505258 4.003824 4.504302
## 273 2.2090679 1.614319 3.143674 6.287347 6.967060 8.326487 9.940805 9.685913
## 274 0.7617674 1.777457 6.601984 5.078449 7.617674 11.680433 12.442201 9.395131
## 275 0.8261955 1.817630 4.186057 6.334165 6.995122 9.583867 9.694027 6.995122
## dthrt23 dthrt24 dthrt25 dthrt26 dthrt27 dthrt28 dthrt29 dthrt30
## 270 10.772670 10.527836 11.507170 9.303669 6.610502 4.651835 7.589835 4.896668
## 271 6.592149 5.466660 5.949012 5.305876 4.984307 4.019603 3.698035 3.215682
## 272 5.505258 9.008603 3.002868 2.502390 2.502390 4.504302 5.004780 3.503346
## 273 10.280662 12.574694 9.515985 6.542240 6.117419 4.757992 4.842957 4.418136
## 274 7.363752 8.887286 6.601984 5.586294 4.824527 6.348062 5.840217 4.316682
## 275 6.609564 6.058767 5.673209 3.470021 3.690340 3.580180 2.864144 1.597311
## dthrt31 dthrt32
## 270 4.651835 1.9586672
## 271 2.572546 0.9647047
## 272 4.003824 0.5004780
## 273 2.803817 1.5293547
## 274 3.047070 0.7617674
## 275 1.597311 0.9363549To calculate the COVID-19 rate, how fast it spreads in municipality and state level, we will be looking at the cumulative data, for it records the total number of cases in a municipality, whether they are recovered, death cases or still active. So, we will be creating a new data.frame that consists cumulative values.
COVID19RATE <- cbind(MXCO[,1:6], MXCO[,27:46])
head(COVID19RATE)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 cumul13
## 270 09002 09 002 Azcapotzalco 414711 408441 19
## 271 09003 09 003 Coyoacán 620416 621952 25
## 272 09004 09 004 Cuajimalpa de Morelos 186391 199809 66
## 273 09005 09 005 Gustavo A. Madero 1185772 1176967 28
## 274 09006 09 006 Iztacalco 384326 393821 11
## 275 09007 09 007 Iztapalapa 1815786 1815551 41
## cumul14 cumul15 cumul16 cumul17 cumul18 cumul19 cumul20 cumul21 cumul22
## 270 46 87 156 254 391 615 931 1267 1590
## 271 44 107 213 375 554 831 1150 1492 1866
## 272 75 93 122 157 202 260 323 414 573
## 273 96 193 397 734 1107 1653 2330 3108 3843
## 274 27 65 188 343 552 806 1125 1430 1706
## 275 100 225 526 1061 1754 2617 3583 4659 5532
## cumul23 cumul24 cumul25 cumul26 cumul27 cumul28 cumul29 cumul30 cumul31
## 270 1939 2249 2594 2922 3238 3539 3839 4171 4442
## 271 2255 2554 2936 3208 3590 3903 4393 4820 5216
## 272 714 851 991 1117 1285 1444 1615 1835 1949
## 273 4674 5366 5977 6502 7075 7587 8212 8875 9380
## 274 1989 2209 2454 2660 2825 2997 3183 3416 3612
## 275 6416 7234 7981 8646 9408 10078 10822 11550 12170
## cumul32
## 270 4466
## 271 5246
## 272 1956
## 273 9439
## 274 3674
## 275 12305Let’s check if COVID19RATE data.frame and MXCO SpatialPolygonsDataFrame have the same number of observations by checking the total number of rows.
nrow(MXCO)## [1] 176nrow(COVID19RATE)## [1] 176Now that the 2 data.frame corresponds to each other, we can do a facet plot to see how the cumulative COVID-19 cases increases from e-week 13 to e-week 32.
COVID19.map <- tm_shape(COVID19RATE)+
tm_fill(col = c("cumul13", "cumul14", "cumul15", "cumul16", "cumul17", "cumul18", "cumul19", "cumul20", "cumul21", "cumul22", "cumul23", "cumul24", "cumul25", "cumul26", "cumul27", "cumul28", "cumul29", "cumul30", "cumul31", "cumul32"),
n =5,
palette = "Blues",
style = "jenks",
title = "Cumulative COVID-19 cases") +
tm_borders(alpha = 0.5) +
tm_layout(legend.height = 0.35,
legend.width = 0.45,
frame = TRUE,
panel.labels = c("w13 cumul cases", "w14 cumul cases", "w15 cumul cases", "w16 cumul cases", "w17 cumul cases", "w18 cumul cases", "w19 cumul cases", "w20 cumul cases", "w21 cumul cases", "w22 cumul cases", "w23 cumul cases", "w24 cumul cases", "w25 cumul cases", "w26 cumul cases", "w27 cumul cases", "w28 cumul cases", "w29 cumul cases", "w30 cumul cases", "w31 cumul cases", "w32 cumul cases"))
COVID19.map
We can see that there’s no unexpected cumulative results since all the plots look similar and its colour more intense as each e-week pass by.
Then, we can create new columns for COVID-19 infection rate for each e-week, number of cases of COVID-19 per 10,000 population.
COVID19RATE$Pop2020 <- as.double(COVID19RATE$Pop2020)
COVID19RATE <- COVID19RATE %>%
mutate(RATE13 = round((COVID19RATE$cumul13/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE14 = round((COVID19RATE$cumul14/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE15 = round((COVID19RATE$cumul15/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE16 = round((COVID19RATE$cumul16/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE17 = round((COVID19RATE$cumul17/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE18 = round((COVID19RATE$cumul18/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE19 = round((COVID19RATE$cumul19/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE20 = round((COVID19RATE$cumul20/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE21 = round((COVID19RATE$cumul21/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE22 = round((COVID19RATE$cumul22/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE23 = round((COVID19RATE$cumul23/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE24 = round((COVID19RATE$cumul24/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE25 = round((COVID19RATE$cumul25/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE26 = round((COVID19RATE$cumul26/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE27 = round((COVID19RATE$cumul27/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE28 = round((COVID19RATE$cumul28/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE29 = round((COVID19RATE$cumul29/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE30 = round((COVID19RATE$cumul30/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE31 = round((COVID19RATE$cumul31/(COVID19RATE$Pop2020/10000)),2))
COVID19RATE <- COVID19RATE %>%
mutate(RATE32 = round((COVID19RATE$cumul32/(COVID19RATE$Pop2020/10000)),2))
head(COVID19RATE)## # A tibble: 6 x 46
## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 cumul13 cumul14 cumul15 cumul16
## <chr> <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 09002 09 002 Azcap~ 414711 408441 19 46 87 156
## 2 09003 09 003 Coyoa~ 620416 621952 25 44 107 213
## 3 09004 09 004 Cuaji~ 186391 199809 66 75 93 122
## 4 09005 09 005 Gusta~ 1185772 1176967 28 96 193 397
## 5 09006 09 006 Iztac~ 384326 393821 11 27 65 188
## 6 09007 09 007 Iztap~ 1815786 1815551 41 100 225 526
## # ... with 36 more variables: cumul17 <dbl>, cumul18 <dbl>, cumul19 <dbl>,
## # cumul20 <dbl>, cumul21 <dbl>, cumul22 <dbl>, cumul23 <dbl>, cumul24 <dbl>,
## # cumul25 <dbl>, cumul26 <dbl>, cumul27 <dbl>, cumul28 <dbl>, cumul29 <dbl>,
## # cumul30 <dbl>, cumul31 <dbl>, cumul32 <dbl>, RATE13 <dbl>, RATE14 <dbl>,
## # RATE15 <dbl>, RATE16 <dbl>, RATE17 <dbl>, RATE18 <dbl>, RATE19 <dbl>,
## # RATE20 <dbl>, RATE21 <dbl>, RATE22 <dbl>, RATE23 <dbl>, RATE24 <dbl>,
## # RATE25 <dbl>, RATE26 <dbl>, RATE27 <dbl>, RATE28 <dbl>, RATE29 <dbl>,
## # RATE30 <dbl>, RATE31 <dbl>, RATE32 <dbl>3. Spatio-temporal distribution of COVID-19 rates at municipality level
Let’s plot the COVID-19 cases per 10,000 population (rate) for each e-week.
COVID19RATE.map <- tm_shape(COVID19RATE)+
tm_fill(col = c("RATE13", "RATE14", "RATE15", "RATE16", "RATE17", "RATE18", "RATE19", "RATE20", "RATE21", "RATE22", "RATE23", "RATE24", "RATE25", "RATE26", "RATE27", "RATE28", "RATE29", "RATE30", "RATE31", "RATE32"),
n =5,
style = "jenks")+
tm_borders(alpha = 0.5) +
tm_layout(legend.height = 0.35,
legend.width = 0.45,
frame = TRUE,
panel.labels = c("w13 rate", "w14 rate", "w15 rate", "w16 rate", "w17 rate", "w18 rate", "w19 rate", "w20 rate", "w21 rate", "w22 rate", "w23 rate", "w24 rate", "w25 rate", "w26 rate", "w27 rate", "w28 rate", "w29 rate", "w30 rate", "w31 rate", "w32 rate"))
COVID19RATE.map
We can see some changes in the darker orange municipalities. In e-week 13, the infection rate per 10,000 population is still little but going into e-week 17, the whole Mexico City is already high in COVID-19 infections. This trend continues on until e-week 32 with increasing infections.
Then, we will plot the start e-week 13 & final e-week 32 COVID-19 infection rate against the population data in 2020 to see the differences.
COVID19RATE13.map <- tm_shape(COVID19RATE) +
tm_fill(col = "RATE13",
n = 5,
style = "jenks",
title = "COVID-19 rate in e-week 13") +
tm_borders(alpha = 0.5) +
tm_layout(legend.height = 0.35,
legend.width = 0.45,
frame = FALSE)
COVID19RATE32.map <- tm_shape(COVID19RATE) +
tm_fill(col = "RATE32",
n = 5,
style = "jenks",
title = "COVID-19 rate in e-week 32") +
tm_borders(alpha = 0.5) +
tm_layout(legend.height = 0.35,
legend.width = 0.45,
frame = FALSE)
POP2020.map <- tm_shape(COVID19RATE)+
tm_fill(col = "Pop2020",
n =5,
style = "jenks",
palette = "Greens",
title = "Population in 2020") +
tm_borders(alpha = 0.5) +
tm_layout(legend.height = 0.35,
legend.width = 0.45,
frame = FALSE) tmap_arrange(COVID19RATE32.map, POP2020.map, asp = NA, ncol = 2)
The spatio-temporal patterns observations are:
Firstly, we can see that there is kind of cluster in the COVID-19 rate (cases per 10,000 population) which is mainly located in Mexico City
Secondly, we can see that the municipalities with higher population does not necessarily mean higher rate of COVID-19 infection (disproportionate)
Thirdly, we can somehow see that the COVID-19 infection rate travels up from Mexico City to Mexico State and not south even though it’s closer to Milpa Alta municipality in Mexico City
There seems to be a new possible cluster near Atizapán municipality in Mexico State
4. Local Moran’s I analysis
Moran’s I is a measure of the degree to which the value at a target site is similar to values at adjacent sites.
Moran’s I is large & positive when the value for a given target (or for all locations in the global case) is similar to adjacent values and negative when the value at a target is dissimilar to adjacent values.
We will be using the Inverse Distance Spatial Weighting method (IDW) and these are the required steps:
Create Queen contiguity based neighbours using poly2nb() function of spdep package
Derive spatial weight matrices using nbdists() function of spdep package
Create spatial weight matrix using nb2listw() function of spdep package which will supplement the neighbours list created with spatial weights matrices from nbdists() function for the chosen coding scheme (Style=“B” means basic binary coding which is a more robust option than Style=“W”)
4.1 Deriving spatial weights matrix
wm_q <- poly2nb(COVID19RATE, queen = TRUE)
summary(wm_q)## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 962
## Percentage nonzero weights: 3.10563
## Average number of links: 5.465909
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 14
## 3 3 18 41 31 31 22 13 9 3 1 1
## 3 least connected regions:
## 721 739 765 with 1 link
## 1 most connected region:
## 761 with 14 linksFrom the summary, we can find out that:
There are 176 municipalities in the three states in Mexico, Mexico City, Mexico State & Morelos State
Most connected area unit has 14 neighbours which is shown by the first row “14” (no. of neighbours) with the next row showing 1 meaning only 1 municipality has 11 neighbours
There are 3 municipalities with only 1 neighbour, which must be the bigger municipalities at the edge of the study area
plot(COVID19RATE, border = 'lightgrey')
plot(wm_q, coordinates(COVID19RATE), pch = 19, cex = 0.6, add = TRUE, col = "red")
We can see from the plot that there are so much more links in the centre of the map compared to the edges, we choose Inverse Distance Weighting (IDW) method because it is the most appropriate method to model continuous data whereby the closer two municipalities are, the more likely they are to influence each other in terms of their COVID-19 infection rate.
Next, we will be applying the function to inverse the value of the distance matrix of our neighbour list.
dist <- nbdists(wm_q, coordinates(COVID19RATE), longlat = TRUE)
ids <- lapply(dist, function(x) 1/(x))
ids[1]## [[1]]
## [1] 0.13780415 0.14476598 0.14784979 0.08542328 0.14933148dist[1]## [[1]]
## [1] 7.256676 6.907700 6.763621 11.706411 6.696512For IDW method, we can see that by inversing the distance, the larger distances get scaled down and the smaller distances get scaled up.
4.2 Row-standardisation spatial weight matrix
We will derive the row standardised distance weight matrix using the Style argument “B” which means basic binary coding. This will give us the final spatial weight matrix!
wm_ids <- nb2listw(wm_q, glist = ids, style = "B", zero.policy = TRUE)
wm_ids## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 962
## Percentage nonzero weights: 3.10563
## Average number of links: 5.465909
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 176 30976 87.66922 21.54584 211.6514To see the final summary of our spatial weight matrix, we unlist the weights variable of wm_ids.
summary(unlist(wm_ids$weights))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02925 0.05687 0.07776 0.09113 0.10860 0.489674.3 Computing Local Moran’s I
We use the localmoran() function of spdep package to compute the Local Moran’s I values as well as the z-value and requires variable to be used & the spatial weights as it’s arguments.
index <- order(COVID19RATE$NOMGEO)
localMI.13 <- localmoran(COVID19RATE$RATE13, wm_ids)
localMI.32 <- localmoran(COVID19RATE$RATE32, wm_ids)
head(localMI.32)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 270 2.617647 -0.003800998 0.08497784 8.992665 1.206530e-19
## 271 1.862400 -0.002954675 0.05483002 7.966220 8.180098e-16
## 272 0.990469 -0.002450652 0.04792994 4.535349 2.875421e-06
## 273 1.806102 -0.004546294 0.07600730 6.567596 2.556709e-11
## 274 3.026496 -0.004730396 0.14042562 8.089007 3.007658e-16
## 275 1.676524 -0.004528051 0.07578346 6.106522 5.091270e-10Using the printCoefmat() function, we can see see the content of the Local Moran matrix.
printCoefmat(data.frame(localMI.32[index,], row.names = COVID19RATE$NOMGEO[index]), check.names = FALSE)## Ii E.Ii Var.Ii Z.Ii
## Ã\201lvaro Obregón 2.38020049 -0.00386484 0.06621313 9.26501418
## Acambay de RuÃz Castañeda 0.06136346 -0.00129141 0.01195643 0.57299919
## Acolman 0.00799447 -0.00460437 0.10550256 0.03878817
## Aculco 0.13360080 -0.00124929 0.01258535 1.20203902
## Almoloya de Alquisiras 0.08185568 -0.00158860 0.01858622 0.61206987
## Almoloya de Juárez 0.04225081 -0.00178784 0.01321903 0.38303124
## Almoloya del RÃo 1.35632113 -0.00620682 0.34551872 2.31797971
## Amacuzac 0.10162782 -0.00189672 0.02649196 0.63604273
## Amanalco -0.00447816 -0.00196824 0.01581055 -0.01996121
## Amatepec 0.08996570 -0.00080847 0.00653661 1.12275782
## Amecameca -0.01755715 -0.00207982 0.03176566 -0.08683946
## Apaxco 0.03977394 -0.00108437 0.01763674 0.30766013
## Atenco 0.03927278 -0.00385999 0.07556961 0.15690391
## Atizapán 1.32194257 -0.00463191 0.27570521 2.52644094
## Atizapán de Zaragoza 0.00412558 -0.00285147 0.04030208 0.03475425
## Atlacomulco -0.03118129 -0.00171410 0.01837432 -0.21738680
## Atlatlahucan 0.13934441 -0.00371067 0.05926392 0.58763560
## Atlautla 0.10764287 -0.00256022 0.04225299 0.53612376
## Axapusco 0.07621078 -0.00229359 0.04668231 0.36334370
## Axochiapan 0.01003688 -0.00102580 0.01010913 0.11002811
## Ayala 0.08384543 -0.00283760 0.02627007 0.53481489
## Ayapango 0.08253021 -0.00340289 0.05870173 0.35467864
## Azcapotzalco 2.61764720 -0.00380100 0.08497784 8.99266537
## Benito Juárez 1.69692774 -0.00459620 0.10548393 5.23895795
## Calimaya 0.02120395 -0.00427997 0.06855305 0.09733137
## Capulhuac -0.00186167 -0.00204543 0.04484618 0.00086773
## Chalco 0.51801542 -0.00475898 0.09686377 1.67970734
## Chapa de Mota 0.12016965 -0.00142458 0.01524842 0.98469232
## Chapultepec 0.25188090 -0.00353073 0.11519140 0.75254183
## Chiautla 0.14707439 -0.00765834 0.34297803 0.26420987
## Chicoloapan 0.02494544 -0.00251684 0.04678299 0.12696752
## Chiconcuac 0.07491139 -0.00379657 0.20606720 0.17338607
## Chimalhuacán 0.00525783 -0.00283469 0.06170883 0.03257689
## Coacalco de Berriozábal 0.01429395 -0.00297344 0.06881632 0.06582349
## Coatepec Harinas 0.20803516 -0.00255703 0.02539654 1.32146303
## Coatetelco 0.35699575 -0.00325456 0.05401282 1.55008657
## Coatlán del RÃo 0.16221413 -0.00312752 0.05614879 0.69776962
## Cocotitlán 0.25336203 -0.00296358 0.09869884 0.81589818
## Coyoacán 1.86240022 -0.00295468 0.05483002 7.96621971
## Coyotepec -0.02580396 -0.00289514 0.07121996 -0.08584239
## Cuajimalpa de Morelos 0.99046900 -0.00245065 0.04792994 4.53534859
## Cuauhtémoc 3.43086489 -0.00525012 0.13253662 9.43842986
## Cuautitlán -0.01877543 -0.00669645 0.22277087 -0.02559180
## Cuautitlán Izcalli -0.02332615 -0.00309760 0.04770386 -0.09261642
## Cuautla -0.02248695 -0.00190254 0.02611848 -0.12736926
## Cuernavaca 0.05454788 -0.00281957 0.03307663 0.31543140
## Donato Guerra 0.14582072 -0.00170481 0.02294497 0.97392021
## Ecatepec de Morelos 0.04045010 -0.00469235 0.06069097 0.18324116
## Ecatzingo 0.19049185 -0.00239743 0.04527688 0.90650426
## El Oro -0.02026502 -0.00117239 0.01035819 -0.18759620
## Emiliano Zapata 0.19033174 -0.00310851 0.04949374 0.86950424
## Gustavo A. Madero 1.80610176 -0.00454629 0.07600730 6.56759558
## Huehuetoca 0.01262405 -0.00231636 0.05133804 0.06593908
## Hueypoxtla 0.01974638 -0.00119875 0.01391152 0.17758068
## Huitzilac -0.18938502 -0.00205173 0.02025181 -1.31638546
## Huixquilucan 0.06984377 -0.00245819 0.03826368 0.36962109
## Isidro Fabela 0.08043917 -0.00239363 0.03729053 0.42894644
## Ixtapaluca -0.00153226 -0.00232361 0.02638016 0.00487229
## Ixtapan de la Sal 0.11623761 -0.00204762 0.02541464 0.74197380
## Ixtapan del Oro 0.12114453 -0.00142692 0.02109019 0.84401284
## Ixtlahuaca 0.06322399 -0.00153356 0.01384202 0.55041563
## Iztacalco 3.02649608 -0.00473040 0.14042562 8.08900689
## Iztapalapa 1.67652404 -0.00452805 0.07578346 6.10652234
## Jaltenco 0.03061434 -0.00436343 0.13386612 0.09559975
## Jantetelco 0.20688160 -0.00231512 0.05098795 0.92644813
## Jilotepec 0.07723975 -0.00155951 0.01737532 0.59779974
## Jilotzingo 0.00820772 -0.00257035 0.04041411 0.05361353
## Jiquipilco 0.11947638 -0.00181311 0.01688505 0.93340904
## Jiutepec 0.10751526 -0.00212571 0.03356214 0.59847780
## Jocotitlán 0.02814105 -0.00218384 0.02198332 0.20452789
## Jojutla 0.00425599 -0.00218993 0.03705342 0.03348661
## Jonacatepec de Leandro Valle 0.12032406 -0.00216015 0.04498405 0.57749849
## Joquicingo 0.03714787 -0.00290584 0.04552752 0.18771793
## Juchitepec -0.18927030 -0.00455987 0.06873750 -0.70452198
## La Magdalena Contreras 1.27411425 -0.00168860 0.02904368 7.48613650
## La Paz -0.00814395 -0.00355034 0.06546356 -0.01795374
## Lerma 0.00986002 -0.00390911 0.04537204 0.06464164
## Luvianos 0.01576115 -0.00081819 0.00642210 0.20688488
## Malinalco 0.10782284 -0.00226254 0.02504132 0.69566635
## Mazatepec 0.25932940 -0.00368770 0.07673494 0.94948349
## Melchor Ocampo 0.09904773 -0.00452924 0.21900168 0.22132960
## Metepec 0.25348635 -0.00517980 0.13710480 0.69857538
## Mexicaltzingo 0.34596878 -0.00374397 0.14795388 0.90917667
## Miacatlán 0.17363485 -0.00287529 0.03468317 0.94778591
## Miguel Hidalgo 2.51404377 -0.00444575 0.08996224 8.39672650
## Milpa Alta 2.43921362 -0.00273028 0.02745364 14.73790343
## Morelos 0.12535994 -0.00197211 0.01896642 0.92458116
## Naucalpan de Juárez 0.20846575 -0.00360058 0.04811617 0.96677673
## Nextlalpan -0.00947141 -0.00666156 0.19146763 -0.00642148
## Nezahualcóyotl 0.45020584 -0.00442395 0.07564713 1.65295745
## Nicolás Romero 0.05647386 -0.00257649 0.02999081 0.34097956
## Nopaltepec 0.07857198 -0.00104454 0.03172790 0.44697426
## Ocoyoacac -0.00171138 -0.00350696 0.04724101 0.00826122
## Ocuilan 0.07538681 -0.00264950 0.02480655 0.49546570
## Ocuituco 0.20460006 -0.00255808 0.04725138 0.95300406
## Otumba -0.00010526 -0.00189753 0.02625912 0.01106024
## Otzoloapan 0.12462386 -0.00207179 0.03177160 0.71079136
## Otzolotepec 0.00082881 -0.00295927 0.04853448 0.01719468
## Ozumba 0.08408000 -0.00424363 0.10738834 0.26952455
## Papalotla 0.04978964 -0.00299141 0.11750898 0.15397225
## Polotitlán 0.07833260 -0.00064691 0.00741448 0.91722138
## Puente de Ixtla -0.04896251 -0.00243815 0.03200957 -0.26004023
## Rayón 0.00529871 -0.00371301 0.11428388 0.02665724
## San Antonio la Isla 0.10192193 -0.00542473 0.15991948 0.26843419
## San Felipe del Progreso 0.07265408 -0.00166130 0.01345954 0.64056546
## San José del Rincón 0.07906399 -0.00105041 0.00808778 0.89083189
## San MartÃn de las Pirámides 0.00052017 -0.00270478 0.04977896 0.01445436
## San Mateo Atenco 0.03595997 -0.00200315 0.04533043 0.17830643
## San Simón de Guerrero 0.00774519 -0.00125368 0.01618593 0.07073255
## Santo Tomás 0.09104815 -0.00163946 0.02832092 0.55076708
## Soyaniquilpan de Juárez 0.02377180 -0.00060928 0.01079509 0.23466037
## Sultepec 0.11626651 -0.00126147 0.00942164 1.21081564
## Tecámac 0.00073049 -0.00387024 0.05569136 0.01949540
## Tejupilco -0.01808392 -0.00165502 0.01135312 -0.15418811
## Temamatla 0.12858372 -0.00400037 0.11138719 0.39725904
## Temascalapa -0.00856599 -0.00153749 0.01739237 -0.05329463
## Temascalcingo 0.00988376 -0.00112626 0.00979216 0.11126256
## Temascaltepec 0.08570288 -0.00245761 0.02103179 0.60790526
## Temixco 0.18398609 -0.00268477 0.04131283 0.91840458
## Temoac 0.28136844 -0.00227448 0.04278327 1.37130782
## Temoaya 0.08579648 -0.00238752 0.02805923 0.52644365
## Tenancingo -0.00771465 -0.00239353 0.02800844 -0.03179498
## Tenango del Aire -0.04349932 -0.00381235 0.09897788 -0.12614756
## Tenango del Valle 0.00150220 -0.00370043 0.04732709 0.02391486
## Teoloyucan -0.00576389 -0.00271537 0.05957713 -0.01248961
## Teotihuacán -0.04056023 -0.00331471 0.06166197 -0.14999105
## Tepalcingo 0.03525654 -0.00137852 0.01397023 0.30995246
## Tepetlaoxtoc 0.03052375 -0.00304781 0.03974729 0.16839055
## Tepetlixpa 0.14931699 -0.00381514 0.10940964 0.46295472
## Tepotzotlán 0.01710402 -0.00293845 0.03560074 0.10622374
## Tepoztlán -0.09557419 -0.00270789 0.03001186 -0.53605789
## Tequixquiac 0.04669538 -0.00188856 0.02669927 0.29733292
## Tetecala 0.17478824 -0.00265308 0.07454476 0.64989951
## Tetela del Volcán 0.19653648 -0.00235167 0.04206356 0.96974129
## Texcaltitlán 0.08579127 -0.00215508 0.02422806 0.56501315
## Texcalyacac -0.01757667 -0.00530230 0.15569499 -0.03110726
## Texcoco 0.02961380 -0.00421520 0.05322196 0.14663706
## Tezoyuca 0.01530760 -0.00322754 0.10328133 0.05767468
## Tianguistenco -0.20321233 -0.00726855 0.12497914 -0.55425894
## Timilpan 0.09390960 -0.00181881 0.01595034 0.75797671
## Tláhuac 2.28965046 -0.00283078 0.05418360 9.84853993
## Tlalmanalco 0.05544390 -0.00304555 0.03871812 0.29724901
## Tlalnepantla -0.08507526 -0.00264343 0.04095356 -0.40733254
## Tlalnepantla de Baz 0.23025950 -0.00402794 0.07035599 0.88328010
## Tlalpan 2.27256887 -0.00325104 0.03477691 12.20372479
## Tlaltizapán de Zapata 0.08780002 -0.00263435 0.02910882 0.53005531
## Tlaquiltenango 0.01633319 -0.00155268 0.01152247 0.16662389
## Tlatlaya 0.06230551 -0.00032514 0.00307416 1.12959788
## Tlayacapan -0.07643449 -0.00285448 0.04835054 -0.33462557
## Toluca 0.00618236 -0.00303853 0.03143444 0.05200799
## Tonanitla 0.01621607 -0.00330746 0.06856852 0.07455830
## Tonatico 0.03763913 -0.00122949 0.02187540 0.26279746
## Totolapan 0.15432419 -0.00337669 0.06558662 0.61578128
## Tultepec -0.00203543 -0.00505462 0.15898588 0.00757201
## Tultitlán -0.01373999 -0.00559929 0.09637321 -0.02622310
## Valle de Bravo -0.00910735 -0.00220108 0.01965712 -0.04925881
## Valle de Chalco Solidaridad -0.35433852 -0.00276965 0.05230162 -1.53727952
## Venustiano Carranza 2.81423859 -0.00394928 0.12464823 7.98227876
## Villa de Allende 0.14552383 -0.00144996 0.01776699 1.10263811
## Villa del Carbón 0.12625206 -0.00169946 0.01699185 0.98157862
## Villa Guerrero 0.09989074 -0.00202697 0.02442637 0.65210897
## Villa Victoria 0.10907166 -0.00141346 0.01167682 1.02244863
## Xalatlaco -0.00717122 -0.00173877 0.02537388 -0.03410376
## Xochimilco 3.98440891 -0.00246553 0.03497573 21.31813352
## Xochitepec 0.21773490 -0.00299211 0.05279423 0.96064348
## Xonacatlán -0.02493190 -0.00222945 0.04000400 -0.11350657
## Xoxocotla 0.12377171 -0.00357876 0.06894761 0.48499909
## Yautepec 0.13332166 -0.00319509 0.03733060 0.70656721
## Yecapixtla 0.20863692 -0.00380613 0.04983446 0.95165085
## Zacatepec -0.08930214 -0.00248969 0.05350497 -0.37530551
## Zacazonapan 0.08472760 -0.00172558 0.02099237 0.59669196
## Zacualpan 0.12742975 -0.00116665 0.00996065 1.28850135
## Zacualpan de Amilpas 0.27041029 -0.00248266 0.04717361 1.25644292
## Zinacantepec 0.00175722 -0.00155826 0.01186718 0.03043488
## Zumpahuacán 0.08387945 -0.00244315 0.03107834 0.48966112
## Zumpango -0.01265014 -0.00424784 0.06359480 -0.03331865
## Pr.z...0.
## Ã\201lvaro Obregón 0.0000
## Acambay de RuÃz Castañeda 0.2833
## Acolman 0.4845
## Aculco 0.1147
## Almoloya de Alquisiras 0.2702
## Almoloya de Juárez 0.3508
## Almoloya del RÃo 0.0102
## Amacuzac 0.2624
## Amanalco 0.5080
## Amatepec 0.1308
## Amecameca 0.5346
## Apaxco 0.3792
## Atenco 0.4377
## Atizapán 0.0058
## Atizapán de Zaragoza 0.4861
## Atlacomulco 0.5860
## Atlatlahucan 0.2784
## Atlautla 0.2959
## Axapusco 0.3582
## Axochiapan 0.4562
## Ayala 0.2964
## Ayapango 0.3614
## Azcapotzalco 0.0000
## Benito Juárez 0.0000
## Calimaya 0.4612
## Capulhuac 0.4997
## Chalco 0.0465
## Chapa de Mota 0.1624
## Chapultepec 0.2259
## Chiautla 0.3958
## Chicoloapan 0.4495
## Chiconcuac 0.4312
## Chimalhuacán 0.4870
## Coacalco de Berriozábal 0.4738
## Coatepec Harinas 0.0932
## Coatetelco 0.0606
## Coatlán del RÃo 0.2427
## Cocotitlán 0.2073
## Coyoacán 0.0000
## Coyotepec 0.5342
## Cuajimalpa de Morelos 0.0000
## Cuauhtémoc 0.0000
## Cuautitlán 0.5102
## Cuautitlán Izcalli 0.5369
## Cuautla 0.5507
## Cuernavaca 0.3762
## Donato Guerra 0.1650
## Ecatepec de Morelos 0.4273
## Ecatzingo 0.1823
## El Oro 0.5744
## Emiliano Zapata 0.1923
## Gustavo A. Madero 0.0000
## Huehuetoca 0.4737
## Hueypoxtla 0.4295
## Huitzilac 0.9060
## Huixquilucan 0.3558
## Isidro Fabela 0.3340
## Ixtapaluca 0.4981
## Ixtapan de la Sal 0.2291
## Ixtapan del Oro 0.1993
## Ixtlahuaca 0.2910
## Iztacalco 0.0000
## Iztapalapa 0.0000
## Jaltenco 0.4619
## Jantetelco 0.1771
## Jilotepec 0.2750
## Jilotzingo 0.4786
## Jiquipilco 0.1753
## Jiutepec 0.2748
## Jocotitlán 0.4190
## Jojutla 0.4866
## Jonacatepec de Leandro Valle 0.2818
## Joquicingo 0.4255
## Juchitepec 0.7594
## La Magdalena Contreras 0.0000
## La Paz 0.5072
## Lerma 0.4742
## Luvianos 0.4180
## Malinalco 0.2433
## Mazatepec 0.1712
## Melchor Ocampo 0.4124
## Metepec 0.2424
## Mexicaltzingo 0.1816
## Miacatlán 0.1716
## Miguel Hidalgo 0.0000
## Milpa Alta 0.0000
## Morelos 0.1776
## Naucalpan de Juárez 0.1668
## Nextlalpan 0.5026
## Nezahualcóyotl 0.0492
## Nicolás Romero 0.3666
## Nopaltepec 0.3274
## Ocoyoacac 0.4967
## Ocuilan 0.3101
## Ocuituco 0.1703
## Otumba 0.4956
## Otzoloapan 0.2386
## Otzolotepec 0.4931
## Ozumba 0.3938
## Papalotla 0.4388
## Polotitlán 0.1795
## Puente de Ixtla 0.6026
## Rayón 0.4894
## San Antonio la Isla 0.3942
## San Felipe del Progreso 0.2609
## San José del Rincón 0.1865
## San MartÃn de las Pirámides 0.4942
## San Mateo Atenco 0.4292
## San Simón de Guerrero 0.4718
## Santo Tomás 0.2909
## Soyaniquilpan de Juárez 0.4072
## Sultepec 0.1130
## Tecámac 0.4922
## Tejupilco 0.5613
## Temamatla 0.3456
## Temascalapa 0.5213
## Temascalcingo 0.4557
## Temascaltepec 0.2716
## Temixco 0.1792
## Temoac 0.0851
## Temoaya 0.2993
## Tenancingo 0.5127
## Tenango del Aire 0.5502
## Tenango del Valle 0.4905
## Teoloyucan 0.5050
## Teotihuacán 0.5596
## Tepalcingo 0.3783
## Tepetlaoxtoc 0.4331
## Tepetlixpa 0.3217
## Tepotzotlán 0.4577
## Tepoztlán 0.7040
## Tequixquiac 0.3831
## Tetecala 0.2579
## Tetela del Volcán 0.1661
## Texcaltitlán 0.2860
## Texcalyacac 0.5124
## Texcoco 0.4417
## Tezoyuca 0.4770
## Tianguistenco 0.7103
## Timilpan 0.2242
## Tláhuac 0.0000
## Tlalmanalco 0.3831
## Tlalnepantla 0.6581
## Tlalnepantla de Baz 0.1885
## Tlalpan 0.0000
## Tlaltizapán de Zapata 0.2980
## Tlaquiltenango 0.4338
## Tlatlaya 0.1293
## Tlayacapan 0.6310
## Toluca 0.4793
## Tonanitla 0.4703
## Tonatico 0.3964
## Totolapan 0.2690
## Tultepec 0.4970
## Tultitlán 0.5105
## Valle de Bravo 0.5196
## Valle de Chalco Solidaridad 0.9379
## Venustiano Carranza 0.0000
## Villa de Allende 0.1351
## Villa del Carbón 0.1632
## Villa Guerrero 0.2572
## Villa Victoria 0.1533
## Xalatlaco 0.5136
## Xochimilco 0.0000
## Xochitepec 0.1684
## Xonacatlán 0.5452
## Xoxocotla 0.3138
## Yautepec 0.2399
## Yecapixtla 0.1706
## Zacatepec 0.6463
## Zacazonapan 0.2754
## Zacualpan 0.0988
## Zacualpan de Amilpas 0.1045
## Zinacantepec 0.4879
## Zumpahuacán 0.3122
## Zumpango 0.5133Let’s bind both the localMI for e-week 13 and e-week 32 into our COVID19RATE SpatialPolygonsDataFrame to create COVID19RATE.localMI.
COVID19RATE.localMI <- cbind(COVID19RATE, localMI.13)
names(COVID19RATE.localMI)[47] <- "Ii_13"
COVID19RATE.localMI <- cbind(COVID19RATE.localMI, localMI.32)
names(COVID19RATE.localMI)[52] <- "Ii_32"
COVID19RATE.localMI <- COVID19RATE.localMI %>%
arrange(desc(Ii_32))
head(COVID19RATE.localMI)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 cumul13 cumul14
## 281 09013 09 013 Xochimilco 415007 418060 11 32
## 283 09015 09 015 Cuauhtémoc 531831 548606 33 62
## 274 09006 09 006 Iztacalco 384326 393821 11 27
## 285 09017 09 017 Venustiano Carranza 430978 433231 9 35
## 270 09002 09 002 Azcapotzalco 414711 408441 19 46
## 284 09016 09 016 Miguel Hidalgo 372889 379624 112 143
## cumul15 cumul16 cumul17 cumul18 cumul19 cumul20 cumul21 cumul22 cumul23
## 281 72 137 288 440 688 1018 1361 1725 2152
## 283 136 233 367 546 787 1090 1436 1744 2093
## 274 65 188 343 552 806 1125 1430 1706 1989
## 285 75 149 291 429 680 990 1309 1598 1858
## 270 87 156 254 391 615 931 1267 1590 1939
## 284 175 238 342 453 592 778 996 1199 1383
## cumul24 cumul25 cumul26 cumul27 cumul28 cumul29 cumul30 cumul31 cumul32
## 281 2588 2911 3252 3540 3880 4321 4698 5034 5086
## 283 2403 2707 2989 3254 3512 3829 4148 4353 4387
## 274 2209 2454 2660 2825 2997 3183 3416 3612 3674
## 285 2153 2363 2589 2858 3117 3396 3724 4027 4048
## 270 2249 2594 2922 3238 3539 3839 4171 4442 4466
## 284 1617 1865 2048 2246 2427 2711 2989 3183 3214
## RATE13 RATE14 RATE15 RATE16 RATE17 RATE18 RATE19 RATE20 RATE21 RATE22
## 281 0.26 0.77 1.72 3.28 6.89 10.52 16.46 24.35 32.56 41.26
## 283 0.60 1.13 2.48 4.25 6.69 9.95 14.35 19.87 26.18 31.79
## 274 0.28 0.69 1.65 4.77 8.71 14.02 20.47 28.57 36.31 43.32
## 285 0.21 0.81 1.73 3.44 6.72 9.90 15.70 22.85 30.21 36.89
## 270 0.47 1.13 2.13 3.82 6.22 9.57 15.06 22.79 31.02 38.93
## 284 2.95 3.77 4.61 6.27 9.01 11.93 15.59 20.49 26.24 31.58
## RATE23 RATE24 RATE25 RATE26 RATE27 RATE28 RATE29 RATE30 RATE31 RATE32
## 281 51.48 61.90 69.63 77.79 84.68 92.81 103.36 112.38 120.41 121.66
## 283 38.15 43.80 49.34 54.48 59.31 64.02 69.80 75.61 79.35 79.97
## 274 50.51 56.09 62.31 67.54 71.73 76.10 80.82 86.74 91.72 93.29
## 285 42.89 49.70 54.54 59.76 65.97 71.95 78.39 85.96 92.95 93.44
## 270 47.47 55.06 63.51 71.54 79.28 86.65 93.99 102.12 108.75 109.34
## 284 36.43 42.59 49.13 53.95 59.16 63.93 71.41 78.74 83.85 84.66
## Ii_13 E.Ii Var.Ii Z.Ii Pr.z...0. Ii_32
## 281 0.03719376 -0.002465534 0.02941878 0.2312239 4.085704e-01 3.984409
## 283 1.75684627 -0.005250121 0.11160827 5.2745022 6.655846e-08 3.430865
## 274 0.13589918 -0.004730396 0.11805691 0.4092898 3.411635e-01 3.026496
## 285 0.05903000 -0.003949281 0.10467113 0.1946634 4.228282e-01 2.814239
## 270 1.03846696 -0.003800998 0.07146737 3.8987508 4.834510e-05 2.617647
## 284 9.04381770 -0.004445751 0.07578673 32.8676495 3.187145e-237 2.514044
## E.Ii.1 Var.Ii.1 Z.Ii.1 Pr.z...0..1
## 281 -0.002465534 0.03497573 21.318134 3.853574e-101
## 283 -0.005250121 0.13253662 9.438430 1.892035e-21
## 274 -0.004730396 0.14042562 8.089007 3.007658e-16
## 285 -0.003949281 0.12464823 7.982279 7.182807e-16
## 270 -0.003800998 0.08497784 8.992665 1.206530e-19
## 284 -0.004445751 0.08996224 8.396726 2.295489e-174.4 Mapping Local Moran’s I statistic & p-value
localMI.I.13 <- tm_shape(COVID19RATE.localMI) +
tm_fill(col = "Ii_13",
style = "pretty",
palette = "RdBu",
title = "L. Moran's I stats (e-w13)") +
tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE)
localMI.p.13 <- tm_shape(COVID19RATE.localMI) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette = "-Blues",
title = "L. Moran's I p-values (e-w13)") +
tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE)
tmap_arrange(localMI.I.13, localMI.p.13, asp = 1, ncol = 2)## Variable(s) "Ii_13" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
localMI.I.32 <- tm_shape(COVID19RATE.localMI) +
tm_fill(col = "Ii_32",
style = "pretty",
palette = "RdBu",
title = "L. Moran's I statistics (e-w 32)") +
tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE)
localMI.p.32 <- tm_shape(COVID19RATE.localMI) +
tm_fill(col = "Pr.z...0..1",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette = "-Blues",
title = "L. Moran's I p-values (e-w 32)") +
tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE)
tmap_arrange(localMI.I.32, localMI.p.32, asp = 1, ncol = 2)## Variable(s) "Ii_32" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
For the Local Moran’s I statistics map, we can see that there are some municipalities in the Mexico City state that are considered high in Local Moran’s I statistics which means that have positive autocorrelation and their COVID-19 rate values are similar to adjacent values.
We can also see from the Local Moran’s I p-values that the high I statistics municipalities are contained in the cluster of statistically significant dark blue area which depicts that these municipalities’ p-value are less than 0.001 or smaller than the confidence interval of 99.9 %.
Although there are municipalities with negative Local Moran’s I statistics indicating outliers and negative autocorrelation, they are all not stastically significant because their p-values are not less than 0.001 if we set the confidence interval to be 99.9%.
spatio-temporal pattern of COVID-19 infection rate by Local Moran’s I:
From arranging the
COVID19RATE.localMISpatialPolygonsDataFrame in descending order according to their Ii values, we can see that there is a municipality that has Local Moran’s I statistics higher than 4. The municipality is Xochimilco, in Mexico CityAll the municipalities with negative spatial autocorrelation are all not statistically significant
We can see that the municipalities with high COVID-19 infection rate in Mexico City are surrounded by municipalities with negative Local Moran’s I statistics (outliers) and this shows that the COVID-19 infection rate in Mexico City did not influence the infection rate in Mexico State and Moleros State.
5. Local Getis-Ord Gi* analysis
Getis-Ord Gi* statistics will look at neighbours within a defined proximity to identify where high or low values cluster spatially/in space.
A significant positive value of G ( G> 0) indicates that there is a " cluster" of high values of the variable analyzed with reference to their average. Also, a significant but negative value of G (G < 0) is shown as a group of low values relative to the average of the variable analyzed.
Hot spot area = Local Gi>0, Cold spot area = Local Gi<0
These are the required steps:
Deriving distance-based neighbour list using either fixed/adaptive distance weighting scheme
Deriving spatial weight matrix by converting nb object to weight list object
Computing Gi statistics using localG() function of spdep package
5.1 Deriving distance-based neighbour list
We will derive the fixed distance weight matrices using dnearneigh() function of spdep package. This function identifies neighbours of region points by Euclidean distance between lower & upper bounds in km.
We will also derive the adaptive distance weight matrices using knearneigh() function of spdep package. This function returns the neighbour points by specifying the k, number of maximum neighbours.
Fixed Distance Weighting scheme
Determining the cut-off distance for the upper boundary to be used for the dnearneigh() function.
knb1 <- knn2nb(knearneigh(coordinates(COVID19RATE)))
k1dists <- unlist(nbdists(knb1, coordinates(COVID19RATE), longlat = TRUE))
summary(k1dists)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.042 5.751 8.099 8.758 11.052 20.545We can see that the maximum distance between municipalities is 20.545 km so we will round it up and use 21 km as our upper boundary.
dnb <- dnearneigh(coordinates(COVID19RATE), 0, 21, longlat = TRUE)
dnb## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 1788
## Percentage nonzero weights: 5.772211
## Average number of links: 10.15909plot(COVID19RATE, border="light grey")
plot(dnb, coordinates(COVID19RATE), add= TRUE, col = "red")
The map looks unbalanced, let’s check the number of disjoint connected subgraphs in the dnb fixed distance weight matrix (whether all the edges are connected in one network or is there a another network in the same map) using n.comp.nb() function of spdep package.
The result variables:
nc value - number of disjoint connected networks of neighbours
comp.id - vector with indices of the disjoint connected networks of neighbours
n_comp <- n.comp.nb(dnb)
n_comp$nc## [1] 2table(n_comp$comp.id)##
## 1 2
## 173 3We can see from the results that there are 2 network of edges, 1 network with 173 municipalities and another one with 3 municipalities. Before we make our decision, let’s also compute using the Adaptive Distance Weighting scheme (fixed number of neighbours).
Adaptive Distance Weighting scheme
Before we create the neighbour list for adaptive distance weight matrix, we can observe the distribution of our variable of interest, RATE, infection rate of COVID-19 in municipalities per 10,000 population.
ggplot(data = as.data.frame(COVID19RATE$RATE32), aes(x=`COVID19RATE$RATE32`))+
geom_histogram(bins = 20, color = "black", fill = "light blue")
We can see that the distribution is very skewed to the left. In this case, the rule of thumb is to evaluate within the context of at least 8 neighbours.
knb <- knn2nb(knearneigh(coordinates(COVID19RATE), k=8, longlat = TRUE), row.names = row.names(COVID19RATE$RATE))plot(COVID19RATE, border="light grey")
plot(knb, coordinates(COVID19RATE), add= TRUE, col = "red")
Since the plot for adaptive distance weighting scheme and having fixed 8 neighbours for each municipality looks more complete, we will go ahead and use the adaptive proximity matrix, knb.
5.2 Deriving Spatial weight matrix
Once again, we will be using nb2listw() function of spdep package to supplement the neighbours list (knb) created with spatial weights matrices from nbdists() function for the chosen coding scheme. We will be using Style=“B” which means basic binary coding.
knb_wm <- nb2listw(knb, style = "B")
summary(knb_wm)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 1408
## Percentage nonzero weights: 4.545455
## Average number of links: 8
## Non-symmetric neighbours list
## Link number distribution:
##
## 8
## 176
## 176 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 8 links
## 176 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 8 links
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 176 30976 1408 2526 461825.3 Computing Gi statistics
We use the localG() function of spdep package to compute the Local Gi statistics of each municipality based on the spatial weights object used. It will return a vector of G/Gstar values with attributes “gstari” set to TRUE/FALSE, “call” set to the function call and class “localG”. The Gi statistics is represented as Z-score.
We will be computing the Local Gi statistics for the COVID-19 infection rate per 10,000 population for e-week 13 and e-week 32.
index2 <- order(COVID19RATE$NOMGEO)
localGi.adaptive.13 <- localG(COVID19RATE$RATE13, knb_wm)
localGi.adaptive.32 <- localG(COVID19RATE$RATE32, knb_wm)
localGi.adaptive.32## [1] 3.60913038 5.86441740 3.79509782 4.37279789 4.75953237 5.22579544
## [7] 3.56576096 3.00047100 4.67032486 4.28310803 6.16151544 6.25691257
## [13] 5.95461008 5.76653324 5.30195784 5.09683702 -1.17196792 -0.32325102
## [19] -1.49779170 -1.81754139 -0.50670280 1.83124946 -1.38577029 -1.39083928
## [25] -0.72530401 -0.54798922 -0.18528317 1.03070939 1.74317212 -1.26521111
## [31] -1.02576533 0.19074473 0.20984256 1.13689089 2.11573378 0.53485945
## [37] -0.94084240 0.12958411 -0.30519241 -0.26518991 0.18497666 -1.66392887
## [43] 2.72529677 -0.55016204 -0.41248774 -0.59335214 0.90450602 -1.62085756
## [49] 0.36321026 -1.43478898 -0.11025055 -0.53559875 3.58749182 -0.45397248
## [55] 0.13198332 -1.63448410 -1.56663496 -1.40534206 3.60224678 0.34386644
## [61] -1.71640354 -0.12096803 -1.15024478 -1.18845195 1.91217706 -0.30928847
## [67] 0.64787766 -1.12888848 0.09101272 2.76651728 2.61171898 -1.35933682
## [73] 3.09867776 3.17078935 0.33275720 -0.65561103 -0.01517555 2.75460586
## [79] -0.88852753 -1.31835777 0.21450092 -1.16864719 -0.42306327 -1.58406651
## [85] -0.61851981 2.54741561 -1.50990002 2.51222012 2.79846031 -1.12467538
## [91] 0.62853655 2.77610590 -1.46276701 -1.84291163 -1.93533448 -1.53083863
## [97] -0.12807248 -1.78560832 0.95803462 0.65731033 -1.22988580 -0.92602654
## [103] -0.47601333 -0.94974674 0.34371328 0.11775363 0.26744212 -0.40029310
## [109] -0.26657030 -1.48389265 -0.21385058 -0.53326578 -1.91904885 2.07956917
## [115] -0.98086930 -0.22503395 2.07075159 -1.41928299 0.49401913 2.38142842
## [121] -1.45213461 1.07955320 -1.71368807 0.18782603 0.56374486 -1.91803230
## [127] -1.62509538 -1.56879838 -0.64295200 -1.49468179 -0.31519969 -1.13534507
## [133] -1.74327811 0.32904720 -1.39810145 -0.01768785 0.34163371 3.60788110
## [139] -1.39293867 -1.49776051 0.38946773 -1.17653535 -0.93642337 -1.80033391
## [145] -1.58687942 -1.89972780 -1.55479687 -1.72154185 -1.80798342 1.83938575
## [151] -1.60127003 -1.67668127 -1.36534893 -1.42696468 -1.72711162 -1.79056572
## [157] -1.84651714 -1.52824584 -1.81267453 -1.32554606 0.83032521 -1.68589230
## [163] -1.81967248 0.65352159 -0.99067133 -0.78718837 -1.45650506 -1.15626122
## [169] -1.59216883 -1.10657001 -1.41488972 -1.49998913 -2.06222955 -1.59190874
## [175] -1.44162416 -1.07479226
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = COVID19RATE$RATE32, listw = knb_wm)
## attr(,"class")
## [1] "localG"Then, we join the Local Gi statistics vector to our COVID19RATE SpatialPolygonsDataFrame.
COVID19RATE.localGi <- cbind(COVID19RATE, as.matrix(localGi.adaptive.13), as.matrix(localGi.adaptive.32))
head(COVID19RATE.localGi)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 cumul13 cumul14
## 1 09002 09 002 Azcapotzalco 414711 408441 19 46
## 2 09003 09 003 Coyoacán 620416 621952 25 44
## 3 09004 09 004 Cuajimalpa de Morelos 186391 199809 66 75
## 4 09005 09 005 Gustavo A. Madero 1185772 1176967 28 96
## 5 09006 09 006 Iztacalco 384326 393821 11 27
## 6 09007 09 007 Iztapalapa 1815786 1815551 41 100
## cumul15 cumul16 cumul17 cumul18 cumul19 cumul20 cumul21 cumul22 cumul23
## 1 87 156 254 391 615 931 1267 1590 1939
## 2 107 213 375 554 831 1150 1492 1866 2255
## 3 93 122 157 202 260 323 414 573 714
## 4 193 397 734 1107 1653 2330 3108 3843 4674
## 5 65 188 343 552 806 1125 1430 1706 1989
## 6 225 526 1061 1754 2617 3583 4659 5532 6416
## cumul24 cumul25 cumul26 cumul27 cumul28 cumul29 cumul30 cumul31 cumul32
## 1 2249 2594 2922 3238 3539 3839 4171 4442 4466
## 2 2554 2936 3208 3590 3903 4393 4820 5216 5246
## 3 851 991 1117 1285 1444 1615 1835 1949 1956
## 4 5366 5977 6502 7075 7587 8212 8875 9380 9439
## 5 2209 2454 2660 2825 2997 3183 3416 3612 3674
## 6 7234 7981 8646 9408 10078 10822 11550 12170 12305
## RATE13 RATE14 RATE15 RATE16 RATE17 RATE18 RATE19 RATE20 RATE21 RATE22 RATE23
## 1 0.47 1.13 2.13 3.82 6.22 9.57 15.06 22.79 31.02 38.93 47.47
## 2 0.40 0.71 1.72 3.42 6.03 8.91 13.36 18.49 23.99 30.00 36.26
## 3 3.30 3.75 4.65 6.11 7.86 10.11 13.01 16.17 20.72 28.68 35.73
## 4 0.24 0.82 1.64 3.37 6.24 9.41 14.04 19.80 26.41 32.65 39.71
## 5 0.28 0.69 1.65 4.77 8.71 14.02 20.47 28.57 36.31 43.32 50.51
## 6 0.23 0.55 1.24 2.90 5.84 9.66 14.41 19.74 25.66 30.47 35.34
## RATE24 RATE25 RATE26 RATE27 RATE28 RATE29 RATE30 RATE31 RATE32
## 1 55.06 63.51 71.54 79.28 86.65 93.99 102.12 108.75 109.34
## 2 41.06 47.21 51.58 57.72 62.75 70.63 77.50 83.86 84.35
## 3 42.59 49.60 55.90 64.31 72.27 80.83 91.84 97.54 97.89
## 4 45.59 50.78 55.24 60.11 64.46 69.77 75.41 79.70 80.20
## 5 56.09 62.31 67.54 71.73 76.10 80.82 86.74 91.72 93.29
## 6 39.84 43.96 47.62 51.82 55.51 59.61 63.62 67.03 67.78
## structure.c.3.79281262858071..4.12737686355648..6.5223927677069..
## 1 3.7928126
## 2 4.1273769
## 3 6.5223928
## 4 3.3852665
## 5 3.7010410
## 6 0.7631226
## structure.c.3.6091303847798..5.86441739551681..3.79509781846556..
## 1 3.609130
## 2 5.864417
## 3 3.795098
## 4 4.372798
## 5 4.759532
## 6 5.225795We can see that the Z-score’s column name is not ideal so we can change the column name.
names(COVID19RATE.localGi)[47] <- "gstat13"
names(COVID19RATE.localGi)[48] <- "gstat32"
COVID19RATE.localGi <- COVID19RATE.localGi %>%
arrange(desc(gstat32))
head(COVID19RATE.localGi)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 cumul13 cumul14
## 12 09013 09 013 Xochimilco 415007 418060 11 32
## 11 09012 09 012 Tlalpan 650567 682234 41 70
## 13 09014 09 014 Benito Juárez 385439 433708 23 59
## 2 09003 09 003 Coyoacán 620416 621952 25 44
## 14 09015 09 015 Cuauhtémoc 531831 548606 33 62
## 15 09016 09 016 Miguel Hidalgo 372889 379624 112 143
## cumul15 cumul16 cumul17 cumul18 cumul19 cumul20 cumul21 cumul22 cumul23
## 12 72 137 288 440 688 1018 1361 1725 2152
## 11 153 330 537 736 1045 1357 1668 2070 2437
## 13 102 178 278 406 551 705 902 1129 1301
## 2 107 213 375 554 831 1150 1492 1866 2255
## 14 136 233 367 546 787 1090 1436 1744 2093
## 15 175 238 342 453 592 778 996 1199 1383
## cumul24 cumul25 cumul26 cumul27 cumul28 cumul29 cumul30 cumul31 cumul32
## 12 2588 2911 3252 3540 3880 4321 4698 5034 5086
## 11 2871 3296 3794 4303 4740 5241 5764 6210 6275
## 13 1481 1678 1802 1977 2112 2278 2432 2617 2638
## 2 2554 2936 3208 3590 3903 4393 4820 5216 5246
## 14 2403 2707 2989 3254 3512 3829 4148 4353 4387
## 15 1617 1865 2048 2246 2427 2711 2989 3183 3214
## RATE13 RATE14 RATE15 RATE16 RATE17 RATE18 RATE19 RATE20 RATE21 RATE22 RATE23
## 12 0.26 0.77 1.72 3.28 6.89 10.52 16.46 24.35 32.56 41.26 51.48
## 11 0.60 1.03 2.24 4.84 7.87 10.79 15.32 19.89 24.45 30.34 35.72
## 13 0.53 1.36 2.35 4.10 6.41 9.36 12.70 16.26 20.80 26.03 30.00
## 2 0.40 0.71 1.72 3.42 6.03 8.91 13.36 18.49 23.99 30.00 36.26
## 14 0.60 1.13 2.48 4.25 6.69 9.95 14.35 19.87 26.18 31.79 38.15
## 15 2.95 3.77 4.61 6.27 9.01 11.93 15.59 20.49 26.24 31.58 36.43
## RATE24 RATE25 RATE26 RATE27 RATE28 RATE29 RATE30 RATE31 RATE32 gstat13
## 12 61.90 69.63 77.79 84.68 92.81 103.36 112.38 120.41 121.66 0.9124837
## 11 42.08 48.31 55.61 63.07 69.48 76.82 84.49 91.02 91.98 4.1455178
## 13 34.15 38.69 41.55 45.58 48.70 52.52 56.07 60.34 60.82 4.2083720
## 2 41.06 47.21 51.58 57.72 62.75 70.63 77.50 83.86 84.35 4.1273769
## 14 43.80 49.34 54.48 59.31 64.02 69.80 75.61 79.35 79.97 3.7172515
## 15 42.59 49.13 53.95 59.16 63.93 71.41 78.74 83.85 84.66 2.5606405
## gstat32
## 12 6.256913
## 11 6.161515
## 13 5.954610
## 2 5.864417
## 14 5.766533
## 15 5.3019585.4 Mapping Local Gi statistic
localGI.map.13 <- tm_shape(COVID19RATE.localGi) +
tm_fill(col = "gstat13",
style = "pretty",
palette = "-RdBu",
title = "Local Gi statistics (e-week 13)") +
tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE)
localGI.map.32 <- tm_shape(COVID19RATE.localGi) +
tm_fill(col = "gstat32",
style = "pretty",
palette = "-RdBu",
title = "Local Gi statistics (e-week 32)") +
tm_borders(alpha = 0.5) +
tm_layout(frame = FALSE)
tmap_arrange(localGI.map.13, localGI.map.32, asp = 1, ncol = 2)## Variable(s) "gstat13" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Variable(s) "gstat32" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
In the choropleth map above, we can see that the municipalities shaded from light pink to red are the hot spot areas and the municipalities shaded in blue are the cold spot areas. The darkness of the colours represents the intensity of the Gi values.
spatio-temporal pattern of COVID-19 infection rate by Local Getis-Ord Gi:
We can see that in e-week 13 map, there are two other hot spot areas/ small clusters at the edge of the study area that have the potential to become big clusters
We can see that the hotspot areas are the municipalities in Mexico City where the COVID-19 infection starts and these municipalities with local Gi statistics over 6 are Xochimilco & Tlalpan in e-week 32 map
We can also see that the distribution of hot spot areas to cold spot areas are balanced, the areas outside of hot spot areas are all cold spot areas with local Gi statistics between -2 to 0 in e-week 32 map
There is 1 outlier municipalities in the southeast region of the study area and this municipality is Zacualpan de Amilpas in Morelos State