IS415_Take-Home_EX02
Sean Lee
9 October 2020
1 Overview
1.1 Background Context
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 (SARS-CoV-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 mid-January 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.
1.2 Objectives
In this take-home exercise, the objective 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.
- Understand the COVID-19 Rate in the study area from E-Week 13 to E-Week 32.
- Show the spatio-temporal distribution of COVID-19 Rates at municipality level using appropriate thematic mapping techniques
- Perform Local Moran’s I analysis and local Getis-Ord Gi* analysis and visualise the results using appropriate thematic mapping techniques. Describe the spatio-temporal patterns revealed by the maps.
2 Import packages & Libraries
packages = c('rgdal', 'spdep', 'raster','spatstat', 'tmap', 'tidyverse','sf')
for (p in packages) {
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}2.1 Import Geospatial shp File
mex <- readOGR(dsn = "data/geospatial", layer = "municipalities_COVID")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\65918\Desktop\NUS\AY_2021_SEM1\IS415\Take Home Assignments\IS415_Take-Home_EX02\data\geospatial", layer: "municipalities_COVID"
## with 2465 features
## It has 198 fields2.2 Extracting out Study Area
mex_sa <- mex[(mex$CVE_ENT == '09' | mex$CVE_ENT == '15' | mex$CVE_ENT == '17'),]
plot(mex_sa, main = "Study Area")
2.3 Filtering out required Columns
mex_sa_tp <- mex_sa[,c(1:6,19:38)]
mex_sa_tp@data$`Total New` <- as.numeric(apply(mex_sa_tp@data[,7:26], 1, sum))
mex_sa_tp@data$`Covid Rate (Per 10000)` <- mex_sa_tp@data$`Total New` / (mex_sa_tp@data$`Pop2020`/10000)
#Further filter out columns to drop redundant columns for analysis
mex_sa_fin <- mex_sa_tp[, c(1:4,6,27,28)]
crs(mex_sa_fin)## CRS arguments:
## +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5
## +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs3 Exploratory Data Analysis on Variables
3.1 Barplots for Categorical Variables
barplot(table(mex_sa_fin$CVE_ENT)[c(9,15,17)],
main = "Barplot of Municipalities in each CVE_ENT",
ylab = "Count")
- We can see that Mexico State (15) has the highest number of municipalities contained within it.
Plot Area by CVE_ENT
qtm(mex_sa_fin, "CVE_ENT") +
tm_legend(show = FALSE)
3.2 Histograms for Continuous Variables
hist(mex_sa_fin$Pop2020, main = "Municipality Population Distribution")
hist(mex_sa_fin$`Total New`, main = "Municipality Covid Cases Distribution")
hist(mex_sa_fin$`Covid Rate (Per 10000)`, main = "Municipality Covid Rate (Per 10000) Distribution")
- The population distribution also shows the presence of outliers - municipalities with much larger populations than other municipalities. This is a similar story that can be seen in the Covid Cases distribution plot as well.
- When the cases are balanced with population count, we can see that the distribution is less skewed than the other 2 previous plots. The tail is shorter, but outliers are still present. It would be interesting to visualise all these variables on a chloropleth map to get a basic understanding of the spatial distribution.
4 Thematic Mapping
4.1 Chloropleth
tm_shape(mex_sa_fin) +
tm_fill("Covid Rate (Per 10000)",
palette = "Blues"
) +
tm_layout(main.title = "Chloropleth of Covid Rate (Per 10000)",
main.title.size = 1.2,
main.title.position = "center",
legend.height = 0.43,
legend.width = 0.35,
legend.position = c("right","bottom"),
legend.outside = TRUE,
frame = FALSE) +
tm_borders(alpha = 0.5)
From the above chloropleth, we can see that the covid rate is predominantly higher in a set of municipalities that are shaded darker blue. This points to the possibility of spatial correlation and clustering. However, it would be more meaningful for us to plot a box map as a chloropleth’s interpretation can be very subjective, depending on how the map maker makes the plot.
4.2 Box Map
Box Map Function
boxbreaks <- function(v,mult=1.5) {
qv <- unname(quantile(v,na.rm = TRUE))
iqr <- qv[4] - qv[2]
upfence <- qv[4] + mult * iqr
lofence <- qv[2] - mult * iqr
# initialize break points vector
bb <- vector(mode="numeric",length=7)
# logic for lower and upper fences
if (lofence < qv[1]) {
# no lower outliers
bb[1] <- lofence
bb[2] <- floor(qv[1])
} else {
bb[2] <- lofence
bb[1] <- qv[1]
}
if(upfence > qv[5]) {
# no upper outliers
bb[7] <- upfence
bb[6] <- ceiling(qv[5])
} else {
bb[6] <- upfence
bb[7] <- qv[5]
}
bb[3:5] <- qv[2:4]
return(bb)
}Create get.var function
get.var <- function(vname,df) {
v <- df[vname] %>% st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}boxmap <- function(vnam,df,titlesize,legx, legy, legtitle=NA,mtitle="Box Map",mult=1.5){
var <- get.var(vnam,df)
bb <- boxbreaks(var)
tm_shape(df) +
tm_fill(vnam,
title=legtitle,
breaks=bb,
palette="-RdBu",
labels = c("lower outlier", "< 25%", "25% - 50%", "50% - 75%","> 75%", "upper outlier")
) +
tm_fill(vnam,
df %>%
filter(is.na(vnam)),
labels = c("NA")
) +
tm_borders(alpha =0.2) +
tm_layout(main.title = paste("Box Map of ",toString(vnam)," by Mexican Municipality"),
main.title.size = titlesize,
main.title.position = "center",
frame = FALSE,
legend.height = legx,
legend.width = legy,
legend.position = c("right","bottom"),
legend.outside = TRUE
)
}Plot Boxmap
mex_sa_finsf <- st_as_sf(mex_sa_fin)
boxmap("Covid Rate (Per 10000)", mex_sa_finsf
, 0.7,0.45,0.35)
The box map shows us a clearer picture of the current covid situation in the study area. This is probably due to the fact that the distribution of the variable is skewed, causing many of the values to appear as low in the initial chloropleth map. From the above box map, we can see that that distribution seems to be centred around the Mexico City, which decreases in general incidence as we get away from the city cente. This points to the possibility that spatial correlation is present.
5 Derive Weights Matrix
5.1 Contiguity Based Neighbours
5.1.1 Create QUEEN contiguity based neighbours
mex_q <- poly2nb(mex_sa_fin, queen = TRUE)
summary(mex_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:
## 721 739 765 with 1 link
## 1 most connected region:
## 761 with 14 links5.1.2 Create ROOK contiguity based neighbours
mex_r <- poly2nb(mex_sa_fin, queen = FALSE)
summary(mex_r)## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 946
## Percentage nonzero weights: 3.053977
## Average number of links: 5.375
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 14
## 3 3 18 47 27 31 22 14 6 4 1
## 3 least connected regions:
## 721 739 765 with 1 link
## 1 most connected region:
## 761 with 14 links5.1.3 Visualise contiguity based neighbour maps
par(mfrow=c(1,2))
plot(mex_sa_fin, border="lightgrey")
plot(mex_q, coordinates(mex_sa_fin), pch = 19, cex = 0.4, add = TRUE, col= "red", main="Queen Contiguity")
plot(mex_sa_fin, border="lightgrey")
plot(mex_r, coordinates(mex_sa_fin), pch = 19, cex = 0.4, add = TRUE, col = "red", main="Rook Contiguity")
Visually, there is little difference from the 2 methods. However, the queen contiguity method does have a higher average number of links.
5.2 Distance Based Neighbours
5.2.1 Fixed Distance Based Neighbours
coords <- coordinates(mex_sa_fin)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat = TRUE))
summary(k1dists)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1267 6953 10703 10470 13830 19827
The summary report shows that the largest first nearest neighbour distance is 19827 km, so using this as the upper threshold gives certainty that all units will have at least one neighbour.
mex_f <- dnearneigh(coords, 0, 19827, longlat = TRUE)
mex_f## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 30784
## Percentage nonzero weights: 99.38017
## Average number of links: 174.9091
The first nearest neighbour distance is highly skewed in distribution. By using the largest one to calculate the fixed distance weight matrix, it causes almost every single municipality to be a neighbour of all the municipality in the study region. This is not a meaningful method.
5.2.2 Adaptive Distance Based Neighbours
knn6 <- knn2nb(knearneigh(coords, k=6))
knn6## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 1056
## Percentage nonzero weights: 3.409091
## Average number of links: 6
## Non-symmetric neighbours list
This is a better method than the fixed distance method given how sparsely distributed the population is in our study area - as evidenced by the skewed distribution in the histogram for study area in the EDA portion of this report.
5.2.3 Visualise distance based neighbours
par(mfrow=c(1,2))
plot(mex_sa_fin, border="lightgrey")
plot(mex_f, coords, add=TRUE)
plot(k1, coords, add=TRUE, col="red", length=0.08)
plot(mex_sa_fin, border="lightgrey")
plot(knn6, coords, pch = 19, cex = 0.6, add = TRUE, col = "red")
Further proves that fixed distance method is not a good approach for this context.
5.3 Row standardised weights matrix
rswm_q <- nb2listw(mex_q, style="W", zero.policy = TRUE)
rswm_q## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 962
## Percentage nonzero weights: 3.10563
## Average number of links: 5.465909
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 176 30976 176 71.10695 731.679
Used style = ‘W’ for simplicity
6 Spatial Analysis
6.1 Spatial lag
COVID.lag <- lag.listw(rswm_q, mex_sa_fin$`Covid Rate (Per 10000)`)
COVID.lag## [1] 65.173501 83.919003 56.334680 58.841201 69.311153 69.831564
## [7] 66.897849 54.460995 78.951095 86.628452 72.977175 105.188880
## [13] 80.881738 86.497243 71.390864 74.573109 18.430921 27.737866
## [19] 14.423635 8.818124 24.422779 60.137774 16.501641 13.762601
## [25] 23.508170 19.893033 27.552122 35.310943 27.892150 16.669582
## [31] 20.322240 38.075314 25.812475 41.208353 29.509486 44.907087
## [37] 17.949159 54.472268 25.413166 27.721732 57.078502 14.590662
## [43] 40.425085 22.993079 26.799172 21.897430 29.690661 19.024676
## [49] 34.577960 14.085890 25.902176 23.436928 56.906521 22.877109
## [55] 31.345833 12.464604 19.622254 16.927449 37.071844 30.400805
## [61] 11.155867 29.270920 18.313267 20.654229 24.118175 45.444376
## [67] 36.069927 21.110282 38.100668 40.429342 42.517254 18.208294
## [73] 49.650651 56.415698 34.063709 16.019590 40.945406 62.080004
## [79] 21.346132 15.835091 32.250376 19.406752 30.280135 17.260040
## [85] 18.893468 34.496693 13.015758 31.259766 54.609969 18.707917
## [91] 31.262303 41.814313 20.893620 15.992568 19.706406 16.481862
## [97] 30.716693 14.663874 35.155219 37.274139 28.293344 22.702941
## [103] 21.290845 18.542540 37.226166 27.497065 36.935359 26.366913
## [109] 29.803270 14.772992 29.605548 24.473185 17.617340 32.969111
## [115] 25.948673 22.300648 46.131139 16.854869 38.624164 49.540459
## [121] 3.519392 29.252209 18.770633 36.949700 37.596190 13.772708
## [127] 14.287846 14.442223 21.209975 15.718976 30.525001 22.537199
## [133] 12.476061 23.599802 18.210045 29.324160 32.713512 58.209694
## [139] 19.037816 17.479400 28.616720 20.669144 22.456103 14.732444
## [145] 16.416916 14.048310 15.693980 14.703250 17.335272 57.362738
## [151] 17.793050 15.452228 25.844834 19.594349 15.897136 14.668824
## [157] 12.370439 12.968487 14.771576 21.108920 42.558303 13.548369
## [163] 12.081496 48.242039 18.675616 27.655608 13.932382 19.863409
## [169] 10.830841 20.794401 16.323098 18.015965 11.669421 13.576173
## [175] 17.352059 23.677386Append back to original mex_sa_fin
lag.list <- list(mex_sa_fin$NOMGEO, lag.listw(rswm_q, mex_sa_fin$`Covid Rate (Per 10000)`))
lag.res <- as.data.frame(lag.list)
colnames(lag.res) <- c("NOMGEO", "Lag COVID Rate")
mex_sa_fin@data <- left_join(mex_sa_fin@data,lag.res)## Joining, by = "NOMGEO"Visualise spatial lag COVID Rate
covid <- qtm(mex_sa_fin, "Covid Rate (Per 10000)")
lag_covid <- qtm(mex_sa_fin, "Lag COVID Rate")
tmap_arrange(covid, lag_covid, asp=1, ncol=2)
6.2 Global Moran’s I
Used to investiage whether there is an overall pattern of clustering in the Covid Rate over the entire study area.
moran.test(mex_sa_fin$`Covid Rate (Per 10000)`, listw=rswm_q, zero.policy = TRUE, na.action=na.omit)##
## Moran I test under randomisation
##
## data: mex_sa_fin$`Covid Rate (Per 10000)`
## weights: rswm_q
##
## Moran I statistic standard deviate = 10.854, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.495940195 -0.005714286 0.002135988
The p value obtained is lesser than 0.05. Thus, we can say at the 5% significance level that we have sufficient evidence to reject the null hypothesis that the distribution of the variable COVID Rate is randomly distributed. The positive Moran I statistic of 0.496 indicates to us that there is a positive spatial autocorrelation (Clustering). While this global moran’s I allows us to see that the overall pattern is clustered, it does not indicate where exactly such patterns occur in the study area. Thus, we need to perform local spatial autocorrelation.
6.3 Monte Carlo Moran’s I
set.seed(1234) #allow reproducible results
bperm= moran.mc(mex_sa_fin$`Covid Rate (Per 10000)`, listw=rswm_q, nsim=999, zero.policy = TRUE, na.action=na.omit)
bperm##
## Monte-Carlo simulation of Moran I
##
## data: mex_sa_fin$`Covid Rate (Per 10000)`
## weights: rswm_q
## number of simulations + 1: 1000
##
## statistic = 0.49594, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greaterVisualise Results
mean(bperm$res[1:999])## [1] -0.006244441var(bperm$res[1:999])## [1] 0.002138681summary(bperm$res[1:999])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.133572 -0.039006 -0.009037 -0.006244 0.023391 0.193374hist(bperm$res, freq=TRUE, breaks=25, xlab="Simulated Moran's I")
abline(v=bperm$statistic, col="red") 
We can see that the probability of obtaining a value as extreme as our test statistic is extremely small. Thus, reject the null hypothesis that the distribution of COVID rate in the study area is random.
6.4 Local Moran’s I
fips <- order(mex_sa_fin$NOMGEO)
localMI <- localmoran(mex_sa_fin$`Covid Rate (Per 10000)`, rswm_q)
head(localMI)## Ii E.Ii Var.Ii Z.Ii Pr(z > 0)
## 270 3.715790 -0.005714286 0.1857254 8.635408 2.925889e-18
## 271 3.907305 -0.005714286 0.1857254 9.079802 5.438625e-20
## 272 2.268934 -0.005714286 0.2334478 4.707815 1.251934e-06
## 273 1.899870 -0.005714286 0.1141417 5.640347 8.485396e-09
## 274 3.310531 -0.005714286 0.1857254 7.695043 7.072335e-15
## 275 1.980378 -0.005714286 0.1141417 5.878644 2.068199e-09printCoefmat(data.frame(localMI[fips,], row.names=mex_sa_fin$NOMGEO[fips]), check.names=FALSE)## Ii E.Ii Var.Ii
## Ã\201lvaro Obregón 3.13141221 -0.00571429 0.13118547
## Acambay de RuÃz Castañeda 0.26390635 -0.00571429 0.23344785
## Acolman 0.01426135 -0.00571429 0.13118547
## Aculco 0.56632612 -0.00571429 0.23344785
## Almoloya de Alquisiras 0.30401140 -0.00571429 0.23344785
## Almoloya de Juárez 0.13581163 -0.00571429 0.13118547
## Almoloya del RÃo 0.67470812 -0.00571429 0.23344785
## Amacuzac 0.32183542 -0.00571429 0.23344785
## Amanalco -0.01484814 -0.00571429 0.13118547
## Amatepec 0.65595393 -0.00571429 0.31298526
## Amecameca -0.05132161 -0.00571429 0.23344785
## Apaxco 0.19887632 -0.00571429 0.47206007
## Atenco 0.05565255 -0.00571429 0.15391044
## Atizapán 0.74051905 -0.00571429 0.31298526
## Atizapán de Zaragoza 0.02399602 -0.00571429 0.15391044
## Atlacomulco -0.11524758 -0.00571429 0.18572540
## Atlatlahucan 0.21459988 -0.00571429 0.13118547
## Atlautla 0.24540192 -0.00571429 0.18572540
## Axapusco 0.09806448 -0.00571429 0.23344785
## Axochiapan 0.05576287 -0.00571429 0.31298526
## Ayala 0.18529881 -0.00571429 0.10088550
## Ayapango 0.12507698 -0.00571429 0.15391044
## Azcapotzalco 3.71578995 -0.00571429 0.18572540
## Benito Juárez 2.04673506 -0.00571429 0.15391044
## Calimaya 0.02624586 -0.00571429 0.11414173
## Capulhuac -0.00302759 -0.00571429 0.31298526
## Chalco 0.92752879 -0.00571429 0.10088550
## Chapa de Mota 0.46452020 -0.00571429 0.23344785
## Chapultepec 0.23978916 -0.00571429 0.23344785
## Chiautla 0.12379509 -0.00571429 0.13118547
## Chicoloapan 0.05317621 -0.00571429 0.23344785
## Chiconcuac 0.09787837 -0.00571429 0.31298526
## Chimalhuacán 0.00927770 -0.00571429 0.23344785
## Coacalco de Berriozábal 0.04804869 -0.00571429 0.23344785
## Coatepec Harinas 0.45174810 -0.00571429 0.11414173
## Coatetelco 0.55220367 -0.00571429 0.15391044
## Coatlán del RÃo 0.30183344 -0.00571429 0.15391044
## Cocotitlán 0.43055795 -0.00571429 0.31298526
## Coyoacán 3.90730541 -0.00571429 0.18572540
## Coyotepec -0.04554039 -0.00571429 0.23344785
## Cuajimalpa de Morelos 2.26893414 -0.00571429 0.23344785
## Cuauhtémoc 3.75433271 -0.00571429 0.15391044
## Cuautitlán -0.13585213 -0.00571429 0.11414173
## Cuautitlán Izcalli -0.02486671 -0.00571429 0.15391044
## Cuautla -0.06787939 -0.00571429 0.23344785
## Cuernavaca 0.11106901 -0.00571429 0.13118547
## Donato Guerra 0.37624778 -0.00571429 0.23344785
## Ecatepec de Morelos 0.04146713 -0.00571429 0.08160370
## Ecatzingo 0.44738058 -0.00571429 0.23344785
## El Oro -0.09683642 -0.00571429 0.23344785
## Emiliano Zapata 0.34151463 -0.00571429 0.15391044
## Gustavo A. Madero 1.89986950 -0.00571429 0.11414173
## Huehuetoca 0.06620238 -0.00571429 0.23344785
## Hueypoxtla 0.08465995 -0.00571429 0.31298526
## Huitzilac -0.64879100 -0.00571429 0.15391044
## Huixquilucan 0.10760397 -0.00571429 0.18572540
## Isidro Fabela 0.21351661 -0.00571429 0.18572540
## Ixtapaluca -0.00068109 -0.00571429 0.15391044
## Ixtapan de la Sal 0.34631182 -0.00571429 0.18572540
## Ixtapan del Oro 0.42327324 -0.00571429 0.31298526
## Ixtlahuaca 0.24341734 -0.00571429 0.18572540
## Iztacalco 3.31053141 -0.00571429 0.18572540
## Iztapalapa 1.98037802 -0.00571429 0.11414173
## Jaltenco 0.00755830 -0.00571429 0.15391044
## Jantetelco 0.42730045 -0.00571429 0.23344785
## Jilotepec 0.31007892 -0.00571429 0.18572540
## Jilotzingo 0.01164928 -0.00571429 0.18572540
## Jiquipilco 0.36485488 -0.00571429 0.15391044
## Jiutepec 0.28604550 -0.00571429 0.23344785
## Jocotitlán 0.08604397 -0.00571429 0.13118547
## Jojutla 0.02376588 -0.00571429 0.23344785
## Jonacatepec de Leandro Valle 0.23448221 -0.00571429 0.23344785
## Joquicingo 0.09431862 -0.00571429 0.15391044
## Juchitepec -0.43747374 -0.00571429 0.10088550
## La Magdalena Contreras 3.80785505 -0.00571429 0.31298526
## La Paz -0.02838077 -0.00571429 0.15391044
## Lerma 0.01477010 -0.00571429 0.09028051
## Luvianos 0.10721642 -0.00571429 0.31298526
## Malinalco 0.27576759 -0.00571429 0.15391044
## Mazatepec 0.35469114 -0.00571429 0.15391044
## Melchor Ocampo 0.07343547 -0.00571429 0.31298526
## Metepec 0.18947611 -0.00571429 0.13118547
## Mexicaltzingo 0.44736022 -0.00571429 0.31298526
## Miacatlán 0.32815993 -0.00571429 0.13118547
## Miguel Hidalgo 2.88083576 -0.00571429 0.13118547
## Milpa Alta 4.54149654 -0.00571429 0.11414173
## Morelos 0.33742780 -0.00571429 0.15391044
## Naucalpan de Juárez 0.30194898 -0.00571429 0.11414173
## Nextlalpan -0.03645473 -0.00571429 0.11414173
## Nezahualcóyotl 0.50115932 -0.00571429 0.11414173
## Nicolás Romero 0.12350691 -0.00571429 0.13118547
## Nopaltepec 0.43630779 -0.00571429 0.94928452
## Ocoyoacac -0.00101117 -0.00571429 0.11414173
## Ocuilan 0.16230965 -0.00571429 0.11414173
## Ocuituco 0.45881920 -0.00571429 0.23344785
## Otumba 0.00078201 -0.00571429 0.23344785
## Otzoloapan 0.30792007 -0.00571429 0.18572540
## Otzolotepec 0.01036922 -0.00571429 0.15391044
## Ozumba 0.11314263 -0.00571429 0.13118547
## Papalotla 0.09712461 -0.00571429 0.31298526
## Polotitlán 0.59015434 -0.00571429 0.47206007
## Puente de Ixtla -0.12795927 -0.00571429 0.15391044
## Rayón 0.00363773 -0.00571429 0.23344785
## San Antonio la Isla 0.09152224 -0.00571429 0.15391044
## San Felipe del Progreso 0.25841641 -0.00571429 0.15391044
## San José del Rincón 0.45265169 -0.00571429 0.23344785
## San MartÃn de las Pirámides 0.00009449 -0.00571429 0.18572540
## San Mateo Atenco 0.09725230 -0.00571429 0.31298526
## San Simón de Guerrero 0.02798668 -0.00571429 0.31298526
## Santo Tomás 0.25920407 -0.00571429 0.31298526
## Soyaniquilpan de Juárez 0.22277735 -0.00571429 0.94928452
## Sultepec 0.51615315 -0.00571429 0.18572540
## Tecámac 0.00022609 -0.00571429 0.11414173
## Tejupilco -0.06464077 -0.00571429 0.13118547
## Temamatla 0.20128425 -0.00571429 0.18572540
## Temascalapa -0.03381126 -0.00571429 0.23344785
## Temascalcingo 0.05679423 -0.00571429 0.23344785
## Temascaltepec 0.23541759 -0.00571429 0.10088550
## Temixco 0.40424972 -0.00571429 0.18572540
## Temoac 0.63410025 -0.00571429 0.23344785
## Temoaya 0.18005254 -0.00571429 0.13118547
## Tenancingo -0.01925520 -0.00571429 0.15391044
## Tenango del Aire -0.05415416 -0.00571429 0.18572540
## Tenango del Valle 0.00271679 -0.00571429 0.10088550
## Teoloyucan -0.00619852 -0.00571429 0.23344785
## Teotihuacán -0.09974196 -0.00571429 0.15391044
## Tepalcingo 0.14441834 -0.00571429 0.23344785
## Tepetlaoxtoc 0.03383955 -0.00571429 0.13118547
## Tepetlixpa 0.27109306 -0.00571429 0.18572540
## Tepotzotlán 0.03150028 -0.00571429 0.13118547
## Tepoztlán -0.32562726 -0.00571429 0.13118547
## Tequixquiac 0.12717342 -0.00571429 0.23344785
## Tetecala 0.40067870 -0.00571429 0.31298526
## Tetela del Volcán 0.47880979 -0.00571429 0.23344785
## Texcaltitlán 0.23391857 -0.00571429 0.15391044
## Texcalyacac -0.00817712 -0.00571429 0.15391044
## Texcoco 0.03057381 -0.00571429 0.09028051
## Tezoyuca 0.02884594 -0.00571429 0.31298526
## Tianguistenco -0.13626139 -0.00571429 0.06301054
## Timilpan 0.31583391 -0.00571429 0.15391044
## Tláhuac 6.24857236 -0.00571429 0.18572540
## Tlalmanalco 0.10702064 -0.00571429 0.13118547
## Tlalnepantla -0.24518797 -0.00571429 0.18572540
## Tlalnepantla de Baz 0.25768033 -0.00571429 0.13118547
## Tlalpan 3.53735045 -0.00571429 0.10088550
## Tlaltizapán de Zapata 0.21236172 -0.00571429 0.13118547
## Tlaquiltenango 0.05996626 -0.00571429 0.15391044
## Tlatlaya 1.09905005 -0.00571429 0.94928452
## Tlayacapan -0.15746976 -0.00571429 0.18572540
## Toluca -0.03825691 -0.00571429 0.10088550
## Tonanitla 0.01782151 -0.00571429 0.18572540
## Tonatico 0.17498718 -0.00571429 0.47206007
## Totolapan 0.27301347 -0.00571429 0.18572540
## Tultepec -0.00134689 -0.00571429 0.18572540
## Tultitlán -0.01468854 -0.00571429 0.09028051
## Valle de Bravo -0.02563639 -0.00571429 0.13118547
## Valle de Chalco Solidaridad -0.53484420 -0.00571429 0.18572540
## Venustiano Carranza 3.77273368 -0.00571429 0.23344785
## Villa de Allende 0.52745224 -0.00571429 0.23344785
## Villa del Carbón 0.43882908 -0.00571429 0.18572540
## Villa Guerrero 0.28384665 -0.00571429 0.18572540
## Villa Victoria 0.43743371 -0.00571429 0.18572540
## Xalatlaco -0.03817267 -0.00571429 0.23344785
## Xochimilco 9.32722299 -0.00571429 0.18572540
## Xochitepec 0.40194830 -0.00571429 0.18572540
## Xonacatlán -0.00951967 -0.00571429 0.23344785
## Xoxocotla 0.23390315 -0.00571429 0.15391044
## Yautepec 0.26680939 -0.00571429 0.11414173
## Yecapixtla 0.30399490 -0.00571429 0.10088550
## Zacatepec -0.18739238 -0.00571429 0.23344785
## Zacazonapan 0.22850884 -0.00571429 0.18572540
## Zacualpan 0.65163149 -0.00571429 0.23344785
## Zacualpan de Amilpas 0.59318036 -0.00571429 0.23344785
## Zinacantepec 0.00901322 -0.00571429 0.15391044
## Zumpahuacán 0.19379577 -0.00571429 0.15391044
## Zumpango -0.01007570 -0.00571429 0.10088550
## Z.Ii Pr.z...0.
## Ã\201lvaro Obregón 8.66142143 0.0000
## Acambay de RuÃz Castañeda 0.55803082 0.2884
## Acolman 0.05515157 0.4780
## Aculco 1.18394565 0.1182
## Almoloya de Alquisiras 0.64103579 0.2607
## Almoloya de Juárez 0.39074473 0.3480
## Almoloya del RÃo 1.40826267 0.0795
## Amacuzac 0.67792598 0.2489
## Amanalco -0.02521803 0.5101
## Amatepec 1.18271017 0.1185
## Amecameca -0.09439296 0.5376
## Apaxco 0.29777420 0.3829
## Atenco 0.15642267 0.4378
## Atizapán 1.33386755 0.0911
## Atizapán de Zaragoza 0.07573090 0.4698
## Atlacomulco -0.25416192 0.6003
## Atlatlahucan 0.60827443 0.2715
## Atlautla 0.58269205 0.2801
## Axapusco 0.21478976 0.4150
## Axochiapan 0.10988839 0.4562
## Ayala 0.60137970 0.2738
## Ayapango 0.33338396 0.3694
## Azcapotzalco 8.63540803 0.0000
## Benito Juárez 5.23164671 0.0000
## Calimaya 0.09459899 0.4623
## Capulhuac 0.00480237 0.4981
## Chalco 2.93819355 0.0017
## Chapa de Mota 0.97323907 0.1652
## Chapultepec 0.50811575 0.3057
## Chiautla 0.35756776 0.3603
## Chicoloapan 0.12188501 0.4515
## Chiconcuac 0.18516847 0.4265
## Chimalhuacán 0.03102874 0.4876
## Coacalco de Berriozábal 0.11127263 0.4557
## Coatepec Harinas 1.35404518 0.0879
## Coatetelco 1.42212017 0.0775
## Coatlán del RÃo 0.78393217 0.2165
## Cocotitlán 0.77982227 0.2177
## Coyoacán 9.07980204 0.0000
## Coyotepec -0.08242763 0.5328
## Cuajimalpa de Morelos 4.70781452 0.0000
## Cuauhtémoc 9.58427430 0.0000
## Cuautitlán -0.38519565 0.6500
## Cuautitlán Izcalli -0.04881909 0.5195
## Cuautla -0.12866242 0.5512
## Cuernavaca 0.32243180 0.3736
## Donato Guerra 0.79054264 0.2146
## Ecatepec de Morelos 0.16516423 0.4344
## Ecatzingo 0.93776541 0.1742
## El Oro -0.18859446 0.5748
## Emiliano Zapata 0.88507861 0.1881
## Gustavo A. Madero 5.64034689 0.0000
## Huehuetoca 0.14884512 0.4408
## Hueypoxtla 0.16154098 0.4358
## Huitzilac -1.63918793 0.9494
## Huixquilucan 0.26294459 0.3963
## Isidro Fabela 0.50870513 0.3055
## Ixtapaluca 0.01282951 0.4949
## Ixtapan de la Sal 0.81684417 0.2070
## Ixtapan del Oro 0.76680110 0.2216
## Ixtlahuaca 0.57808702 0.2816
## Iztacalco 7.69504290 0.0000
## Iztapalapa 5.87864445 0.0000
## Jaltenco 0.03383153 0.4865
## Jantetelco 0.89620578 0.1851
## Jilotepec 0.73276907 0.2318
## Jilotzingo 0.04029055 0.4839
## Jiquipilco 0.94457237 0.1724
## Jiutepec 0.60385198 0.2730
## Jocotitlán 0.25333913 0.4000
## Jojutla 0.06101477 0.4757
## Jonacatepec de Leandro Valle 0.49713202 0.3095
## Joquicingo 0.25498160 0.3994
## Juchitepec -1.35933808 0.9130
## La Magdalena Contreras 6.81662978 0.0000
## La Paz -0.05777635 0.5230
## Lerma 0.06817512 0.4728
## Luvianos 0.20185992 0.4200
## Malinalco 0.71749089 0.2365
## Mazatepec 0.91866524 0.1791
## Melchor Ocampo 0.14147758 0.4437
## Metepec 0.53890919 0.2950
## Mexicaltzingo 0.80985578 0.2090
## Miacatlán 0.92180705 0.1783
## Miguel Hidalgo 7.96959463 0.0000
## Milpa Alta 13.45931183 0.0000
## Morelos 0.87466137 0.1909
## Naucalpan de Juárez 0.91065402 0.1812
## Nextlalpan -0.09098878 0.5362
## Nezahualcóyotl 1.50029770 0.0668
## Nicolás Romero 0.35677210 0.3606
## Nopaltepec 0.45367595 0.3250
## Ocoyoacac 0.01392076 0.4944
## Ocuilan 0.49733488 0.3095
## Ocuituco 0.96143979 0.1682
## Otumba 0.01344531 0.4946
## Otzoloapan 0.72775965 0.2334
## Otzolotepec 0.04099650 0.4836
## Ozumba 0.32815694 0.3714
## Papalotla 0.18382114 0.4271
## Polotitlán 0.86726512 0.1929
## Puente de Ixtla -0.31159969 0.6223
## Rayón 0.01935576 0.4923
## San Antonio la Isla 0.24785370 0.4021
## San Felipe del Progreso 0.67326313 0.2504
## San José del Rincón 0.94867496 0.1714
## San MartÃn de las Pirámides 0.01347873 0.4946
## San Mateo Atenco 0.18404938 0.4270
## San Simón de Guerrero 0.06023936 0.4760
## Santo Tomás 0.47353284 0.3179
## Soyaniquilpan de Juárez 0.23451580 0.4073
## Sultepec 1.21094536 0.1130
## Tecámac 0.01758296 0.4930
## Tejupilco -0.16269256 0.5646
## Temamatla 0.48032104 0.3155
## Temascalapa -0.05815200 0.5232
## Temascalcingo 0.12937319 0.4485
## Temascaltepec 0.75917213 0.2239
## Temixco 0.95128373 0.1707
## Temoac 1.32421702 0.0927
## Temoaya 0.51289127 0.3040
## Tenancingo -0.03451549 0.5138
## Tenango del Aire -0.11240026 0.5447
## Tenango del Valle 0.02654414 0.4894
## Teoloyucan -0.00100221 0.5004
## Teotihuacán -0.23967441 0.5947
## Tepalcingo 0.31072782 0.3780
## Tepetlaoxtoc 0.10920582 0.4565
## Tepetlixpa 0.64230596 0.2603
## Tepotzotlán 0.10274724 0.4591
## Tepoztlán -0.88326088 0.8115
## Tequixquiac 0.27503620 0.3916
## Tetecala 0.72641410 0.2338
## Tetela del Volcán 1.00281408 0.1580
## Texcaltitlán 0.61081869 0.2707
## Texcalyacac -0.00627771 0.5025
## Texcoco 0.12077226 0.4519
## Tezoyuca 0.06177527 0.4754
## Tianguistenco -0.52006860 0.6985
## Timilpan 0.81961905 0.2062
## Tláhuac 14.51249650 0.0000
## Tlalmanalco 0.31125448 0.3778
## Tlalnepantla -0.55567664 0.7108
## Tlalnepantla de Baz 0.72721701 0.2335
## Tlalpan 11.15487509 0.0000
## Tlaltizapán de Zapata 0.60209502 0.2736
## Tlaquiltenango 0.16741823 0.4335
## Tlatlaya 1.13389136 0.1284
## Tlayacapan -0.35213461 0.6376
## Toluca -0.10245618 0.5408
## Tonanitla 0.05461266 0.4782
## Tonatico 0.26300441 0.3963
## Totolapan 0.64676210 0.2589
## Tultepec 0.01013414 0.4960
## Tultitlán -0.02986768 0.5119
## Valle de Bravo -0.05500377 0.5219
## Valle de Chalco Solidaridad -1.22779726 0.8902
## Venustiano Carranza 7.82021168 0.0000
## Villa de Allende 1.10348881 0.1349
## Villa del Carbón 1.03152197 0.1511
## Villa Guerrero 0.67189951 0.2508
## Villa Victoria 1.02828412 0.1519
## Xalatlaco -0.06717876 0.5268
## Xochimilco 21.65622193 0.0000
## Xochitepec 0.94594350 0.1721
## Xonacatlán -0.00787596 0.5031
## Xoxocotla 0.61077938 0.2707
## Yautepec 0.80664418 0.2099
## Yecapixtla 0.97507879 0.1648
## Zacatepec -0.37601712 0.6465
## Zacazonapan 0.54349321 0.2934
## Zacualpan 1.36050123 0.0868
## Zacualpan de Amilpas 1.23952559 0.1076
## Zinacantepec 0.03754009 0.4850
## Zumpahuacán 0.50854659 0.3055
## Zumpango -0.01373133 0.5055Map local Moran’s I values
mex_sa_fin.localMI <- cbind(mex_sa_fin,localMI)Visualise Local Moran’s I Test-Statistics & P-value
localMI.map <- tm_shape(mex_sa_fin.localMI) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map <- tm_shape(mex_sa_fin.localMI) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Reds",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map, pvalue.map, asp=1, ncol=2)## 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.
From the above plots, we can see that the Local Moran’s I p value is statistically significant mostly in Mexico City (9). The rest of the study area is mostly statistically insignificant (p>0.100). We can also see that for Mexico City (9), it exhibits clustering in terms of COVID-19 Rates. This is evidenced by the positive and high local moran statistics in the region. Interestingly enough, the area South of Mexico City also has a dispersed pattern, as evidenced by the negative local moran’s I. However, this value is not statistically significant so it does not warrant further study.
6.5 Create LISA Cluster Map
6.5.1 Morans Scatterplot
moran.plot(mex_sa_fin$`Covid Rate (Per 10000)`, rswm_q, labels=as.character(mex_sa_fin$NOMGEO), xlab="COVID Rate (Per 10000)", ylab="Spatially Lag COVID Rate")
From the Moran’s Plot, we can see that most of the municipalities with high covid rates are surrounded by those with high covid rates as well. Those with lower covid rates are also surrounded by those with lower covid rates. This is expected given how Mexico City was in lockdown for some period of time and not exposed to the other Study Areas. Some areas such as Ocovoacac seem to defy this logic though, but these are generally in the minority. Let us ascertain whether those that seem to have high COVID rates surrounded by high covid rate municipalities are in Mexico City only.
moran.plot(mex_sa_fin$`Covid Rate (Per 10000)`, rswm_q, labels=as.character(mex_sa_fin$CVE_ENT), xlab="COVID Rate (Per 10000)", ylab="Spatially Lag COVID Rate")
Yes! Indeed, those municipalities with high COVID rates that are surrounded by other high COVID rate municipalities are found in Mexico City. This could point to some effectiveness of the government lockdown instilled in Mexico City.
6.5.2 Morans Scatterplot with Standardised Variable
mex_sa_fin$`Z.Covid Rate` <- scale(mex_sa_fin$`Covid Rate (Per 10000)`) %>% as.vector moran.plot(mex_sa_fin$`Z.Covid Rate`, rswm_q, labels=as.character(mex_sa_fin$NOMGEO), xlab="z-Covid Rate", ylab="Spatially Lag z-Covid Rate")
Prepare LISA Cluster Map
quadrant <- vector(mode="numeric",length=nrow(localMI))
DV <- mex_sa_fin$`Covid Rate (Per 10000)`- mean(mex_sa_fin$`Covid Rate (Per 10000)`)
C_mI <- localMI[,1] - mean(localMI[,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[,5]>signif] <- 0Build LISA Cluster Map
mex_sa_fin.localMI$quadrant <- quadrant
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(mex_sa_fin.localMI) +
tm_fill(col = "quadrant", style = "cat", palette = colors[c(sort(unique(quadrant)))+1], labels = clusters[c(sort(unique(quadrant)))+1]) +
tm_view(set.zoom.limits = c(11,17)) +
tm_borders(alpha=0.5)
The area with significant local Moran I’s p values are concentrated mostly in Mexico City. Only one other municipality also has a significant p value other than regions in Mexico City. This is in line with the other analysis that we have done above.
6.6 Getis-Ord Gi* Analysis
Useful to conduct this since the Local Moran’s I does not help us identify the degree of clustering that is present. We use the adaptive distance weights matrix as the fixed distance one we created is not meaningful.
knn6_lw <- nb2listw(knn6, style = 'B') # Use Adaptive Weights neighbours we got in previous sectionVisualise Centroids
plot(mex_sa_fin, border="lightgrey")
plot(knn6, coords, pch = 19, cex = 0.6, add = TRUE, col = "red")
Compute Gi Statistics
fips <- order(mex_sa_fin$NOMGEO)
gi.adaptive <- localG(mex_sa_fin$`Covid Rate (Per 10000)`, knn6_lw)
mex_sa_fin.gi <- cbind(mex_sa_fin, as.matrix(gi.adaptive))
names(mex_sa_fin.gi)[10] <- "gstat_adaptive"Visualise local Gi
tm_shape(mex_sa_fin.gi) +
tm_fill(col = "gstat_adaptive",
style = "pretty",
palette="-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)## 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.
The chloropleth shows us the municipalities that are hot spot areas (red), and those that are cold spot areas (blue). The intensity of the colours give us how intense the Gi values are. This further supports the local Moran’s I analysis that we did above, as it shows that the hot spots are mostly concentrated in the Mexico City area, and its surrounding municipalities. Mexico State has a mixture of mostly cool states with some semi-hot states, while Morelos State is mostly consisting of cool spots. This could be because it is less densely populated, at 364 people/km^{2} than the other 2 regions (Mexico City is at 5920 people/km^{2} while Mexico State is at 679 people/km^{2}) and hence the spread of COVID is slower.
7 Conclusions
In conclusion, we do see a significant degree in clustering in the study area, with a much higher covid rate in Mexico City municipalities as compared to Mexico State and Morelos State. This could be a result of the former’s higher population density, and the swift lockdown of the Mexico City when the COVID-19 outbreak got out of control.