Setting the Scene
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.
The COVID-19 pandemic in Mexico is part of the ongoing worldwide pandemic of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARSCoV-2). The virus was confirmed to have reached Mexico in February 2020. However, the National Council of Science and Technology (CONACYT) reported two cases of COVID-19 in midJanuary 2020 in the states of Nayarit and Tabasco, one case per state.
As of September 26, there had been 726,431 confirmed cases of COVID-19 in Mexico and 76,243 reported deaths, although the Secretariat of Health, through the “Programa Centinela” (Spanish for “Sentinel Program”) estimated in mid July 2020 that there were more than 2,875,734 cases in Mexico, because they were considering the total number of cases confirmed as a statistical sample. (Note: This paragraph is extracted from wiki: https://en.wikipedia.org/wiki/COVID19_pandemic_in_Mexico).
In this take-home exercise, you are task 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. The specific tasks are as following tasks:
Firstly, I will start with installing all the R packages that required for this take home assignment exercise 02.
This code chunk installs R packages such as rgal, sf, spdep, tmap, tidyverse
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/THIN SU YEE/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/THIN SU YEE/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()
This code chunk uses st_read() of sf package to import the municipalities_COVID shapefile into R. The imported shapefile will be in SimpleFeaturesDataFrame object of SF.
## Reading layer `municipalities_COVID' from data source `C:\IS415_Geospatial\TakeHomeAssignment_02\Take Home\Takehome\TakeHomeAssignment_02\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_LCC
This code chunk checks the first 4 records of sf_covid dataframe contents
## Simple feature collection with 4 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 new4
## 1 01001 01 001 Aguascalientes 797010 961977 0 0 0 0
## 2 01002 01 002 Asientos 45492 50864 0 0 0 0
## 3 01003 01 003 Calvillo 54136 60760 0 0 0 0
## 4 01004 01 004 Cosío 15042 16918 0 0 0 0
## new5 new6 new7 new8 new9 new10 new11 new12 new13 new14 new15 new16 new17
## 1 0 0 0 0 0 0 1 7 30 8 10 51 60
## 2 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 8
## 4 0 0 0 0 0 0 0 0 0 0 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## 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
## actvrt8 actvrt9 actvr10 actvr11 actvr12 actvr13 actvr14 actvr15
## 1 0 0 0.1039526 0.4158104 1.663241 4.573914 5.093677 6.133203
## 2 0 0 0.0000000 0.0000000 1.966027 1.966027 1.966027 0.000000
## 3 0 0 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000
## 4 0 0 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000
## actvr16 actvr17 actvr18 actvr19 actvr20 actvr21 actvr22
## 1 9.771543 15.59289 24.324906 26.403958 30.873919 32.537160 47.29843
## 2 0.000000 0.00000 1.966027 1.966027 5.898081 13.762189 15.72822
## 3 0.000000 14.81238 14.812377 18.104016 3.291639 3.291639 1.64582
## 4 0.000000 0.00000 0.000000 5.910864 41.376049 82.752098 100.48469
## 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
## 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
## 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
## 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
## dthrt24 dthrt25 dthrt26 dthrt27 dthrt28 dthrt29 dthrt30 dthrt31
## 1 2.80672 2.598815 2.494862 1.975099 1.975099 1.975099 2.079052 1.351384
## 2 0.00000 0.000000 3.932054 0.000000 3.932054 3.932054 0.000000 0.000000
## 3 0.00000 0.000000 1.645820 0.000000 0.000000 0.000000 1.645820 0.000000
## 4 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 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...
Observation: Based on the observation, it can see that there are many other variables that are of not interest to this take home assignment. Therefore, data wrangling needs to be performed on this data before conducting the analysis.
This code chunk will check the projection of sf_covid dataframe
## Coordinate Reference System:
## User input: MEXICO_ITRF_2008_LCC
## wkt:
## PROJCRS["MEXICO_ITRF_2008_LCC",
## BASEGEOGCRS["ITRF2008",
## DATUM["International Terrestrial Reference Frame 2008",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]],
## ID["EPSG",1061]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",12,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-102,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",17.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",29.5,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",2500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
Observation: Based on the observation, the projection of the imported data is MEXICO_ITRF_2008_LCC. Therefore, this needs to be converted to EPSG 6372 which is the corresponding projection system for Mexico ITRF2008 / LCC
This code chunk will transform the projection of the imported dataframe.
## Coordinate Reference System:
## User input: EPSG:6372
## wkt:
## PROJCRS["Mexico ITRF2008 / LCC",
## BASEGEOGCRS["Mexico ITRF2008",
## DATUM["Mexico ITRF2008",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",6365]],
## CONVERSION["Mexico LCC",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",12,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-102,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",17.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",29.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",2500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Mexico"],
## BBOX[12.1,-122.19,32.72,-84.64]],
## ID["EPSG",6372]]
Observation: Based on the output, it can see that the projection system is changed to EPSG 6372
In this part, we will prepare the study data by:
* Extracting Mexico City (9), Mexico State (15) and Morelos State (17) * Selecting the CVE_ENT, Cumulative Cases from e-week 13 to 32, Population 2020 data and the city names * Calculate the Covid-19 rate (cases per 10000 population) from e-week 13 to 32 by using the formula ((Cumulative/Population)10000) Saving the calculation in new fields using mutate
This code chunk will extract and prepare the data to fit the interest of this assignment.
mexico_covid_rate <- sf_covid%>%
filter(CVE_ENT %in% c('09', '15', '17'))%>%
select(`CVE_ENT`, `CVEGEO`, `NOMGEO`, `Pop2020`, 51:70, )%>%
mutate_at(vars(contains("cumul")), funs(covid_rate=((./Pop2020)*10000)))## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Simple feature collection with 4 features and 44 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2776218 ymin: 806588.9 xmax: 2808743 ymax: 846749.5
## projected CRS: MEXICO_ITRF_2008_LCC
## CVE_ENT CVEGEO NOMGEO Pop2020 cumul13 cumul14 cumul15 cumul16
## 1 09 09002 Azcapotzalco 408441 19 46 87 156
## 2 09 09003 Coyoacán 621952 25 44 107 213
## 3 09 09004 Cuajimalpa de Morelos 199809 66 75 93 122
## 4 09 09005 Gustavo A. Madero 1176967 28 96 193 397
## cumul17 cumul18 cumul19 cumul20 cumul21 cumul22 cumul23 cumul24 cumul25
## 1 254 391 615 931 1267 1590 1939 2249 2594
## 2 375 554 831 1150 1492 1866 2255 2554 2936
## 3 157 202 260 323 414 573 714 851 991
## 4 734 1107 1653 2330 3108 3843 4674 5366 5977
## cumul26 cumul27 cumul28 cumul29 cumul30 cumul31 cumul32
## 1 2922 3238 3539 3839 4171 4442 4466
## 2 3208 3590 3903 4393 4820 5216 5246
## 3 1117 1285 1444 1615 1835 1949 1956
## 4 6502 7075 7587 8212 8875 9380 9439
## geometry cumul13_covid_rate cumul14_covid_rate
## 1 MULTIPOLYGON (((2794860 837... 0.4651835 1.1262337
## 2 MULTIPOLYGON (((2799635 820... 0.4019603 0.7074501
## 3 MULTIPOLYGON (((2787230 825... 3.3031545 3.7535847
## 4 MULTIPOLYGON (((2802176 843... 0.2378996 0.8156558
## cumul15_covid_rate cumul16_covid_rate cumul17_covid_rate cumul18_covid_rate
## 1 2.130051 3.819401 6.218768 9.572986
## 2 1.720390 3.424702 6.029404 8.907440
## 3 4.654445 6.105831 7.857504 10.109655
## 4 1.639808 3.373077 6.236369 9.405531
## cumul19_covid_rate cumul20_covid_rate cumul21_covid_rate cumul22_covid_rate
## 1 15.05725 22.79399 31.02039 38.92851
## 2 13.36116 18.49017 23.98899 30.00232
## 3 13.01243 16.16544 20.71979 28.67739
## 4 14.04457 19.79665 26.40686 32.65172
## cumul23_covid_rate cumul24_covid_rate cumul25_covid_rate cumul26_covid_rate
## 1 47.47320 55.06303 63.50978 71.54032
## 2 36.25682 41.06426 47.20622 51.57954
## 3 35.73413 42.59067 49.59737 55.90339
## 4 39.71224 45.59176 50.78307 55.24369
## cumul27_covid_rate cumul28_covid_rate cumul29_covid_rate cumul30_covid_rate
## 1 79.27706 86.64654 93.99154 102.12001
## 2 57.72150 62.75404 70.63246 77.49794
## 3 64.31142 72.26902 80.82719 91.83771
## 4 60.11214 64.46230 69.77256 75.40568
## cumul31_covid_rate cumul32_covid_rate
## 1 108.75500 109.34260
## 2 83.86499 84.34735
## 3 97.54315 97.89349
## 4 79.69637 80.19766
Observation: After extracting the 3 central mexico city and state, 176 observations derieved. Moreover, the contents of the newly created mexico_covid_rate dataframe only contains 45 columns such as, CVE_ENT to differentiate the 3 central mexico city and state, CVEGEO to uniquely identify each municipal in the 3 central mexico area, NOMGEO to show their names, POP20 to see the population area of each municipal, the cumul cases across e-week 13 to 32, the cumulative covid rate as per 10000 population for 19 weeks.
## tmap mode set to interactive viewing
This code chunk will create a choropleth map showing the location of municipalities within the study area by using tm_shape() and tm_polygon() of tmap package.
tm_shape(mexico_covid_rate) +
tm_fill("CVE_ENT",
popup.vars=c("CVE_ENT"="CVE_ENT",
"Name"="NOMGEO",
"CVEGEO"="CVEGEO",
"Population" = "Pop2020"))+
tm_text("CVEGEO", size=1)## tmap mode set to plotting
This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week 13.
e_week13<- tm_shape(mexico_covid_rate) +
tm_fill("cumul13_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week13",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week 14.
e_week14 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul14_covid_rate", n = 5, style="jenks", palette = "Reds")+
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week14",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week15.
e_week15 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul15_covid_rate", n = 5, style="jenks", palette = "Reds")+
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week15",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week16.
e_week16 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul16_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week16",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week17.
e_week17 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul17_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week17",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week18.
e_week18 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul18_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week18",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week19.
e_week19 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul19_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week19",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week20.
e_week20 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul20_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week20",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week21.
e_week21 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul21_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week21",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week22.
e_week22 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul22_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week22",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week23.
e_week23 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul23_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week23",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week24.
e_week24 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul24_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week24",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week25.
e_week25 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul25_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week25",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week26.
e_week26 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul26_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week26",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week27.
e_week27 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul27_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week27",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week28.
e_week28 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul28_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week28",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week29.
e_week29 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul29_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week29",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week30.
e_week30 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul30_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week30",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week31.
e_week31 <- tm_shape(mexico_covid_rate) +
tm_fill("cumul31_covid_rate", n = 5, style="jenks", palette = "Reds") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = "e-Week31",
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE) This code chunk will use the thematic mapping techniques of tmap package to create a map that show the area that are affected by the Covid19 virus in e-week32.
This code chunk will show the spatio-temporal distribution of covid 19 rate in mexico cities across 20 weeks.
Observation: Based on the spatial temporal analysis of Covid-19 rate in central mexico at municipal level, it can see that the Covid-19 rates are increasing over 20 weeks. In the e-week 13, the covid-19 rate range from 0 to 3.3303 which can be considered as very low and there are many municipals which are yet to be affected by this contiguous virus. As of e-week 13, Mexico city (9) is seemed to be the mostly affected with virus while the spread in mexico state (15) and Morelos State (17) is still relatively low. Within the 4 weeks from e-week 13 to e-week 16, the covid 19 rate increased by double to 6.269 and the virus has seemed to have slowly spreading to the neighboring municipalities that are sharing boundaries with the mexico city as the covid-19 rates in the surrounding areas are increasing. Within e-week 17 to e-week 20, the covid-19 rate in the central mexico range from 0 to 30.14 and many area of mexico city(9) and Morelos State(17) has largely affected. Within e-week 21 to e-week 24, the covid-19 rate range from 0 to 84.95 and multiple area has the very high covid rate as illustrated by the darker color on the map. Within e-week 25 and e-week28, the rate range from 0.53 to 133.5. This number explained that the virus has spread to every municipal in the central mexico since there is no more rate of 0. The map also shows that the spread of virus is mostly clustered around mexico city(9).Within e-week 29 to e-week 32, the covid-19 rate range from 1.1 to 170.4. Overall, we can see from the map that the infection of covid-19 is quite bad in mexico city(9) across 20 weeks and the regions in the mexico city has the most highest covid-19 rate throughout the 20 weeks.
Since the study aims to find out the spread of covid cases across the central mexico, Queen contiguity weight matrices is used to build a neighbor list based on regions with contiguous boundaries. The reason is because it is highly possible for contiguous boundaries to be accessed by another municipal that shared boundary. Therefore, it is more likely for the virus to be spread quckily.
This code chunk will compute Queen contiguity weight matrix using poly2nd() of spdep package.
## 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 links
Observation: Based on the summary report, the study area includes 60 municipals from Central Mexico. The summary report shows that there are 3 least connected area with 1 neighbor whereas there are 1 most connected area with 14 neighbors. In the case of this study, those most connected municipalities can be the most vulnerable in spreading the virus easily compared to least connected area because most connected area have more spatial interaction as they share the boundaries with more neighbors. In another word, they are most connected and sharing the boundaries with 14 other neighbor. Therefore, it is a super spreader and has the higher chance of spreading virus faster compared to the the least connected area. The most connected municipalities is 117 in the wm_q list.
This code chunk checks the name and CVEGEO number of the 1 most connected municipal.
## [1] "15"
## [1] "15101"
## [1] "Tianguistenco"
Observation: Based on the findings, “Tianguistenco” is the most connected municipal with 14 neighbors. Based on its CVE_ENT number, it can also see that it belongs to CVE_ENT 15, Mexico State, which is connecting to many other municipalities. Therefore, once the virus reached to the area, there is a higher chance that the virus will spread to the nearby regions as well since they are located very nearby and connected to each other.
This code chunk will find out the area that wm_q[117] is connected to.
## [[1]]
## [1] 11 22 28 34 35 43 59 65 67 70 78 79 89 114
This code chunk will find out the name and cvegeo of neighbouring municipalities of the queen matrix’s wm_q[117].
## [1] "Tlalpan" "Almoloya del Río" "Atizapán"
## [4] "Calimaya" "Capulhuac" "Chapultepec"
## [7] "Xalatlaco" "Joquicingo" "Lerma"
## [10] "Metepec" "Ocoyoacac" "Ocoyoacac"
## [13] "San Antonio la Isla" "Texcalyacac"
## [1] "09002" "15006" "15012" "15018" "15019" "15027" "15043" "15049" "15051"
## [10] "15054" "15062" "15062" "15073" "15098"
Observation: Based on the output, the neighbouring municipalities are “Tlalpan”, “Almoloya del Río”, “Atizapán”, “Calimaya”, Capulhua“, Chapultepec”, “Xalatlaco”, “Joquicingo”, “Lerma”, “Metepec”, “Ocoyoacac”, “Ocoyoacac”, “San Antonio la Isla” and “Texcalyacac”. And they belong to Mexico City (9) and Mexico State (15). Therefore, once the virus is reached to the Tianguistenco, the most connected region, it is likely for the virus to be spread to these neighboring regions in mexico city and mexico state.
This code chunk will find out the names and it details of the least connected regions which are wmq[77,95, 121]
## [1] "Nopaltepec"
## [1] "15"
## [1] "15061"
## [[1]]
## [1] 32
## [1] "Soyaniquilpan de Juárez"
## [1] "15"
## [1] "15079"
## [[1]]
## [1] 61
## [1] "Tlatlaya"
## [1] "15"
## [1] "15105"
## [[1]]
## [1] 24
Observation: Based on the output, we can see that “Nopaltepec”, “Soyaniquilpan de Juárez”, and “Tlatlaya” are the least connected regions and they belong to the the city code of 15 which is the mexico state. When we observe their location on the map, we can also see that they are located at the each corner of the mexico state. That could be the reason why they are least connected. It also meant that these regions have the lesser chance of spreading or getting infected compared to those very connected regions.
This code chunk will plot the Queen Contiguity Based Neighbours Map using the plot() function.
coords <- st_centroid(mexico_covid_rate$geometry)
plot(mexico_covid_rate$geometry, border="grey", main="Queen Contiguity-Based Neighbours", asp=0.8)
plot(wm_q, st_coordinates(coords), pch = 5, cex = 0.6, add = TRUE, col= "hotpink2" )In this assignment, adaptive distance weight matrix will be used to control the numbers of neighbor by making each municipal have the exact number of neighbors. Since it is said by researches that choosing smaller value for K can create noises that can have influence on the result, while the increasing value of K can provide smoother boundary. Therefore, we will use ‘k=sqrt(N)’ to determine the k value in this case.
This code chunk calculate the value of k to use in adaptive distance weight matrix
## [1] 13.2665
Observation: By using the sqrt() function to calculate the value of k to use in adaptive distance weight matrix, it is found that 13 should be used as k value.
This code chunk will compute adaptive distance weight matrix.
## Warning in knearneigh(coords, k = 13, longlat = FALSE): dnearneigh: longlat
## argument overrides object
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 2288
## Percentage nonzero weights: 7.386364
## Average number of links: 13
## Non-symmetric neighbours list
Observation: Based on the output, it can see that every municipal have the same number of neighbors. By using adaptive distance weight matrix, it ensure that every municipal has a neighbor and avoid the nonzero neighbor problem. Moreover, it also ensures that every region has the equal number of neighbors so that when we analyze the auto spatial correlation, the result will be more accurate.
This code chunk will display the content of the matrix by using str()
## List of 176
## $ : int [1:13] 4 5 9 13 14 15 16 29 36 49 ...
## $ : int [1:13] 3 5 6 7 9 10 11 12 13 14 ...
## $ : int [1:13] 2 7 9 11 13 15 35 53 59 67 ...
## $ : int [1:13] 1 5 13 14 15 16 27 29 36 49 ...
## $ : int [1:13] 1 2 4 6 10 12 13 14 15 16 ...
## $ : int [1:13] 2 4 5 10 12 13 14 15 16 47 ...
## $ : int [1:13] 2 3 9 11 12 13 14 15 35 53 ...
## $ : int [1:13] 10 11 12 38 66 99 105 138 150 161 ...
## $ : int [1:13] 1 2 3 5 7 11 13 14 15 16 ...
## $ : int [1:13] 2 5 6 8 12 38 41 45 47 74 ...
## $ : int [1:13] 2 3 6 7 8 9 10 12 13 59 ...
## $ : int [1:13] 2 5 6 7 8 9 10 11 13 16 ...
## $ : int [1:13] 1 2 3 4 5 6 7 9 12 14 ...
## $ : int [1:13] 1 2 4 5 6 9 13 15 16 49 ...
## $ : int [1:13] 1 2 3 4 5 9 13 14 16 29 ...
## $ : int [1:13] 1 2 4 5 6 13 14 15 47 49 ...
## $ : int [1:13] 19 30 42 58 61 64 72 80 87 90 ...
## $ : int [1:13] 27 36 44 46 49 81 85 91 97 108 ...
## $ : int [1:13] 17 30 42 61 64 72 80 87 90 95 ...
## $ : int [1:13] 37 56 93 96 98 102 104 106 113 123 ...
## $ : int [1:13] 23 34 58 63 70 71 83 90 103 122 ...
## $ : int [1:13] 28 35 43 59 65 70 71 78 88 89 ...
## $ : int [1:13] 21 48 57 93 94 102 113 122 126 127 ...
## $ : int [1:13] 20 82 93 94 96 98 102 113 121 126 ...
## $ : int [1:13] 31 33 38 41 50 55 66 84 99 105 ...
## $ : int [1:13] 39 40 51 52 60 69 75 100 107 111 ...
## $ : int [1:13] 4 18 36 44 45 46 47 49 74 85 ...
## $ : int [1:13] 22 35 43 59 65 70 71 78 88 89 ...
## $ : int [1:13] 1 4 14 15 36 54 62 73 76 111 ...
## $ : int [1:13] 17 19 42 58 61 63 64 72 80 90 ...
## $ : int [1:13] 25 33 50 66 84 105 110 119 143 157 ...
## $ : int [1:13] 18 44 46 77 81 85 91 97 100 108 ...
## $ : int [1:13] 25 31 38 41 50 66 84 99 105 110 ...
## $ : int [1:13] 22 28 43 65 70 71 88 89 92 106 ...
## $ : int [1:13] 22 28 43 59 67 70 71 78 88 89 ...
## $ : int [1:13] 4 40 49 60 69 75 97 107 120 124 ...
## $ : int [1:13] 20 34 56 93 102 104 106 113 123 129 ...
## $ : int [1:13] 8 10 25 33 41 45 55 66 86 99 ...
## $ : int [1:13] 40 51 60 69 75 107 111 112 124 125 ...
## $ : int [1:13] 36 39 51 60 69 75 107 111 124 125 ...
## $ : int [1:13] 10 25 33 38 45 47 55 66 86 99 ...
## $ : int [1:13] 30 39 51 54 61 63 64 72 76 95 ...
## $ : int [1:13] 22 28 34 35 65 70 71 88 89 92 ...
## $ : int [1:13] 18 27 45 46 49 81 85 91 97 108 ...
## $ : int [1:13] 6 10 27 41 44 46 47 55 74 85 ...
## $ : int [1:13] 18 27 44 45 47 49 85 91 97 108 ...
## $ : int [1:13] 5 6 10 16 27 44 45 46 74 85 ...
## $ : int [1:13] 21 23 57 82 93 94 102 126 127 130 ...
## $ : int [1:13] 4 16 18 27 36 44 46 97 116 120 ...
## $ : int [1:13] 25 31 33 66 84 110 143 157 163 168 ...
## $ : int [1:13] 26 39 40 60 69 75 107 111 112 124 ...
## $ : int [1:13] 26 39 40 51 60 69 75 97 100 107 ...
## $ : int [1:13] 1 2 3 7 9 13 14 15 62 67 ...
## $ : int [1:13] 29 53 62 63 67 73 76 83 103 111 ...
## $ : int [1:13] 10 25 33 38 41 45 47 86 99 105 ...
## $ : int [1:13] 20 37 65 68 96 104 106 113 123 129 ...
## $ : int [1:13] 23 48 82 93 94 98 102 126 127 130 ...
## $ : int [1:13] 21 30 63 64 72 80 83 90 103 118 ...
## $ : int [1:13] 3 7 11 22 28 35 65 78 79 89 ...
## $ : int [1:13] 36 39 40 49 69 75 97 107 124 125 ...
## $ : int [1:13] 17 19 30 39 42 51 64 72 87 95 ...
## $ : int [1:13] 1 3 15 29 53 54 67 73 76 83 ...
## $ : int [1:13] 21 42 54 58 62 64 67 72 76 83 ...
## $ : int [1:13] 17 21 30 42 58 63 72 80 90 101 ...
## $ : int [1:13] 22 28 34 35 43 71 79 88 89 104 ...
## $ : int [1:13] 8 25 33 38 41 84 99 105 110 143 ...
## $ : int [1:13] 3 28 35 43 53 62 70 71 78 83 ...
## $ : int [1:13] 56 65 79 104 123 129 135 146 148 156 ...
## $ : int [1:13] 36 39 40 51 60 75 97 107 124 125 ...
## $ : int [1:13] 22 28 34 35 43 67 71 88 89 92 ...
## $ : int [1:13] 22 28 34 35 43 70 88 89 92 106 ...
## $ : int [1:13] 17 30 42 54 58 61 63 64 76 83 ...
## $ : int [1:13] 1 3 4 9 13 14 15 29 53 54 ...
## $ : int [1:13] 2 4 5 6 10 13 14 16 27 45 ...
## $ : int [1:13] 36 39 40 49 60 69 97 107 124 125 ...
## $ : int [1:13] 29 39 40 54 62 73 83 111 120 125 ...
## $ : int [1:13] 18 32 44 46 52 81 85 91 97 100 ...
## $ : int [1:13] 3 7 22 28 35 43 53 59 67 89 ...
## $ : int [1:13] 22 28 59 65 68 88 104 114 117 148 ...
## $ : int [1:13] 17 19 21 30 58 64 72 90 101 118 ...
## $ : int [1:13] 18 27 32 44 46 77 85 91 100 108 ...
## $ : int [1:13] 23 24 48 57 93 94 98 102 113 126 ...
## $ : int [1:13] 53 54 62 63 67 70 73 76 78 92 ...
## $ : int [1:13] 25 31 33 50 66 105 110 143 157 163 ...
## $ : int [1:13] 18 27 44 45 46 47 49 81 91 108 ...
## $ : int [1:13] 5 6 10 12 16 38 41 45 47 55 ...
## $ : int [1:13] 17 19 30 42 61 64 72 80 90 95 ...
## $ : int [1:13] 22 28 34 35 43 65 70 71 89 92 ...
## $ : int [1:13] 22 28 34 35 43 65 70 71 88 92 ...
## $ : int [1:13] 17 21 30 58 63 64 72 80 101 118 ...
## $ : int [1:13] 18 27 32 44 46 77 81 85 97 100 ...
## $ : int [1:13] 22 28 34 35 43 67 70 71 78 88 ...
## $ : int [1:13] 20 23 37 82 94 96 98 102 113 126 ...
## $ : int [1:13] 23 48 57 82 93 98 102 113 126 127 ...
## $ : int [1:13] 17 19 26 39 42 51 61 72 87 111 ...
## $ : int [1:13] 20 24 37 56 93 98 102 113 121 123 ...
## $ : int [1:13] 18 36 49 60 69 75 91 100 108 116 ...
## $ : int [1:13] 20 24 57 82 93 94 96 102 113 121 ...
## $ : int [1:13] 8 10 25 33 38 41 66 84 86 105 ...
## [list output truncated]
## - attr(*, "region.id")= chr [1:176] "1" "2" "3" "4" ...
## - attr(*, "call")= language knearneigh(x = coords, k = 13, longlat = FALSE)
## - attr(*, "sym")= logi FALSE
## - attr(*, "type")= chr "knn"
## - attr(*, "knn-k")= num 13
## - attr(*, "class")= chr "nb"
Observation: Based on the output, we can see that every municipal has the exact number of neighbor and there is no region without a neighbor. Therefore, we are ready to find out if there is spatial autocorrelation between each region and its neighbors.
plot(mexico_covid_rate$geometry, border="grey", main="K nearest neighbours, k=13", asp=0.8)
plot(wm_knn, coords, pch = 5, cex = 0.6, add = TRUE, col= "hotpink2")Observation: Based on the K nearest neighbor map, we can notice that every municipal is connected together and there is no single municipal that does not have a link to another.
Both the Global Moran’s I and Local’s Moran’s I are used to detect the spatial autocorrelation. However, Local Moran I is used to gain a deeper understanding on which objects are similar or different to the objects in their neighborhood area based on the Moran’s I statistic. It identifies concentrations of high values, concentrations of low values, and spatial outliers in the neighborhoods. Significant and negative values will be considered as outliers if the location is associated with relatively low values in the surrounding location whereas significant and positive values will be considered as cluster if the location is associated with relatively high values of the surrounding locations.
In this step, we will perform row-standardization by assigning weight to each neighboring polygon. Row standardization is used to create proportional weights in cases where features have an unequal number of neighbors. Row standardization involves dividing each neighbor weight for a feature by the sum of all neighbor weights for that feature, and is recommended whenever the distribution of your features is potentially biased due to sampling design or an imposed aggregation scheme.
Spatial weights are often row standardized with the use of binary weighting strategies. Therefore, this study will also be using Binary weighting strategies to assign a weight value of 1 to neighboring features and 0 to non-neighbors polygon features with the use of ‘style=“B”’.
This code chunk will be performing row-standardised weights matrix
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 962
## Percentage nonzero weights: 3.10563
## Average number of links: 5.465909
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 176 30976 962 1924 23920
This code chunk will see the weight of the first municipal’s four neighbours
## [[1]]
## [1] 1 1 1 1 1
Observation: Based on the output, the neighbors of the first municipal in the list are given the weight value of 1 to show that they are neighbors.
In this Local Moran’s I analysis, the following steps will be performed: * 1. Computing local Moran’s I : To compute local Moran’s I, the localmoran() function of spdep will be used. It computes Ii values, given a set of zi values and a listw object providing neighbour weighting information for the polygon associated with the zi values.
This code chunk will be used to create a function to compute local Moran’s I of covid 19 rates at central mexico’s municipalities level, append the local moran I statistic into initial SpatialPolygonDataFrame and then create the map content for covid-19 rate in municipalities for each e-week using local Moran’s I values.
localMImap_covid_rate <- function(weekly_covid_rate, i)
{
fips <- order(mexico_covid_rate$NOMGEO)
localMI <- localmoran(weekly_covid_rate, rswm_q)
mexico_covid.localMI <- cbind(mexico_covid_rate,localMI)
localMI.map <- tm_shape(mexico_covid.localMI) +
tm_fill(col = "Ii",
style = "pretty",
title="Local Moran Statistics") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = sprintf("e-week %d",i),
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE)
}This code chunks will create a function to compute local Moran’s I of covid 19 rates at central mexico’s municipalities level and create choropleth map using local moran’s I p-values.
pvalue_covid_rate <- function(weekly_covid_rate, i)
{
fips <- order(mexico_covid_rate$NOMGEO)
localMI <- localmoran(weekly_covid_rate, rswm_q)
mexico_covid.localMI <- cbind(mexico_covid_rate,localMI)
pvalue.map <- tm_shape(mexico_covid.localMI) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title="Local Moran I's p-values") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = sprintf("e-week %d",i),
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE)
}This code chunk will call the localMImap_covid_rate and pvalue_covid_rate functions to create choropleth maps that shows temporal spatial pattens across 20 e-weeks.
week_13_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul13_covid_rate, 13)
week_14_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul14_covid_rate, 14)
week_15_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul15_covid_rate, 15)
week_16_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul16_covid_rate, 16)
week_17_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul17_covid_rate, 17)
week_18_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul18_covid_rate, 18)
week_19_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul19_covid_rate, 19)
week_20_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul20_covid_rate, 20)
week_21_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul21_covid_rate, 21)
week_22_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul22_covid_rate, 22)
week_23_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul23_covid_rate, 23)
week_24_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul24_covid_rate, 24)
week_25_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul25_covid_rate, 25)
week_26_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul26_covid_rate, 26)
week_27_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul27_covid_rate, 27)
week_28_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul28_covid_rate, 28)
week_29_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul29_covid_rate, 29)
week_30_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul30_covid_rate, 30)
week_31_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul31_covid_rate, 31)
week_32_localMI <- localMImap_covid_rate(mexico_covid_rate$cumul32_covid_rate, 32)
week_13_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul13_covid_rate, 13)
week_14_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul14_covid_rate, 14)
week_15_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul15_covid_rate, 15)
week_16_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul16_covid_rate, 16)
week_17_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul17_covid_rate, 17)
week_18_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul18_covid_rate, 18)
week_19_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul19_covid_rate, 19)
week_20_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul20_covid_rate, 20)
week_21_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul21_covid_rate, 21)
week_22_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul22_covid_rate, 22)
week_23_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul23_covid_rate, 23)
week_24_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul24_covid_rate, 24)
week_25_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul25_covid_rate, 25)
week_26_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul26_covid_rate, 26)
week_27_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul27_covid_rate, 27)
week_28_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul28_covid_rate, 28)
week_29_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul29_covid_rate, 29)
week_30_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul30_covid_rate, 30)
week_31_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul31_covid_rate, 31)
week_32_pvalue <- pvalue_covid_rate(mexico_covid_rate$cumul32_covid_rate, 32)
tmap_arrange(week_13_localMI, week_13_pvalue, week_14_localMI, week_14_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_15_localMI, week_15_pvalue, week_16_localMI, week_16_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_17_localMI, week_17_pvalue, week_18_localMI, week_18_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_19_localMI, week_19_pvalue, week_20_localMI, week_20_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_21_localMI, week_21_pvalue, week_22_localMI, week_22_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_23_localMI, week_23_pvalue, week_24_localMI, week_24_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_25_localMI, week_25_pvalue, week_26_localMI, week_26_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_27_localMI, week_27_pvalue, week_28_localMI, week_28_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_29_localMI, week_29_pvalue, week_30_localMI, week_30_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
tmap_arrange(week_31_localMI, week_31_pvalue, week_32_localMI, week_32_pvalue, asp=0, ncol=4, sync= TRUE)## 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.
## 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.
The P values from week 13 to 32.
Define Hypotheses and critical values
Ho = The distribution of COVID-19 are randomly distributed.
H1= The distribution of COVID-19 are not randomly distributed.
Confident level= 99%
critical value 0.001
From the temporal spatial patterns of Local Moran’s I statistics map plotted for e-week 13 to e-week 32, we can see that there are significant areas which has the positive values. Based on the interpretation of local Moran I, the positive statistic values indicate that there is clustering as the location is surrounded by similar high values neighbors. However, there are also some municipals with the significant negative values. Based on the interpretation of local Moran I, the negative statistics values indicate that those areas are outlier as they have the dissimilar low values and surrounded by relatively low values neighbors.
Based on that, we can derive that Mexico city has more positive spatial autocorrelation as it has the highest statistic values across 20 weeks and the covid-19 rates are similar to surrounding neighborhood.
Moreover, we can also see from the Local Moran’s I p-values across week 13 to week 32 that the municipalities in the mexico city has the p-value less than 0.001 which is the less than the confidence interval of 99.9%. Therefore, we reject the null hypothesis that the Covid-19 rate are randomly distributed the since there is a significant evidence that there is a cluster of covid 19 in specific areas across 20 weeks and the covid-19 is not randomly distributed. For the rest of the area, they are of more than 0.001. Therefore, we cannot reject the null hypothesis. Thus, conclude that those area with more than p-value of 0.001 is considered as randomly distributed and they are outliers.
Beside detecting cluster and outliers, localised spatial statistics like Getis and Ord’s G-statistics can be also used to detect hot spot and/or cold spot areas. The term ‘hot spot’ has been used generically across disciplines to describe a region or value that is higher relative to its surroundings (Lepers et al 2005, Aben et al 2012, Isobe et al 2015).
Therefore, the alternative way to detect spatial anomalies is to use Getis and Ord’s G-statistics. This technique look at neighbors within a defined proximity to identify features with either high (i.e. hot spots) or low values (cold spots) are clustered spatially.
In this Local Getis Ord’s G* Analysis analysis, the following steps will be performed:
In this step, a fixed number of neighbors will be used to derive the proximity matrix instead of using a fixed distance. This technique is called adaptive distance method.
This code chunk will compute adaptive distance weight matrix using knn2nb() of spdep package.
## Warning in knearneigh(coords, k = 13, longlat = FALSE): dnearneigh: longlat
## argument overrides object
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 2288
## Percentage nonzero weights: 7.386364
## Average number of links: 13
## Non-symmetric neighbours list
## Link number distribution:
##
## 13
## 176
## 176 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 13 links
## 176 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 13 links
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 176 30976 2288 4156 121174
This code chunk will use to plot the K nearest neighbor map.
plot(mexico_covid_rate$geometry, border="grey", main="K nearest neighbours k=13", asp=0.8)
plot(knb, coords, pch = 19, cex = 0.6, add = TRUE, col = "hotpink2")In this Local Gi analysis, the following steps will be performed: * 1. Computing Gi analysis: the localmoran() function of spdep will be used to compute Gi analysis. It will return the Ii values, given a set of zi values and a list object providing neighbor weighting information for the polygon associated with the zi values.
This code chunk will create a function to compute local Gi statistics on the covid-19 rate across 20 weeks in central mexico’s municipalities and create the map contents for spatial temporal analysis.
gi_covid_rate <- function(weekly_covid_rate, i)
{
fips <- order(mexico_covid_rate$NOMGEO)
gi.adaptive <- localG(weekly_covid_rate, knb_lw)
mexico_covid_rate.gi <- cbind(mexico_covid_rate, as.matrix(gi.adaptive))
names(mexico_covid_rate.gi)[45] <- "gstat_adaptive"
tm_shape(mexico_covid_rate.gi) +
tm_fill(col = "gstat_adaptive",
style = "pretty",
palette = "-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)+
tm_layout(main.title = sprintf("e-week %d",i),
main.title.position = "center",
legend.show = TRUE,
legend.text.size = 0.4,
frame = TRUE)
}This code chunk will call the gi_covid_rate function to create choropleth maps that shows temporal spatial pattens of across 20 e-weeks using local Gi analysis.
week_13_gi <- gi_covid_rate(mexico_covid_rate$cumul13_covid_rate, 13)
week_14_gi <- gi_covid_rate(mexico_covid_rate$cumul14_covid_rate, 14)
week_15_gi <- gi_covid_rate(mexico_covid_rate$cumul15_covid_rate, 15)
week_16_gi <- gi_covid_rate(mexico_covid_rate$cumul16_covid_rate, 16)
week_17_gi <- gi_covid_rate(mexico_covid_rate$cumul17_covid_rate, 17)
week_18_gi <- gi_covid_rate(mexico_covid_rate$cumul18_covid_rate, 18)
week_19_gi <- gi_covid_rate(mexico_covid_rate$cumul19_covid_rate, 19)
week_20_gi <- gi_covid_rate(mexico_covid_rate$cumul20_covid_rate, 20)
week_21_gi <- gi_covid_rate(mexico_covid_rate$cumul21_covid_rate, 21)
week_22_gi <- gi_covid_rate(mexico_covid_rate$cumul22_covid_rate, 22)
week_23_gi <- gi_covid_rate(mexico_covid_rate$cumul23_covid_rate, 23)
week_24_gi <- gi_covid_rate(mexico_covid_rate$cumul24_covid_rate, 24)
week_25_gi <- gi_covid_rate(mexico_covid_rate$cumul25_covid_rate, 25)
week_26_gi <- gi_covid_rate(mexico_covid_rate$cumul26_covid_rate, 26)
week_27_gi <- gi_covid_rate(mexico_covid_rate$cumul27_covid_rate, 27)
week_28_gi <- gi_covid_rate(mexico_covid_rate$cumul28_covid_rate, 28)
week_29_gi <- gi_covid_rate(mexico_covid_rate$cumul29_covid_rate, 29)
week_30_gi <- gi_covid_rate(mexico_covid_rate$cumul30_covid_rate, 30)
week_31_gi <- gi_covid_rate(mexico_covid_rate$cumul31_covid_rate, 31)
week_32_gi <- gi_covid_rate(mexico_covid_rate$cumul32_covid_rate, 32)
tmap_arrange(week_13_gi, week_14_gi, week_15_gi, week_16_gi, asp = 0, ncol=4) ## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "gstat_adaptive" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
Based on the observation from temporal analysis patterns of choropleth map, we can see that the the darkness of the colors on the map indicate the intensity of the Local Gi values. In this case, municipalities that shaded in the light red to dark red are the hotspot areas whereas the municipalities that shaded in the light blue to dark blues are the cold spot areas. Across 20 weeks, we can see that mexico city has the hotspot areas where the spread of covid-19 started. Moreover, the hotspot area remains the same across 20 weeks. Based on the map, we can also see that areas outside of the hotspot are the cold spot areas. From this, we can conclude that the hotspot and cold spot areas are distributed equally across 20 weeks.