Background

Since late December 2019, an outbreak of a novel coronavirus disease (COVID-19; previously known as 2019-nCoV) was reported in Wuhan, China, which has subsequently affected 210 countries worldwide. In general, COVID-19 is an acute resolved disease but it can also be deadly, with a 2% case fatality rate.

In this take home exercise the aim is 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.

Packages Installation

Install the following packages for this study.

packages = c('rgdal', 'maptools', 'raster','spatstat', 'tmap','sf', 'dplyr', 'tidyverse','ggplot2','spdep','tidyselect','gganimate','cowplot','RColorBrewer','plotly','sp','FRK','animation','ggpubr')

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-17, (SVN revision 1070)
## 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/sihua/OneDrive/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/sihua/OneDrive/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Loading required package: maptools
## Checking rgeos availability: TRUE
## Loading required package: raster
## Loading required package: spatstat
## Loading required package: spatstat.data
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## Loading required package: rpart
## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.64-1 is out of date by more than 5 months; we recommend upgrading to the latest version.
## 
## Attaching package: 'spatstat'
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
## Loading required package: tmap
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: tidyverse
## Registered S3 method overwritten by 'cli':
##   method     from    
##   print.boxx spatstat
## -- Attaching packages ------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v stringr 1.4.0
## v tidyr   1.1.2     v forcats 0.5.0
## v readr   1.4.0
## -- Conflicts ---------------------------------- tidyverse_conflicts() --
## x dplyr::collapse() masks nlme::collapse()
## x tidyr::extract()  masks raster::extract()
## x dplyr::filter()   masks stats::filter()
## x dplyr::lag()      masks stats::lag()
## x dplyr::select()   masks raster::select()
## 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: tidyselect
## Loading required package: gganimate
## 
## Attaching package: 'gganimate'
## The following object is masked from 'package:raster':
## 
##     animate
## Loading required package: cowplot
## Loading required package: RColorBrewer
## Loading required package: plotly
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:raster':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
## Loading required package: FRK
## 
## Attaching package: 'FRK'
## The following object is masked from 'package:raster':
## 
##     distance
## Loading required package: animation
## Loading required package: ggpubr
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
## 
##     get_legend
## The following objects are masked from 'package:spatstat':
## 
##     border, rotate
## The following object is masked from 'package:raster':
## 
##     rotate

Data Preparation

This section is to prepar the data in order to perform analysis.

Data Sources

This study will be using the following datasets:

  • Reported COVID-19 cases and other related data at the municipality level.
    • Mexico City (9)
    • Mexico State (15)
    • Morelos State (17)

Data Import

we will be using st_read() function to read in the municipalities_COIVD shape file.

Reported COVID-19 cases:

mcovid<- st_read(dsn =  "data/geospatial",layer="municipalities_COVID")
## Reading layer `municipalities_COVID' from data source `C:\GIS415\TakeHome\IS415_Take-home_Ex02_PengSzuHua\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

Geospatial data wrangling

The data preparation is preformed to:

  1. Extracting only 3 municipality level
  • Mexico City (9)
  • Mexico State (15)
  • Morelos State (17)
  1. Extract only 19 epidemiological period
    • Week 13 to Week 32
  2. Calculate the COVID-19 rate from e-week 13 until e-week 32
extractmcovid <- mcovid %>%
  filter(`CVE_ENT` %in% c('09','15','17'))%>%
  select(-num_range("cumul", 1:12, width = NULL, vars = NULL),
         -contains(c("death","activ","new","actvr","dthrt","dethrt","Pop2010","CVE_MUN")))%>%
  mutate_at(vars(contains("cumul")),funs(covidRate = ./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.

Spatio-temporal distribution of COVID-19 rates at municipality level

To better presenting spatio-temporal distribution of COVID-19 rate, animated map will be used.

the code below is to select location ID and geometry.

extractmcovidlocation <- extractmcovid%>%
  select("CVEGEO","geometry")

the code below is to select the area name, region code, as well as the covid Rate columns.

DFextractmcovid<-st_set_geometry(extractmcovid, NULL)
covidRatetable<- DFextractmcovid %>%
  select("CVEGEO","CVE_ENT","NOMGEO",contains("Rate"))

The code below is to convert column to row on the 19 covid Rate column so that animated map can read the week as a single column.

covidRatetable<- gather(covidRatetable, "week", "covidRate", 4:23)

The code below is to use left join function to map the geometry to the Covid rate data.

DFextractmcovidlocation1<-left_join(covidRatetable,extractmcovidlocation,by="CVEGEO")

The code below is to transformed the data frame into sf format in order to plot on map.

covid19_shapes<-st_as_sf(DFextractmcovidlocation1)

The code below is to Plot a map to show the following 3 region:

  • Mexico City (9)
  • Mexico State (15)
  • Morelos State (17)
region<-tm_shape(extractmcovid) + 
  tm_fill("CVE_ENT", 
  popup.vars=c("CVE_ENT"="CVE_ENT",
  "Name"="NOMGEO",
  "CVEGEO"="CVEGEO"))

The code below is to plot the animated map by weeks on the covid Rate. Thus ,we can perform the spatio-temporal distribution analysis.

covidRate_map<-
  ggplot()+
  geom_sf(data = covid19_shapes,fill = "white")+
  geom_sf(data = covid19_shapes,aes(fill=`covidRate`))+
  ggtitle("COVID-19 Rate at municipality level")+
  xlab("")+
  ylab("")+
  labs(subtitle = "Weeks: {current_frame}",
       caption = "week 13 - week 32")+
  cowplot::background_grid(major = "none", minor = "none") +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        axis.line = element_blank(),
        legend.background = element_blank(),
        legend.position="bottom",
        plot.background = element_blank(),
        panel.background = element_blank(),
        legend.text = element_text(size=12),
        legend.title = element_text(colour="black", size=12, face="bold"),
        
        plot.title=element_text(size=20, face="bold",hjust =0.5),
        plot.subtitle = element_text(hjust = 0.5,size=12),
        plot.caption = element_text(size = 11,
                                    hjust = .5,
                                    color = "black",
                                    face = "bold"))+
  scale_fill_distiller("Rate",palette ="Reds",type = "div",
                       direction = 1)+
  transition_manual(week)
  covidspread <-animate(covidRate_map, nframe=27,fps = 2, end_pause = 20,height = 500, width =500)

The code below is to show the region map and the covid-spread map.

par(mfrow=c(1,2))
region

covidspread

Analysis:

From the animated map above shows that:

  • The COVID-19 rate distribution are increasing from week 13 to week 32

  • Mexico City is the most serious area.

  • The spread of the COVID-19 are spread outwards from Mexico City to Mexico State and Morelos State .

  • Even the south part of the Mexico city is the most serious area, however compare to Morelos State, Mexico State seems to be more affected even Morelos State is near to the most serious area.

Local Moran’s I analysis

Creating (QUEEN) contiguity based neighbours

Queen contiguity based neighbors method is created as queens method covered 360 degrees around the object.

sp_CovidRate<-as(extractmcovid, 'Spatial')
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum International Terrestrial Reference Frame 2008 in
## CRS definition

Creating (QUEEN) contiguity based neighbours

The code below is to create queen contiguity based neighbours and using row-standardized weights.

covid19_q <- poly2nb(sp_CovidRate, queen=TRUE)

summary(covid19_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
rcovid19_q <- nb2listw(covid19_q, style="B", zero.policy = TRUE)
summary(rcovid19_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 
## 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
## 
## Weights style: B 
## Weights constants summary:
##     n    nn  S0   S1    S2
## B 176 30976 962 1924 23920

Analysis:

The summary report above shows that there are total of 176 area units. The most connected area unit has 14 neighbours. There are 3 area units with only one neighbours.

Local Monte Carlo Moran’s I

In this section is to create a function to

  • Calculate Local Monte Carlo Moran’s I [ p-values and local moran statistics]

  • Plot the results in to 2 different map

The code below is to create a function to calculate local monte carlo moran.

GenerateCovidRateMap <- function(covidRateweek) 
{

fips <- order(extractmcovid$NOMGEO)
localMI <- localmoran(covidRateweek, rcovid19_q)

printCoefmat(data.frame(localMI[fips,], row.names=extractmcovid$NOMGEO[fips]), check.names=FALSE)

covid.localMI <- cbind(extractmcovid,localMI)

localMI.map <- tm_shape(covid.localMI) +
  tm_fill(col = "Ii", 
          style = "pretty", 
          title = "local moran statistics") +
  tm_borders(alpha = 0.5)

pvalue.map <- tm_shape(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's I p-values") +
  tm_borders(alpha = 0.5)

tmap_arrange(localMI.map, pvalue.map, asp=1, ncol=2)

}

the code below is to create a for loop to loop throught the function and create the individual local map and p-value map for each e-week.

weeklist <- list(extractmcovid$cumul13_covidRate,extractmcovid$cumul14_covidRate,extractmcovid$cumul15_covidRate,extractmcovid$cumul16_covidRate,extractmcovid$cumul17_covidRate,extractmcovid$cumul18_covidRate,extractmcovid$cumul19_covidRate,extractmcovid$cumul20_covidRate,extractmcovid$cumul21_covidRate,extractmcovid$cumul22_covidRate,extractmcovid$cumul23_covidRate,extractmcovid$cumul24_covidRate,extractmcovid$cumul25_covidRate,extractmcovid$cumul26_covidRate,extractmcovid$cumul27_covidRate,extractmcovid$cumul28_covidRate,extractmcovid$cumul29_covidRate,extractmcovid$cumul30_covidRate,extractmcovid$cumul31_covidRate,extractmcovid$cumul32_covidRate)

for(i in 1:length(weeklist)) 
  {
   assign(paste("week", i, sep = ""), GenerateCovidRateMap(weeklist[[i]]))
}

the code below is to combined all the generated chart in to animated formation by using walk() function and fig.show=animated so that we can see the changes from week 13 to week 32.

weeksmc <- list (week1,week2,week3,week4,week5,week6,week7,week8,week9,week10,week11,week12,week13,week14,week15,week16,week17,week18,week19,week20)
walk(weeksmc, print)
## 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.
## 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.
## 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.
## 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.
## 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.
## 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.
## 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.
## 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.
## 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.
## 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 above chart showing I and 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= 95%

critical value 0.05

Analysis:

  • The green color showing positive auto-correlation, brown is negative auto-correlation.

  • Intermix of high value with low value

  • The higher the number shows a cluster (green), the lower the color shows outlier (brown).

  • In the I value chart ,from week 13, the alternative hypothesis on statistic of bigger than > 0 (In green color) so there is positive auto-correlation. In conclusion, the spatial pattern resemble cluster distribution in the mexico city area.

  • In the local moran satistic chart,the brown color which is around Morelos State area are having statistic of smaller than > 0 so there is negative auto-correlation. In conclusion, the spatial pattern not cluster distribution outside of the mexico city area.

  • From the P value chart shows that, from week 13 the covid19 in the mexico city are small enough to be considered statiscally significantnot and the covid cases are randomly distributed as p-value are smaller than the alpha value, we failed to accept the null hypothesis.

  • However when the time past by, from week 13 to week 32 we can see that there are more area having a p-value lesser than the alph values. Which means that at the start of week 13, only the COVID-19 cases in mexico city are not randomly distributed and the rest of the area are randomly distributed. slowly, the the rest of the area changed from randomly distributed to not randomly distributed.

  • The COVID-19 cases spread outwards from mexico city. from week 13 slowly, there are increase of the cluster area for the COVID-19 places and they are not randomly distributed.

Local Getis-Ord Gi* analysis

The Getis-Ord Gi* analysis is used to detect spatial anomalies as well as hot spot and/or cold spot areas. By looking at the neighbor within a defined proximity to identify where either high or low values cluster spatially. it is statistically significant hot-spots are recognized as areas of high values where other areas within a neighborhood range also share high values too.

The analysis consists of three steps:

  • Deriving spatial weight matrix
  • Computing Gi statistics
  • Mapping Gi statistics

Creating adaptive proximity matrix

The code below is to set EPSG as 4488(Mexico).

spcovid <- spTransform(sp_CovidRate, CRS("+init=epsg:4488"))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Mexico_ITRF92 in CRS definition

The code chunk below will be used to derive an adaptive distance weight matrix with k = 8.

coords<- st_centroid(extractmcovid$geometry)
knb <- knn2nb(knearneigh(coords, k=8, longlat = FALSE))
## Warning in knearneigh(coords, k = 8, 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: 1408 
## Percentage nonzero weights: 4.545455 
## Average number of links: 8 
## Non-symmetric neighbours list
## Link number distribution:
## 
##   8 
## 176 
## 176 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 8 links
## 176 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 8 links
## 
## Weights style: B 
## Weights constants summary:
##     n    nn   S0   S1    S2
## B 176 30976 1408 2526 46164

The below plot shows the centroids of the counties and neighbors by using the code chunk below.

plot(extractmcovid$geometry, border="lightgrey")
plot(knb, coords, pch = 19, cex = 0.6, add = TRUE, col = "hotpink2")

Computing Gi statistics & Visualising local Gi

The code below is to create a function to

  • Calculate the GI* statistic

  • Map Gi* values with adaptive distance weights

  • Plot them on to the geospatial map.

GenerateCovidRateGiMap <- function(covidRateweek) 
{

fips <- order(spcovid$NOMGEO)
gi.adaptive <- localG(covidRateweek, knb_lw)
spcovid.gi <- cbind(spcovid, as.matrix(gi.adaptive))
names(spcovid.gi)[45] <- "gstat_adaptive"

tm_shape(spcovid.gi) +
  tm_fill(col = "gstat_adaptive",
          style = "pretty",
          palette = "-RdBu",
          title = "local Gi") +
  tm_borders(alpha = 0.5)

}

The code below is to use for loop to loop throught all the 20 weeks and using the above function to calculate GI* statistic for each e-week.

weeklist <- list(spcovid$cumul13_covidRate,spcovid$cumul14_covidRate,spcovid$cumul15_covidRate,spcovid$cumul16_covidRate,spcovid$cumul17_covidRate,spcovid$cumul18_covidRate,spcovid$cumul19_covidRate,spcovid$cumul20_covidRate,spcovid$cumul21_covidRate,spcovid$cumul22_covidRate,spcovid$cumul23_covidRate,spcovid$cumul24_covidRate,spcovid$cumul25_covidRate,spcovid$cumul26_covidRate,spcovid$cumul27_covidRate,spcovid$cumul28_covidRate,spcovid$cumul29_covidRate,spcovid$cumul30_covidRate,spcovid$cumul31_covidRate,spcovid$cumul32_covidRate)

for(i in 1:length(weeklist)) {
   assign(paste("weekgi", i, sep = ""), GenerateCovidRateGiMap(weeklist[[i]]))
}

The below code is to put all the GI* chart into animate form so that it will be able to show the transition from week to week.

weeksgi <- list (weekgi1,weekgi2,weekgi3,weekgi4,weekgi5,weekgi6,weekgi7,weekgi8,weekgi9,weekgi10,weekgi11,weekgi12,weekgi13,weekgi14,weekgi15,weekgi16,weekgi17,weekgi18,weekgi19,weekgi20)
walk(weeksgi, print)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## 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.

Analysis:

In the choropleth map above, region shaded in red are the hot spot areas and area shaded in blue are the cold spot areas. The darkness of the colours representing the intensity of the Gi values.

From the map, we can see that the hot spot area is always the Mexico city area, however, even if the Morelos State is the neighbour of the Mexico city, it remained as cold spot area for the covid spread. On the other hand, Mexico State also the neighbour of Mexico city but it is more affected and some of the sub area that nearer to the Mexici city actually become hot spot area for the Covid-19. Which this revel the true situation as compared to Morelos state, Mexico State are more developed. Hence there are more people in the Mexico state and higher chance to have a higher spread of COVID-19.