1. Overview

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

1.1. The Tasks

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:

    1. Using appropriate geospatial data wrangling methods,
    • extract municipalities located with the study area, and
    • calculate the COVID-19 rate (i.e. cases per 10000 population) from e-week 13 until e-week 32.
    1. 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. SMU Classification: Restricted
    1. 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.
    1. Perform local Getis-Ord Gi* analysis and display the results using appropriate thematic mapping techniques. Describe the spatio-temporal patterns reveal by the maps.

2. Install R Packages

Firstly, I will start with installing all the R packages that required for this take home assignment exercise 02.

  • rgal, sf, spdep, tmap, tidyverse

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()

3. Loading the data

3.1. Importing the data

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.

sf_covid <- st_read(dsn = "data/geospatial", layer="municipalities_COVID")
## 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

3.2. Checking the content

This code chunk checks the first 4 records of sf_covid dataframe contents

head(sf_covid, n=4)
## 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.

3.3. Checking the projection

This code chunk will check the projection of sf_covid dataframe

st_crs(sf_covid)
## 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

3.4. Checking the projection

This code chunk will transform the projection of the imported dataframe.

sf_covid_6372 = st_transform(sf_covid, 6372)
st_crs(sf_covid_6372)
## 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

4. Data Wrangling Methods

4.1. Extract and Load Data

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.
head(mexico_covid_rate, n = 4)
## 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.

4.2. Visualing municipalities located within the study area of Central Mexico

tmap_mode('view')
## 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)

Observation: Based on the map, it can see that the municipalities are grouped into CVE_ENT such as Mexico City(9), Mexico State(15), and Morelos State(17) based on their locations. The municipalities within each CVE_ENT are very close or mostly connected to each other.

tmap_mode('plot')
## tmap mode set to plotting

5. Spatio-Temporal Distribution

5.1. Creating map contents 20 weeks

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.

e_week32 <- tm_shape(mexico_covid_rate) + 
  tm_fill("cumul32_covid_rate", n = 5, style="jenks", palette = "Reds") +
  tm_borders(alpha = 0.5)+
  tm_layout(main.title = "e-Week32",
            main.title.position = "center",
            legend.show =  TRUE,
            legend.text.size = 0.4,
            frame = TRUE) 

5.2. Visualizing the spatio-temporal distribution of covid 19 rate in central mexico

This code chunk will show the spatio-temporal distribution of covid 19 rate in mexico cities across 20 weeks.

tmap_arrange(e_week13, e_week14, e_week15, e_week16, asp = 0, ncol=4) 

tmap_arrange(e_week17, e_week18, e_week19, e_week20, asp = 0, ncol=4)

tmap_arrange(e_week21, e_week22, e_week23, e_week24, asp = 0, ncol=4)

tmap_arrange(e_week25, e_week26, e_week27, e_week28, asp = 0, ncol=4)

tmap_arrange(e_week29, e_week30, e_week31, e_week32, asp = 0, ncol=4)

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.

6. Modelling Spatial Neighbours

6.1. Creating Queen contiguity based neighbors

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.

wm_q <- poly2nb(mexico_covid_rate, queen=TRUE)
summary(wm_q)
## Neighbour list object:
## Number of regions: 176 
## Number of nonzero links: 962 
## Percentage nonzero weights: 3.10563 
## Average number of links: 5.465909 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 14 
##  3  3 18 41 31 31 22 13  9  3  1  1 
## 3 least connected regions:
## 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.

6.2. Nearest neighbours analysis based on Queen matrices

6.2.1. The most connected region in central mexico

This code chunk checks the name and CVEGEO number of the 1 most connected municipal.

mexico_covid_rate$CVE_ENT[117]
## [1] "15"
mexico_covid_rate$CVEGEO[117]
## [1] "15101"
mexico_covid_rate$NOMGEO[117]
## [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.

wm_q[117]
## [[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].

mexico_covid_rate$NOMGEO[c(11, 22, 28, 34, 35, 43, 59, 65, 67, 70, 78, 78, 89, 114)]
##  [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"
mexico_covid_rate$CVEGEO[c(1, 22, 28, 34, 35, 43, 59, 65, 67, 70, 78, 78, 89, 114)]
##  [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.

6.2.2. The least connected regions in central mexico

This code chunk will find out the names and it details of the least connected regions which are wmq[77,95, 121]

mexico_covid_rate$NOMGEO[77]
## [1] "Nopaltepec"
mexico_covid_rate$CVE_ENT[77]
## [1] "15"
mexico_covid_rate$CVEGEO[77]
## [1] "15061"
wm_q[77]
## [[1]]
## [1] 32
mexico_covid_rate$NOMGEO[95]
## [1] "Soyaniquilpan de Juárez"
mexico_covid_rate$CVE_ENT[95]
## [1] "15"
mexico_covid_rate$CVEGEO[95]
## [1] "15079"
wm_q[95]
## [[1]]
## [1] 61
mexico_covid_rate$NOMGEO[121]
## [1] "Tlatlaya"
mexico_covid_rate$CVE_ENT[121]
## [1] "15"
mexico_covid_rate$CVEGEO[121]
## [1] "15105"
wm_q[121]
## [[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.

6.3. Plotting Queen Contiguity Based Neighbours

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" )

6.4. Computing adaptive distance weight matrix

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.

6.4.1. Determining the k value

This code chunk calculate the value of k to use in adaptive distance weight matrix

k <- sqrt(176)
k
## [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.

6.4.2. Compuging the adaptive distance weight matrix

This code chunk will compute adaptive distance weight matrix.

wm_knn <- knn2nb(knearneigh(coords, k=13, longlat = FALSE))
## Warning in knearneigh(coords, k = 13, longlat = FALSE): dnearneigh: longlat
## argument overrides object
wm_knn
## 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()

str(wm_knn)
## 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.

6.5. Plotting adaptive distance based neighbors map

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.

7. Performing Local Moran’s I analysis

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.

7.1. Row-standardised weights matrix

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

rswm_q <- nb2listw(wm_q, style="B", 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: 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

rswm_q$weights[1]
## [[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.

7.2. Computing Local Moran’s I

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.

    1. Appending the local Moran’s I: Before mapping the local Moran’s I map, it is wise to append the local Moran’s I dataframe (i.e. localMI) onto mexico_covid_rate SpatialPolygonDataFrame. The output SpatialPolygonDataFrame will be called as mexico_covid.localMI
    1. Mapping local Moran’s I values: Using choropleth mapping functions of tmap package to plot the local Moran’s I values.
    1. Mapping local Moran’s I p-values: Using choropleth mapping function of tmap package to plot the Moran’s I p-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) 
}

7.3. Visualizating Local Moran’s I Spatial Temporal Analysis

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

7.4. Geocomunication on Moran I’s Analysis

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.

8. Performing Getis Ord’s G* Analysis

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:

    1. Deriving spatial weight matrix
    1. Computing Gi statistics
    1. Visualizing Gi statistics

8.1. Deriving Spatial weight Matrix

8.1.1. Creating adaptive proximity matrix

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.

knb <- knn2nb(knearneigh(coords, k=13, longlat = FALSE))
## Warning in knearneigh(coords, k = 13, longlat = FALSE): dnearneigh: longlat
## argument overrides object
knb_lw <- nb2listw(knb, style = 'B')
summary(knb_lw)
## 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")

8.2. Computing Gi statistics

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.

    1. Appending the Local Gi: cbind() will be used to append the local G value to the existing SpatialPolygonDataFrame, mexico_covid_rate and the output SpatialPolygonDataFrame will be called as mexico_covid.rate.gi.
    1. Mapping local gi values: Using choropleth mapping functions of tmap package to plot the local gi’s 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) 
}

8.3. Visualizing Gi statistics

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.

tmap_arrange(week_17_gi, week_18_gi, week_19_gi, week_20_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.

tmap_arrange(week_21_gi, week_22_gi, week_23_gi, week_24_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.

tmap_arrange(week_25_gi, week_26_gi, week_27_gi, week_28_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.

tmap_arrange(week_29_gi, week_30_gi, week_31_gi, week_32_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.

8.4. Geocommunication on Local G’s Analysis

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.