Take-Home Exercise 2
Choy Jing Wei
10/11/2020
1 Introduction
Background
Since late December 2019, an outbreak of a novel coronavirus disease (COVID-19; previously known as 2019-nCoV) was reported in Wuhan, China, which has subsequently affected 210 countries worldwide. In general, COVID-19 is an acute resolved disease but it can also be deadly, with a 2% case fatality rate.
Country of Interest
The country of interest is Mexico, we will be focusing on the Central Mexico Area: (i.e. Mexico City (9), Mexico State (15) and Morelos State (17)).
as of September 26 2020, there had been 726,431 confirmed cases of COVID-19 in Mexico and 76,243 reported deaths.It is also estimated in mid July 2020 that there were more than 2,875,734 cases in Mexico.
(Extracted from wiki: https://en.wikipedia.org/wiki/COVID-19_pandemic_in_Mexico)
1.1 Objectives of Study
To analyse the spatio-temporal patterns of COVID-19 case at Central Mexico (i.e. Mexico City (9), Mexico State (15) and Morelos State (17)) by using localized spatial statistics methods.
1.2 Tasks
Using appropriate geospatial data wrangling methods:
o Extract municipalities located with the study area, and
o Calculate the COVID-19 rate (i.e. cases per 10000 population) from e-week 13 until e-week 32.
Show the spatio-temporal distribution of COVID-19 rates at municipality level by using appropriate thematic mapping technique and describe the spatio-temporal patterns reveals.
Perform local Moran’s I analysis and display the results by using appropriate thematic mapping techniques. Describe the spatio-temporal patterns reveal by the maps.
Perform local Getis-Ord Gi* analysis and display the results using appropriate thematic mapping techniques. Describe the spatio-temporal patterns reveal by the maps.
1.3 About the data
For the purpose of this study, a GIS data called municipalities_COVID is given.
It is in ESRI shapefile format. The shapefile consists of reported COVID-19 cases and other related data at the municipality level. The data was extracted from Secretary of Health (Spanish: Secretaría de Salud), Government of Mexico homepage’s (http://datosabiertos.salud.gob.mx/gobmx/salud/datos_abiertos/).
2 Getting Started
2.1 Installing R Packages
These are the packages installed in this assignment.
packages = c('rgdal', 'sf', 'spdep', 'tmap', 'tidyverse')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}## Loading required package: rgdal## Loading required package: sp## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/jingwei/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/jingwei/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.## Loading required package: sf## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1## Loading required package: spdep## 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')`## Loading required package: tmap## Loading required package: tidyverse## -- Attaching packages --------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.2
## v tidyr 1.1.1 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ------------------------------------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()3 Geospatial Data Wrangling
Import municipalities_COVID shapefile into R.
municipalities_COVID <- st_read(dsn = "data/geospatial", layer = "municipalities_COVID")## Reading layer `municipalities_COVID' from data source `C:\IS415\Assignment2\IS415_Assignment_2\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 2465 features and 198 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 911292 ymin: 319149.1 xmax: 4082997 ymax: 2349615
## projected CRS: MEXICO_ITRF_2008_LCC3.1 General overview
Look at the general details of the dataframe.
head(municipalities_COVID)## Simple feature collection with 6 features and 198 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2410092 ymin: 1067540 xmax: 2514978 ymax: 1159778
## projected CRS: MEXICO_ITRF_2008_LCC
## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 new1 new2 new3
## 1 01001 01 001 Aguascalientes 797010 961977 0 0 0
## 2 01002 01 002 Asientos 45492 50864 0 0 0
## 3 01003 01 003 Calvillo 54136 60760 0 0 0
## 4 01004 01 004 Cosío 15042 16918 0 0 0
## 5 01005 01 005 Jesús María 99590 130184 0 0 0
## 6 01006 01 006 Pabellón de Arteaga 41862 50032 0 0 0
## new4 new5 new6 new7 new8 new9 new10 new11 new12 new13 new14 new15 new16 new17
## 1 0 0 0 0 0 0 0 1 7 30 8 10 51 60
## 2 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 8
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 1 0 0 6 3
## 6 0 0 0 0 0 0 0 0 0 4 1 0 0 0
## new18 new19 new20 new21 new22 new23 new24 new25 new26 new27 new28 new29 new30
## 1 87 101 101 99 185 257 312 280 258 234 295 307 267
## 2 1 0 0 5 2 1 8 6 4 9 11 6 7
## 3 1 1 1 0 0 3 3 8 11 8 20 35 13
## 4 0 0 6 7 3 7 6 12 1 0 2 10 9
## 5 3 4 5 5 6 5 11 10 15 14 15 13 10
## 6 6 3 2 14 6 17 8 19 13 7 8 21 14
## new31 new32 cumul1 cumul2 cumul3 cumul4 cumul5 cumul6 cumul7 cumul8 cumul9
## 1 219 23 0 0 0 0 0 0 0 0 0
## 2 11 0 0 0 0 0 0 0 0 0 0
## 3 13 3 0 0 0 0 0 0 0 0 0
## 4 4 0 0 0 0 0 0 0 0 0 0
## 5 15 2 0 0 0 0 0 0 0 0 0
## 6 19 0 0 0 0 0 0 0 0 0 0
## cumul10 cumul11 cumul12 cumul13 cumul14 cumul15 cumul16 cumul17 cumul18
## 1 0 1 8 38 46 56 107 167 254
## 2 0 0 1 1 1 1 1 1 2
## 3 0 0 0 0 0 0 0 8 9
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 1 1 1 7 10 13
## 6 0 0 0 4 5 5 5 5 11
## cumul19 cumul20 cumul21 cumul22 cumul23 cumul24 cumul25 cumul26 cumul27
## 1 355 456 555 740 997 1309 1589 1847 2081
## 2 2 2 7 9 10 18 24 28 37
## 3 10 11 11 11 14 17 25 36 44
## 4 0 6 13 16 23 29 41 42 42
## 5 17 22 27 33 38 49 59 74 88
## 6 14 16 30 36 53 61 80 93 100
## cumul28 cumul29 cumul30 cumul31 cumul32 activ1 activ2 activ3 activ4 activ5
## 1 2376 2683 2950 3169 3192 0 0 0 0 0
## 2 48 54 61 72 72 0 0 0 0 0
## 3 64 99 112 125 128 0 0 0 0 0
## 4 44 54 63 67 67 0 0 0 0 0
## 5 103 116 126 141 143 0 0 0 0 0
## 6 108 129 143 162 162 0 0 0 0 0
## activ6 activ7 activ8 activ9 activ10 activ11 activ12 activ13 activ14 activ15
## 1 0 0 0 0 1 4 16 45 53 75
## 2 0 0 0 0 0 0 1 1 1 1
## 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 1 1 3
## 6 0 0 0 0 0 0 0 5 5 5
## activ16 activ17 activ18 activ19 activ20 activ21 activ22 activ23 activ24
## 1 139 203 309 393 500 622 848 1133 1424
## 2 1 1 2 2 4 9 10 14 19
## 3 0 9 9 11 11 11 12 15 18
## 4 0 0 0 1 7 14 18 27 34
## 5 7 12 14 19 23 28 34 44 53
## 6 5 5 13 14 22 34 41 60 70
## activ25 activ26 activ27 activ28 activ29 activ30 activ31 activ32 death1 death2
## 1 1713 1972 2248 2544 2845 3079 3189 3192 0 0
## 2 27 31 41 51 56 70 72 72 0 0
## 3 28 39 53 90 108 124 127 128 0 0
## 4 42 42 43 50 56 66 67 67 0 0
## 5 69 83 96 109 122 135 143 143 0 0
## 6 88 95 104 119 135 148 162 162 0 0
## death3 death4 death5 death6 death7 death8 death9 death10 death11 death12
## 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
## death13 death14 death15 death16 death17 death18 death19 death20 death21
## 1 0 0 1 2 2 3 7 6 4
## 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 2
## 6 0 0 0 0 0 0 0 0 1
## death22 death23 death24 death25 death26 death27 death28 death29 death30
## 1 17 18 27 25 24 19 19 19 20
## 2 0 0 0 0 2 0 2 2 0
## 3 0 0 0 0 1 0 0 0 1
## 4 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 1 3
## 6 1 0 1 1 1 0 1 0 1
## death31 death32 actvrt1 actvrt2 actvrt3 actvrt4 actvrt5 actvrt6 actvrt7
## 1 13 7 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
## actvrt8 actvrt9 actvr10 actvr11 actvr12 actvr13 actvr14 actvr15
## 1 0 0 0.1039526 0.4158104 1.663241 4.5739139 5.0936769 6.133203
## 2 0 0 0.0000000 0.0000000 1.966027 1.9660271 1.9660271 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.0000000 0.0000000 0.000000
## 5 0 0 0.0000000 0.0000000 0.000000 0.7681436 0.7681436 2.304431
## 6 0 0 0.0000000 0.0000000 0.000000 9.9936041 9.9936041 9.993604
## actvr16 actvr17 actvr18 actvr19 actvr20 actvr21 actvr22
## 1 9.771543 15.592888 24.324906 26.403958 30.873919 32.537160 47.29843
## 2 0.000000 0.000000 1.966027 1.966027 5.898081 13.762189 15.72822
## 3 0.000000 14.812377 14.812377 18.104016 3.291639 3.291639 1.64582
## 4 0.000000 0.000000 0.000000 5.910864 41.376049 82.752098 100.48469
## 5 4.608861 8.449579 8.449579 9.217723 8.449579 10.754010 11.52215
## 6 0.000000 0.000000 15.989767 17.988487 33.978254 41.973137 53.96546
## actvr23 actvr24 actvr25 actvr26 actvr27 actvr28 actvr29
## 1 65.801989 83.36998 89.91899 87.21622 85.65693 86.38460 90.75061
## 2 19.660271 19.66027 33.42246 33.42246 43.25260 47.18465 49.15068
## 3 6.583278 11.52074 26.33311 39.49967 57.60369 102.04082 113.56155
## 4 118.217283 118.21728 141.86074 88.66296 53.19778 47.28691 82.75210
## 5 16.131015 19.20359 26.88502 29.95760 33.03017 30.72574 29.95760
## 6 75.951391 71.95395 93.93988 69.95523 67.95651 61.96035 79.94883
## actvr30 actvr31 actvr32 dethrt1 dethrt2 dethrt3 dethrt4 dethrt5 dethrt6
## 1 86.38460 67.04942 36.07155 0 0 0 0 0 0
## 2 57.01478 41.28657 31.45643 0 0 0 0 0 0
## 3 116.85319 60.89533 32.91639 0 0 0 0 0 0
## 4 135.94988 100.48469 65.01951 0 0 0 0 0 0
## 5 29.95760 26.11688 16.13101 0 0 0 0 0 0
## 6 87.94372 85.94500 53.96546 0 0 0 0 0 0
## dethrt7 dethrt8 dethrt9 dthrt10 dthrt11 dthrt12 dthrt13 dthrt14 dthrt15
## 1 0 0 0 0 0 0 0 0 0.1039526
## 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
## dthrt16 dthrt17 dthrt18 dthrt19 dthrt20 dthrt21 dthrt22 dthrt23
## 1 0.2079052 0.2079052 0.3118578 0.7276681 0.6237155 0.4158104 1.767194 1.871147
## 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 0.0000000 0.000000 0.000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.5362871 0.000000 0.000000
## 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.9987208 1.998721 0.000000
## dthrt24 dthrt25 dthrt26 dthrt27 dthrt28 dthrt29 dthrt30 dthrt31
## 1 2.806720 2.598815 2.494862 1.975099 1.975099 1.9750992 2.079052 1.351384
## 2 0.000000 0.000000 3.932054 0.000000 3.932054 3.9320541 0.000000 0.000000
## 3 0.000000 0.000000 1.645820 0.000000 0.000000 0.0000000 1.645820 0.000000
## 4 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000 0.000000
## 5 0.000000 0.000000 0.000000 0.000000 0.000000 0.7681436 2.304431 0.000000
## 6 1.998721 1.998721 1.998721 0.000000 1.998721 0.0000000 1.998721 0.000000
## dthrt32 geometry
## 1 0.7276681 MULTIPOLYGON (((2489073 111...
## 2 0.0000000 MULTIPOLYGON (((2494680 114...
## 3 0.0000000 MULTIPOLYGON (((2429607 112...
## 4 0.0000000 MULTIPOLYGON (((2470518 115...
## 5 0.0000000 MULTIPOLYGON (((2465527 111...
## 6 0.0000000 MULTIPOLYGON (((2473997 112...3.2 Check if there is any NA values
Check if there is any NA values in the dataframe.
sum(is.na(municipalities_COVID))## [1] 03.3 Filter out the states and city
Filter out the states and city that we are interested in: Mexico City, Mexico State and Morelos State.
municipalities_COVID_filtered_by_states_and_city_only <- municipalities_COVID %>%
filter(CVE_ENT %in% c("09", "15","17"))3.4 Calculate the Covid Rate from week 13 until week 32
Calculate the total cumulative numbers for week 1 to 12 and week 13 to 32
municipalities_COVID_filtered_by_states_and_city_only <- municipalities_COVID_filtered_by_states_and_city_only %>%
mutate(`cumulative_w1_to_w12` = `new1`+`new2`+`new3`+`new4`+`new5`+`new6`+`new7`+`new8`+`new9`+`new10`+`new11`+`new12`) %>%
mutate(`cumulative_w13_to_w32` = `new13`+`new14`+`new15`+`new16`+`new17`+`new18`+`new19`+`new20`+`new21`+`new22`+`new23`+`new24`+
`new25`+`new26`+`new27`+`new28`+`new29`+`new30`+`new31`+`new32`)Calculate the COVID-19 rate (i.e. cases per 10000 population) from e-week 13 until e-week 32
We will be using the formula: (accumulated cases / local population) x 10,000, to calculate the Covid cases per 10,000 population.
municipalities_COVID_rates_w13_to_w32 <- municipalities_COVID_filtered_by_states_and_city_only %>%
dplyr::select(CVEGEO,CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13, cumul14, cumul15, cumul16, cumul17, cumul18, cumul19, cumul20, cumul21, cumul22,cumul23, cumul24, cumul25, cumul26, cumul27, cumul28, cumul29, cumul30, cumul31, cumul32) %>%
mutate(across(starts_with("cumul"), list(rate = ~./ `Pop2020` * 10000))) %>%
dplyr::select(CVEGEO,CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate,cumul23_rate, cumul24_rate, cumul25_rate, cumul26_rate, cumul27_rate, cumul28_rate, cumul29_rate, cumul30_rate, cumul31_rate, cumul32_rate)
municipalities_COVID_rates_w13_to_w32## Simple feature collection with 176 features and 25 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2645804 ymin: 707708.1 xmax: 2855437 ymax: 921773.9
## projected CRS: MEXICO_ITRF_2008_LCC
## First 10 features:
## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2020 cumul13_rate
## 1 09002 09 002 Azcapotzalco 408441 0.46518347
## 2 09003 09 003 Coyoacán 621952 0.40196028
## 3 09004 09 004 Cuajimalpa de Morelos 199809 3.30315451
## 4 09005 09 005 Gustavo A. Madero 1176967 0.23789962
## 5 09006 09 006 Iztacalco 393821 0.27931471
## 6 09007 09 007 Iztapalapa 1815551 0.22582676
## 7 09008 09 008 La Magdalena Contreras 245147 0.61187777
## 8 09009 09 009 Milpa Alta 139371 0.00000000
## 9 09010 09 010 Álvaro Obregón 755537 0.74119467
## 10 09011 09 011 Tláhuac 366586 0.08183619
## cumul14_rate cumul15_rate cumul16_rate cumul17_rate cumul18_rate
## 1 1.12623365 2.1300506 3.819401 6.218768 9.572986
## 2 0.70745009 1.7203900 3.424702 6.029404 8.907440
## 3 3.75358467 4.6544450 6.105831 7.857504 10.109655
## 4 0.81565583 1.6398081 3.373077 6.236369 9.405531
## 5 0.68559066 1.6504960 4.773742 8.709541 14.016520
## 6 0.55079698 1.2392932 2.897192 5.843956 9.660979
## 7 1.10137999 1.8356333 3.385724 6.893823 8.974207
## 8 0.07175094 1.6502716 4.233305 9.829879 14.852444
## 9 1.19120573 1.9324004 3.573617 5.453075 8.060492
## 10 0.30006601 0.8183619 2.155020 4.773778 8.456406
## cumul19_rate cumul20_rate cumul21_rate cumul22_rate cumul23_rate
## 1 15.05725 22.79399 31.02039 38.92851 47.47320
## 2 13.36116 18.49017 23.98899 30.00232 36.25682
## 3 13.01243 16.16544 20.71979 28.67739 35.73413
## 4 14.04457 19.79665 26.40686 32.65172 39.71224
## 5 20.46615 28.56628 36.31091 43.31917 50.50518
## 6 14.41436 19.73506 25.66163 30.47009 35.33913
## 7 11.99280 17.41812 23.29215 29.08459 34.22436
## 8 21.59703 30.13539 44.12683 55.82223 68.95265
## 9 11.21057 15.47244 20.56815 25.69034 30.75958
## 10 13.99399 22.25944 31.86155 38.21750 46.48295
## cumul24_rate cumul25_rate cumul26_rate cumul27_rate cumul28_rate
## 1 55.06303 63.50978 71.54032 79.27706 86.64654
## 2 41.06426 47.20622 51.57954 57.72150 62.75404
## 3 42.59067 49.59737 55.90339 64.31142 72.26902
## 4 45.59176 50.78307 55.24369 60.11214 64.46230
## 5 56.09147 62.31257 67.54338 71.73310 76.10056
## 6 39.84465 43.95911 47.62191 51.81898 55.50932
## 7 39.20097 46.58429 53.88604 62.37074 71.71207
## 8 84.95311 97.15077 106.98065 119.60881 130.44321
## 9 36.82149 42.57899 47.43646 52.57188 58.11760
## 10 52.94801 60.55878 68.11499 76.16221 84.67317
## cumul29_rate cumul30_rate cumul31_rate cumul32_rate
## 1 93.99154 102.12001 108.75500 109.34260
## 2 70.63246 77.49794 83.86499 84.34735
## 3 80.82719 91.83771 97.54315 97.89349
## 4 69.77256 75.40568 79.69637 80.19766
## 5 80.82352 86.73991 91.71680 93.29111
## 6 59.60725 63.61705 67.03199 67.77557
## 7 85.78526 95.86085 107.16019 107.64970
## 8 145.08040 159.21533 169.61922 170.40848
## 9 65.58249 72.09442 77.71956 78.54016
## 10 94.24801 103.14087 110.91531 111.95190
## geometry
## 1 MULTIPOLYGON (((2794860 837...
## 2 MULTIPOLYGON (((2799635 820...
## 3 MULTIPOLYGON (((2787230 825...
## 4 MULTIPOLYGON (((2802176 843...
## 5 MULTIPOLYGON (((2808291 828...
## 6 MULTIPOLYGON (((2812453 823...
## 7 MULTIPOLYGON (((2792518 818...
## 8 MULTIPOLYGON (((2814877 806...
## 9 MULTIPOLYGON (((2794396 824...
## 10 MULTIPOLYGON (((2816579 817...4 Show the spatio-temporal distribution of COVID-19 rates at municipality level
Since there are 20 weeks and the Covid rates are increasing exponentially.
If we use only 1 scale which that ranges very widely for all of the 20 weeks, the maps at the beginning will not show much changes.
Thus, we will be using a few sets of scales so as to better allow viewers to understand how has the Covid rates change from 1 week to another, especially in the earlier weeks.
For Week 13 to 16 we will be using range between < 1 to 10 <
For Week 17 to 20 we will be using range between < 5 to 20 <
For Week 21 to 24 we will be using range between < 10 to 40 <
For Week 25 to 28 we will be using range between < 15 to 60 <
For Week 29 to 32 we will be using range between < 20 to 80 <
Additionally at the end of this section, we will be comparing the change from week 12 and week 3.
4.1 Week 13 to 16 spatio-temporal distribution of COVID-19 rates at municipality level
covid_map_wk13_14 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul13_rate", "cumul14_rate"),
breaks=c(-Inf, 1, 4, 7, 10, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 13", "Covid Rate for Week 14"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk15_16 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul15_rate", "cumul16_rate"),
breaks=c(-Inf, 1, 4, 7, 10, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 15", "Covid Rate for Week 16"),
legend.title.size = 1,
legend.text.size = 0.6)covid_map_wk13_14
covid_map_wk15_16
Based on the maps shown for the 4 weeks (13,14,15,16), it is obvious that there are 3 clusters where covid cases are between 1 to 4 per 10,000 population.
From week 13 to 14, the cluster at the north east region began to spread to its neighbouring areas at quite a quick rate compared to the initial 3 clusters. The areas that now has 1 to 4 cases per 10,000 is now increased from 3 to 9. However, the other 2 initial cluster does not seem to be spreading as quickly.
From week 15 to 16, we can now tell that only 1 out of the initial 3 clusters we suspected to be the source, the north east cluster’s situation has deteriorated and spread further north east, including some surrounding areas. Certain areas that initially has 1 to 4 cases per 10,000 has increased to 4 to 7 cases per 10,000.
4.2 Week 17 to 20 spatio-temporal distribution of COVID-19 rates at municipality level
covid_map_wk17_18 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul17_rate", "cumul18_rate"),
breaks=c(-Inf, 5, 10, 15, 20, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 17", "Covid Rate for Week 18"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk19_20 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul19_rate", "cumul20_rate"),
breaks=c(-Inf, 5, 10, 15, 20, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 19", "Covid Rate for Week 20"),
legend.title.size = 1,
legend.text.size = 0.6)covid_map_wk17_18
covid_map_wk19_20
Based on the maps shown for the 4 weeks (17,18,19,20), we can see that the situation at north east region still stands out even after the scale has been adjusted.
From week 17 to 18, we can see that the north east cluster is beginning to spread to neighbouring regions, the cases per 10,000 is now reaching 10 to 15 for the initial regions.
From week 18 to 20, the situation is now more obvious as it can be seen visually that there is a significant jump from week 18 to 19, many other regions are now affected and the number of cases has now intensified, some regions now has 20 or more cases.
From week 19 to 20, almost 6 regions have reached 20 or more cases per 10,000 and a region in the south has also reached 15 to 20 cases per 10,000 which could be the next potential cluster.
4.3 Week 21 to 24 spatio-temporal distribution of COVID-19 rates at municipality level
covid_map_wk21_22 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul21_rate", "cumul22_rate"),
breaks=c(-Inf, 10, 20, 30, 40, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 21", "Covid Rate for Week 22"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk23_24 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul23_rate", "cumul24_rate"),
breaks=c(-Inf, 10, 20, 30, 40, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 23", "Covid Rate for Week 24"),
legend.title.size = 1,
legend.text.size = 0.6)covid_map_wk21_22
covid_map_wk23_24
Based on the maps shown for the 4 weeks (21,22,23,24), the north east cluster still remains at the upper bound of the scale, reaching 40 or more cases.
From week 21 to 22, we can see that most of the cases are still at the north east and east side of Mexico.
From week 22 to 24, we can see that the cases at the south is gradually increasing from less than 10 cases to 20 cases in many regions. We can also tell that the Covid situation is beginning to spread westwards. More regions in the west are now affected with 10 to 20 cases.
4.4 Week 25 to 28 spatio-temporal distribution of COVID-19 rates at municipality level
covid_map_wk25_26 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul25_rate", "cumul26_rate"),
breaks=c(-Inf, 15, 25, 45, 60, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 25", "Covid Rate for Week 26"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk27_28 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul27_rate", "cumul28_rate"),
breaks=c(-Inf, 15, 25, 45, 60, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 27", "Covid Rate for Week 28"),
legend.title.size = 1,
legend.text.size = 0.6)covid_map_wk25_26
covid_map_wk27_28
Based on the maps shown for the 4 weeks (25,26,27,28), the north east region has reached 60 or more cases per 10,000. It can be concluded that most regions are now affected, many have up to 15 to 25 cases per 10,000.
By week 28, parts that have previously had less than 15 cases per 10,000 are now between 15 to 25 cases per 10,000.
4.5 Week 29 to 32 spatio-temporal distribution of COVID-19 rates at municipality level
covid_map_wk29_30 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul29_rate", "cumul30_rate"),
breaks=c(-Inf, 20, 40, 60, 80, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 29", "Covid Rate for Week 30"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk31_32 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul31_rate", "cumul32_rate"),
breaks=c(-Inf, 20, 40, 60, 80, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 31", "Covid Rate for Week 32"),
legend.title.size = 1,
legend.text.size = 0.6)covid_map_wk29_30
covid_map_wk31_32
Based on the maps shown for the 4 weeks (29,30,31,32), at week 29, we can see that many regions has reached 20 to 40 cases per 10,000.
Also, the situation at the north east region is still serious as they yet again hit the upper bound of the scale, 80 or more cases per 10,000.
4.6 Week 13 to 32 spatio-temporal distribution of COVID-19 rates at municipality level
We will now compare the changes from week 13, 20, 27, 32 using the same scale to see the overall changes.
We will be using the scale: < 15 to 120 <
covid_map_wk13_20 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul13_rate", "cumul20_rate"),
breaks=c(-Inf, 15, 30, 60, 120, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 13", "Covid Rate for Week 20"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk27_32 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul27_rate", "cumul32_rate"),
breaks=c(-Inf, 15, 30, 60, 120, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 27", "Covid Rate for Week 32"),
legend.title.size = 1,
legend.text.size = 0.6)covid_map_wk13_20
covid_map_wk27_32
In conclusion, we can see that there is significant difference visually from week 13 to 32. Within the 20 weeks, the north eastern part has been hit the hardest by Covid cases, from less than 15 cases per 10,000 to even 120 or more at certain areas. We can also see that Covid usually spread to the neighbouring regions first and areas that have been hit with Covid cases are very likely to have a greater rate per 10,000 by week 32.
We can also clearly see that the covid rate has been increasing exponentially. Between week 13 and week 20, it can be seen that a cluster of cases has been appearing at the north east region within a week. However, given another week, at week 27, the number of regions that has 15 to 30 cases per 10,000 has increased drastically and more areas have seen that increase.
Between week 27 and week 32, some regions that had 15 to 30 cases previously now has reached 30 to 60 cases. The original north east cluster now has 120 or more cases per 10,000.
5 Using Queen Contiguity Neighbours
Based on the graphs we previously plotted, there is some evidence that the covid cases usually spread to neighbouring regions. Thus, we can conduct queen contiguity based neighbours to identify the regions linked to high number of neighbours.
5.1 Converting to SpatialPolygonsDataFrame
spdf_municipalities_COVID_rates_w13_to_w32 <- as_Spatial(municipalities_COVID_rates_w13_to_w32)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown based on GRS80 ellipsoid in CRS definition(spdf_municipalities_COVID_rates_w13_to_w32)## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 25
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 3.30315451255949, 3.76688512844288, 4.6544449949702, 6.26936126272312, 9.82987852566172, 14.8524441957079, 21.5970323811984, 30.1353940202768, 44.1268269582625, 55.8222298756556, ...5.2 Compute Queen contiguity weight matrix
spdf_municipalities_COVID_rates_w13_to_w32_queen <- poly2nb(spdf_municipalities_COVID_rates_w13_to_w32, queen=TRUE)
summary(spdf_municipalities_COVID_rates_w13_to_w32_queen)## 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:
## 77 95 121 with 1 link
## 1 most connected region:
## 117 with 14 linksFrom the data generated, we can see that Municipality Code 117 has the most links (14) and many other areas have 4 to 6 links based on the frequency table. The average number of links is also 5.465909 which is about 5.
5.3 Display the complete weight matrix
str(spdf_municipalities_COVID_rates_w13_to_w32_queen)## List of 176
## $ : int [1:5] 4 14 15 73 120
## $ : int [1:5] 6 9 11 12 13
## $ : int [1:4] 9 15 53 78
## $ : int [1:8] 1 14 16 36 49 74 120 125
## $ : int [1:5] 6 13 14 16 74
## $ : int [1:8] 2 5 10 12 13 74 86 138
## $ : int [1:3] 9 11 78
## $ : int [1:8] 10 11 12 41 66 150 161 164
## $ : int [1:7] 2 3 7 11 13 15 78
## $ : int [1:5] 6 8 12 41 138
## $ : int [1:9] 2 7 8 9 12 59 78 117 150
## $ : int [1:5] 2 6 8 10 11
## $ : int [1:6] 2 5 6 9 14 15
## $ : int [1:6] 1 4 5 13 15 16
## $ : int [1:7] 1 3 9 13 14 53 73
## $ : int [1:4] 4 5 14 74
## $ : int [1:4] 19 30 101 118
## $ : int [1:7] 27 44 49 97 108 109 116
## $ : int [1:4] 17 61 87 118
## $ : int [1:4] 37 96 113 133
## $ : int [1:7] 23 58 90 103 122 130 134
## $ : int [1:4] 28 89 114 117
## $ : int [1:7] 21 48 102 126 127 130 134
## $ : int [1:3] 96 98 121
## $ : int [1:4] 31 33 84 119
## $ : int [1:2] 52 112
## $ : int [1:6] 18 44 46 49 115 116
## $ : int [1:3] 22 89 117
## $ : int [1:6] 54 62 73 76 120 137
## $ : int [1:5] 17 64 72 101 118
## $ : int [1:5] 25 50 84 163 171
## $ : int [1:4] 77 81 91 100
## $ : int [1:6] 25 66 84 105 110 119
## $ : int [1:8] 43 70 71 88 89 106 117 122
## $ : int [1:3] 67 78 117
## $ : int [1:4] 4 49 124 125
## $ : int [1:8] 20 56 102 106 113 129 133 134
## $ : int [1:3] 41 99 119
## $ : int [1:4] 51 107 111 136
## $ : int [1:8] 69 75 107 111 124 125 136 137
## $ : int [1:9] 8 10 38 55 66 99 105 119 138
## $ : int [1:4] 61 72 118 128
## $ : int [1:4] 34 70 71 117
## $ : int [1:7] 18 27 46 85 109 115 116
## $ : int [1:4] 47 55 86 115
## $ : int [1:3] 27 44 115
## $ : int [1:4] 45 74 86 115
## $ : int [1:4] 23 57 126 127
## $ : int [1:11] 4 18 27 36 60 74 97 115 120 125 ...
## $ : int [1:4] 31 157 163 171
## $ : int [1:4] 39 111 112 136
## $ : int [1:3] 26 112 136
## $ : int [1:5] 3 15 67 73 78
## $ : int [1:5] 29 62 76 83 103
## $ : int [1:6] 41 45 86 115 119 138
## $ : int [1:5] 37 123 129 133 135
## $ : int [1:3] 48 94 126
## $ : int [1:5] 21 63 64 90 103
## $ : int [1:4] 11 79 117 150
## $ : int [1:6] 49 75 97 125 136 141
## $ : int [1:5] 19 42 87 95 118
## $ : int [1:5] 29 54 73 83 131
## $ : int [1:6] 58 64 72 76 103 128
## $ : int [1:7] 30 58 63 72 80 90 101
## $ : int [1:6] 68 79 104 106 114 117
## $ : int [1:9] 8 33 41 84 99 105 110 164 168
## $ : int [1:10] 35 53 70 73 78 83 92 117 122 131
## $ : int [1:6] 65 79 104 135 146 156
## $ : int [1:3] 40 75 124
## $ : int [1:7] 34 43 67 71 92 117 122
## $ : int [1:3] 34 43 70
## $ : int [1:6] 30 42 63 64 118 128
## $ : int [1:8] 1 15 29 53 62 67 120 131
## $ : int [1:8] 4 5 6 16 47 49 86 115
## $ : int [1:8] 40 60 69 97 124 125 136 141
## $ : int [1:7] 29 54 63 103 111 128 137
## $ : int 32
## $ : int [1:8] 3 7 9 11 35 53 67 117
## $ : int [1:8] 59 65 68 104 117 148 150 156
## $ : int [1:4] 64 90 101 140
## $ : int [1:4] 32 91 108 109
## $ : int [1:5] 94 102 126 132 139
## $ : int [1:6] 54 62 67 103 122 131
## $ : int [1:7] 25 31 33 66 110 143 171
## $ : int [1:3] 44 109 115
## $ : int [1:6] 6 45 47 55 74 138
## $ : int [1:2] 19 61
## $ : int [1:4] 34 89 106 114
## $ : int [1:6] 22 28 34 88 114 117
## $ : int [1:6] 21 58 64 80 130 140
## $ : int [1:5] 32 81 100 108 109
## $ : int [1:3] 67 70 122
## $ : int [1:3] 98 102 113
## $ : int [1:3] 57 82 126
## $ : int 61
## $ : int [1:5] 20 24 98 113 133
## $ : int [1:8] 18 49 60 75 100 108 136 141
## $ : int [1:7] 24 93 96 102 113 132 139
## $ : int [1:5] 38 41 66 105 119
## [list output truncated]
## - attr(*, "class")= chr "nb"
## - attr(*, "region.id")= chr [1:176] "1" "2" "3" "4" ...
## - attr(*, "call")= language poly2nb(pl = spdf_municipalities_COVID_rates_w13_to_w32, queen = TRUE)
## - attr(*, "type")= chr "queen"
## - attr(*, "sym")= logi TRUE5.4 Finding out if highest links means highest covid rate
municipalities_COVID_rates_w13_to_w32$NOMGEO[117]## [1] "Tianguistenco"municipalities_COVID_rates_w13_to_w32$cumul32[117]## [1] 24.49762municipalities_COVID_rates_w13_to_w32$Pop2020[117]## [1] 82457We can also tell that even though Tianguistenco has the highest number of links, its covid rate as of week 32 is only 24.49762 cases per 10,000 population.
5.5 Finding out more details on the area with the highest covid rate
which.max(municipalities_COVID_rates_w13_to_w32$cumul32)## [1] 8municipalities_COVID_rates_w13_to_w32$cumul32[8]## [1] 170.4085municipalities_COVID_rates_w13_to_w32$NOMGEO[8]## [1] "Milpa Alta"spdf_municipalities_COVID_rates_w13_to_w32_queen[8]## [[1]]
## [1] 10 11 12 41 66 150 161 164municipalities_COVID_rates_w13_to_w32$Pop2020[8]## [1] 139371We can tell that Milpa Alta is the municipality with the highest covid rate at week 32. It has almost 170 cases per 10,000 population which is about 1.7% and has 8 links connected to it which can be considered highly connected.
5.6 Plotting Queen contiguity based neighbours map
Plotting the Queen contiguity based neighbours map comparing it with the covid rate map for week 32 to identify visually if there is any similarity.
plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(spdf_municipalities_COVID_rates_w13_to_w32_queen, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col= "red")
covid_map_wk32 <- tm_shape(municipalities_COVID_rates_w13_to_w32) +
tm_fill(col = c("cumul32_rate"),
breaks=c(-Inf, 20, 40, 60, 80, Inf),
palette="Greens",
title = "Rate Per 10000") +
tm_borders(alpha = 0.5) +
tm_layout(c("Covid Rate for Week 32"),
legend.title.size = 1,
legend.text.size = 0.6)
covid_map_wk32
qtm(spdf_municipalities_COVID_rates_w13_to_w32, "cumul32_rate")
Based on the queen contiguity map, we can see that the areas with the most connected links are situated towards near the north east of Mexico and similarly most of the covid cases are also situated there too. So there could be a possibility that the links may have accelerated the spread of covid and the ones with the most links are potentially open to being clusters.
6 Performing Global Spatial Autocorrelation
6.1 Moran’s I
Computing row-standardised weights matrix
rswm_q <- nb2listw(spdf_municipalities_COVID_rates_w13_to_w32_queen, style="W", zero.policy = TRUE)
rswm_q## 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: W
## Weights constants summary:
## n nn S0 S1 S2
## W 176 30976 176 71.10695 731.6796.1.1 Moron’s I test for Week 13
Performs Moran’s I statistical testing using moran.test() of spdep
moran.test(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate, listw=rswm_q, zero.policy = TRUE, na.action=na.omit)##
## Moran I test under randomisation
##
## data: spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate
## weights: rswm_q
##
## Moran I statistic standard deviate = 8.3371, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.346519030 -0.005714286 0.0017849716.1.2 Computing Monte Carlo Moran’s I for Week 13
Performs permutation test for Moran’s I statistic by using moran.mc() of spdep
H0 = The distributions are randomly distributed.
H1= The distributions are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05
set.seed(1234)
bperm= moran.mc(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate, listw=rswm_q, nsim=99, zero.policy = TRUE, na.action=na.omit)
bperm##
## Monte-Carlo simulation of Moran I
##
## data: spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate
## weights: rswm_q
## number of simulations + 1: 100
##
## statistic = 0.34652, observed rank = 100, p-value = 0.01
## alternative hypothesis: greaterSince the p-value is smaller than alpha value of 0.05, we can reject the null hypothesis.
Our global Moran’I statistic is 0.34652 which is greater than 0, tells us that it is positive autocorrelation and clustered.
6.1.3 Computing Monte Carlo Moran’s I for Week 32
Performs permutation test for Moran’s I statistic by using moran.mc() of spdep
H0 = The distributions are randomly distributed.
H1= The distributions are not randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05
set.seed(1234)
bperm= moran.mc(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate, listw=rswm_q, nsim=99, zero.policy = TRUE, na.action=na.omit)
bperm##
## Monte-Carlo simulation of Moran I
##
## data: spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate
## weights: rswm_q
## number of simulations + 1: 100
##
## statistic = 0.49659, observed rank = 100, p-value = 0.01
## alternative hypothesis: greaterSince the p-value is smaller than alpha value of 0.05, we can reject the null hypothesis.
Our global Moran’I statistic is 0.49659 which is greater than 0, tells us that it is positive autocorrelation and clustered.
7 Performing Local Spatial Autocorrelation
7.1 Local Moran’s I for Week 13
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate)
localMI_w13 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate, rswm_q)
head(localMI_w13)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 1.42938461 -0.005714286 0.15624522 3.6306005 0.0001413813
## 2 0.52791359 -0.005714286 0.15624522 1.3500042 0.0885073115
## 3 21.14786143 -0.005714286 0.19615364 47.7623404 0.0000000000
## 4 0.06784766 -0.005714286 0.09638259 0.2369485 0.4063483824
## 5 0.17785066 -0.005714286 0.15624522 0.4643938 0.3211828265
## 6 0.04471045 -0.005714286 0.09638259 0.1624218 0.43548684837.2 Mapping both local Moran’s I values and p-values for week 13
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w13 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w13)
localMI.map_w13 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w13) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w13 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w13) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w13, pvalue.map_w13, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and 3 regions that have greater auto correlation than the rest.
For the local moran p value map, we can see only a few regions are statistically significant while the rest are most statistically insignificant.
7.3 Plotting Moran scatterplot for week 13 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate) %>% as.vector
nci13 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 13 Rate", ylab="Spatially Lag z-Week 13 Rate")
7.4 prepare a LISA cluster map for week 13
quadrant_13 <- vector(mode="numeric",length=nrow(localMI_w13))
DV_13 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate)
C_mI_13 <- localMI_w13[,1] - mean(localMI_w13[,1])
signif <- 0.05
quadrant_13[DV_13 >0 & C_mI_13>0] <- 4
quadrant_13[DV_13 <0 & C_mI_13<0] <- 1
quadrant_13[DV_13 <0 & C_mI_13>0] <- 2
quadrant_13[DV_13 >0 & C_mI_13<0] <- 3
quadrant_13[localMI_w13[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w13$quadrant13 <- quadrant_13
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w13) +
tm_fill(col = "quadrant13", style = "cat", palette = colors[c(sort(unique(quadrant_13)))+1], labels = clusters[c(sort(unique(quadrant_13)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant. We can conclude that that is a cluster formed with high values surrounded by high values.
7.5 Local Moran’s I for Week 14
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate)
localMI_w14 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate, rswm_q)
head(localMI_w14)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 3.6806637 -0.005714286 0.1689328 8.968977 1.496405e-19
## 2 1.2456693 -0.005714286 0.1689328 3.044623 1.164861e-03
## 3 18.8476095 -0.005714286 0.2122042 40.927101 0.000000e+00
## 4 0.8569830 -0.005714286 0.1040257 2.674780 3.738911e-03
## 5 0.9906167 -0.005714286 0.1689328 2.424078 7.673644e-03
## 6 0.3498874 -0.005714286 0.1040257 1.102538 1.351140e-017.6 Mapping both local Moran’s I values and p-values for week 14
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w14 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w14)
localMI.map_w14 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w14) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w14 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w14) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w14, pvalue.map_w14, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and 3 regions that have greater auto correlation than the rest.
For the local moran p value map, we can see only a few regions are statistically significant while the rest are most statistically insignificant.
7.7 Plotting Moran scatterplot for week 14 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate) %>% as.vector
nci14 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 14 Rate", ylab="Spatially Lag z-Week 14 Rate")
7.8 prepare a LISA cluster map for week 14
quadrant_14 <- vector(mode="numeric",length=nrow(localMI_w14))
DV_14 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate)
C_mI_14 <- localMI_w14[,1] - mean(localMI_w14[,1])
signif <- 0.05
quadrant_14[DV_14 >0 & C_mI_14>0] <- 4
quadrant_14[DV_14 <0 & C_mI_14<0] <- 1
quadrant_14[DV_14 <0 & C_mI_14>0] <- 2
quadrant_14[DV_14 >0 & C_mI_14<0] <- 3
quadrant_14[localMI_w14[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w14$quadrant14 <- quadrant_14
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w14) +
tm_fill(col = "quadrant14", style = "cat", palette = colors[c(sort(unique(quadrant_14)))+1], labels = clusters[c(sort(unique(quadrant_14)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant including a new area at south of the cluster. We can conclude that that is a cluster formed with high values surrounded by high values.
7.9 Local Moran’s I for Week 15
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate)
localMI_w15 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate, rswm_q)
head(localMI_w15)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 5.234822 -0.005714286 0.1821012 12.280595 5.757426e-35
## 2 3.301115 -0.005714286 0.1821012 7.749175 4.624574e-15
## 3 13.426999 -0.005714286 0.2288630 28.078636 8.932469e-174
## 4 1.944503 -0.005714286 0.1119585 5.828471 2.796872e-09
## 5 2.823879 -0.005714286 0.1821012 6.630828 1.669049e-11
## 6 1.142206 -0.005714286 0.1119585 3.430705 3.010071e-047.10 Mapping both local Moran’s I values and p-values for week 15
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w15 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w15)
localMI.map_w15 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w15) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w15 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w15) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w15, pvalue.map_w15, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and a few regions that have greater auto correlation than the rest.
For the local moran p value map, we can see only more regions are statistically significant compared to the previous week while the rest are still most statistically insignificant.
7.11 Plotting Moran scatterplot for week 15 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate) %>% as.vector
nci15 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 15 Rate", ylab="Spatially Lag z-Week 15 Rate")
7.12 prepare a LISA cluster map for week 15
quadrant_15 <- vector(mode="numeric",length=nrow(localMI_w15))
DV_15 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate)
C_mI_15 <- localMI_w15[,1] - mean(localMI_w15[,1])
signif <- 0.05
quadrant_15[DV_15 >0 & C_mI_15>0] <- 4
quadrant_15[DV_15 <0 & C_mI_15<0] <- 1
quadrant_15[DV_15 <0 & C_mI_15>0] <- 2
quadrant_15[DV_15 >0 & C_mI_15<0] <- 3
quadrant_15[localMI_w15[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w15$quadrant15 <- quadrant_15
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w15) +
tm_fill(col = "quadrant15", style = "cat", palette = colors[c(sort(unique(quadrant_15)))+1], labels = clusters[c(sort(unique(quadrant_15)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant and expanded slightly south as well. We can conclude that that is a cluster formed with high values surrounded by high values.
7.13 Local Moran’s I for Week 16
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate)
localMI_w16 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate, rswm_q)
head(localMI_w16)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 5.352148 -0.005714286 0.1882250 12.349601 2.447695e-35
## 2 4.923713 -0.005714286 0.1882250 11.362080 3.229426e-30
## 3 8.477413 -0.005714286 0.2366100 17.439718 2.060512e-68
## 4 2.806128 -0.005714286 0.1156475 8.268427 6.786917e-17
## 5 6.554371 -0.005714286 0.1882250 15.120663 5.917243e-52
## 6 2.586863 -0.005714286 0.1156475 7.623660 1.232912e-147.14 Mapping both local Moran’s I values and p-values for week 16
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w16 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w16)
localMI.map_w16 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w16) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w16 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w16) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w16, pvalue.map_w16, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and a few regions have greater auto correlation than the rest.
For the local moran p value map, we can see some regions are statistically significant while the rest still remained as statistically insignificant.
7.15 Plotting Moran scatterplot for week 16 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate) %>% as.vector
nci16 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 16 Rate", ylab="Spatially Lag z-Week 16 Rate")
7.16 prepare a LISA cluster map for week 16
quadrant_16 <- vector(mode="numeric",length=nrow(localMI_w16))
DV_16 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate)
C_mI_16 <- localMI_w16[,1] - mean(localMI_w16[,1])
signif <- 0.05
quadrant_16[DV_16 >0 & C_mI_16>0] <- 4
quadrant_16[DV_16 <0 & C_mI_16<0] <- 1
quadrant_16[DV_16 <0 & C_mI_16>0] <- 2
quadrant_16[DV_16 >0 & C_mI_16<0] <- 3
quadrant_16[localMI_w16[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w16$quadrant16 <- quadrant_16
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w16) +
tm_fill(col = "quadrant16", style = "cat", palette = colors[c(sort(unique(quadrant_16)))+1], labels = clusters[c(sort(unique(quadrant_16)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant and expanded to a new area at the east as well. We can conclude that that is a cluster formed with high values surrounded by high values.
7.17 Local Moran’s I for Week 17
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate)
localMI_w17 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate, rswm_q)
head(localMI_w17)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.617384 -0.005714286 0.1893433 10.624492 1.146289e-26
## 2 5.250207 -0.005714286 0.1893433 12.078803 6.834785e-34
## 3 5.089952 -0.005714286 0.2380247 10.444556 7.758962e-26
## 4 3.048399 -0.005714286 0.1163212 8.954796 1.701815e-19
## 5 7.548933 -0.005714286 0.1893433 17.361580 8.061205e-68
## 6 3.721528 -0.005714286 0.1163212 10.928443 4.214122e-287.18 Mapping both local Moran’s I values and p-values for week 17
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w17 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w17)
localMI.map_w17 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w17) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w17 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w17) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w17, pvalue.map_w17, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and some regions have greater auto correlation than the rest.
For the local moran p value map, we can see only a few regions are statistically significant while the rest are most statistically insignificant, however, there is some signs of clustering at the south west area.
7.19 Plotting Moran scatterplot for week 17 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate) %>% as.vector
nci17 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 17 Rate", ylab="Spatially Lag z-Week 17 Rate")
7.20 prepare a LISA cluster map for week 17
quadrant_17 <- vector(mode="numeric",length=nrow(localMI_w17))
DV_17 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate)
C_mI_17 <- localMI_w17[,1] - mean(localMI_w17[,1])
signif <- 0.05
quadrant_17[DV_17 >0 & C_mI_17>0] <- 4
quadrant_17[DV_17 <0 & C_mI_17<0] <- 1
quadrant_17[DV_17 <0 & C_mI_17>0] <- 2
quadrant_17[DV_17 >0 & C_mI_17<0] <- 3
quadrant_17[localMI_w17[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w17$quadrant17 <- quadrant_17
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w17) +
tm_fill(col = "quadrant17", style = "cat", palette = colors[c(sort(unique(quadrant_17)))+1], labels = clusters[c(sort(unique(quadrant_17)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant and expanded to a new area at the east as well. We can conclude that that is a cluster formed with high values surrounded by high values.
7.21 Local Moran’s I for Week 18
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate)
localMI_w18 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate, rswm_q)
head(localMI_w18)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.576061 -0.005714286 0.1895715 10.523187 3.377458e-26
## 2 5.187772 -0.005714286 0.1895715 11.928134 4.222249e-33
## 3 3.563690 -0.005714286 0.2383133 7.311754 1.318387e-13
## 4 3.061775 -0.005714286 0.1164586 8.988707 1.250779e-19
## 5 8.639745 -0.005714286 0.1895715 19.856451 4.846153e-88
## 6 4.581027 -0.005714286 0.1164586 13.440593 1.748110e-417.22 Mapping both local Moran’s I values and p-values for week 18
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w18 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w18)
localMI.map_w18 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w18) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w18 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w18) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w18, pvalue.map_w18, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and some regions that have greater auto correlation than the rest.
For the local moran p value map, we can see only a few regions are statistically significant while the rest are most statistically insignificant.
7.23 Plotting Moran scatterplot for week 18 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate) %>% as.vector
nci18 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 18 Rate", ylab="Spatially Lag z-Week 18 Rate")
7.24 prepare a LISA cluster map for week 18
quadrant_18 <- vector(mode="numeric",length=nrow(localMI_w18))
DV_18 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate)
C_mI_18 <- localMI_w18[,1] - mean(localMI_w18[,1])
signif <- 0.05
quadrant_18[DV_18 >0 & C_mI_18>0] <- 4
quadrant_18[DV_18 <0 & C_mI_18<0] <- 1
quadrant_18[DV_18 <0 & C_mI_18>0] <- 2
quadrant_18[DV_18 >0 & C_mI_18<0] <- 3
quadrant_18[localMI_w18[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w18$quadrant18 <- quadrant_18
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w18) +
tm_fill(col = "quadrant18", style = "cat", palette = colors[c(sort(unique(quadrant_18)))+1], labels = clusters[c(sort(unique(quadrant_18)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant. We can conclude that that is a cluster formed with high values surrounded by high values.
Also, there is an area at the south west region being deemed as low high significance, which could potentially be a cluster.
7.25 Local Moran’s I for Week 19
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate)
localMI_w19 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate, rswm_q)
head(localMI_w19)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.873789 -0.005714286 0.1894634 11.210189 1.817073e-29
## 2 5.347326 -0.005714286 0.1894634 12.298094 4.636816e-35
## 3 2.501142 -0.005714286 0.2381766 5.136649 1.398406e-07
## 4 3.249307 -0.005714286 0.1163935 9.540901 7.079846e-22
## 5 8.786938 -0.005714286 0.1894634 20.200272 4.867990e-91
## 6 4.838074 -0.005714286 0.1163935 14.197788 4.727534e-467.26 Mapping both local Moran’s I values and p-values for week 19
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w19 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w19)
localMI.map_w19 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w19) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w19 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w19) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w19, pvalue.map_w19, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and some regions that have greater auto correlation than the rest.
For the local moran p value map, we can see only a few regions are statistically significant while the rest are most statistically insignificant.
7.27 Plotting Moran scatterplot for week 19 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate) %>% as.vector
nci19 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 19 Rate", ylab="Spatially Lag z-Week 19 Rate")
7.28 prepare a LISA cluster map for week 19
quadrant_19 <- vector(mode="numeric",length=nrow(localMI_w19))
DV_19 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate)
C_mI_19 <- localMI_w19[,1] - mean(localMI_w19[,1])
signif <- 0.05
quadrant_19[DV_19 >0 & C_mI_19>0] <- 4
quadrant_19[DV_19 <0 & C_mI_19<0] <- 1
quadrant_19[DV_19 <0 & C_mI_19>0] <- 2
quadrant_19[DV_19 >0 & C_mI_19<0] <- 3
quadrant_19[localMI_w18[,5]>signif] <- 0spdf_municipalities_COVID_rates_w13_to_w32.localMI_w19$quadrant19 <- quadrant_19
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w19) +
tm_fill(col = "quadrant19", style = "cat", palette = colors[c(sort(unique(quadrant_19)))+1], labels = clusters[c(sort(unique(quadrant_19)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant is also shown as highly significant. We can conclude that that is a cluster formed with high values surrounded by high values.
However, the map seems slightly different from the previous week’s map, we can see that 2 areas beside the original cluster identified are now deemed as high-low.
The low high cluster still remains the same as the previous week.
7.29 Local Moran’s I for Week 20
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate)
localMI_w20 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate, rswm_q)
head(localMI_w20)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 5.206203 -0.005714286 0.1892798 11.979683 2.270276e-33
## 2 5.026026 -0.005714286 0.1892798 11.565543 3.080027e-31
## 3 1.918684 -0.005714286 0.2379444 3.945093 3.988444e-05
## 4 3.316441 -0.005714286 0.1162830 9.742310 9.948965e-23
## 5 8.502413 -0.005714286 0.1892798 19.556081 1.830798e-85
## 6 4.709147 -0.005714286 0.1162830 13.826458 8.825011e-447.30 Mapping both local Moran’s I values and p-values for week 20
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w20 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w20)
localMI.map_w20 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w20) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w20 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w20) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w20, pvalue.map_w20, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest when it was still at week 19.
For the local moran p value map, we can see more regions are statistically significant compared to when it was in week 19 while the rest are still most statistically insignificant.
7.31 Plotting Moran scatterplot for week 20 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate) %>% as.vector
nci20 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 20 Rate", ylab="Spatially Lag z-Week 20 Rate")
7.32 prepare a LISA cluster map for week 20
quadrant_20 <- vector(mode="numeric",length=nrow(localMI_w20))
DV_20 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate)
C_mI_20 <- localMI_w20[,1] - mean(localMI_w20[,1])
signif <- 0.05
quadrant_20[DV_20 >0 & C_mI_20>0] <- 4
quadrant_20[DV_20 <0 & C_mI_20<0] <- 1
quadrant_20[DV_20 <0 & C_mI_20>0] <- 2
quadrant_20[DV_20 >0 & C_mI_20<0] <- 3
quadrant_20[localMI_w20[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w20## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, 0, 0, 0, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 30.1353940202768, 44.1268269582625, 55.8222298756556, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w20$quadrant20 <- quadrant_20
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w20) +
tm_fill(col = "quadrant20", style = "cat", palette = colors[c(sort(unique(quadrant_20)))+1], labels = clusters[c(sort(unique(quadrant_20)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 19 has now expanded while there is also a cluster forming at the south west region with low-high statistically significant.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.33 Local Moran’s I for Week 21
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate)
localMI_w21 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate, rswm_q)
head(localMI_w21)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 5.222450 -0.005714286 0.1884320 12.044032 1.042554e-33
## 2 4.658494 -0.005714286 0.1884320 10.744857 3.132936e-27
## 3 1.741859 -0.005714286 0.2368718 3.590696 1.648983e-04
## 4 3.274386 -0.005714286 0.1157722 9.640176 2.704748e-22
## 5 7.776664 -0.005714286 0.1884320 17.928132 3.556450e-72
## 6 4.444203 -0.005714286 0.1157722 13.078254 2.191978e-397.34 Mapping both local Moran’s I values and p-values for week 21
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w21 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w21)
localMI.map_w21 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w21) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w21 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w21) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w21, pvalue.map_w21, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest.
For the local moran p value map, we can see more regions are statistically significant while the rest are still most statistically insignificant.
7.35 Plotting Moran scatterplot for week 21 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate) %>% as.vector
nci21 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 21 Rate", ylab="Spatially Lag z-Week 21 Rate")
7.36 prepare a LISA cluster map for week 21
quadrant_21 <- vector(mode="numeric",length=nrow(localMI_w21))
DV_21 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate)
C_mI_21 <- localMI_w21[,1] - mean(localMI_w21[,1])
signif <- 0.05
quadrant_21[DV_21 >0 & C_mI_21>0] <- 4
quadrant_21[DV_21 <0 & C_mI_21<0] <- 1
quadrant_21[DV_21 <0 & C_mI_21>0] <- 2
quadrant_21[DV_21 >0 & C_mI_21<0] <- 3
quadrant_21[localMI_w21[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w21## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, 0, 0, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 44.1268269582625, 55.8222298756556, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w21$quadrant21 <- quadrant_21
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w21) +
tm_fill(col = "quadrant21", style = "cat", palette = colors[c(sort(unique(quadrant_21)))+1], labels = clusters[c(sort(unique(quadrant_21)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 20 still remains the same while there is also a cluster forming at the south west region with low-high statistically significant along with a region with low low.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.37 Local Moran’s I for Week 22
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate)
localMI_w22 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate, rswm_q)
head(localMI_w22)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.942213 -0.005714286 0.1884614 11.397566 2.149853e-30
## 2 4.484690 -0.005714286 0.1884614 10.343661 2.235604e-25
## 3 1.886702 -0.005714286 0.2369090 3.887996 5.053757e-05
## 4 3.046531 -0.005714286 0.1157899 8.969826 1.484913e-19
## 5 6.893008 -0.005714286 0.1884614 15.891230 3.644180e-57
## 6 3.876947 -0.005714286 0.1157899 11.410224 1.858862e-307.38 Mapping both local Moran’s I values and p-values for week 22
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w22 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w22)
localMI.map_w22 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w22) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w22 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w22) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w22, pvalue.map_w22, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest.
For the local moran p value map, we can see more regions are statistically significant while the rest are still most statistically insignificant.
7.39 Plotting Moran scatterplot for week 22 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate) %>% as.vector
nci22 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 22 Rate", ylab="Spatially Lag z-Week 22 Rate")
7.40 prepare a LISA cluster map for week 22
quadrant_22 <- vector(mode="numeric",length=nrow(localMI_w22))
DV_22 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate)
C_mI_22 <- localMI_w22[,1] - mean(localMI_w22[,1])
signif <- 0.05
quadrant_22[DV_22 >0 & C_mI_22>0] <- 4
quadrant_22[DV_22 <0 & C_mI_22<0] <- 1
quadrant_22[DV_22 <0 & C_mI_22>0] <- 2
quadrant_22[DV_22 >0 & C_mI_22<0] <- 3
quadrant_22[localMI_w22[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w22## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, 0, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 55.8222298756556, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w22$quadrant22 <- quadrant_22
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w22) +
tm_fill(col = "quadrant22", style = "cat", palette = colors[c(sort(unique(quadrant_22)))+1], labels = clusters[c(sort(unique(quadrant_22)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 21 still remains the same with a region at the west deemed as high low.
Also, there is also a cluster forming at the south west region with low-high statistically significant.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.41 Local Moran’s I for Week 23
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate)
localMI_w23 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate, rswm_q)
head(localMI_w23)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.752828 -0.005714286 0.1877510 10.982034 2.331512e-28
## 2 4.225626 -0.005714286 0.1877510 9.765327 7.929952e-23
## 3 1.839598 -0.005714286 0.2360104 3.798433 7.280701e-05
## 4 2.962513 -0.005714286 0.1153620 8.739082 1.175060e-18
## 5 6.067700 -0.005714286 0.1877510 14.016570 6.171857e-45
## 6 3.429137 -0.005714286 0.1153620 10.112921 2.420945e-247.42 Mapping both local Moran’s I values and p-values for week 23
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w23 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w23)
localMI.map_w23 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w23) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w23 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w23) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w23, pvalue.map_w23, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest.
For the local moran p value map, we can see more regions are statistically significant while the rest are still most statistically insignificant.
7.43 Plotting Moran scatterplot for week 23 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate) %>% as.vector
nci23 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 23 Rate", ylab="Spatially Lag z-Week 23 Rate")
7.44 prepare a LISA cluster map for week 23
quadrant_23 <- vector(mode="numeric",length=nrow(localMI_w23))
DV_23 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate)
C_mI_23 <- localMI_w23[,1] - mean(localMI_w23[,1])
signif <- 0.05
quadrant_23[DV_23 >0 & C_mI_23>0] <- 4
quadrant_23[DV_23 <0 & C_mI_23<0] <- 1
quadrant_23[DV_23 <0 & C_mI_23>0] <- 2
quadrant_23[DV_23 >0 & C_mI_23<0] <- 3
quadrant_23[localMI_w23[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w23## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w23$quadrant23 <- quadrant_23
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w23) +
tm_fill(col = "quadrant23", style = "cat", palette = colors[c(sort(unique(quadrant_23)))+1], labels = clusters[c(sort(unique(quadrant_23)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 22 still remains the same while there is also a cluster forming at the south west region with low-high statistically significant along with a region with low low.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.45 Local Moran’s I for Week 24
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate)
localMI_w24 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate, rswm_q)
head(localMI_w24)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.319319 -0.005714286 0.1869899 10.001849 7.478887e-24
## 2 3.752512 -0.005714286 0.1869899 8.691081 1.795031e-18
## 3 1.844210 -0.005714286 0.2350476 3.815717 6.789401e-05
## 4 2.624207 -0.005714286 0.1149035 7.758470 4.297989e-15
## 5 5.090588 -0.005714286 0.1869899 11.785448 2.319740e-32
## 6 2.878158 -0.005714286 0.1149035 8.507646 8.875021e-187.46 Mapping both local Moran’s I values and p-values for week 24
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w24 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w24)
localMI.map_w24 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w24) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w24 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w24) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w24, pvalue.map_w24, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest.
For the local moran p value map, we can see more regions are statistically significant while the rest are still most statistically insignificant.
7.47 Plotting Moran scatterplot for week 24 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate) %>% as.vector
nci24 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 24 Rate", ylab="Spatially Lag z-Week 24 Rate")
7.48 prepare a LISA cluster map for week 24
quadrant_24 <- vector(mode="numeric",length=nrow(localMI_w24))
DV_24 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate)
C_mI_24 <- localMI_w24[,1] - mean(localMI_w24[,1])
signif <- 0.05
quadrant_24[DV_24 >0 & C_mI_24>0] <- 4
quadrant_24[DV_24 <0 & C_mI_24<0] <- 1
quadrant_24[DV_24 <0 & C_mI_24>0] <- 2
quadrant_24[DV_24 >0 & C_mI_24<0] <- 3
quadrant_24[localMI_w24[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w24## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w24$quadrant24 <- quadrant_24
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w24) +
tm_fill(col = "quadrant24", style = "cat", palette = colors[c(sort(unique(quadrant_24)))+1], labels = clusters[c(sort(unique(quadrant_24)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 23 still remains the same while there is also a cluster forming at the south west region with low-high statistically significant.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.49 Local Moran’s I for Week 25
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate)
localMI_w25 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate, rswm_q)
head(localMI_w25)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.221236 -0.005714286 0.1862921 9.793321 6.013972e-23
## 2 3.644553 -0.005714286 0.1862921 8.457217 1.369165e-17
## 3 1.898852 -0.005714286 0.2341647 3.935821 4.145644e-05
## 4 2.423037 -0.005714286 0.1144831 7.178146 3.533156e-13
## 5 4.579407 -0.005714286 0.1862921 10.623158 1.162794e-26
## 6 2.576927 -0.005714286 0.1144831 7.632968 1.147052e-147.50 Mapping both local Moran’s I values and p-values for week 25
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w25 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w25)
localMI.map_w25 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w25) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w25 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w25) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w25, pvalue.map_w25, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest.
For the local moran p value map, we can see more regions are statistically significant while the rest are still most statistically insignificant.
7.51 Plotting Moran scatterplot for week 25 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate) %>% as.vector
nci25 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 25 Rate", ylab="Spatially Lag z-Week 25 Rate")
7.52 prepare a LISA cluster map for week 25
quadrant_25 <- vector(mode="numeric",length=nrow(localMI_w25))
DV_25 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate)
C_mI_25 <- localMI_w25[,1] - mean(localMI_w25[,1])
signif <- 0.05
quadrant_25[DV_25 >0 & C_mI_25>0] <- 4
quadrant_25[DV_25 <0 & C_mI_25<0] <- 1
quadrant_25[DV_25 <0 & C_mI_25>0] <- 2
quadrant_25[DV_25 >0 & C_mI_25<0] <- 3
quadrant_25[localMI_w25[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w25## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w25$quadrant25 <- quadrant_25
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w25) +
tm_fill(col = "quadrant25", style = "cat", palette = colors[c(sort(unique(quadrant_25)))+1], labels = clusters[c(sort(unique(quadrant_25)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 24 still remains the same while there is also a cluster forming at the south west region with low-high statistically significant.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.53 Local Moran’s I for Week 26
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate)
localMI_w26 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate, rswm_q)
head(localMI_w26)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.171369 -0.005714286 0.1863513 9.676245 1.902155e-22
## 2 3.507135 -0.005714286 0.1863513 8.137542 2.016922e-16
## 3 1.874779 -0.005714286 0.2342397 3.885452 5.106979e-05
## 4 2.287963 -0.005714286 0.1145188 6.777882 6.097506e-12
## 5 4.236309 -0.005714286 0.1863513 9.826679 4.321084e-23
## 6 2.389023 -0.005714286 0.1145188 7.076517 7.391118e-137.54 Mapping both local Moran’s I values and p-values for week 26
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w26 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w26)
localMI.map_w26 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w26) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w26 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w26) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w26, pvalue.map_w26, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions have positive autocorrelation and more regions that have greater auto correlation than the rest.
For the local moran p value map, we can see more regions are statistically significant while the rest are still most statistically insignificant.
7.55 Plotting Moran scatterplot for week 26 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate) %>% as.vector
nci26 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 26 Rate", ylab="Spatially Lag z-Week 26 Rate")
7.56 prepare a LISA cluster map for week 26
quadrant_26 <- vector(mode="numeric",length=nrow(localMI_w26))
DV_26 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate)
C_mI_26 <- localMI_w26[,1] - mean(localMI_w26[,1])
signif <- 0.05
quadrant_26[DV_26 >0 & C_mI_26>0] <- 4
quadrant_26[DV_26 <0 & C_mI_26<0] <- 1
quadrant_26[DV_26 <0 & C_mI_26>0] <- 2
quadrant_26[DV_26 >0 & C_mI_26<0] <- 3
quadrant_26[localMI_w26[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w26## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w26$quadrant26 <- quadrant_26
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w26) +
tm_fill(col = "quadrant26", style = "cat", palette = colors[c(sort(unique(quadrant_26)))+1], labels = clusters[c(sort(unique(quadrant_26)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 25 still remains the same while there is also a cluster forming at the south west region with low-high statistically significant.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.57 Local Moran’s I for Week 27
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate)
localMI_w27 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate, rswm_q)
head(localMI_w27)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 4.058873 -0.005714286 0.1859911 9.424762 2.155417e-21
## 2 3.531871 -0.005714286 0.1859911 8.202776 1.174493e-16
## 3 1.996319 -0.005714286 0.2337840 4.140606 1.731948e-05
## 4 2.184905 -0.005714286 0.1143018 6.479485 4.601815e-11
## 5 3.895230 -0.005714286 0.1859911 9.045315 7.461743e-20
## 6 2.254400 -0.005714286 0.1143018 6.685037 1.154329e-117.58 Mapping both local Moran’s I values and p-values for week 27
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w27 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w27)
localMI.map_w27 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w27) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w27 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w27) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w27, pvalue.map_w27, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions still have positive autocorrelation and the regions that have greater auto correlation, remains relative similar now as well.
For the local moran p value map, we can see the regions are statistically significant remains relatively similar while the rest are still most statistically insignificant.
7.59 Plotting Moran scatterplot for week 27 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate) %>% as.vector
nci27 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 27 Rate", ylab="Spatially Lag z-Week 27 Rate")
7.60 prepare a LISA cluster map for week 27
quadrant_27 <- vector(mode="numeric",length=nrow(localMI_w27))
DV_27 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate)
C_mI_27 <- localMI_w27[,1] - mean(localMI_w27[,1])
signif <- 0.05
quadrant_27[DV_27 >0 & C_mI_27>0] <- 4
quadrant_27[DV_27 <0 & C_mI_27<0] <- 1
quadrant_27[DV_27 <0 & C_mI_27>0] <- 2
quadrant_27[DV_27 >0 & C_mI_27<0] <- 3
quadrant_27[localMI_w27[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w27## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w27$quadrant27 <- quadrant_27
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w27) +
tm_fill(col = "quadrant27", style = "cat", palette = colors[c(sort(unique(quadrant_27)))+1], labels = clusters[c(sort(unique(quadrant_27)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 26 still remains the same while there is also a cluster forming at the south west region with low-high statistically significant along with a region with low low.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.61 Local Moran’s I for Week 28
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate)
localMI_w28 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate, rswm_q)
head(localMI_w28)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 3.882421 -0.005714286 0.1855013 9.027508 8.781276e-20
## 2 3.430963 -0.005714286 0.1855013 7.979309 7.357729e-16
## 3 2.037457 -0.005714286 0.2331644 4.231298 1.161734e-05
## 4 2.035897 -0.005714286 0.1140068 6.046552 7.398942e-10
## 5 3.568091 -0.005714286 0.1855013 8.297692 5.307671e-17
## 6 2.098522 -0.005714286 0.1140068 6.232024 2.302230e-107.62 Mapping both local Moran’s I values and p-values for week 28
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w28 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w28)
localMI.map_w28 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w28) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w28 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w28) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w28, pvalue.map_w28, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions still have positive autocorrelation and the regions that have greater auto correlation, remains relative similar now as well.
For the local moran p value map, we can see the regions are statistically significant remains relatively similar while the rest are still most statistically insignificant.
7.63 Plotting Moran scatterplot for week 28 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate) %>% as.vector
nci28 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 28 Rate", ylab="Spatially Lag z-Week 28 Rate")
7.64 prepare a LISA cluster map for week 28
quadrant_28 <- vector(mode="numeric",length=nrow(localMI_w28))
DV_28 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate)
C_mI_28 <- localMI_w28[,1] - mean(localMI_w28[,1])
signif <- 0.05
quadrant_28[DV_28 >0 & C_mI_28>0] <- 4
quadrant_28[DV_28 <0 & C_mI_28<0] <- 1
quadrant_28[DV_28 <0 & C_mI_28>0] <- 2
quadrant_28[DV_28 >0 & C_mI_28<0] <- 3
quadrant_28[localMI_w28[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w28## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w28$quadrant28 <- quadrant_28
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w28) +
tm_fill(col = "quadrant28", style = "cat", palette = colors[c(sort(unique(quadrant_28)))+1], labels = clusters[c(sort(unique(quadrant_28)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 27 still remains the same while the cluster forming at the south west region with low-high statistically significant.
We can conclude that there is a cluster formed with high values surrounded by high values and also a cluster formed with low-high values surrounded by low-high values.
7.65 Local Moran’s I for Week 29
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate)
localMI_w29 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate, rswm_q)
head(localMI_w29)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 3.797299 -0.005714286 0.1857596 8.823731 5.535406e-19
## 2 3.622823 -0.005714286 0.1857596 8.418912 1.899986e-17
## 3 2.167046 -0.005714286 0.2334911 4.496521 3.453724e-06
## 4 1.946904 -0.005714286 0.1141623 5.779044 3.756318e-09
## 5 3.323899 -0.005714286 0.1857596 7.725350 5.577332e-15
## 6 2.015667 -0.005714286 0.1141623 5.982558 1.098302e-097.66 Mapping both local Moran’s I values and p-values for week 29
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w29 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w29)
localMI.map_w29 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w29) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w29 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w29) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w29, pvalue.map_w29, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions still have positive autocorrelation and the regions that have greater auto correlation, remains relative similar now as well.
For the local moran p value map, we can see the regions are statistically significant remains relatively similar while the rest are still most statistically insignificant.
7.67 Plotting Moran scatterplot for week 29 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate) %>% as.vector
nci29 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 29 Rate", ylab="Spatially Lag z-Week 29 Rate")
7.68 prepare a LISA cluster map for week 29
quadrant_29 <- vector(mode="numeric",length=nrow(localMI_w29))
DV_29 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate)
C_mI_29 <- localMI_w29[,1] - mean(localMI_w29[,1])
signif <- 0.05
quadrant_29[DV_29 >0 & C_mI_29>0] <- 4
quadrant_29[DV_29 <0 & C_mI_29<0] <- 1
quadrant_29[DV_29 <0 & C_mI_29>0] <- 2
quadrant_29[DV_29 >0 & C_mI_29<0] <- 3
quadrant_29[localMI_w29[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w29## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w29$quadrant29 <- quadrant_29
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w29) +
tm_fill(col = "quadrant29", style = "cat", palette = colors[c(sort(unique(quadrant_29)))+1], labels = clusters[c(sort(unique(quadrant_29)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 28 still remains the same while the cluster forming at the south west region with low-high statistically significant has become smaller this week and only 1 region remains in the low high range.
We can conclude that there is a cluster formed with high values surrounded by high values.
7.69 Local Moran’s I for Week 30
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate)
localMI_w30 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate, rswm_q)
head(localMI_w30)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 3.726912 -0.005714286 0.1857550 8.660526 2.347963e-18
## 2 3.655502 -0.005714286 0.1857550 8.494839 9.910364e-18
## 3 2.355058 -0.005714286 0.2334853 4.885671 5.153849e-07
## 4 1.887906 -0.005714286 0.1141596 5.604498 1.044294e-08
## 5 3.194746 -0.005714286 0.1857550 7.425781 5.605826e-14
## 6 1.914873 -0.005714286 0.1141596 5.684310 6.567080e-097.70 Mapping both local Moran’s I values and p-values for week 30
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w30 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w30)
localMI.map_w30 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w30) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w30 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w30) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w30, pvalue.map_w30, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions still have positive autocorrelation and the regions that have greater auto correlation, remains relative similar now as well.
For the local moran p value map, we can see the regions are statistically significant remains relatively similar while the rest are still most statistically insignificant.
7.71 Plotting Moran scatterplot for week 30 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate) %>% as.vector
nci30 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 30 Rate", ylab="Spatially Lag z-Week 30 Rate")
7.72 prepare a LISA cluster map for week 30
quadrant_30 <- vector(mode="numeric",length=nrow(localMI_w30))
DV_30 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate)
C_mI_30 <- localMI_w30[,1] - mean(localMI_w30[,1])
signif <- 0.05
quadrant_30[DV_30 >0 & C_mI_30>0] <- 4
quadrant_30[DV_30 <0 & C_mI_30<0] <- 1
quadrant_30[DV_30 <0 & C_mI_30>0] <- 2
quadrant_30[DV_30 >0 & C_mI_30<0] <- 3
quadrant_30[localMI_w30[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w30## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w30$quadrant30 <- quadrant_30
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w30) +
tm_fill(col = "quadrant30", style = "cat", palette = colors[c(sort(unique(quadrant_30)))+1], labels = clusters[c(sort(unique(quadrant_30)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 29 still remains the same while the cluster forming at the south west region with low-high statistically significant has disappeared and is now deemed at insignifcant.
We can conclude that there is a cluster formed with high values surrounded by high values.
7.73 Local Moran’s I for Week 31
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate)
localMI_w31 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate, rswm_q)
head(localMI_w31)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 3.728797 -0.005714286 0.1857269 8.665556 2.246581e-18
## 2 3.869239 -0.005714286 0.1857269 8.991438 1.220088e-19
## 3 2.399819 -0.005714286 0.2334497 4.978687 3.200860e-07
## 4 1.880550 -0.005714286 0.1141426 5.583142 1.181057e-08
## 5 3.223644 -0.005714286 0.1857269 7.493399 3.355619e-14
## 6 1.933110 -0.005714286 0.1141426 5.738714 4.769910e-097.74 Mapping both local Moran’s I values and p-values for week 31
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w31 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w31)
localMI.map_w31 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w31) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w31 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w31) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w31, pvalue.map_w31, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions still have positive autocorrelation and the regions that have greater auto correlation, remains relative similar now as well.
For the local moran p value map, we can see the regions are statistically significant remains relatively similar while the rest are still most statistically insignificant.
7.75 Plotting Moran scatterplot for week 31 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate) %>% as.vector
nci31 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 31 Rate", ylab="Spatially Lag z-Week 31 Rate")
7.76 prepare a LISA cluster map for week 31
quadrant_31 <- vector(mode="numeric",length=nrow(localMI_w31))
DV_31 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate)
C_mI_31 <- localMI_w31[,1] - mean(localMI_w31[,1])
signif <- 0.05
quadrant_31[DV_31 >0 & C_mI_31>0] <- 4
quadrant_31[DV_31 <0 & C_mI_31<0] <- 1
quadrant_31[DV_31 <0 & C_mI_31>0] <- 2
quadrant_31[DV_31 >0 & C_mI_31<0] <- 3
quadrant_31[localMI_w31[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w31## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w31$quadrant31 <- quadrant_31
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w31) +
tm_fill(col = "quadrant31", style = "cat", palette = colors[c(sort(unique(quadrant_31)))+1], labels = clusters[c(sort(unique(quadrant_31)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 30 still remains the same.
We can conclude that there is a cluster formed with high values surrounded by high values.
7.77 Local Moran’s I for Week 32
fips <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate)
localMI_w32 <- localmoran(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate, rswm_q)
head(localMI_w32)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 1 3.740132 -0.005714286 0.1857746 8.690739 1.800447e-18
## 2 3.897909 -0.005714286 0.1857746 9.056799 6.716768e-20
## 3 2.407429 -0.005714286 0.2335101 4.993789 2.960310e-07
## 4 1.885342 -0.005714286 0.1141714 5.596621 1.092849e-08
## 5 3.291571 -0.005714286 0.1857746 7.650033 1.004635e-14
## 6 1.966097 -0.005714286 0.1141714 5.835616 2.679608e-097.78 Mapping both local Moran’s I values and p-values for week 32
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w32 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32,localMI_w32)
localMI.map_w32 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w32) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map_w32 <- tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w32) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map_w32, pvalue.map_w32, asp=1, ncol=2)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "Ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
For the local moran map, we can see that most regions still have positive autocorrelation and the regions that have greater auto correlation, remains relative similar now as well.
For the local moran p value map, we can see the regions are statistically significant remains relatively similar while the rest are still most statistically insignificant.
7.79 Plotting Moran scatterplot for week 32 with standardised variable
spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate <- scale(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate) %>% as.vector
nci32 <- moran.plot(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate, rswm_q, labels=as.character(spdf_municipalities_COVID_rates_w13_to_w32$NOMGEO), xlab="z-Week 32 Rate", ylab="Spatially Lag z-Week 32 Rate")
7.80 prepare a LISA cluster map for week 32
quadrant_32 <- vector(mode="numeric",length=nrow(localMI_w32))
DV_32 <- spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate - mean(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate)
C_mI_32 <- localMI_w32[,1] - mean(localMI_w32[,1])
signif <- 0.05
quadrant_32[DV_32 >0 & C_mI_32>0] <- 4
quadrant_32[DV_32 <0 & C_mI_32<0] <- 1
quadrant_32[DV_32 <0 & C_mI_32>0] <- 2
quadrant_32[DV_32 >0 & C_mI_32<0] <- 3
quadrant_32[localMI_w32[,5]>signif] <- 0
spdf_municipalities_COVID_rates_w13_to_w32.localMI_w32## class : SpatialPolygonsDataFrame
## features : 176
## extent : 2645804, 2855437, 707708.1, 921773.9 (xmin, xmax, ymin, ymax)
## crs : +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_defs
## variables : 30
## names : CVEGEO, CVE_ENT, CVE_MUN, NOMGEO, Pop2020, cumul13_rate, cumul14_rate, cumul15_rate, cumul16_rate, cumul17_rate, cumul18_rate, cumul19_rate, cumul20_rate, cumul21_rate, cumul22_rate, ...
## min values : 09002, 09, 001, Acambay de Ruíz Castañeda, 4236, -0.324549207597564, -0.452111871422612, -0.586901880299524, -0.735406597084209, -0.831681807833266, -0.896690434436167, -0.943636514653789, -0.963701234626591, -0.973811433117046, -1.02939048515341, ...
## max values : 17035, 17, 125, Zumpango, 1815551, 7.64012216045284, 6.67058117734819, 5.594342182035, 4.43202286955697, 4.06554058011044, 4.17228848410956, 4.1844231806437, 4.16133154367836, 4.68123238078507, 4.74276398609257, ...spdf_municipalities_COVID_rates_w13_to_w32.localMI_w32$quadrant32 <- quadrant_32
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.localMI_w32) +
tm_fill(col = "quadrant32", style = "cat", palette = colors[c(sort(unique(quadrant_32)))+1], labels = clusters[c(sort(unique(quadrant_32)))+1], popup.vars = c("Postal.Code")) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
From the LISA map above, we can see that the area that was previously deemed as statistically significant in week 31 still remains the same.
We can conclude that there is a cluster formed with high values surrounded by high values.
In conclusion, we can see that throughout the 20 weeks, there is some fluctuations in terms of being significant or insignificant. However, it is certain that the area in the north east is a cluster as it has appeared as high high for each of the 20 weeks.
During the 20 weeks, we can also see some clusters appearing at the south west region, but it has disappeared in recent weeks and is now deemed as insignificant.
8 Perform local Getis-Ord Gi* analysis
coords <- coordinates(spdf_municipalities_COVID_rates_w13_to_w32)8.1 Creating adaptive proximity matrix for week 13
knb13 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate))
knb13 <- include.self(knb13)
knb_lw13 <- nb2listw(knb13, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb13, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.2 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 13
fips13 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate)
gi.adaptive13 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul13_rate, knb_lw13)
spdf_municipalities_COVID_rates_w13_to_w32.gi13 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive13))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi13)[6] <- "gstat_adaptive13"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi13) +
tm_fill(col = "gstat_adaptive13",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive13" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
We could see that the ones that are in blue are cold spots while the ones that are slightly redder in colour are hotspot areas.
Currently at week 13, we could see some hotspots around Mexico but generally, there are more coldspots than hotspots.
8.3 Creating adaptive proximity matrix for week 14
knb14 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate))
knb14 <- include.self(knb14)
knb_lw14 <- nb2listw(knb14, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb14, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.4 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 14
fips14 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate)
gi.adaptive14 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul14_rate, knb_lw14)
spdf_municipalities_COVID_rates_w13_to_w32.gi14 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive14))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi14)[7] <- "gstat_adaptive14"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi14) +
tm_fill(col = "gstat_adaptive14",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive14" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 14, we could see that there are more hotspots appearing compared to the previous week.
8.5 Creating adaptive proximity matrix for week 15
knb15 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate))
knb15 <- include.self(knb15)
knb_lw15 <- nb2listw(knb15, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb15, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.6 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 15
fips15 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate)
gi.adaptive15 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul15_rate, knb_lw15)
spdf_municipalities_COVID_rates_w13_to_w32.gi15 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive15))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi15)[8] <- "gstat_adaptive15"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi15) +
tm_fill(col = "gstat_adaptive15",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive15" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 15, we could see that there are slightly lesser hotspots appearing compared to the previous week.
However, the hotspots that are previously identified, some of which are intensified especially around the north east region.
8.7 Creating adaptive proximity matrix for week 16
knb16 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate))
knb16 <- include.self(knb16)
knb_lw16 <- nb2listw(knb16, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb16, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.8 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 16
fips16 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate)
gi.adaptive16 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul16_rate, knb_lw16)
spdf_municipalities_COVID_rates_w13_to_w32.gi16 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive16))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi16)[9] <- "gstat_adaptive16"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi16) +
tm_fill(col = "gstat_adaptive16",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive16" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 16, the hotspots identified are somewhat similar to the previous week with some areas becoming coldspots and some becoming hotspots. However, the main cluster of hotspot still remains at the north east region.
8.9 Creating adaptive proximity matrix for week 17
knb17 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate))
knb17 <- include.self(knb17)
knb_lw17 <- nb2listw(knb17, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb17, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.10 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 17
fips17 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate)
gi.adaptive17 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul17_rate, knb_lw17)
spdf_municipalities_COVID_rates_w13_to_w32.gi17 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive17))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi17)[10] <- "gstat_adaptive17"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi17) +
tm_fill(col = "gstat_adaptive17",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive17" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 17, the hotspots identified are somewhat similar to the previous week with some areas becoming coldspots and some becoming hotspots. However, the main cluster of hotspot still remains at the north east region and has intensified in colour.
8.11 Creating adaptive proximity matrix for week 18
knb18 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate))
knb18 <- include.self(knb18)
knb_lw18 <- nb2listw(knb18, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb18, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.12 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 18
fips18 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate)
gi.adaptive18 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul18_rate, knb_lw18)
spdf_municipalities_COVID_rates_w13_to_w32.gi18 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive18))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi18)[11] <- "gstat_adaptive18"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi18) +
tm_fill(col = "gstat_adaptive18",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive18" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 18, the hotspots identified are somewhat similar to the previous week with some areas becoming coldspots and some becoming hotspots. However, the main cluster of hotspot still remains at the north east region.
8.13 Creating adaptive proximity matrix for week 19
knb19 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate))
knb19 <- include.self(knb19)
knb_lw19 <- nb2listw(knb19, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb19, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.14 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 19
fips19 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate)
gi.adaptive19 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul19_rate, knb_lw19)
spdf_municipalities_COVID_rates_w13_to_w32.gi19 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive19))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi19)[12] <- "gstat_adaptive19"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi19) +
tm_fill(col = "gstat_adaptive19",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive19" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 19, the hotspots identified are somewhat similar to the previous week with some areas becoming coldspots and some becoming hotspots. However, the main cluster of hotspot still remains at the north east region, 2 areas at the south region have intensified as well.
8.15 Creating adaptive proximity matrix for week 20
knb20 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate))
knb20 <- include.self(knb20)
knb_lw20 <- nb2listw(knb20, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb20, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.16 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 20
fips20 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate)
gi.adaptive20 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul20_rate, knb_lw20)
spdf_municipalities_COVID_rates_w13_to_w32.gi20 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive20))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi20)[13] <- "gstat_adaptive20"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi20) +
tm_fill(col = "gstat_adaptive20",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive20" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 20, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspot still remains at the north east region.
8.17 Creating adaptive proximity matrix for week 21
knb21 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate))
knb21 <- include.self(knb21)
knb_lw21 <- nb2listw(knb21, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb21, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.18 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 21
fips21 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate)
gi.adaptive21 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul21_rate, knb_lw21)
spdf_municipalities_COVID_rates_w13_to_w32.gi21 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive21))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi21)[14] <- "gstat_adaptive21"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi21) +
tm_fill(col = "gstat_adaptive21",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive21" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 21, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspot still remains at the north east region.
8.19 Creating adaptive proximity matrix for week 22
knb22 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate))
knb22 <- include.self(knb22)
knb_lw22 <- nb2listw(knb22, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb22, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.20 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 22
fips22 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate)
gi.adaptive22 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul22_rate, knb_lw22)
spdf_municipalities_COVID_rates_w13_to_w32.gi22 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive22))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi22)[15] <- "gstat_adaptive22"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi22) +
tm_fill(col = "gstat_adaptive22",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive22" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 22, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspot still remains at the north east region.
Additionally, 3 areas have dropped below -1 and deemed as coldspots.
8.21 Creating adaptive proximity matrix for week 23
knb23 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate))
knb23 <- include.self(knb23)
knb_lw23 <- nb2listw(knb23, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb23, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.22 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 23
fips23 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate)
gi.adaptive23 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul23_rate, knb_lw23)
spdf_municipalities_COVID_rates_w13_to_w32.gi23 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive23))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi23)[16] <- "gstat_adaptive23"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi23) +
tm_fill(col = "gstat_adaptive23",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive23" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 23, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspot still remains at the north east region. A few areas are now at 0 to 1 from the previous week’s -1 to 0.
8.23 Creating adaptive proximity matrix for week 24
knb24 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate))
knb24 <- include.self(knb24)
knb_lw24 <- nb2listw(knb24, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb24, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.24 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 24
fips24 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate)
gi.adaptive24 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul24_rate, knb_lw24)
spdf_municipalities_COVID_rates_w13_to_w32.gi24 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive24))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi24)[17] <- "gstat_adaptive24"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi24) +
tm_fill(col = "gstat_adaptive24",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive24" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 24, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
Areas near the south have intensified while an area at the west is now between 0 to 1.
8.25 Creating adaptive proximity matrix for week 25
knb25 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate))
knb25 <- include.self(knb25)
knb_lw25 <- nb2listw(knb25, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb25, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.26 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 25
fips25 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate)
gi.adaptive25 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul25_rate, knb_lw25)
spdf_municipalities_COVID_rates_w13_to_w32.gi25 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive25))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi25)[18] <- "gstat_adaptive25"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi25) +
tm_fill(col = "gstat_adaptive25",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive25" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 25, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
8.27 Creating adaptive proximity matrix for week 26
knb26 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate))
knb26 <- include.self(knb26)
knb_lw26 <- nb2listw(knb26, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb26, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.28 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 26
fips26 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate)
gi.adaptive26 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul26_rate, knb_lw26)
spdf_municipalities_COVID_rates_w13_to_w32.gi26 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive26))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi26)[19] <- "gstat_adaptive26"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi26) +
tm_fill(col = "gstat_adaptive26",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive26" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 26, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
Also, the hotspots are beginning to spread towards the west, with a few regions now between 0 to 1.
8.29 Creating adaptive proximity matrix for week 27
knb27 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate))
knb27 <- include.self(knb27)
knb_lw27 <- nb2listw(knb27, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb27, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.30 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 27
fips27 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate)
gi.adaptive27 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul27_rate, knb_lw27)
spdf_municipalities_COVID_rates_w13_to_w32.gi27 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive27))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi27)[20] <- "gstat_adaptive27"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi27) +
tm_fill(col = "gstat_adaptive27",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive27" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 27, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
Also, the hotspots are beginning to spread towards the north, with a few regions now between 0 to 1.
8.31 Creating adaptive proximity matrix for week 28
knb28 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate))
knb28 <- include.self(knb28)
knb_lw28 <- nb2listw(knb28, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb28, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.32 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 28
fips28 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate)
gi.adaptive28 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul28_rate, knb_lw28)
spdf_municipalities_COVID_rates_w13_to_w32.gi28 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive28))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi28)[21] <- "gstat_adaptive28"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi28) +
tm_fill(col = "gstat_adaptive28",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive28" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 28, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
Also, the hotspots are beginning to spread towards the north, with a few regions now between 0 to 1.
8.33 Creating adaptive proximity matrix for week 29
knb29 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate))
knb29 <- include.self(knb29)
knb_lw29 <- nb2listw(knb29, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb29, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.34 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 29
fips29 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate)
gi.adaptive29 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul29_rate, knb_lw29)
spdf_municipalities_COVID_rates_w13_to_w32.gi29 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive29))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi29)[22] <- "gstat_adaptive29"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi29) +
tm_fill(col = "gstat_adaptive29",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive29" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 29, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region while some hotspots are now within the 1 to 2 range.
8.35 Creating adaptive proximity matrix for week 30
knb30 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate))
knb30 <- include.self(knb30)
knb_lw30 <- nb2listw(knb30, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb30, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.36 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 30
fips30 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate)
gi.adaptive30 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul30_rate, knb_lw30)
spdf_municipalities_COVID_rates_w13_to_w32.gi30 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive30))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi30)[23] <- "gstat_adaptive30"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi30) +
tm_fill(col = "gstat_adaptive30",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive30" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 30, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
8.37 Creating adaptive proximity matrix for week 31
knb31 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate))
knb31 <- include.self(knb31)
knb_lw31 <- nb2listw(knb31, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb31, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.38 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 30
fips31 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate)
gi.adaptive31 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul31_rate, knb_lw31)
spdf_municipalities_COVID_rates_w13_to_w32.gi31 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive31))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi31)[24] <- "gstat_adaptive31"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi31) +
tm_fill(col = "gstat_adaptive31",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive31" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 31, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
8.39 Creating adaptive proximity matrix for week 32
knb32 <- knn2nb(knearneigh(coords, k=5, longlat = FALSE), row.names=row.names(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate))
knb32 <- include.self(knb32)
knb_lw32 <- nb2listw(knb32, style = 'B')plot(spdf_municipalities_COVID_rates_w13_to_w32, border="lightgrey")
plot(knb32, coordinates(spdf_municipalities_COVID_rates_w13_to_w32), pch = 19, cex = 0.6, add = TRUE, col = "red")
8.40 Gi* statistics using adaptive distance and Mapping Gi values with adaptive distance weights for week 32
fips32 <- order(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate)
gi.adaptive32 <- localG(spdf_municipalities_COVID_rates_w13_to_w32$cumul32_rate, knb_lw32)
spdf_municipalities_COVID_rates_w13_to_w32.gi32 <- cbind(spdf_municipalities_COVID_rates_w13_to_w32, as.matrix(gi.adaptive32))
names(spdf_municipalities_COVID_rates_w13_to_w32.gi32)[25] <- "gstat_adaptive32"
tm_shape(spdf_municipalities_COVID_rates_w13_to_w32.gi32) +
tm_fill(col = "gstat_adaptive32",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output## Variable(s) "gstat_adaptive32" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
At week 32, the hotspots identified are somewhat similar to the previous week. However, the main cluster of hotspots still remains at the north east region.
In conclusion, we can see that in the last 20 weeks (13 to 32 weeks), we can see that the hotspots cluster in the north east region remains consistently at the 3 to 6 range. While certain parts in the west fluctuates between -1 to 1.
We could also tell that the next potential cluster could be in the north as there are areas that are intensified for the last few weeks. Also, there are signs that the spread could also be towards the west of Mexico as well. Thus authorities should take note of these potential areas and ramp up safety measures to isolate/contain the spread as well as informing locals.