Detecting Spatiotemporal Patterns of COVID-19 in Central Mexico

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

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. (Note: This paragraph is extracted from wiki: https://en.wikipedia.org/wiki/COVID-19_pandemic_in_Mexico).

2 Task and Objectives

In this take-home exercise, you are task to analyse the spatio-temporal patterns of COVID-19 case at Central Mexico (i.e. Mexico City (9), Mexico State (15) and Morelos State (17) by using localized spatial statistics methods.

3 Data Preparation

For the purpose of this study, a GIS data called municipalities_COVID is given. It is in ESRI shapefile format. The shapefile consists of reported COVID-19 cases and other related data at the municipality level. The data was extracted from Secretary of Health (Spanish: Secretaría de Salud), Government of Mexico homepage’s (http://datosabiertos.salud.gob.mx/gobmx/salud/datos_abiertos/).

The original data are stored in csv file format. The cases are recorded daily from 1st January 2020 until 6th August 2020. They had been aggregated up to epidemiological week (e-week) and municipality levels. There are a total of 32 weeks of records in the municipalities_COVID GIS data.

In this study, you are required to focus on Mexico City, Mexico State and Morelos State. Their CVE_ENT codes are 9, 15 and 17 respectively. In terms of epidemiological period, you are required to focus on e-week 13 to e-week 32.

3.1 Loading libraries

The code chunk below will check if the R packages in the packaging list have been installed. if not, install the library. After the installation, it will also load the R packages in R.

packages <- c('rgdal', 'spdep',  'tmap', 'tidyverse', 'prettydoc', 'sf', 'magick')
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p, character.only = T)
}

3.2 Read data

4 Understanding data

4.1 Extract municipalities located with the study area

Filter out Mexico City, Mexico State and Morelos State. As mentioned above, their CVE_ENT codes are 9, 15, and 17 respectively. And created another combined object which containing the three areas.

mexi_city <- muni_COVID[(muni_COVID$CVE_ENT == '09'),]
mexi_state <- muni_COVID[(muni_COVID$CVE_ENT == '15'),]
mor_state <- muni_COVID[(muni_COVID$CVE_ENT == '17'),]

mexi_combined <- muni_COVID[(muni_COVID$CVE_ENT == '09' | muni_COVID$CVE_ENT == '15' | muni_COVID$CVE_ENT == '17'),]

Plot to check the map.

plot(mexi_city, main = "Mexico City")

plot(mexi_state, main = "Mexico State")

plot(mor_state, main = "Morelos State")

plot(mexi_combined, main = "Combined Study Area")

qtm(mexi_combined, "CVE_ENT",  title = "Mexico City, Mexico State and Morelos State") +
  tm_legend(show = TRUE)

4.2 Calculate the COVID-19 rate (i.e. cases per 10000 population) from e-week 13 until e-week 32

The code chunk below loop through from e-week 13 to e-week 32. For each loop, it takes the number in columns started with “cumul”(cumulative) follow by i and divide by population / 10000. Then, store the result in a new column called “c19_rate_{i}”, with i representing week number.

Note: The “paste0” function concatenate vectors after converting to character.

for (n in c(13:32)) {
  mexi_city[[paste0("c19_rate_", n)]] <- mexi_city[[paste0("cumul",n)]] / (mexi_city$Pop2020 / 10000)
  mexi_state[[paste0("c19_rate_", n)]] <- mexi_state[[paste0("cumul",n)]] / (mexi_state$Pop2020 / 10000)
  mor_state[[paste0("c19_rate_", n)]] <- mor_state[[paste0("cumul",n)]] / (mor_state$Pop2020 / 10000) 
}

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

Use appropriate thematic mapping technique to show and describe the spatio-temporal patterns reveals.

5.1 Convert to sf object and merge into one

After coverting the three objects, we merge them into one Simple Feature DataFrame by using “rbind” function for better visualization in the following analysis.

mexi_city_sf <- st_as_sf(mexi_city)
mexi_state_sf <- st_as_sf(mexi_state)
mor_state_sf <- st_as_sf(mor_state)
central_mexi_sf <- rbind(mexi_city_sf, mexi_state_sf, mor_state_sf)

5.2 Filtered Columns

The following code chunk stores the columns name we want into selected_columns object. And create a subset from central_mexi_sf object by selecting the columns we need.

selected_columns = c("NOMGEO")
for (n in c(13:32)){
  selected_columns <- c(selected_columns, paste0("c19_rate_", n))
}
filtered_central_mexi_sf <- subset(central_mexi_sf, select = selected_columns)

5.3 Adding Covid-19 Rate to Each Borough

For preparing the facet map, we need to add COVID-19 rate(from week 13 to 32) to each borough.

Firstly, we need to create a central_mexi_weeks Simple Feature DataFrome to store the result. And extract the unique boroughs from the central_mexi_sf.

In the for loop, we loop through the boroughs to get their geometry value. Then in each borough, we get the borough’s COVID-19 rate from week 13 to 32 and assign it to c19_rate object. And use rbind function to bind the borough name, borough geometry, COVID-19 rate, and week number to central_mexi_weeks object.

central_mexi_weeks <- data.frame()
uinique_boroughs <- unique(filtered_central_mexi_sf$NOMGEO)

for (borough in uinique_boroughs) {
  geometry <- subset(filtered_central_mexi_sf, NOMGEO == borough)$geometry
  for (week in c(13:32)) {
      c19_rate <- subset(filtered_central_mexi_sf, NOMGEO == borough)[[paste0("c19_rate_",week)]]
      central_mexi_weeks <- rbind(central_mexi_weeks, data.frame(borough, geometry, c19_rate, week))
  }
}

central_mexi_weeks_sf <- st_as_sf(central_mexi_weeks)

5.4 Saving the result as RData and load

As the calculation above is time consuming, we will save the object as .rdata for saving time.

save(central_mexi_weeks_sf, file="../data/rdata/central_mexi_weeks_sf.RData")
load("../data/rdata/central_mexi_weeks_sf.RData")

5.5 Plotting the choropleth map

tm_shape(central_mexi_weeks_sf) +
    tm_fill("c19_rate", 
            style="jenks") +
    tm_borders() +
    tm_facets(by = "week", nrow=1, ncol=1, as.layers = TRUE)
    tm_layout(title.position = c("right","bottom"))

Choropleth map Based on COVID-19 Rate

As the choropleth map showing above, we can see that COVID-19 rate increase from week 17. And the virus origniated from borough - Milpa Alta which located in Mexico City. And spread out to the neighbour boroughs. Follow by Morelos State and Mexio State. And the COVID-19 rate raised from 0 to 170 in 5 months.

6 Local Moran’s I analysis

Perform local Moran’s I analysis and display the results by using appropriate thematic mapping techniques. Describe the spatio-temporal patterns reveal by the maps.

The following code check will

• loop from week 13 to week 32
• create a subset for specific week
• create (queen) continuity based neighbours of the subset using ploy2nd fucntion
• create the row-standardized weights matrix with nb2listw function
• calculate local Moran’s I value

moran_central_mexi <- data.frame()

for (week in c(13:32)) {
  
  subset <- central_mexi_weeks_sf %>% filter(week == week)
  
  wm_q <- poly2nb(subset, queen=T)
  rswm_q <- nb2listw(wm_q, zero.policy = TRUE)
  
  fips <- order(subset$borough)
  localMI <- localmoran(subset$c19_rate, rswm_q)
  
  quadrant <- vector(mode="numeric",length=nrow(localMI))
  
  DV <- subset$c19_rate - mean(subset$c19_rate)     
  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] <- 0
  
  subset$quadrant <- quadrant
  
  moran_central_mexi <- rbind(moran_central_mexi, cbind(subset, localMI))
}

moran_central_mexi_sf <- st_as_sf(moran_central_mexi)

6.1 Saving as .RData object and load it

As the calculation above is extremely time consuming, we will save the object as .rdata for saving time.

save(moran_central_mexi_sf, file="../data/rdata/moran_central_mexi_sf.RData")
load("../data/rdata/moran_central_mexi_sf.RData")

6.2 Mapping local Moran’s I values and p-values

tm_shape(moran_central_mexi_sf) +
  tm_fill(col = "Ii", 
          style = "pretty",
          title = "local moran statistics") +
  tm_facets("week", nrow=1, ncol=1) +
  tm_borders(alpha = 0.5)

tm_shape(moran_central_mexi_sf) +
  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_facets("week", nrow=1, ncol=1) +
  tm_borders(alpha = 0.5)

Local Moran’s I Statistics

Local Moran’s I p-values

From the local Moran statistics graph we can see that Mexico City has a trend of increase in local Moran’s value. And other boroughs in the study area shows a slightly decrease of the value. Based on the local Moran’s I p-values graph, we also can understand the most of the p-values in the west part of the study area is less than 0.05 which is small enough for the results to be considered statistically significant. However, the p-values kept increasing from week 13 to 32 in west part of Mexico State. Apart from that, the p-values started to decrease from week 20 to week 32 in Mexico City.

7 Local Getis-Ord Gi* analysis

7.1 Deriving distance-based weight matrix and calculating Gi statisitcs

central_mexi_coor <- rbind(coordinates(mexi_city), coordinates(mexi_state), coordinates(mor_state))

The code chunk below will

• create a Simple Feature DataFrame object for storing the result
• loop from week 13 to week 32
• create a subset and assin the specific week to it
• create adaptive proximity matrix of sebst using adaptive distance function (knn2nb)
• calculate Gi statistics based on adaptive distance matrix
• bind to the Simple Feature DataFrame object to combine the result

central_mexi_gi <- data.frame()
fips <- order(central_mexi_weeks_sf$borough)

for (i in c(13:32)) {
  
  subset <- central_mexi_weeks_sf %>% filter(week == i)
  
  knb <- knn2nb(
    knearneigh(
      central_mexi_coor, 
      k=8, 
      longlat = FALSE
    ), 
    row.names=row.names(subset$m))
  
  knb_lw <- nb2listw(knb, style = 'B')
  
  gi.adaptive <- localG(subset$c19_rate, knb_lw)
  
  gi <- cbind(subset, as.matrix(gi.adaptive))
  names(gi)[4] <- "gstat_adaptive"
  
  central_mexi_gi <- rbind(central_mexi_gi, gi)
}

head(central_mexi_gi)

7.2 Mapping Gi values based on adaptive distance weights

tm_shape(central_mexi_gi) +
  tm_fill(col = "gstat_adaptive",
          style = "pretty",
          palette = "-RdYlBu",
          title = "Local Gi Based on \nAdaptive Distance Weights") +
  tm_facets("week", ncol=1, nrow=1) +
  tm_borders(alpha = 0.5)

Local Gi Based on Adaptive Distance Weights

From the local Gi graph we can see that the hot spot areas which represents high COVID-19 rate is located at Mexico City. Furthermore, the neighbor boroughs surround Mexico City also have an increase of local Gi value through the weeks we study. The cold spot areas, on the other hand, mainly comprise of counties located on the south-western part of central Mexico.

8 References

• https://rpubs.com/wynlo/IS415_Take-home_Ex02