Harvey Lauw

Note: Recommend to use Navigation bar on the left as code chunkz may be too long for scrolling.

Overview

Objectives

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.

Installing and Launching R Packages

packages = c('rgdal', 'maptools', 'raster','spatstat', 'tmap', 'sf', 'spdep', 'tidyverse','ggplot2','gridExtra')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}

Exploratory Spatial Data Analysis

Importing data

data <- readOGR(dsn = "data/geospatial",
            layer="municipalities_COVID")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Harvey\Desktop\IS415-Geospatial Analytics and Applications\IS415_Take-home_Ex02_Harvey\data\geospatial", layer: "municipalities_COVID"
## with 2465 features
## It has 198 fields

Data wrangling

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.

municipalities <- data

#Remove population2010 columns and Week 1 to 12 columns, Convert all cases to COVID rates: Cases per 10,000 population
options(scipen=999)
municipalities@data <- municipalities@data %>%
  select(-c(5,7:18,39:50,71:82,103:114,135:146,167:178)) %>% 
  rename(`Location ID`= CVEGEO, `City/State ID`= CVE_ENT, `Town ID` = CVE_MUN, `Town name` = NOMGEO) %>%
  mutate_if(is.double, funs(./Pop2020*10000))

#Ensure every dataset in each city/state is not left out
mexico_city <- municipalities[municipalities$`City/State ID` == c("09"),]
mexico_state <- municipalities[municipalities$`City/State ID` == c("15"),]
morelos_state <- municipalities[municipalities$`City/State ID` == c("17"),]

Population Histogram

Understanding the population demographics by City/States. Mexico City’s Towns have lesser population than that in Mexico State and Morelos State.

mec1 <- ggplot(mexico_city@data, aes(x=Pop2020)) + 
  geom_histogram(color="blue", fill="lightblue") +
  labs(title = "Mexico City")

mes1 <- ggplot(mexico_state@data, aes(x=Pop2020)) + 
  geom_histogram(color="green", fill="lightblue") +
  labs(title = "Mexico State")

mos1 <- ggplot(morelos_state@data, aes(x=Pop2020)) + 
  geom_histogram(color="brown", fill="lightblue") +
  labs(title = "Morelos State")

grid.arrange(mec1, mes1, mos1, nrow = 1)

COVID-19 cumul Cases Rate by municipal

Show the spatio-temporal distribution of COVID-19 rates at municipality level - using appropriate thematic mapping technique and describe the spatio-temporal patterns reveals.

Visualize in Thematic maps instead of tables for cumul cases rate in week 13 to week 32.

Weekly study of COVID cumul cases rate in Mexico City, Mexico State and Morelas State. Using cumul cases rate as part of this research instead of

#Combine Mexico City, Mexico State & Morelos State Spatial polygon data frames together
combined_2020 <- bind(mexico_city,mexico_state)
combined_2020 <- bind(combined_2020,morelos_state)

cumul13 <- tm_shape(combined_2020) +
  tm_fill("cumul13",
          style = "cont",
          breaks = c(0,3.3),
          palette="BuGn") +
  tm_borders() 

cumul14 <- tm_shape(combined_2020) +
  tm_fill("cumul14",
          style = "cont",
          breaks = c(0,3.7),
          palette="BuGn") +
  tm_borders() 

cumul15 <- tm_shape(combined_2020) +
  tm_fill("cumul15",
          style = "cont",
          breaks = c(0,4.6),
          palette="BuGn") +
  tm_borders() 

cumul16 <- tm_shape(combined_2020) +
  tm_fill("cumul16",
          style = "cont",
          breaks = c(0,6.2),
          palette="BuGn") +
  tm_borders() 

cumul17 <- tm_shape(combined_2020) +
  tm_fill("cumul17",
          style = "cont",
          breaks = c(0,9.8),
          palette="BuGn") +
  tm_borders() 

cumul18 <- tm_shape(combined_2020) +
  tm_fill("cumul18",
          style = "cont",
          breaks = c(0,14.9),
          palette="BuGn") +
  tm_borders() 

cumul19 <- tm_shape(combined_2020) +
  tm_fill("cumul19",
          style = "cont",
          breaks = c(0,21.6),
          palette="BuGn") +
  tm_borders() 

cumul20 <- tm_shape(combined_2020) +
  tm_fill("cumul20",
          style = "cont",
          breaks = c(0,30.1),
          palette="BuGn") +
  tm_borders() 

cumul21 <- tm_shape(combined_2020) +
  tm_fill("cumul21",
          style = "cont",
          breaks = c(0,44.1),
          palette="BuGn") +
  tm_borders() 

cumul22 <- tm_shape(combined_2020) +
  tm_fill("cumul22",
          style = "cont",
          breaks = c(0,55.8),
          palette="BuGn") +
  tm_borders() 

cumul23 <- tm_shape(combined_2020) +
  tm_fill("cumul23",
          style = "cont",
          breaks = c(0,68.9),
          palette="BuGn") +
  tm_borders() 

cumul24 <- tm_shape(combined_2020) +
  tm_fill("cumul24",
          style = "cont",
          breaks = c(0,85),
          palette="BuGn") +
  tm_borders() 

cumul25 <- tm_shape(combined_2020) +
  tm_fill("cumul25",
          style = "cont",
          breaks = c(0,97.1),
          palette="BuGn") +
  tm_borders() 

cumul26 <- tm_shape(combined_2020) +
  tm_fill("cumul26",
          style = "cont",
          breaks = c(0,107),
          palette="BuGn") +
  tm_borders() 

cumul27 <- tm_shape(combined_2020) +
  tm_fill("cumul27",
          style = "cont",
          breaks = c(0,119.6),
          palette="BuGn") +
  tm_borders() 

cumul28 <- tm_shape(combined_2020) +
  tm_fill("cumul28",
          style = "cont",
          breaks = c(0,133.5),
          palette="BuGn") +
  tm_borders() 

cumul29 <- tm_shape(combined_2020) +
  tm_fill("cumul29",
          style = "cont",
          breaks = c(0,145.1),
          palette="BuGn") +
  tm_borders() 

cumul30 <- tm_shape(combined_2020) +
  tm_fill("cumul30",
          style = "cont",
          breaks = c(0,159.2),
          palette="BuGn") +
  tm_borders() 

cumul31 <- tm_shape(combined_2020) +
  tm_fill("cumul31",
          style = "cont",
          breaks = c(0,169.6),
          palette="BuGn") +
  tm_borders() 

cumul32 <- tm_shape(combined_2020) +
  tm_fill("cumul32",
          style = "cont",
          breaks = c(0,170.4),
          palette="BuGn") +
  tm_borders()

Week 13 to Week 16

As seen from the scale of the legends the maximum COVID rate of cumul cases are increasing drastically from approximately 3.3 cumul Cases/10,000 people to 6.2 cumul Cases/10,000 people, numbers may seem to be low but the rate has increased by 2.9 cases/10,000 people since the past 4 weeks. Most of the cases are clustered within the Mexico City and also at the borders of Mexico State and Morelos State.

tmap_mode("plot")
tmap_arrange(cumul13,cumul14,cumul15,cumul16,ncol = 2, nrow = 2, sync = TRUE)

Week 17 to Week 20

As seen from the scale of the legends the maximum COVID rate of cumul cases are increasing drastically from approximately 3.3 cumul Cases/10,000 people in week 13 to 30.1 cumul Cases/10,000 people, numbers may seem to be low but the rate has increased by 26.8 cases/10,000 people since the past 8 weeks. Most of the cases are still clustered within the Mexico City and also at the borders of Mexico State and Morelos State.

tmap_arrange(cumul17,cumul18,cumul19,cumul20,ncol = 4, nrow = 1, sync = TRUE)

Week 21 to Week 24

As seen from the scale of the legends the maximum COVID rate of cumul cases are stablizing and starting to be controlled from approximately 3.3 cumul Cases/10,000 people in week 13 to 85 cumul Cases/10,000 people, numbers are still high at the rate that has increased 7 times since the past 12 weeks. Most of the cases are clustered within the Mexico City and also at the borders of Mexico State and Morelos State.

tmap_arrange(cumul21,cumul22,cumul23,cumul24,ncol = 4, nrow = 1, sync = TRUE)

Week 25 to Week 28

As seen from the scale of the legends the maximum COVID rate of cumul cases are stablizing and decreasing to 0.07 cumul Cases/10,000 people, numbers are still high at the rate that has increased by 81.7 cases/10,000 people since the past 13 weeks. Most of the cases are clustered within the Mexico City and also at the borders of Mexico State and Morelos State.

tmap_arrange(cumul25,cumul26,cumul27,cumul28,ncol = 4, nrow = 1, sync = TRUE)

Week 29 to Week 32

As seen from the scale of the legends the maximum COVID rate of cumul cases are stablizing and decreasing tremendously to 2 cumul Cases/10,000 people, numbers are still at a higher rate that has increased by 166.7 cases/10,000 people the past 20 weeks. at the final week there is a sudden drop from 169.6 cumul Cases/10,000 people to 170.4 cumul Cases/10,000 people which is only an increase of 0.8 Cases/10,000 people.

tmap_arrange(cumul29,cumul30,cumul31,cumul32,ncol = 4, nrow = 1, sync = TRUE)

Weekly choropleth comparison with the same scale

Setting 0 as the lower bound & max value of cumul32 column of 175 as the upper bound of the scale, in the increment of 35 to view the color gradient change in 5 color shades.

Week 13 to 22 COVID-19 rate is still not as prominent in the first 10 weeks.

tm_shape(combined_2020) +
  tm_fill(c("cumul13", "cumul14","cumul15", "cumul16","cumul17", "cumul18","cumul19", "cumul20","cumul21", "cumul22"),
          title = c("Week 13", "Week 14","Week 15", "Week 16","Week 17", "Week 18","Week 19", "Week 20","Week 21", "Week 22"),
          breaks = c(0, 35, 70, 105, 140, 175),
          palette="BuGn") +
  tm_borders() +
  tm_layout(legend.title.size = 1,
          legend.text.size = 0.5,
          legend.position = c("LEFT", "TOP"))

Week 23 to 32

COVID-19 rate is starting to become prominent in the last 10 weeks, and we are able to see a clear distribution change across the different weeks.

tm_shape(combined_2020) +
  tm_fill(c("cumul23", "cumul24","cumul25", "cumul26","cumul27", "cumul28","cumul29", "cumul30","cumul31", "cumul32"),
          title = c("Week 23", "Week 24","Week 25", "Week 26","Week 27", "Week 28","Week 29", "Week 30","Week 31", "Week 32"),
          breaks = c(0, 35, 70, 105, 140, 175),
          palette="BuGn") +
  tm_borders() +
  tm_layout(legend.title.size = 1,
          legend.text.size = 0.5,
          legend.position = c("LEFT", "TOP"))

Perform local Moran’s I analysis

Display the results by using appropriate thematic mapping techniques
Describe the spatio-temporal patterns reveal by the maps.

Rook vs Queen Contiguity Matrices vs Distance based weight matrix

Distance based weight matrix is has a generated an average number of links of 175 which means every town is linked to one another, and unfortunately does seem to be a good set of neighbours to be used. In this case, we are left with the Rook & Queen contiguity and we will choose based on the highest number of nonzero links and the highest average number of links which is the Queen Contiguity Matrix for both Local Moran’s and Getis-Ord Gi* analysis.

#Queens Matrix
combined_2020_q <- poly2nb(combined_2020, queen=TRUE)
combined_2020_rswm_q <- nb2listw(combined_2020_q, zero.policy = TRUE)
combined_2020_rswm_q

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176 
## Number of nonzero links: 962 
## Percentage nonzero weights: 3.10563 
## Average number of links: 5.465909 
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1      S2
## W 176 30976 176 71.10695 731.679

#Rook Matrix
combined_2020_r <- poly2nb(combined_2020, queen=FALSE)
combined_2020_rswm_r <- nb2listw(combined_2020_r, zero.policy = TRUE)
combined_2020_rswm_r

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176 
## Number of nonzero links: 946 
## Percentage nonzero weights: 3.053977 
## Average number of links: 5.375 
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0      S1       S2
## W 176 30976 176 72.1476 731.4063

#Distance Matrix
k1 <- knn2nb(knearneigh(coordinates(combined_2020)))
k1dists <- unlist(nbdists(k1, coordinates(combined_2020), longlat = TRUE))
all.linked <- max(unlist(nbdists(k1, coordinates(combined_2020))))
combined_2020_d <- dnearneigh(coordinates(combined_2020), 0, all.linked, longlat = TRUE)
combined_2020_d

## Neighbour list object:
## Number of regions: 176 
## Number of nonzero links: 30800 
## Percentage nonzero weights: 99.43182 
## Average number of links: 175

par(mfrow=c(1,3))
plot(combined_2020, border="lightgrey",main="Queen")
plot(combined_2020_rswm_q, coordinates(combined_2020), pch = 19, cex = 0.6, add = TRUE, col= "red")

plot(combined_2020, border="lightgrey",main="Rook")
plot(combined_2020_rswm_r, coordinates(combined_2020), pch = 19, cex = 0.6, add = TRUE, col= "red")

plot(combined_2020, border="lightgrey",main="Distance")
plot(combined_2020_d, coordinates(combined_2020), add=TRUE)
plot(k1, coordinates(combined_2020), add=TRUE, col="red", length=0.08)

Choosing Weeks of Study

From the previous weekly analysis on the COVID cumul cases rate, big change occus from week 13 to week 20 and then week 31 to week 32, Performing local Moran’s I analysis on Weeks 13, 20, 31 & 32 which are the crucial weeks.

fips <- order(combined_2020$`Town name`)
localMI_13 <- localmoran(combined_2020$cumul13, combined_2020_rswm_q)
combined_2020.localMI_13 <- cbind(combined_2020,localMI_13)

morans_13 <- tm_shape(combined_2020.localMI_13) +
  tm_fill(col = "Ii", 
          style = "pretty",
          palette = "RdBu",
          title = "local moran statistics") +
  tm_borders(alpha = 0.5)

localMI_14 <- localmoran(combined_2020$cumul14, combined_2020_rswm_q)
combined_2020.localMI_14 <- cbind(combined_2020,localMI_14)

localMI_15 <- localmoran(combined_2020$cumul15, combined_2020_rswm_q)
combined_2020.localMI_15 <- cbind(combined_2020,localMI_15)

localMI_16 <- localmoran(combined_2020$cumul16, combined_2020_rswm_q)
combined_2020.localMI_16 <- cbind(combined_2020,localMI_16)

localMI_17 <- localmoran(combined_2020$cumul17, combined_2020_rswm_q)
combined_2020.localMI_17 <- cbind(combined_2020,localMI_17)

localMI_18 <- localmoran(combined_2020$cumul18, combined_2020_rswm_q)
combined_2020.localMI_18 <- cbind(combined_2020,localMI_18)

localMI_19 <- localmoran(combined_2020$cumul19, combined_2020_rswm_q)
combined_2020.localMI_19 <- cbind(combined_2020,localMI_19)

localMI_20 <- localmoran(combined_2020$cumul20, combined_2020_rswm_q)
combined_2020.localMI_20 <- cbind(combined_2020,localMI_20)

morans_20 <- tm_shape(combined_2020.localMI_20) +
  tm_fill(col = "Ii", 
          style = "pretty",
          palette = "RdBu",
          title = "local moran statistics") +
  tm_borders(alpha = 0.5)

localMI_21 <- localmoran(combined_2020$cumul21, combined_2020_rswm_q)
combined_2020.localMI_21 <- cbind(combined_2020,localMI_21)

localMI_22 <- localmoran(combined_2020$cumul22, combined_2020_rswm_q)
combined_2020.localMI_22 <- cbind(combined_2020,localMI_22)

localMI_23 <- localmoran(combined_2020$cumul23, combined_2020_rswm_q)
combined_2020.localMI_23 <- cbind(combined_2020,localMI_23)

localMI_24 <- localmoran(combined_2020$cumul24, combined_2020_rswm_q)
combined_2020.localMI_24 <- cbind(combined_2020,localMI_24)

localMI_25 <- localmoran(combined_2020$cumul25, combined_2020_rswm_q)
combined_2020.localMI_25 <- cbind(combined_2020,localMI_25)

localMI_26 <- localmoran(combined_2020$cumul26, combined_2020_rswm_q)
combined_2020.localMI_26 <- cbind(combined_2020,localMI_26)

localMI_27 <- localmoran(combined_2020$cumul27, combined_2020_rswm_q)
combined_2020.localMI_27 <- cbind(combined_2020,localMI_27)

localMI_28 <- localmoran(combined_2020$cumul28, combined_2020_rswm_q)
combined_2020.localMI_28 <- cbind(combined_2020,localMI_28)

localMI_29 <- localmoran(combined_2020$cumul29, combined_2020_rswm_q)
combined_2020.localMI_29 <- cbind(combined_2020,localMI_29)

localMI_30 <- localmoran(combined_2020$cumul30, combined_2020_rswm_q)
combined_2020.localMI_30 <- cbind(combined_2020,localMI_30)

localMI_31 <- localmoran(combined_2020$cumul31, combined_2020_rswm_q)
combined_2020.localMI_31 <- cbind(combined_2020,localMI_31)

morans_31 <- tm_shape(combined_2020.localMI_31) +
  tm_fill(col = "Ii", 
          style = "pretty",
          palette = "RdBu",
          title = "local moran statistics") +
  tm_borders(alpha = 0.5)

localMI_32 <- localmoran(combined_2020$cumul32, combined_2020_rswm_q)
combined_2020.localMI_32 <- cbind(combined_2020,localMI_32)

morans_32 <- tm_shape(combined_2020.localMI_32) +
  tm_fill(col = "Ii", 
          style = "pretty",
          palette = "RdBu",
          title = "local moran statistics") +
  tm_borders(alpha = 0.5)

morans_13_pvalue <- tm_shape(combined_2020.localMI_13) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 13") +
  tm_borders(alpha = 0.5)

morans_14_pvalue <- tm_shape(combined_2020.localMI_14) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 14") +
  tm_borders(alpha = 0.5)

morans_15_pvalue <- tm_shape(combined_2020.localMI_15) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 15") +
  tm_borders(alpha = 0.5)

morans_16_pvalue <- tm_shape(combined_2020.localMI_16) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 16") +
  tm_borders(alpha = 0.5)

morans_17_pvalue <- tm_shape(combined_2020.localMI_17) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 17") +
  tm_borders(alpha = 0.5)

morans_18_pvalue <- tm_shape(combined_2020.localMI_18) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 18") +
  tm_borders(alpha = 0.5)

morans_19_pvalue <- tm_shape(combined_2020.localMI_19) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 19") +
  tm_borders(alpha = 0.5)

morans_20_pvalue <- tm_shape(combined_2020.localMI_20) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 20") +
  tm_borders(alpha = 0.5)

morans_21_pvalue <- tm_shape(combined_2020.localMI_21) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 21") +
  tm_borders(alpha = 0.5)

morans_22_pvalue <- tm_shape(combined_2020.localMI_22) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 22") +
  tm_borders(alpha = 0.5)

morans_23_pvalue <- tm_shape(combined_2020.localMI_23) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 23") +
  tm_borders(alpha = 0.5)

morans_24_pvalue <- tm_shape(combined_2020.localMI_24) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 24") +
  tm_borders(alpha = 0.5)

morans_25_pvalue <- tm_shape(combined_2020.localMI_25) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 25") +
  tm_borders(alpha = 0.5)

morans_26_pvalue <- tm_shape(combined_2020.localMI_26) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 26") +
  tm_borders(alpha = 0.5)

morans_27_pvalue <- tm_shape(combined_2020.localMI_27) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 27") +
  tm_borders(alpha = 0.5)

morans_28_pvalue <- tm_shape(combined_2020.localMI_28) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 28") +
  tm_borders(alpha = 0.5)

morans_29_pvalue <- tm_shape(combined_2020.localMI_29) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 29") +
  tm_borders(alpha = 0.5)

morans_30_pvalue <- tm_shape(combined_2020.localMI_30) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 30") +
  tm_borders(alpha = 0.5)

morans_31_pvalue <- tm_shape(combined_2020.localMI_31) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 31") +
  tm_borders(alpha = 0.5)

morans_32_pvalue <- tm_shape(combined_2020.localMI_32) +
  tm_fill(col = "Pr.z...0.", 
          style = "pretty",
          breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
          palette = "RdBu",
          title = "Week 32") +
  tm_borders(alpha = 0.5)

Moran’s P-value

Week 13 to 16

tmap_mode("plot")
tmap_arrange(morans_13_pvalue,morans_14_pvalue,morans_15_pvalue,morans_16_pvalue
             ,nrow = 1, ncol = 4)

Week 17 to 20

tmap_arrange(morans_17_pvalue,morans_18_pvalue,morans_19_pvalue,morans_20_pvalue
             ,nrow = 1, ncol = 4)

Week 21 to 24

tmap_arrange(morans_21_pvalue,morans_22_pvalue,morans_23_pvalue,morans_24_pvalue
             ,nrow = 1, ncol = 4)

Week 25 to 28

tmap_arrange(morans_25_pvalue,morans_26_pvalue,morans_27_pvalue,morans_28_pvalue
             ,nrow = 1, ncol = 4)

Week 29 to 32

tmap_arrange(morans_29_pvalue,morans_30_pvalue,morans_31_pvalue,morans_32_pvalue
             ,nrow = 1, ncol = 4)

Local Indicators for Spatial Association (LISA) generation

quadrant <- vector(mode="numeric",length=nrow(localMI_13))
DV <- combined_2020$cumul13 - mean(combined_2020$cumul13)   
C_mI <- localMI_13[,1] - mean(localMI_13[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_13[,5]>signif] <- 0

combined_2020.localMI_13$quadrant <- quadrant
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")

final_morans_13 <- tm_shape(combined_2020.localMI_13) +
  tm_fill(col = "quadrant", title="Week13", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_14))
DV <- combined_2020$cumul14 - mean(combined_2020$cumul14)   
C_mI <- localMI_14[,1] - mean(localMI_14[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_14[,5]>signif] <- 0

combined_2020.localMI_14$quadrant <- quadrant

final_morans_14 <- tm_shape(combined_2020.localMI_14) +
  tm_fill(col = "quadrant", title=" Week14", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_15))
DV <- combined_2020$cumul15 - mean(combined_2020$cumul15)   
C_mI <- localMI_15[,1] - mean(localMI_15[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_15[,5]>signif] <- 0

combined_2020.localMI_15$quadrant <- quadrant

final_morans_15 <- tm_shape(combined_2020.localMI_15) +
  tm_fill(col = "quadrant", title=" Week15", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_16))
DV <- combined_2020$cumul16 - mean(combined_2020$cumul16)   
C_mI <- localMI_16[,1] - mean(localMI_16[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_16[,5]>signif] <- 0

combined_2020.localMI_16$quadrant <- quadrant

final_morans_16 <- tm_shape(combined_2020.localMI_16) +
  tm_fill(col = "quadrant", title=" Week16", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_17))
DV <- combined_2020$cumul17 - mean(combined_2020$cumul17)   
C_mI <- localMI_17[,1] - mean(localMI_17[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_17[,5]>signif] <- 0

combined_2020.localMI_17$quadrant <- quadrant

final_morans_17 <- tm_shape(combined_2020.localMI_17) +
  tm_fill(col = "quadrant", title=" Week17", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_18))
DV <- combined_2020$cumul18 - mean(combined_2020$cumul18)   
C_mI <- localMI_18[,1] - mean(localMI_18[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_18[,5]>signif] <- 0

combined_2020.localMI_18$quadrant <- quadrant

final_morans_18 <- tm_shape(combined_2020.localMI_18) +
  tm_fill(col = "quadrant", title=" Week18", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_19))
DV <- combined_2020$cumul19 - mean(combined_2020$cumul19)   
C_mI <- localMI_19[,1] - mean(localMI_19[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_19[,5]>signif] <- 0

combined_2020.localMI_19$quadrant <- quadrant

final_morans_19 <- tm_shape(combined_2020.localMI_19) +
  tm_fill(col = "quadrant", title=" Week19", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_20))
DV <- combined_2020$cumul20 - mean(combined_2020$cumul20)   
C_mI <- localMI_20[,1] - mean(localMI_20[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_20[,5]>signif] <- 0

combined_2020.localMI_20$quadrant <- quadrant

final_morans_20 <- tm_shape(combined_2020.localMI_20) +
  tm_fill(col = "quadrant", title=" Week20", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_21))
DV <- combined_2020$cumul21 - mean(combined_2020$cumul21)   
C_mI <- localMI_21[,1] - mean(localMI_21[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_21[,5]>signif] <- 0

combined_2020.localMI_21$quadrant <- quadrant

final_morans_21 <- tm_shape(combined_2020.localMI_21) +
  tm_fill(col = "quadrant", title=" Week21", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_22))
DV <- combined_2020$cumul22 - mean(combined_2020$cumul22)   
C_mI <- localMI_22[,1] - mean(localMI_22[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_22[,5]>signif] <- 0

combined_2020.localMI_22$quadrant <- quadrant

final_morans_22 <- tm_shape(combined_2020.localMI_22) +
  tm_fill(col = "quadrant", title=" Week22", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_23))
DV <- combined_2020$cumul23 - mean(combined_2020$cumul23)   
C_mI <- localMI_23[,1] - mean(localMI_23[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_23[,5]>signif] <- 0

combined_2020.localMI_23$quadrant <- quadrant

final_morans_23 <- tm_shape(combined_2020.localMI_23) +
  tm_fill(col = "quadrant", title=" Week23", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_24))
DV <- combined_2020$cumul24 - mean(combined_2020$cumul24)   
C_mI <- localMI_24[,1] - mean(localMI_24[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_24[,5]>signif] <- 0

combined_2020.localMI_24$quadrant <- quadrant

final_morans_24 <- tm_shape(combined_2020.localMI_24) +
  tm_fill(col = "quadrant", title=" Week24", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_25))
DV <- combined_2020$cumul25 - mean(combined_2020$cumul25)   
C_mI <- localMI_25[,1] - mean(localMI_25[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_25[,5]>signif] <- 0

combined_2020.localMI_25$quadrant <- quadrant

final_morans_25 <- tm_shape(combined_2020.localMI_25) +
  tm_fill(col = "quadrant", title=" Week25", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_26))
DV <- combined_2020$cumul26 - mean(combined_2020$cumul26)   
C_mI <- localMI_26[,1] - mean(localMI_26[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_26[,5]>signif] <- 0

combined_2020.localMI_26$quadrant <- quadrant

final_morans_26 <- tm_shape(combined_2020.localMI_26) +
  tm_fill(col = "quadrant", title=" Week26", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_27))
DV <- combined_2020$cumul27 - mean(combined_2020$cumul27)   
C_mI <- localMI_27[,1] - mean(localMI_27[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_27[,5]>signif] <- 0

combined_2020.localMI_27$quadrant <- quadrant

final_morans_27 <- tm_shape(combined_2020.localMI_27) +
  tm_fill(col = "quadrant", title=" Week27", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_28))
DV <- combined_2020$cumul28 - mean(combined_2020$cumul28)   
C_mI <- localMI_28[,1] - mean(localMI_28[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_28[,5]>signif] <- 0

combined_2020.localMI_28$quadrant <- quadrant

final_morans_28 <- tm_shape(combined_2020.localMI_28) +
  tm_fill(col = "quadrant", title=" Week28", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_29))
DV <- combined_2020$cumul29 - mean(combined_2020$cumul29)   
C_mI <- localMI_29[,1] - mean(localMI_29[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_29[,5]>signif] <- 0

combined_2020.localMI_29$quadrant <- quadrant

final_morans_29 <- tm_shape(combined_2020.localMI_29) +
  tm_fill(col = "quadrant", title=" Week29", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_30))
DV <- combined_2020$cumul30 - mean(combined_2020$cumul30)   
C_mI <- localMI_30[,1] - mean(localMI_30[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_30[,5]>signif] <- 0

combined_2020.localMI_30$quadrant <- quadrant

final_morans_30 <- tm_shape(combined_2020.localMI_30) +
  tm_fill(col = "quadrant", title=" Week30", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_31))
DV <- combined_2020$cumul31 - mean(combined_2020$cumul31)   
C_mI <- localMI_31[,1] - mean(localMI_31[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_31[,5]>signif] <- 0

combined_2020.localMI_31$quadrant <- quadrant

final_morans_31 <- tm_shape(combined_2020.localMI_31) +
  tm_fill(col = "quadrant", title=" Week31", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

quadrant <- vector(mode="numeric",length=nrow(localMI_32))
DV <- combined_2020$cumul32 - mean(combined_2020$cumul32)   
C_mI <- localMI_32[,1] - mean(localMI_32[,1]) 
signif <- 0.05   
quadrant[DV >0 & C_mI>0] <- 4      
quadrant[DV <0 & C_mI<0] <- 1      
quadrant[DV <0 & C_mI>0] <- 2
quadrant[DV >0 & C_mI<0] <- 3
quadrant[localMI_32[,5]>signif] <- 0

combined_2020.localMI_32$quadrant <- quadrant

final_morans_32 <- tm_shape(combined_2020.localMI_32) +
  tm_fill(col = "quadrant", title=" Week32", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
  tm_borders(alpha=0.5)

LISA analysis

Week 13 to 16 In the choropleth map below, we can see a distinct high-high quadrant reading for towns located at or near Mexico city.

tmap_mode("plot")
tmap_arrange(final_morans_13,final_morans_14,final_morans_15,final_morans_16
             ,nrow = 1, ncol = 4)

Week 17 to 20 In the choropleth map below, we can see a distinct high-high quadrant reading for towns located at or near Mexico city. With the addition of low-highs in the southwest of Mexico state by end of Week20.

tmap_arrange(final_morans_17,final_morans_18,final_morans_19,final_morans_20
             ,nrow = 1, ncol = 4)

Week 21 to 24 In the choropleth map below, we can see a distinct high-high quadrant reading for towns located at or near Mexico city. With the addition of low-highs in the southwest of Mexico state by end of Week20, even with the fluctuating change from low-low to low-highs over this 4 weeks.

tmap_arrange(final_morans_21,final_morans_22,final_morans_23,final_morans_24
             ,nrow = 1, ncol = 4)

Week 25 to 28 In the choropleth map below, we can see a distinct high-high quadrant reading for towns located at or near Mexico city. With the addition of low-highs in West and East of Morelos state.

tmap_arrange(final_morans_25,final_morans_26,final_morans_27,final_morans_28
             ,nrow = 1, ncol = 4)

Week 29 to 32 In the choropleth map below, we can see a distinct high-high quadrant reading for towns located at or near Mexico city.

tmap_arrange(final_morans_29,final_morans_30,final_morans_31,final_morans_32
             ,nrow = 1, ncol = 4)

Further Geocommunication on Local Moran’s I analysis

Week 13

The test hypotheses at 95% confident interval are:

H0 = The distribution of COVID-19 cumul case rates are randomly distributed. H1 = The distribution of COVID-19 cumul case rates are not randomly distributed.

P-value for COVID-19 cumul case rates in week 13 is 0.025 which is less than the significance level of 0.05. Hence, we reject H0 for COVID-19 cumul case rates in week 13. And that the COVID-19 cumul case rates in week 13 are not randomly distributed in overall.

tmap_arrange(morans_13,morans_13_pvalue,ncol = 2,sync = TRUE)

par(mfrow=c(1,2))
M13 <- moran.mc(combined_2020$cumul13, combined_2020_rswm_q, nsim=39, alternative="greater")
plot(M13)
M13

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  combined_2020$cumul13 
## weights: combined_2020_rswm_q  
## number of simulations + 1: 40 
## 
## statistic = 0.34652, observed rank = 40, p-value = 0.025
## alternative hypothesis: greater

combined_2020$cumul13.z <- scale(combined_2020$cumul13) %>% as.vector 
nci2 <- moran.plot(combined_2020$cumul13.z, combined_2020_rswm_q, labels=as.character(combined_2020$`Town name`), xlab="z-cumul13", ylab="Spatially Lag z-cumul13")

Week 20

P-value for COVID-19 cumul case rates in week 20 is 0.025 which is less than the significance level of 0.05. Hence, we reject H0 for COVID-19 cumul case rates in week 20. And that the COVID-19 cumul case rates in week 20 are not randomly distributed.

However as plotted on the Thematic map, we can tell that there is a cluster that are not randomly distributed within Mexico City. But The towns at the border of all 3 City/State has shown an increase in P-value and some Towns have cases where we do not reject H0 and the COVID-19 cumul case rates are randomly distributed in between each City/State. p-values are starting to be within the range of 0 to 0.2 which means covid-19 case rates in North, West & Sout of Mexico stown is starting to show a pattern that is not randomly distributed.

tmap_arrange(morans_20,morans_20_pvalue,ncol = 2,sync = TRUE)

par(mfrow=c(1,2))
M20 <- moran.mc(combined_2020$cumul20, combined_2020_rswm_q, nsim=39, alternative="greater")
plot(M20)
M20

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  combined_2020$cumul20 
## weights: combined_2020_rswm_q  
## number of simulations + 1: 40 
## 
## statistic = 0.62878, observed rank = 40, p-value = 0.025
## alternative hypothesis: greater

combined_2020$cumul20.z <- scale(combined_2020$cumul20) %>% as.vector 
nci2 <- moran.plot(combined_2020$cumul20.z, combined_2020_rswm_q, labels=as.character(combined_2020$`Town name`), xlab="z-cumul20", ylab="Spatially Lag z-cumul20")

Week 31

P-value for COVID-19 cumul case rates in week 31 is 0.025 which is less than the significance level of 0.05. Hence, we reject H0 for COVID-19 cumul case rates in week 31. And that the COVID-19 cumul case rates in week 31 are not randomly distributed.

tmap_arrange(morans_31,morans_31_pvalue,ncol = 2,sync = TRUE)

par(mfrow=c(1,2))
M31 <- moran.mc(combined_2020$cumul31, combined_2020_rswm_q, nsim=39, alternative="greater")
plot(M31)
M31

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  combined_2020$cumul31 
## weights: combined_2020_rswm_q  
## number of simulations + 1: 40 
## 
## statistic = 0.49458, observed rank = 40, p-value = 0.025
## alternative hypothesis: greater

combined_2020$cumul31.z <- scale(combined_2020$cumul31) %>% as.vector 
nci2 <- moran.plot(combined_2020$cumul31.z, combined_2020_rswm_q, labels=as.character(combined_2020$`Town name`), xlab="z-cumul31", ylab="Spatially Lag z-cumul31")

Week 32

P-value for COVID-19 cumul case rates in week 32 is 0.025 which is less than the significance level of 0.05. Hence, we reject H0 for COVID-19 cumul case rates in week 32. And that the COVID-19 cumul case rates in week 32 are not randomly distributed.

tmap_arrange(morans_32,morans_32_pvalue,ncol = 2,sync = TRUE)

par(mfrow=c(1,2))
M32 <- moran.mc(combined_2020$cumul32, combined_2020_rswm_q, nsim=39, alternative="greater")
plot(M32)
M32

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  combined_2020$cumul32 
## weights: combined_2020_rswm_q  
## number of simulations + 1: 40 
## 
## statistic = 0.49659, observed rank = 40, p-value = 0.025
## alternative hypothesis: greater

combined_2020$cumul32.z <- scale(combined_2020$cumul32) %>% as.vector 
nci2 <- moran.plot(combined_2020$cumul32.z, combined_2020_rswm_q, labels=as.character(combined_2020$`Town name`), xlab="z-cumul32", ylab="Spatially Lag z-cumul32")

Conclusion of Local Moran’s I analysis on the top 4 Crucial weeks

In week 13, the LISA map has shown explicitly that mexico’s city’s cluster is a High-High while encroaching towards the mexico state. In week 20, the LISA map has shown explicitly that mexico’s city’s cluster is a High-High while encroaching towards the mexico state. In week 31, the LISA map has shown explicitly that mexico’s city’s cluster is a High-High while retracting away from the mexico state and expanding south towards Towns within the mexico city. In week 32, the LISA map has shown explicitly that mexico’s city’s cluster is a High-High while staying within the mexico city and moving towards Morelos state.

Perform local Getis-Ord Gi* analysis

Display the results using appropriate thematic mapping techniques.
Describe the spatio-temporal patterns reveal by the maps.

localG_13 <- localG(combined_2020$cumul13, combined_2020_rswm_q)
combined_2020.localG_13 <- cbind(combined_2020,as.matrix(localG_13))
names(combined_2020.localG_13)[130] <- "gstat_adaptive"

g_13 <- tm_shape(combined_2020.localG_13) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 13") +
  tm_borders(alpha = 0.5)

localG_14 <- localG(combined_2020$cumul14, combined_2020_rswm_q)
combined_2020.localG_14 <- cbind(combined_2020,as.matrix(localG_14))
names(combined_2020.localG_14)[130] <- "gstat_adaptive"

g_14 <- tm_shape(combined_2020.localG_14) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 14") +
  tm_borders(alpha = 0.5)

localG_15 <- localG(combined_2020$cumul15, combined_2020_rswm_q)
combined_2020.localG_15 <- cbind(combined_2020,as.matrix(localG_15))
names(combined_2020.localG_15)[130] <- "gstat_adaptive"

g_15 <- tm_shape(combined_2020.localG_15) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 15") +
  tm_borders(alpha = 0.5)

localG_16 <- localG(combined_2020$cumul16, combined_2020_rswm_q)
combined_2020.localG_16 <- cbind(combined_2020,as.matrix(localG_16))
names(combined_2020.localG_16)[130] <- "gstat_adaptive"

g_16 <- tm_shape(combined_2020.localG_16) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 16") +
  tm_borders(alpha = 0.5)

localG_17 <- localG(combined_2020$cumul17, combined_2020_rswm_q)
combined_2020.localG_17 <- cbind(combined_2020,as.matrix(localG_17))
names(combined_2020.localG_17)[130] <- "gstat_adaptive"

g_17 <- tm_shape(combined_2020.localG_17) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 17") +
  tm_borders(alpha = 0.5)

localG_18 <- localG(combined_2020$cumul18, combined_2020_rswm_q)
combined_2020.localG_18 <- cbind(combined_2020,as.matrix(localG_18))
names(combined_2020.localG_18)[130] <- "gstat_adaptive"

g_18 <- tm_shape(combined_2020.localG_18) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 18") +
  tm_borders(alpha = 0.5)

localG_19 <- localG(combined_2020$cumul19, combined_2020_rswm_q)
combined_2020.localG_19 <- cbind(combined_2020,as.matrix(localG_19))
names(combined_2020.localG_19)[130] <- "gstat_adaptive"

g_19 <- tm_shape(combined_2020.localG_19) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 19") +
  tm_borders(alpha = 0.5)

localG_20 <- localG(combined_2020$cumul20, combined_2020_rswm_q)
combined_2020.localG_20 <- cbind(combined_2020,as.matrix(localG_20))
names(combined_2020.localG_20)[130] <- "gstat_adaptive"

g_20 <- tm_shape(combined_2020.localG_20) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 20") +
  tm_borders(alpha = 0.5)

localG_20 <- localG(combined_2020$cumul20, combined_2020_rswm_q)
combined_2020.localG_20 <- cbind(combined_2020,as.matrix(localG_20))
names(combined_2020.localG_20)[130] <- "gstat_adaptive"

g_2 <- tm_shape(combined_2020.localG_20) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 20") +
  tm_borders(alpha = 0.5)

localG_21 <- localG(combined_2020$cumul21, combined_2020_rswm_q)
combined_2020.localG_21 <- cbind(combined_2020,as.matrix(localG_21))
names(combined_2020.localG_21)[130] <- "gstat_adaptive"

g_21 <- tm_shape(combined_2020.localG_21) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 21") +
  tm_borders(alpha = 0.5)

localG_22 <- localG(combined_2020$cumul22, combined_2020_rswm_q)
combined_2020.localG_22 <- cbind(combined_2020,as.matrix(localG_22))
names(combined_2020.localG_22)[130] <- "gstat_adaptive"

g_22 <- tm_shape(combined_2020.localG_22) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 22") +
  tm_borders(alpha = 0.5)

localG_23 <- localG(combined_2020$cumul23, combined_2020_rswm_q)
combined_2020.localG_23 <- cbind(combined_2020,as.matrix(localG_23))
names(combined_2020.localG_23)[130] <- "gstat_adaptive"

g_23 <- tm_shape(combined_2020.localG_23) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 23") +
  tm_borders(alpha = 0.5)

localG_24 <- localG(combined_2020$cumul24, combined_2020_rswm_q)
combined_2020.localG_24 <- cbind(combined_2020,as.matrix(localG_24))
names(combined_2020.localG_24)[130] <- "gstat_adaptive"

g_24 <- tm_shape(combined_2020.localG_24) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 24") +
  tm_borders(alpha = 0.5)

localG_25 <- localG(combined_2020$cumul25, combined_2020_rswm_q)
combined_2020.localG_25 <- cbind(combined_2020,as.matrix(localG_25))
names(combined_2020.localG_25)[130] <- "gstat_adaptive"

g_25 <- tm_shape(combined_2020.localG_25) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 25") +
  tm_borders(alpha = 0.5)

localG_26 <- localG(combined_2020$cumul26, combined_2020_rswm_q)
combined_2020.localG_26 <- cbind(combined_2020,as.matrix(localG_26))
names(combined_2020.localG_26)[130] <- "gstat_adaptive"

g_26 <- tm_shape(combined_2020.localG_26) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 26") +
  tm_borders(alpha = 0.5)

localG_27 <- localG(combined_2020$cumul27, combined_2020_rswm_q)
combined_2020.localG_27 <- cbind(combined_2020,as.matrix(localG_27))
names(combined_2020.localG_27)[130] <- "gstat_adaptive"

g_27 <- tm_shape(combined_2020.localG_27) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 27") +
  tm_borders(alpha = 0.5)

localG_28 <- localG(combined_2020$cumul28, combined_2020_rswm_q)
combined_2020.localG_28 <- cbind(combined_2020,as.matrix(localG_28))
names(combined_2020.localG_28)[130] <- "gstat_adaptive"

g_28 <- tm_shape(combined_2020.localG_28) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 28") +
  tm_borders(alpha = 0.5)

localG_29 <- localG(combined_2020$cumul29, combined_2020_rswm_q)
combined_2020.localG_29 <- cbind(combined_2020,as.matrix(localG_29))
names(combined_2020.localG_29)[130] <- "gstat_adaptive"

g_29 <- tm_shape(combined_2020.localG_29) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 29") +
  tm_borders(alpha = 0.5)

localG_30 <- localG(combined_2020$cumul30, combined_2020_rswm_q)
combined_2020.localG_30 <- cbind(combined_2020,as.matrix(localG_30))
names(combined_2020.localG_30)[130] <- "gstat_adaptive"

g_30 <- tm_shape(combined_2020.localG_30) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 30") +
  tm_borders(alpha = 0.5)

localG_31 <- localG(combined_2020$cumul31, combined_2020_rswm_q)
combined_2020.localG_31 <- cbind(combined_2020,as.matrix(localG_31))
names(combined_2020.localG_31)[130] <- "gstat_adaptive"

g_31 <- tm_shape(combined_2020.localG_31) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 31") +
  tm_borders(alpha = 0.5)

localG_32 <- localG(combined_2020$cumul32, combined_2020_rswm_q)
combined_2020.localG_32 <- cbind(combined_2020,as.matrix(localG_32))
names(combined_2020.localG_32)[130] <- "gstat_adaptive"

g_32 <- tm_shape(combined_2020.localG_32) +
  tm_fill(col = "gstat_adaptive", 
          style = "pretty",
          palette="-RdBu",
          title = "Week 32") +
  tm_borders(alpha = 0.5)

Local Getis-Ord Gi* analysis

Week 13 to 16

tmap_mode("plot")
tmap_arrange(g_13,g_14,g_15,g_16
             ,nrow = 1, ncol = 4)

Week 17 to 20

tmap_arrange(g_17,g_18,g_19,g_20
             ,nrow = 1, ncol = 4)

Week 21 to 24

tmap_arrange(g_21,g_22,g_23,g_24
             ,nrow = 1, ncol = 4)

Week 25 to 28

tmap_arrange(g_25,g_26,g_27,g_28
             ,nrow = 1, ncol = 4)

Week 29 to 32

tmap_arrange(g_29,g_30,g_31,g_32
             ,nrow = 1, ncol = 4)

Conclusion on Local Getis-Ord Gi* analysis

In conclusion, the choropleth over the 20 weeks above shows a clear sign of hot spot areas in the COVID-19 cumul cases rate per 10,000 population are backtracking back to Mexico City where the Hottest spot are originated and centered from, and the hot areas can sometime move towards the west of Mexico state and back to Mexico city again. The cold spot areas mainly comprise of towns located on the West & North of Mexico state as well as the South of Morelos State.

IS415-Geospatial Analytics and Applications: Take Home Ex2