Since late December 2019, an outbreak of a novel coronavirus disease (COVID-19; previously known as 2019-nCoV) was reported in Wuhan, China, which has subsequently affected 210 countries worldwide. In general, COVID-19 is an acute resolved disease but it can also be deadly, with a 2% case fatality rate.
The COVID-19 pandemic in Mexico is part of the ongoing worldwide pandemic of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARSCoV-2). The virus was confirmed to have reached Mexico in February 2020. However, the National Council of Science and Technology (CONACYT) reported two cases of COVID-19 in midJanuary 2020 in the states of Nayarit and Tabasco, one case per state. As of September 26, there had been 726,431 confirmed cases of COVID-19 in Mexico and 76,243 reported deaths, although the Secretariat of Health, through the “Programa Centinela” (Spanish for “Sentinel Program”) estimated in mid July 2020 that there were more than 2,875,734 cases in Mexico, because they were considering the total number of cases confirmed as a statistical sample. (Note: This paragraph is extracted from wiki: https://en.wikipedia.org/wiki/COVID19_pandemic_in_Mexico).
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.
The analysis techniques and methods will be follow along using the Prof’s Kam Hands On Exercise Materials as follow: - IS415 In-class Exercise 4 - https://rpubs.com/tskam/IS415-Hands-on_Ex06 - https://rpubs.com/tskam/IS415-Hands-on_Ex07
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, only 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, and only focus on e-week 13 to e-week 32.
This code chunk will check if the R packages in the packaging list have been installed. if no, go ahead to install the missing one. After the installation, it will also load the R packages in R.
#c is creating a list to a packages
#"<-" && "=" serves as the same function
packages <- c('sp', 'spdep', 'rgdal', 'rgeos', 'raster', 'sf', 'tmap', 'tidyverse', 'readr', 'ggplot2', 'spatstat', 'dplyr')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p, character.only = T)
}The code chunk below uses st_read() function of sf package to import the datas into R as simple feature dataframes. And also import data by using read_csv() of readr package. Read_excel is being used for the excel format of “xls” & “xlsx”. While “csv” uses read_csv.
## OGR data source with driver: ESRI Shapefile
## Source: "D:\Google Drive\SMU\Year 4\Term 1\IS415-Geospatial Analytics and Applications\Take Home Exercises\Take-home_Ex02\data\geospatial", layer: "municipalities_COVID"
## with 2465 features
## It has 198 fields
# mexico_city <- municipalities_COVID_sf %>%
# filter(CVE_ENT == "09")
#
# mexico_state <- municipalities_COVID_sf %>%
# filter(CVE_ENT == "15")
#
# morelos_state <- municipalities_COVID_sf %>%
# filter(CVE_ENT == "17")
mexico_city <- municipalities_COVID[(municipalities_COVID$CVE_ENT == '09'),]
mexico_state <- municipalities_COVID[(municipalities_COVID$CVE_ENT == '15'),]
morelos_state <- municipalities_COVID[(municipalities_COVID$CVE_ENT == '17'),]
mexi_combined <- municipalities_COVID[(municipalities_COVID$CVE_ENT == '09' | municipalities_COVID$CVE_ENT == '15' | municipalities_COVID$CVE_ENT == '17'),]
#
# mexi_combined <- st_as_sf(mexi_combined)
#
# mexi_combined <- mexi_combined %>%
# filter("CVE_ENT" == "09") %>%
# mutate("area" = "Mexico City")
#
# mexi_combined <- mexi_combined %>%
# filter("CVE_ENT" == "15") %>%
# mutate("area" = "Mexico State")
#
# mexi_combined <- mexi_combined %>%
# filter("CVE_ENT" == "17") %>%
# mutate("area" = "Morelos State")Plot based on the data analyses on Cumulative, Active and Death rate of the Covid 19.
Add the sum of week 13 to 32 of new cases to achieve the cumulative cases Similarly we can perform the same method for active and death cases but using the respective column accordingly.
#c(1:6) --> dplyr::select("CVEGEO", "CVE_ENT", "CVE_MUN", "NOMGEO", "Pop2020", "wk13_32_total_count")
mexico_city_filtered <- mexico_city[,c(1:6, 19:38, 51:70, 83:102, 115:134)]
mexico_city_filtered_cumu <- mexico_city[,c(1:6,19:38)]
mexico_state_filtered <- mexico_state[,c(1:6, 19:38, 51:70, 83:102, 115:134)]
morelos_state_filtered <- morelos_state[,c(1:6, 19:38, 51:70, 83:102, 115:134)]
mexi_combined_filtered <- mexi_combined[,c(1:6, 19:38, 51:70, 83:102, 115:134)]
mexico_city_filtered@data$`wk13_32_cumul_count` <- as.numeric(apply(mexico_city_filtered@data[,7:26], 1, sum))
#mexico_city_filtered@data$`wk13_32_cumulative` <- mexico_city_filtered@data$`cumul32` - mexico_city_filtered@data$`cumul12`
mexico_city_filtered@data$`wk13_32_activ_count` <- as.numeric(apply(mexico_city_filtered@data[,47:66], 1, sum))
mexico_city_filtered@data$`wk13_32_death_count` <- as.numeric(apply(mexico_city_filtered@data[,67:86], 1, sum))
mexico_state_filtered@data$`wk13_32_cumul_count` <- as.numeric(apply(mexico_state_filtered@data[,7:26], 1, sum))
#mexico_state_filtered@data$`wk13_32_cumulative` <- mexico_state_filtered@data$`cumul32` - mexico_state_filtered@data$`cumul12`
mexico_state_filtered@data$`wk13_32_activ_count` <- as.numeric(apply(mexico_state_filtered@data[,47:66], 1, sum))
mexico_state_filtered@data$`wk13_32_death_count` <- as.numeric(apply(mexico_state_filtered@data[,67:86], 1, sum))
morelos_state_filtered@data$`wk13_32_cumul_count` <- as.numeric(apply(morelos_state_filtered@data[,7:26], 1, sum))
#morelos_state_filtered@data$`wk13_32_cumulative` <- morelos_state_filtered@data$`cumul32` - morelos_state_filtered@data$`cumul12`
morelos_state_filtered@data$`wk13_32_activ_count` <- as.numeric(apply(morelos_state_filtered@data[,47:66], 1, sum))
morelos_state_filtered@data$`wk13_32_death_count` <- as.numeric(apply(morelos_state_filtered@data[,67:86], 1, sum))
mexi_combined_filtered@data$`wk13_32_cumul_count` <- as.numeric(apply(mexi_combined_filtered@data[,7:26], 1, sum))
#mexi_combined_filtered@data$`wk13_32_cumulative` <- mexi_combined_filtered@data$`cumul32` - mexi_combined_filtered@data$`cumul12`
mexi_combined_filtered@data$`wk13_32_activ_count` <- as.numeric(apply(mexi_combined_filtered@data[,47:66], 1, sum))
mexi_combined_filtered@data$`wk13_32_death_count` <- as.numeric(apply(mexi_combined_filtered@data[,67:86], 1, sum))For effective comparison, calculate the COVID-19 rate from e-week 13 till 32 based on cumulative cases Rate by per case per 10,000 population
mexico_city_filtered@data$`rate_of_covid` <- (mexico_city_filtered@data$`wk13_32_cumul_count` * 10000)/ mexico_city_filtered@data$`Pop2020`
mexico_state_filtered@data$`rate_of_covid` <- (mexico_state_filtered@data$`wk13_32_cumul_count` * 10000)/ mexico_state_filtered@data$`Pop2020`
morelos_state_filtered@data$`rate_of_covid` <- (morelos_state_filtered@data$`wk13_32_cumul_count` * 10000)/ morelos_state_filtered@data$`Pop2020`
mexi_combined_filtered@data$`rate_of_covid` <- (mexi_combined_filtered@data$`wk13_32_cumul_count` * 10000)/ mexi_combined_filtered@data$`Pop2020`
mexico_city_filtered_roc_sf <- st_as_sf(mexico_city_filtered) %>%
select("CVE_ENT", "NOMGEO", "rate_of_covid")
mexico_state_filtered_roc_sf <- st_as_sf(mexico_state_filtered) %>%
select("CVE_ENT", "NOMGEO", "rate_of_covid")
morelos_state_filtered_roc_sf <- st_as_sf(morelos_state_filtered) %>%
select("CVE_ENT", "NOMGEO", "rate_of_covid")
mexi_combined_filtered_roc_sf <- st_as_sf(mexi_combined_filtered) %>%
select("CVE_ENT", "NOMGEO", "rate_of_covid")Calculating the COVID-19 rate from e-week 13 till 32 based on active cases Rate by per case per 10,000 population
mexico_city_filtered@data$`active_rate_of_covid` <- (mexico_city_filtered@data$`wk13_32_activ_count` * 10000)/ mexico_city_filtered@data$`Pop2020`
mexico_state_filtered@data$`active_rate_of_covid` <- (mexico_state_filtered@data$`wk13_32_activ_count` * 10000)/ mexico_state_filtered@data$`Pop2020`
morelos_state_filtered@data$`active_rate_of_covid` <- (morelos_state_filtered@data$`wk13_32_activ_count` * 10000)/ morelos_state_filtered@data$`Pop2020`
mexi_combined_filtered@data$`active_rate_of_covid` <- (mexi_combined_filtered@data$`wk13_32_activ_count` * 10000)/ mexi_combined_filtered@data$`Pop2020`
#print(mexico_city_filtered)
#print(mexico_state_filtered)
#print(morelos_state_filtered)
#print(mexi_combined_filtered)Calculating the COVID-19 rate from e-week 13 till 32 based on death cases Rate by per case per 10,000 population
mexico_city_filtered@data$`death_rate_of_covid` <- (mexico_city_filtered@data$`wk13_32_death_count` * 10000)/ mexico_city_filtered@data$`Pop2020`
mexico_state_filtered@data$`death_rate_of_covid` <- (mexico_state_filtered@data$`wk13_32_death_count` * 10000)/ mexico_state_filtered@data$`Pop2020`
morelos_state_filtered@data$`death_rate_of_covid` <- (morelos_state_filtered@data$`wk13_32_death_count` * 10000)/ morelos_state_filtered@data$`Pop2020`
mexi_combined_filtered@data$`death_rate_of_covid` <- (mexi_combined_filtered@data$`wk13_32_death_count` * 10000)/ mexi_combined_filtered@data$`Pop2020`
mexico_city_filtered_sf <- st_as_sf(mexico_city_filtered)
mexico_state_filtered_sf <- st_as_sf(mexico_state_filtered)
morelos_state_filtered_sf <- st_as_sf(morelos_state_filtered)
mexi_combined_filtered_sf <- st_as_sf(mexi_combined_filtered)
#print(mexico_city_filtered_sf)
#print(mexico_state_filtered_sf)
#print(morelos_state_filtered_sf)
#print(mexi_combined_filtered_sf)#Size bottom left top right
par(mar=c(7, 5, 2, 1))
mexico_city_filtered_bp <- barplot(mexico_city_filtered$wk13_32_cumul_count, ylim = c(0, 16000), main = "Total No. of COVID 19 cases in each Municipalities of Mexico City", ylab="Cumulative Count")
axis(1, at=mexico_city_filtered_bp, labels=mexico_city_filtered$NOMGEO, las=3, cex.lab=4, cex.axis=0.7)#Size bottom left top right
par(mar=c(7, 5, 2, 1))
mexico_state_filtered_bp <- barplot(mexico_state_filtered$wk13_32_cumul_count, ylim = c(0, 8000), main = "Total No. of COVID 19 cases in each Municipalities of Mexico State", ylab="Cumulative Count")
axis(1, at=mexico_state_filtered_bp, labels=mexico_state_filtered$NOMGEO, las=3, cex.lab=4, cex.axis=0.7)#Size bottom left top right
par(mar=c(7, 5, 2, 1))
morelos_state_filtered_bp <- barplot(morelos_state_filtered$wk13_32_cumul_count, ylim = c(0, 1200), main = "Total No. of COVID 19 cases in each Municipalities of Morelos State", ylab="Cumulative Count")
axis(1, at=morelos_state_filtered_bp, labels=morelos_state_filtered$NOMGEO, las=3, cex.lab=4, cex.axis=0.7)#Size bottom left top right
par(mar=c(7, 5, 2, 1))
mexico_city_filtered_roc_bp <- barplot(mexico_city_filtered_roc_sf$rate_of_covid, ylim = c(0, 200), main = "Rate of COVID 19 in each Municipalities of Mexico City", ylab="Rate of COVID 19")
axis(1, at=mexico_city_filtered_roc_bp, labels=mexico_city_filtered_roc_sf$NOMGEO, las=3, cex.lab=4, cex.axis=0.7)#Size bottom left top right
par(mar=c(7, 5, 2, 1))
mexico_state_filtered_roc_bp <- barplot(mexico_state_filtered_roc_sf$rate_of_covid, ylim = c(0, 200), main = "Rate of COVID 19 in each Municipalities of Mexico State", ylab="Rate of COVID 19")
axis(1, at=mexico_state_filtered_roc_bp, labels=mexico_state_filtered_roc_sf$NOMGEO, las=3, cex.lab=4, cex.axis=0.7)#Size bottom left top right
par(mar=c(7, 5, 2, 1))
morelos_state_filtered_roc_bp <- barplot(morelos_state_filtered_roc_sf$rate_of_covid, ylim = c(0, 60), main = "Rate of COVID 19 in each Municipalities of Morelos State", ylab="Rate of COVID 19")
axis(1, at=morelos_state_filtered_roc_bp, labels=morelos_state_filtered_roc_sf$NOMGEO, las=3, cex.lab=4, cex.axis=0.7) ## Histogram to compare Frequency for Municipalities of population and COVID 19 Rate
hist(mexi_combined_filtered$`rate_of_covid`, main = "Municipality COVID Rate Distribution (Per 10,000)")COVID 19 Rate of Mexico City (Per 10,000)
tm_shape(mexico_city_filtered) +
tm_fill("rate_of_covid",
style = "quantile",
thres.poly = 0,
n = 5,
palette = "Oranges",
title = "COVID 19 Rate of Mexico City (Per 10,000)") +
tm_layout(legend.show = TRUE, title.position = c("center", "center"),
title.size = 20) +
tm_borders(alpha = 0.5) + tm_text("NOMGEO", size=0.2) COVID 19 Rate of Mexico State (Per 10,000)
tm_shape(mexico_state_filtered) +
tm_fill("rate_of_covid",
style = "quantile",
thres.poly = 0,
n = 5,
palette = "Oranges",
title = "COVID 19 Rate of Mexico State (Per 10,000)") +
tm_layout(legend.show = TRUE, title.position = c("center", "center"),
title.size = 20) +
tm_borders(alpha = 0.5) + tm_text("NOMGEO", size=0.2) COVID 19 Rate of Morelos State (Per 10,000)
tm_shape(morelos_state_filtered) +
tm_fill("rate_of_covid",
style = "quantile",
thres.poly = 0,
n = 5,
palette = "Oranges",
title = "COVID 19 Rate of Morelos State (Per 10,000)") +
tm_layout(legend.show = TRUE, title.position = c("center", "center"),
title.size = 20) +
tm_borders(alpha = 0.5) + tm_text("NOMGEO", size=0.2) COVID 19 Rate of Three Selected Study Area (Per 10,000)
tm_shape(mexi_combined_filtered) +
tm_fill("rate_of_covid",
style = "quantile",
thres.poly = 0,
n = 5,
palette = "Oranges",
title = "COVID 19 Rate of Three Selected Study Area (Per 10,000)") +
tm_layout(legend.show = TRUE, title.position = c("center", "center"),
title.size = 20) +
tm_borders(alpha = 0.5) + tm_text("NOMGEO", size=0.2) Based on the above result, Mexico City have the most highest concentration of COVID 19 Rate.
Creating Boxbreaks function
boxbreaks <- function(v, mult=1.5) {
qv <- unname(quantile(v))
iqr <- qv[4] - qv[2]
upfence <- qv[4] + mult * iqr
lofence <- qv[2] - mult * iqr
bb <- vector(mode="numeric",length=7)
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)
}Creating the get.var function
get.var <- function(vname,df) {
v <- df[vname] %>% st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}Box Map Function
boxmap <- function(vnam, df, 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_borders() +
tm_layout(main.title = "Box Map of COVID-19 Rate in our Study Area", main.title.position = "center", main.title.size = 1.3, legend.position = c("right", "top"), legend.outside = TRUE , frame = FALSE)
} + tm_text("NOMGEO", size=0.15) Mapping out the result
The result above allows us to better identify the Outlier in the Mexico city.
Using poly2nb() of spdep package to compute contiguity weight matrices for the study area. This function builds a neighbours list based on regions with contiguous boundaries.
wm_q_mexi_combined_filtered <- poly2nb(mexi_combined_filtered, queen=TRUE)
summary(wm_q_mexi_combined_filtered)## 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 links
The summary report above shows that there are 176 area units in Mexico City, Mexico State and Morelos State. The most connect area unit has 14 neighbours. There are 1 area unit with only 11 neighbours.
List complete weight matrix
## List of 176
## $ : int [1:5] 4 14 15 73 120
## $ : int [1:5] 6 9 11 12 13
## $ : int [1:4] 9 15 53 78
## $ : int [1:8] 1 14 16 36 49 74 120 125
## $ : int [1:5] 6 13 14 16 74
## $ : int [1:8] 2 5 10 12 13 74 86 138
## $ : int [1:3] 9 11 78
## $ : int [1:8] 10 11 12 41 66 150 161 164
## $ : int [1:7] 2 3 7 11 13 15 78
## $ : int [1:5] 6 8 12 41 138
## $ : int [1:9] 2 7 8 9 12 59 78 117 150
## $ : int [1:5] 2 6 8 10 11
## $ : int [1:6] 2 5 6 9 14 15
## $ : int [1:6] 1 4 5 13 15 16
## $ : int [1:7] 1 3 9 13 14 53 73
## $ : int [1:4] 4 5 14 74
## $ : int [1:4] 19 30 101 118
## $ : int [1:7] 27 44 49 97 108 109 116
## $ : int [1:4] 17 61 87 118
## $ : int [1:4] 37 96 113 133
## $ : int [1:7] 23 58 90 103 122 130 134
## $ : int [1:4] 28 89 114 117
## $ : int [1:7] 21 48 102 126 127 130 134
## $ : int [1:3] 96 98 121
## $ : int [1:4] 31 33 84 119
## $ : int [1:2] 52 112
## $ : int [1:6] 18 44 46 49 115 116
## $ : int [1:3] 22 89 117
## $ : int [1:6] 54 62 73 76 120 137
## $ : int [1:5] 17 64 72 101 118
## $ : int [1:5] 25 50 84 163 171
## $ : int [1:4] 77 81 91 100
## $ : int [1:6] 25 66 84 105 110 119
## $ : int [1:8] 43 70 71 88 89 106 117 122
## $ : int [1:3] 67 78 117
## $ : int [1:4] 4 49 124 125
## $ : int [1:8] 20 56 102 106 113 129 133 134
## $ : int [1:3] 41 99 119
## $ : int [1:4] 51 107 111 136
## $ : int [1:8] 69 75 107 111 124 125 136 137
## $ : int [1:9] 8 10 38 55 66 99 105 119 138
## $ : int [1:4] 61 72 118 128
## $ : int [1:4] 34 70 71 117
## $ : int [1:7] 18 27 46 85 109 115 116
## $ : int [1:4] 47 55 86 115
## $ : int [1:3] 27 44 115
## $ : int [1:4] 45 74 86 115
## $ : int [1:4] 23 57 126 127
## $ : int [1:11] 4 18 27 36 60 74 97 115 120 125 ...
## $ : int [1:4] 31 157 163 171
## $ : int [1:4] 39 111 112 136
## $ : int [1:3] 26 112 136
## $ : int [1:5] 3 15 67 73 78
## $ : int [1:5] 29 62 76 83 103
## $ : int [1:6] 41 45 86 115 119 138
## $ : int [1:5] 37 123 129 133 135
## $ : int [1:3] 48 94 126
## $ : int [1:5] 21 63 64 90 103
## $ : int [1:4] 11 79 117 150
## $ : int [1:6] 49 75 97 125 136 141
## $ : int [1:5] 19 42 87 95 118
## $ : int [1:5] 29 54 73 83 131
## $ : int [1:6] 58 64 72 76 103 128
## $ : int [1:7] 30 58 63 72 80 90 101
## $ : int [1:6] 68 79 104 106 114 117
## $ : int [1:9] 8 33 41 84 99 105 110 164 168
## $ : int [1:10] 35 53 70 73 78 83 92 117 122 131
## $ : int [1:6] 65 79 104 135 146 156
## $ : int [1:3] 40 75 124
## $ : int [1:7] 34 43 67 71 92 117 122
## $ : int [1:3] 34 43 70
## $ : int [1:6] 30 42 63 64 118 128
## $ : int [1:8] 1 15 29 53 62 67 120 131
## $ : int [1:8] 4 5 6 16 47 49 86 115
## $ : int [1:8] 40 60 69 97 124 125 136 141
## $ : int [1:7] 29 54 63 103 111 128 137
## $ : int 32
## $ : int [1:8] 3 7 9 11 35 53 67 117
## $ : int [1:8] 59 65 68 104 117 148 150 156
## $ : int [1:4] 64 90 101 140
## $ : int [1:4] 32 91 108 109
## $ : int [1:5] 94 102 126 132 139
## $ : int [1:6] 54 62 67 103 122 131
## $ : int [1:7] 25 31 33 66 110 143 171
## $ : int [1:3] 44 109 115
## $ : int [1:6] 6 45 47 55 74 138
## $ : int [1:2] 19 61
## $ : int [1:4] 34 89 106 114
## $ : int [1:6] 22 28 34 88 114 117
## $ : int [1:6] 21 58 64 80 130 140
## $ : int [1:5] 32 81 100 108 109
## $ : int [1:3] 67 70 122
## $ : int [1:3] 98 102 113
## $ : int [1:3] 57 82 126
## $ : int 61
## $ : int [1:5] 20 24 98 113 133
## $ : int [1:8] 18 49 60 75 100 108 136 141
## $ : int [1:7] 24 93 96 102 113 132 139
## $ : int [1:5] 38 41 66 105 119
## [list output truncated]
## - attr(*, "class")= chr "nb"
## - attr(*, "region.id")= chr [1:176] "270" "271" "272" "273" ...
## - attr(*, "call")= language poly2nb(pl = mexi_combined_filtered, queen = TRUE)
## - attr(*, "type")= chr "queen"
## - attr(*, "sym")= logi TRUE
wm_r_mexi_combined_filtered <- poly2nb(mexi_combined_filtered, queen=FALSE)
summary(wm_r_mexi_combined_filtered)## 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 links
The summary report above shows that there are 176 area units in Mexico City, Mexico State and Morelos State. The most connect area unit has 14 neighbours. There are 4 area units with 10 neighbours.
List complete weight matrix
## List of 176
## $ : int [1:5] 4 14 15 73 120
## $ : int [1:5] 6 9 11 12 13
## $ : int [1:4] 9 15 53 78
## $ : int [1:8] 1 14 16 36 49 74 120 125
## $ : int [1:5] 6 13 14 16 74
## $ : int [1:8] 2 5 10 12 13 74 86 138
## $ : int [1:3] 9 11 78
## $ : int [1:8] 10 11 12 41 66 150 161 164
## $ : int [1:7] 2 3 7 11 13 15 78
## $ : int [1:5] 6 8 12 41 138
## $ : int [1:9] 2 7 8 9 12 59 78 117 150
## $ : int [1:5] 2 6 8 10 11
## $ : int [1:6] 2 5 6 9 14 15
## $ : int [1:6] 1 4 5 13 15 16
## $ : int [1:7] 1 3 9 13 14 53 73
## $ : int [1:4] 4 5 14 74
## $ : int [1:4] 19 30 101 118
## $ : int [1:7] 27 44 49 97 108 109 116
## $ : int [1:4] 17 61 87 118
## $ : int [1:4] 37 96 113 133
## $ : int [1:7] 23 58 90 103 122 130 134
## $ : int [1:4] 28 89 114 117
## $ : int [1:7] 21 48 102 126 127 130 134
## $ : int [1:3] 96 98 121
## $ : int [1:4] 31 33 84 119
## $ : int [1:2] 52 112
## $ : int [1:6] 18 44 46 49 115 116
## $ : int [1:3] 22 89 117
## $ : int [1:6] 54 62 73 76 120 137
## $ : int [1:5] 17 64 72 101 118
## $ : int [1:5] 25 50 84 163 171
## $ : int [1:4] 77 81 91 100
## $ : int [1:5] 25 66 84 105 119
## $ : int [1:8] 43 70 71 88 89 106 117 122
## $ : int [1:3] 67 78 117
## $ : int [1:4] 4 49 124 125
## $ : int [1:8] 20 56 102 106 113 129 133 134
## $ : int [1:3] 41 99 119
## $ : int [1:4] 51 107 111 136
## $ : int [1:8] 69 75 107 111 124 125 136 137
## $ : int [1:8] 8 10 38 55 66 99 119 138
## $ : int [1:4] 61 72 118 128
## $ : int [1:4] 34 70 71 117
## $ : int [1:7] 18 27 46 85 109 115 116
## $ : int [1:4] 47 55 86 115
## $ : int [1:3] 27 44 115
## $ : int [1:4] 45 74 86 115
## $ : int [1:4] 23 57 126 127
## $ : int [1:10] 4 18 27 36 60 74 97 115 120 125
## $ : int [1:4] 31 157 163 171
## $ : int [1:4] 39 111 112 136
## $ : int [1:3] 26 112 136
## $ : int [1:5] 3 15 67 73 78
## $ : int [1:5] 29 62 76 83 103
## $ : int [1:6] 41 45 86 115 119 138
## $ : int [1:5] 37 123 129 133 135
## $ : int [1:3] 48 94 126
## $ : int [1:5] 21 63 64 90 103
## $ : int [1:4] 11 79 117 150
## $ : int [1:5] 49 75 125 136 141
## $ : int [1:5] 19 42 87 95 118
## $ : int [1:5] 29 54 73 83 131
## $ : int [1:6] 58 64 72 76 103 128
## $ : int [1:7] 30 58 63 72 80 90 101
## $ : int [1:6] 68 79 104 106 114 117
## $ : int [1:7] 8 33 41 105 110 164 168
## $ : int [1:10] 35 53 70 73 78 83 92 117 122 131
## $ : int [1:6] 65 79 104 135 146 156
## $ : int [1:3] 40 75 124
## $ : int [1:7] 34 43 67 71 92 117 122
## $ : int [1:3] 34 43 70
## $ : int [1:6] 30 42 63 64 118 128
## $ : int [1:8] 1 15 29 53 62 67 120 131
## $ : int [1:8] 4 5 6 16 47 49 86 115
## $ : int [1:8] 40 60 69 97 124 125 136 141
## $ : int [1:7] 29 54 63 103 111 128 137
## $ : int 32
## $ : int [1:8] 3 7 9 11 35 53 67 117
## $ : int [1:8] 59 65 68 104 117 148 150 156
## $ : int [1:4] 64 90 101 140
## $ : int [1:4] 32 91 108 109
## $ : int [1:4] 94 126 132 139
## $ : int [1:6] 54 62 67 103 122 131
## $ : int [1:6] 25 31 33 110 143 171
## $ : int [1:3] 44 109 115
## $ : int [1:6] 6 45 47 55 74 138
## $ : int [1:2] 19 61
## $ : int [1:4] 34 89 106 114
## $ : int [1:6] 22 28 34 88 114 117
## $ : int [1:6] 21 58 64 80 130 140
## $ : int [1:5] 32 81 100 108 109
## $ : int [1:3] 67 70 122
## $ : int [1:3] 98 102 113
## $ : int [1:3] 57 82 126
## $ : int 61
## $ : int [1:5] 20 24 98 113 133
## $ : int [1:7] 18 49 75 100 108 136 141
## $ : int [1:7] 24 93 96 102 113 132 139
## $ : int [1:4] 38 41 105 119
## [list output truncated]
## - attr(*, "class")= chr "nb"
## - attr(*, "region.id")= chr [1:176] "270" "271" "272" "273" ...
## - attr(*, "call")= language poly2nb(pl = mexi_combined_filtered, queen = FALSE)
## - attr(*, "type")= chr "rook"
## - attr(*, "sym")= logi TRUE
par(mfrow=c(1,2))
plot(mexi_combined_filtered, border="lightgrey")
plot(wm_q_mexi_combined_filtered, coordinates(mexi_combined_filtered), pch = 19, cex = 0.6, add = TRUE, col = "red", ain="Queen Contiguity")
plot(mexi_combined_filtered, border="lightgrey")
plot(wm_r_mexi_combined_filtered, coordinates(mexi_combined_filtered), pch = 19, cex = 0.6, add = TRUE, col = "red", ain="ROOK Contiguity") ## Computing distance based neighbours
To derive distance-based weight matrices by using dnearneigh() of spdep package.
The function identifies neighbours of region points by Euclidean distance with a distance band with lower d1= and upper d2= bounds controlled by the bounds= argument. If unprojected coordinates are used and either specified in the coordinates object x or with x as a two column matrix and longlat=TRUE, great circle distances in km will be calculated assuming the WGS84 reference ellipsoid.
coords <- coordinates(mexi_combined_filtered)
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.
Now, we will compute the distance weight matrix by using dnearneigh() as shown in the code chunk below.
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 30784
## Percentage nonzero weights: 99.38017
## Average number of links: 174.9091
Next, we will use str() to display the content of wm_d62 weight matrix.
## List of 176
## $ : int [1:175] 2 3 4 5 6 7 8 9 10 11 ...
## $ : int [1:175] 1 3 4 5 6 7 8 9 10 11 ...
## $ : int [1:173] 1 2 4 5 6 7 8 9 10 11 ...
## $ : int [1:175] 1 2 3 5 6 7 8 9 10 11 ...
## $ : int [1:174] 1 2 3 4 6 7 8 10 11 12 ...
## $ : int [1:175] 1 2 3 4 5 7 8 9 10 11 ...
## $ : int [1:175] 1 2 3 4 5 6 8 9 10 11 ...
## $ : int [1:175] 1 2 3 4 5 6 7 9 10 11 ...
## $ : int [1:174] 1 2 3 4 6 7 8 10 11 12 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 11 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:174] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:174] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:174] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:174] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:174] 1 2 3 4 5 6 7 8 9 10 ...
## $ : int [1:175] 1 2 3 4 5 6 7 8 9 10 ...
## [list output truncated]
## - attr(*, "class")= chr "nb"
## - attr(*, "nbtype")= chr "distance"
## - attr(*, "distances")= num [1:2] 0 19827
## - attr(*, "region.id")= chr [1:176] "1" "2" "3" "4" ...
## - attr(*, "call")= language dnearneigh(x = coords, d1 = 0, d2 = 19827, longlat = TRUE)
## - attr(*, "dnn")= num [1:2] 0 19827
## - attr(*, "bounds")= chr [1:2] "GT" "LE"
## - attr(*, "sym")= logi TRUE
Using table() and card() of spdep to display the structure of the weight matrix
##
## 173 174 175
## 09 1 2 13
## 15 0 8 117
## 17 1 2 32
## [1] 1
##
## 1
## 176
plot(mexi_combined_filtered, border="lightgrey")
plot(wm_d19827, coords, add=TRUE)
plot(k1, coords, add=TRUE, col="red", length=0.08) The red lines show the links of 1st nearest neighbours and the black lines show the links of neighbours within the cut-off distance of 19827km.
par(mfrow=c(1,2))
plot(mexi_combined_filtered, border="lightgrey")
plot(k1, coords, add=TRUE, col="red", length=0.08, main="1st nearest neighbours")
plot(mexi_combined_filtered, border="lightgrey")
plot(wm_d19827, coords, add=TRUE, pch = 19, cex = 0.6, main="Distance link") ### Computing adaptive distance weight matrix
One of the characteristics of fixed distance weight matrix is that more densely settled areas (usually the urban areas) tend to have more neighbours and the less densely settled areas (usually the rural counties) tend to have lesser neighbours. Having many neighbours smoothes the neighbour relationship across more neighbours.
It is possible to control the numbers of neighbours directly using k-nearest neighbours, either accepting asymmetric neighbours or imposing symmetry as shown in the code chunk below.
## 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
## List of 176
## $ : int [1:6] 4 14 15 16 73 120
## $ : int [1:6] 5 6 9 12 13 14
## $ : int [1:6] 7 9 15 53 73 78
## $ : int [1:6] 1 5 14 16 49 120
## $ : int [1:6] 2 6 13 14 16 74
## $ : int [1:6] 2 5 10 13 16 74
## $ : int [1:6] 2 3 9 11 53 78
## $ : int [1:6] 10 11 12 66 161 164
## $ : int [1:6] 2 3 7 13 15 53
## $ : int [1:6] 6 8 12 86 99 138
## $ : int [1:6] 2 7 9 12 59 150
## $ : int [1:6] 2 6 8 10 11 138
## $ : int [1:6] 2 5 9 14 15 16
## $ : int [1:6] 1 4 5 13 15 16
## $ : int [1:6] 1 9 13 14 16 73
## $ : int [1:6] 4 5 6 13 14 74
## $ : int [1:6] 19 30 64 87 101 118
## $ : int [1:6] 44 46 85 97 108 116
## $ : int [1:6] 17 30 61 87 101 118
## $ : int [1:6] 37 56 93 96 113 133
## $ : int [1:6] 23 58 103 122 130 134
## $ : int [1:6] 28 35 43 89 114 117
## $ : int [1:6] 48 102 126 127 130 134
## $ : int [1:6] 93 96 98 121 132 139
## $ : int [1:6] 31 33 50 84 105 119
## $ : int [1:6] 39 51 52 107 112 136
## $ : int [1:6] 18 44 46 49 85 116
## $ : int [1:6] 22 35 43 89 114 117
## $ : int [1:6] 1 62 73 120 125 137
## $ : int [1:6] 17 64 72 90 101 118
## $ : int [1:6] 25 33 50 84 110 163
## $ : int [1:6] 77 81 91 100 108 109
## $ : int [1:6] 25 38 66 84 99 105
## $ : int [1:6] 43 70 71 88 89 106
## $ : int [1:6] 22 28 43 78 114 117
## $ : int [1:6] 40 49 69 124 125 141
## $ : int [1:6] 20 56 104 106 113 129
## $ : int [1:6] 33 41 99 105 119 138
## $ : int [1:6] 40 51 60 69 107 111
## $ : int [1:6] 60 69 75 107 124 137
## $ : int [1:6] 38 55 99 105 119 138
## $ : int [1:6] 61 72 76 95 118 128
## $ : int [1:6] 28 70 71 88 89 92
## $ : int [1:6] 18 27 46 85 109 116
## $ : int [1:6] 47 55 74 86 115 138
## $ : int [1:6] 18 27 44 85 115 116
## $ : int [1:6] 6 27 45 74 86 115
## $ : int [1:6] 23 57 94 126 127 130
## $ : int [1:6] 4 27 36 116 125 141
## $ : int [1:6] 31 84 110 143 157 163
## $ : int [1:6] 39 40 69 107 111 112
## $ : int [1:6] 26 60 75 100 112 136
## $ : int [1:6] 3 7 9 62 73 131
## $ : int [1:6] 29 62 73 76 83 131
## $ : int [1:6] 38 41 45 86 115 119
## $ : int [1:6] 20 37 104 123 129 135
## $ : int [1:6] 48 82 94 126 127 132
## $ : int [1:6] 21 30 63 64 90 103
## $ : int [1:6] 7 22 28 35 78 117
## $ : int [1:6] 40 69 75 107 124 136
## $ : int [1:6] 17 19 42 72 95 118
## $ : int [1:6] 29 53 54 73 76 131
## $ : int [1:6] 54 58 72 83 103 128
## $ : int [1:6] 30 58 63 72 90 118
## $ : int [1:6] 22 28 88 89 104 114
## $ : int [1:6] 33 99 105 110 164 168
## $ : int [1:6] 35 53 78 83 92 131
## $ : int [1:6] 79 104 123 135 146 156
## $ : int [1:6] 40 60 75 107 124 141
## $ : int [1:6] 34 43 71 89 92 122
## $ : int [1:6] 34 43 70 88 89 92
## $ : int [1:6] 30 42 63 64 118 128
## $ : int [1:6] 1 15 29 53 62 120
## $ : int [1:6] 4 5 6 16 47 86
## $ : int [1:6] 40 60 69 124 136 141
## $ : int [1:6] 29 54 62 111 128 137
## $ : int [1:6] 32 81 91 100 108 109
## $ : int [1:6] 3 28 35 59 67 117
## $ : int [1:6] 65 68 114 117 148 150
## $ : int [1:6] 17 30 64 90 101 140
## $ : int [1:6] 32 77 85 91 108 109
## $ : int [1:6] 48 57 94 126 132 139
## $ : int [1:6] 54 62 63 67 103 131
## $ : int [1:6] 31 33 50 66 110 168
## $ : int [1:6] 18 44 46 109 115 116
## $ : int [1:6] 6 10 45 47 74 138
## $ : int [1:6] 17 19 61 95 101 118
## $ : int [1:6] 22 28 43 71 89 114
## $ : int [1:6] 22 28 43 71 88 114
## $ : int [1:6] 30 58 64 80 130 140
## $ : int [1:6] 18 32 77 81 100 108
## $ : int [1:6] 35 43 67 70 71 89
## $ : int [1:6] 20 98 102 113 126 132
## $ : int [1:6] 48 57 82 126 127 132
## $ : int [1:6] 19 42 61 72 118 128
## $ : int [1:6] 20 24 93 98 113 133
## $ : int [1:6] 18 60 75 108 116 141
## $ : int [1:6] 24 82 93 96 132 139
## $ : int [1:6] 33 38 41 66 105 138
## [list output truncated]
## - attr(*, "region.id")= chr [1:176] "1" "2" "3" "4" ...
## - attr(*, "call")= language knearneigh(x = coords, k = 6)
## - attr(*, "sym")= logi FALSE
## - attr(*, "type")= chr "knn"
## - attr(*, "knn-k")= num 6
## - attr(*, "class")= chr "nb"
Derive a spatial weight matrix based on Inversed Distance method. First, we will compute the distances between areas by using nbdists() of spdep.
dist <- nbdists(wm_q_mexi_combined_filtered, coords, longlat = TRUE)
ids <- lapply(dist, function(x) 1/(x))
ids## [[1]]
## [1] 1.625282e-04 7.286620e-05 7.107408e-05 3.382463e-04 7.435495e-05
##
## [[2]]
## [1] 5.811475e-05 6.758899e-05 5.642515e-05 7.017576e-05 8.656476e-05
##
## [[3]]
## [1] 7.892491e-04 6.549298e-05 1.741204e-04 6.317576e-05
##
## [[4]]
## [1] 1.625282e-04 8.303279e-05 6.389969e-05 7.233752e-05 1.165162e-04
## [6] 2.448648e-04 7.937698e-05 9.351276e-05
##
## [[5]]
## [1] 5.822100e-05 6.080187e-05 1.057847e-04 8.574530e-05 9.734377e-05
##
## [[6]]
## [1] 5.811475e-05 5.822100e-05 1.084275e-04 2.913919e-04 1.575909e-04
## [6] 1.447864e-04 7.795105e-05 7.082767e-05
##
## [[7]]
## [1] 9.168834e-05 1.253962e-04 7.378774e-05
##
## [[8]]
## [1] 8.051543e-05 9.717618e-05 1.392697e-04 8.476316e-05 1.178655e-04
## [6] 6.381720e-05 1.259769e-04 1.883358e-04
##
## [[9]]
## [1] 6.758899e-05 7.892491e-04 9.168834e-05 3.408606e-04 2.768554e-04
## [6] 6.487048e-05 6.829824e-05
##
## [[10]]
## [1] 1.084275e-04 8.051543e-05 9.885437e-05 2.440929e-04 1.367894e-04
##
## [[11]]
## [1] 5.642515e-05 1.253962e-04 9.717618e-05 3.408606e-04 2.877919e-04
## [6] 2.639822e-04 5.920533e-05 6.716088e-05 7.704263e-05
##
## [[12]]
## [1] 7.017576e-05 2.913919e-04 1.392697e-04 9.885437e-05 2.877919e-04
##
## [[13]]
## [1] 8.656476e-05 6.080187e-05 1.575909e-04 2.768554e-04 8.288005e-05
## [6] 8.457247e-05
##
## [[14]]
## [1] 7.286620e-05 8.303279e-05 1.057847e-04 8.288005e-05 8.401945e-05
## [6] 8.893802e-05
##
## [[15]]
## [1] 7.107408e-05 6.549298e-05 6.487048e-05 8.457247e-05 8.401945e-05
## [6] 5.500956e-05 8.881606e-05
##
## [[16]]
## [1] 6.389969e-05 8.574530e-05 8.893802e-05 5.586925e-05
##
## [[17]]
## [1] 7.073320e-04 1.435673e-04 5.669883e-05 1.063840e-04
##
## [[18]]
## [1] 1.356883e-04 7.528173e-05 1.276375e-04 9.455630e-05 6.162470e-05
## [6] 7.825048e-05 7.484700e-05
##
## [[19]]
## [1] 7.073320e-04 5.074741e-05 1.147090e-04 1.232539e-04
##
## [[20]]
## [1] 6.207069e-05 1.989998e-04 2.572725e-04 1.955575e-04
##
## [[21]]
## [1] 6.151062e-05 1.274074e-04 1.794158e-04 9.416121e-05 1.114368e-04
## [6] 7.500310e-05 1.043238e-04
##
## [[22]]
## [1] 7.501270e-05 7.591507e-05 2.928491e-04 2.055733e-04
##
## [[23]]
## [1] 6.151062e-05 1.547248e-04 7.573288e-05 6.228594e-05 1.206285e-04
## [6] 1.150039e-04 8.822855e-05
##
## [[24]]
## [1] 9.839306e-05 6.215702e-05 7.026734e-05
##
## [[25]]
## [1] 1.772273e-04 1.597936e-04 7.456176e-05 7.246207e-05
##
## [[26]]
## [1] 1.728270e-04 5.240886e-05
##
## [[27]]
## [1] 1.356883e-04 5.745815e-05 9.279165e-05 8.762809e-05 1.489506e-04
## [6] 6.348349e-05
##
## [[28]]
## [1] 7.501270e-05 8.074777e-05 1.046676e-04
##
## [[29]]
## [1] 1.092489e-04 1.355575e-04 1.246360e-04 6.797692e-05 1.155047e-04
## [6] 6.979176e-05
##
## [[30]]
## [1] 1.435673e-04 9.289630e-05 1.570112e-04 7.459317e-05 6.827966e-05
##
## [[31]]
## [1] 1.772273e-04 1.016815e-04 5.614833e-05 3.042400e-04 9.680763e-05
##
## [[32]]
## [1] 1.918575e-04 6.965538e-05 1.148537e-04 1.461145e-04
##
## [[33]]
## [1] 1.597936e-04 1.126621e-04 1.142536e-04 2.320960e-04 3.757265e-03
## [6] 5.013834e-05
##
## [[34]]
## [1] 5.359027e-04 5.982625e-05 6.981223e-05 7.253707e-05 6.435784e-05
## [6] 8.143763e-05 9.297839e-05 1.617230e-04
##
## [[35]]
## [1] 6.333917e-05 8.473366e-05 1.036850e-04
##
## [[36]]
## [1] 7.233752e-05 6.065400e-05 7.768842e-05 1.265698e-04
##
## [[37]]
## [1] 6.207069e-05 6.851713e-05 1.162834e-04 8.780240e-05 6.920613e-05
## [6] 1.108079e-04 6.175376e-05 6.794225e-05
##
## [[38]]
## [1] 7.531720e-05 5.348017e-05 5.986237e-05
##
## [[39]]
## [1] 1.161859e-04 1.205144e-04 8.003154e-05 1.749139e-04
##
## [[40]]
## [1] 6.310945e-05 1.502051e-04 5.043640e-05 6.274283e-05 6.583885e-05
## [6] 6.063476e-05 5.785814e-05 1.370508e-04
##
## [[41]]
## [1] 8.476316e-05 2.440929e-04 7.531720e-05 1.414978e-04 5.261504e-05
## [6] 1.544491e-04 1.341811e-04 1.053480e-04 9.422543e-05
##
## [[42]]
## [1] 2.958042e-04 7.717293e-05 1.130843e-04 1.658591e-04
##
## [[43]]
## [1] 5.359027e-04 6.367085e-05 6.277917e-05 9.400085e-05
##
## [[44]]
## [1] 7.528173e-05 5.745815e-05 9.571071e-05 7.616479e-05 1.400602e-04
## [6] 8.314604e-05 6.054047e-04
##
## [[45]]
## [1] 1.473610e-04 1.292529e-04 1.071633e-04 6.943647e-05
##
## [[46]]
## [1] 9.279165e-05 9.571071e-05 5.725829e-05
##
## [[47]]
## [1] 1.473610e-04 5.328770e-05 7.801567e-05 7.605560e-05
##
## [[48]]
## [1] 1.547248e-04 6.415608e-05 9.502806e-05 5.052029e-04
##
## [[49]]
## [1] 1.165162e-04 1.276375e-04 8.762809e-05 6.065400e-05 6.157729e-05
## [6] 2.138680e-04 6.349008e-05 8.240884e-05 1.846241e-04 6.681049e-05
## [11] 7.612655e-05
##
## [[50]]
## [1] 1.016815e-04 7.644401e-05 9.641867e-05 1.096279e-04
##
## [[51]]
## [1] 1.161859e-04 2.417046e-04 9.240346e-05 3.270855e-04
##
## [[52]]
## [1] 1.728270e-04 7.494096e-05 1.923145e-04
##
## [[53]]
## [1] 1.741204e-04 5.500956e-05 6.797303e-05 9.671140e-05 7.848568e-05
##
## [[54]]
## [1] 1.092489e-04 8.780628e-05 9.145785e-05 1.069529e-04 1.119878e-04
##
## [[55]]
## [1] 1.414978e-04 1.292529e-04 5.954553e-05 1.312002e-04 2.425602e-04
## [6] 6.763797e-05
##
## [[56]]
## [1] 6.851713e-05 1.426785e-04 1.322068e-04 4.494805e-04 1.397649e-04
##
## [[57]]
## [1] 6.415608e-05 1.996281e-04 8.649248e-05
##
## [[58]]
## [1] 1.274074e-04 7.998431e-05 1.014577e-04 7.556778e-05 2.063834e-04
##
## [[59]]
## [1] 2.639822e-04 4.687696e-04 7.735722e-05 6.138432e-05
##
## [[60]]
## [1] 6.157729e-05 5.644105e-05 1.253675e-04 1.859244e-04 1.274989e-04
## [6] 1.228858e-04
##
## [[61]]
## [1] 5.074741e-05 2.958042e-04 9.027376e-05 6.813486e-05 8.333363e-05
##
## [[62]]
## [1] 1.355575e-04 8.780628e-05 7.087537e-05 1.078882e-04 6.211540e-05
##
## [[63]]
## [1] 7.998431e-05 5.845233e-05 1.097591e-04 1.028617e-04 5.773333e-05
## [6] 1.110776e-04
##
## [[64]]
## [1] 9.289630e-05 1.014577e-04 5.845233e-05 7.342780e-05 6.392337e-05
## [6] 8.015172e-05 6.542939e-05
##
## [[65]]
## [1] 6.064992e-05 2.347901e-04 3.899482e-04 1.669993e-04 9.817634e-05
## [6] 7.966890e-05
##
## [[66]]
## [1] 1.178655e-04 1.126621e-04 5.261504e-05 5.840214e-05 6.909403e-05
## [6] 8.609568e-05 1.094156e-04 1.008985e-04 6.612251e-05
##
## [[67]]
## [1] 6.333917e-05 6.797303e-05 5.479950e-05 6.819981e-05 2.501014e-04
## [6] 3.845427e-04 8.021070e-05 1.622446e-04 4.930344e-04 1.382934e-04
##
## [[68]]
## [1] 6.064992e-05 7.694168e-05 6.040944e-05 7.272001e-05 4.471832e-04
## [6] 1.427255e-04
##
## [[69]]
## [1] 6.310945e-05 5.713328e-05 1.525564e-04
##
## [[70]]
## [1] 5.982625e-05 6.367085e-05 5.479950e-05 1.653097e-04 1.252556e-04
## [6] 8.165058e-05 5.846018e-05
##
## [[71]]
## [1] 6.981223e-05 6.277917e-05 1.653097e-04
##
## [[72]]
## [1] 1.570112e-04 7.717293e-05 1.097591e-04 7.342780e-05 1.147446e-04
## [6] 7.717540e-05
##
## [[73]]
## [1] 3.382463e-04 8.881606e-05 1.246360e-04 9.671140e-05 7.087537e-05
## [6] 6.819981e-05 6.120162e-05 1.023066e-04
##
## [[74]]
## [1] 2.448648e-04 9.734377e-05 1.447864e-04 5.586925e-05 5.328770e-05
## [6] 2.138680e-04 1.676493e-04 1.235205e-04
##
## [[75]]
## [1] 1.502051e-04 5.644105e-05 5.713328e-05 7.671582e-05 8.538918e-05
## [6] 5.975063e-05 7.195585e-05 7.104666e-05
##
## [[76]]
## [1] 6.797692e-05 9.145785e-05 1.028617e-04 8.037875e-05 8.662185e-05
## [6] 8.293435e-05 1.404467e-04
##
## [[77]]
## [1] 0.0001918575
##
## [[78]]
## [1] 6.317576e-05 7.378774e-05 6.829824e-05 5.920533e-05 8.473366e-05
## [6] 7.848568e-05 2.501014e-04 4.617344e-04
##
## [[79]]
## [1] 4.687696e-04 2.347901e-04 7.694168e-05 1.490246e-04 7.638188e-05
## [6] 3.000689e-04 6.315298e-05 8.856767e-05
##
## [[80]]
## [1] 6.392337e-05 8.794944e-05 1.108574e-04 6.711318e-05
##
## [[81]]
## [1] 6.965538e-05 9.734742e-05 2.243529e-04 1.687366e-04
##
## [[82]]
## [1] 1.529304e-04 1.173296e-04 8.839082e-05 4.549339e-04 6.245430e-05
##
## [[83]]
## [1] 1.069529e-04 1.078882e-04 3.845427e-04 5.643471e-05 2.420899e-04
## [6] 1.193623e-04
##
## [[84]]
## [1] 7.456176e-05 5.614833e-05 1.142536e-04 5.840214e-05 1.178126e-04
## [6] 1.427576e-04 9.904001e-05
##
## [[85]]
## [1] 7.616479e-05 7.671596e-05 2.505499e-04
##
## [[86]]
## [1] 7.795105e-05 1.071633e-04 7.801567e-05 5.954553e-05 1.676493e-04
## [6] 4.374973e-04
##
## [[87]]
## [1] 1.147090e-04 9.027376e-05
##
## [[88]]
## [1] 7.253707e-05 1.956579e-04 1.990647e-04 1.073083e-04
##
## [[89]]
## [1] 7.591507e-05 8.074777e-05 6.435784e-05 1.956579e-04 8.824106e-05
## [6] 7.768063e-05
##
## [[90]]
## [1] 1.794158e-04 7.556778e-05 8.015172e-05 8.794944e-05 5.870988e-05
## [6] 7.304769e-05
##
## [[91]]
## [1] 1.148537e-04 9.734742e-05 7.392934e-05 9.734426e-05 6.873991e-05
##
## [[92]]
## [1] 8.021070e-05 1.252556e-04 6.898998e-05
##
## [[93]]
## [1] 9.521408e-05 1.103750e-04 2.486247e-04
##
## [[94]]
## [1] 0.0001996281 0.0001529304 0.0001467189
##
## [[95]]
## [1] 6.813486e-05
##
## [[96]]
## [1] 1.989998e-04 9.839306e-05 1.654548e-04 1.710716e-04 1.119231e-04
##
## [[97]]
## [1] 9.455630e-05 6.349008e-05 1.253675e-04 7.671582e-05 8.285562e-05
## [6] 8.432010e-05 9.803931e-05 1.269201e-04
##
## [[98]]
## [1] 6.215702e-05 9.521408e-05 1.654548e-04 1.534867e-04 9.459651e-05
## [6] 9.267810e-05 9.259669e-05
##
## [[99]]
## [1] 5.348017e-05 1.544491e-04 6.909403e-05 7.184076e-05 2.280988e-04
##
## [[100]]
## [1] 1.461145e-04 7.392934e-05 8.285562e-05 6.810813e-05
##
## [[101]]
## [1] 5.669883e-05 7.459317e-05 6.542939e-05 1.108574e-04
##
## [[102]]
## [1] 7.573288e-05 1.162834e-04 1.173296e-04 1.103750e-04 1.534867e-04
## [6] 1.684752e-04 2.260939e-04 1.430287e-04 9.710852e-05
##
## [[103]]
## [1] 9.416121e-05 1.119878e-04 2.063834e-04 5.773333e-05 8.037875e-05
## [6] 5.643471e-05 5.936774e-05
##
## [[104]]
## [1] 3.899482e-04 6.040944e-05 1.490246e-04 2.629296e-04 6.504732e-05
## [6] 2.120856e-04
##
## [[105]]
## [1] 2.320960e-04 1.341811e-04 8.609568e-05 7.184076e-05 6.360556e-05
##
## [[106]]
## [1] 8.143763e-05 8.780240e-05 1.669993e-04 1.990647e-04 2.629296e-04
## [6] 7.205777e-05 7.535666e-05 5.727673e-05 2.571497e-04
##
## [[107]]
## [1] 0.0001205144 0.0000504364 0.0002379534 0.0003872600
##
## [[108]]
## [1] 6.162470e-05 2.243529e-04 9.734426e-05 8.432010e-05 6.810813e-05
## [6] 2.338874e-04
##
## [[109]]
## [1] 7.825048e-05 1.400602e-04 1.687366e-04 7.671596e-05 6.873991e-05
## [6] 2.338874e-04 6.080557e-05
##
## [[110]]
## [1] 3.757265e-03 1.094156e-04 1.178126e-04 8.142763e-05 6.656400e-05
##
## [[111]]
## [1] 8.003154e-05 6.274283e-05 2.417046e-04 8.662185e-05 2.379534e-04
## [6] 1.299429e-04 1.105673e-04
##
## [[112]]
## [1] 5.240886e-05 9.240346e-05 7.494096e-05 1.227695e-04
##
## [[113]]
## [1] 2.572725e-04 6.920613e-05 2.486247e-04 1.710716e-04 9.459651e-05
## [6] 1.684752e-04
##
## [[114]]
## [1] 2.928491e-04 9.817634e-05 1.073083e-04 8.824106e-05 7.205777e-05
## [6] 3.641001e-04
##
## [[115]]
## [1] 1.489506e-04 8.314604e-05 6.943647e-05 5.725829e-05 7.605560e-05
## [6] 8.240884e-05 1.312002e-04 1.235205e-04 2.505499e-04 6.080557e-05
##
## [[116]]
## [1] 7.484700e-05 6.348349e-05 6.054047e-04
##
## [[117]]
## [1] 6.716088e-05 2.055733e-04 1.046676e-04 9.297839e-05 1.036850e-04
## [6] 9.400085e-05 7.735722e-05 7.966890e-05 1.622446e-04 8.165058e-05
## [11] 4.617344e-04 7.638188e-05 7.768063e-05 3.641001e-04
##
## [[118]]
## [1] 1.063840e-04 1.232539e-04 6.827966e-05 1.130843e-04 8.333363e-05
## [6] 1.147446e-04
##
## [[119]]
## [1] 7.246207e-05 5.013834e-05 5.986237e-05 1.053480e-04 2.425602e-04
## [6] 2.280988e-04 6.360556e-05
##
## [[120]]
## [1] 7.435495e-05 7.937698e-05 1.155047e-04 1.846241e-04 6.120162e-05
## [6] 9.217043e-05 6.066269e-05
##
## [[121]]
## [1] 7.026734e-05
##
## [[122]]
## [1] 1.114368e-04 1.617230e-04 4.930344e-04 5.846018e-05 2.420899e-04
## [6] 6.898998e-05 5.936774e-05 7.535666e-05 6.060526e-05
##
## [[123]]
## [1] 0.0001426785 0.0000709919
##
## [[124]]
## [1] 7.768842e-05 6.583885e-05 1.525564e-04 8.538918e-05 1.768949e-04
##
## [[125]]
## [1] 9.351276e-05 1.265698e-04 6.063476e-05 6.681049e-05 1.859244e-04
## [6] 5.975063e-05 9.217043e-05 1.768949e-04 8.113619e-05 3.582431e-04
##
## [[126]]
## [1] 6.228594e-05 9.502806e-05 8.649248e-05 8.839082e-05 1.467189e-04
## [6] 2.260939e-04 9.371913e-05
##
## [[127]]
## [1] 1.206285e-04 5.052029e-04 5.969137e-05 6.386635e-05
##
## [[128]]
## [1] 1.658591e-04 1.110776e-04 7.717540e-05 8.293435e-05 1.299429e-04
##
## [[129]]
## [1] 1.108079e-04 1.322068e-04 6.504732e-05 5.727673e-05 6.794786e-05
##
## [[130]]
## [1] 7.500310e-05 1.150039e-04 5.870988e-05 5.969137e-05 2.783276e-04
##
## [[131]]
## [1] 0.0000621154 0.0001382934 0.0001023066 0.0001193623
##
## [[132]]
## [1] 4.549339e-04 9.267810e-05 1.430287e-04 9.371913e-05 5.634864e-05
##
## [[133]]
## [1] 1.955575e-04 6.175376e-05 4.494805e-04 1.119231e-04
##
## [[134]]
## [1] 1.043238e-04 8.822855e-05 6.794225e-05 9.710852e-05 2.571497e-04
## [6] 6.060526e-05
##
## [[135]]
## [1] 1.397649e-04 7.272001e-05 2.120856e-04 7.099190e-05 6.794786e-05
## [6] 6.353800e-05
##
## [[136]]
## [1] 1.749139e-04 5.785814e-05 3.270855e-04 1.923145e-04 1.274989e-04
## [6] 7.195585e-05 9.803931e-05 3.872600e-04 1.227695e-04
##
## [[137]]
## [1] 6.979176e-05 1.370508e-04 1.404467e-04 1.105673e-04 6.066269e-05
## [6] 8.113619e-05
##
## [[138]]
## [1] 7.082767e-05 1.367894e-04 9.422543e-05 6.763797e-05 4.374973e-04
##
## [[139]]
## [1] 6.245430e-05 9.259669e-05 5.634864e-05
##
## [[140]]
## [1] 6.711318e-05 7.304769e-05 6.386635e-05 2.783276e-04
##
## [[141]]
## [1] 7.612655e-05 1.228858e-04 7.104666e-05 1.269201e-04 3.582431e-04
##
## [[142]]
## [1] 1.241413e-04 8.539523e-05 6.672324e-05 6.983264e-05
##
## [[143]]
## [1] 1.427576e-04 8.142763e-05 5.544160e-05 5.829138e-05 3.362260e-04
## [6] 1.781684e-04 1.924510e-04
##
## [[144]]
## [1] 1.243936e-04 9.994883e-05 9.129165e-05
##
## [[145]]
## [1] 3.281936e-04 1.928098e-04 5.898352e-04 3.083782e-04 7.644762e-05
## [6] 7.580708e-05 9.554232e-05 6.357697e-05 7.440387e-05
##
## [[146]]
## [1] 4.471832e-04 6.353800e-05 1.241413e-04 6.544388e-05 1.824225e-04
## [6] 1.365450e-04
##
## [[147]]
## [1] 5.544160e-05 3.281936e-04 7.409936e-05 6.044252e-05
##
## [[148]]
## [1] 3.000689e-04 1.656710e-04 5.674852e-05 1.314863e-04 8.146334e-05
## [6] 6.177564e-05 3.450862e-04
##
## [[149]]
## [1] 1.656710e-04 1.640892e-04 7.444357e-05 2.661668e-04 3.303072e-04
## [6] 1.041609e-04
##
## [[150]]
## [1] 6.381720e-05 7.704263e-05 6.138432e-05 6.315298e-05 5.674852e-05
## [6] 6.730039e-05
##
## [[151]]
## [1] 1.243936e-04 1.928098e-04 2.863397e-04 5.493352e-05
##
## [[152]]
## [1] 0.0001314863 0.0001640892 0.0001973662 0.0002709325
##
## [[153]]
## [1] 7.283566e-05 8.508572e-05 2.643966e-04 8.697488e-05
##
## [[154]]
## [1] 9.994883e-05 5.898352e-04 2.863397e-04 2.126969e-04
##
## [[155]]
## [1] 8.539523e-05 6.544388e-05 9.720004e-05 2.421669e-04 1.030413e-04
## [6] 1.613814e-04
##
## [[156]]
## [1] 1.427255e-04 8.856767e-05 1.824225e-04 8.146334e-05 9.720004e-05
## [6] 1.918939e-04 1.342323e-04
##
## [[157]]
## [1] 7.644401e-05 1.385325e-04 7.110111e-05 2.916675e-04
##
## [[158]]
## [1] 6.672324e-05 7.283566e-05 2.421669e-04 7.108129e-05 1.008078e-04
## [6] 6.897220e-05
##
## [[159]]
## [1] 6.177564e-05 7.444357e-05 1.918939e-04 6.121227e-05 1.094021e-04
##
## [[160]]
## [1] 9.129165e-05 3.083782e-04 2.126969e-04 9.969477e-05
##
## [[161]]
## [1] 1.259769e-04 3.450862e-04 6.730039e-05 1.973662e-04 8.760308e-05
## [6] 6.111340e-05 1.766843e-04
##
## [[162]]
## [1] 6.983264e-05 1.365450e-04 1.030413e-04
##
## [[163]]
## [1] 3.042400e-04 9.641867e-05 1.385325e-04 1.668645e-04
##
## [[164]]
## [1] 1.883358e-04 1.008985e-04 8.760308e-05 2.016158e-04 9.397873e-05
##
## [[165]]
## [1] 7.644762e-05 2.661668e-04 2.636170e-04 2.315066e-04 1.279681e-04
## [6] 8.041173e-05 2.298781e-04
##
## [[166]]
## [1] 7.580708e-05 8.508572e-05 7.108129e-05 9.969477e-05 2.636170e-04
## [6] 9.518302e-05
##
## [[167]]
## [1] 5.829138e-05 6.111340e-05 2.016158e-04 6.464456e-05 5.803879e-05
##
## [[168]]
## [1] 6.612251e-05 6.656400e-05 3.362260e-04 9.397873e-05 6.464456e-05
##
## [[169]]
## [1] 3.303072e-04 6.121227e-05 2.315066e-04 1.018157e-04 1.452857e-04
##
## [[170]]
## [1] 1.781684e-04 9.554232e-05 7.409936e-05 1.041609e-04 2.709325e-04
## [6] 1.766843e-04 1.279681e-04 5.803879e-05
##
## [[171]]
## [1] 9.680763e-05 1.096279e-04 9.904001e-05 1.924510e-04 6.357697e-05
## [6] 6.044252e-05 7.110111e-05 8.893700e-05 1.067492e-04
##
## [[172]]
## [1] 2.643966e-04 8.041173e-05 9.518302e-05 9.999673e-05
##
## [[173]]
## [1] 0.0002916675 0.0001668645 0.0000889370 0.0000619717
##
## [[174]]
## [1] 7.440387e-05 5.493352e-05 1.067492e-04 6.197170e-05
##
## [[175]]
## [1] 0.0001613814 0.0001342323 0.0001008078 0.0001094021 0.0001018157
## [6] 0.0001244095
##
## [[176]]
## [1] 8.697488e-05 6.897220e-05 2.298781e-04 1.452857e-04 9.999673e-05
## [6] 1.244095e-04
Next, we need to assign weights to each neighboring polygon. In our case, each neighboring polygon will be assigned equal weight (style=“W”). This is accomplished by assigning the fraction 1/(#ofneighbors) to each neighboring county then summing the weighted income values. While this is the most intuitive way to summaries the neighbors’ values it has one drawback in that polygons along the edges of the study area will base their lagged values on fewer polygons thus potentially over- or under-estimating the true nature of the spatial autocorrelation in the data. For this example, we’ll stick with the style=“W” option for simplicity’s sake but note that other more robust options are available, notably style=“B”.
rswm_q_mexi_combined_filtered <- nb2listw(wm_q_mexi_combined_filtered, style="W", zero.policy = TRUE)
rswm_q_mexi_combined_filtered## 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
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 176
## Number of nonzero links: 962
## Percentage nonzero weights: 3.10563
## Average number of links: 5.465909
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 176 30976 0.1314899 0.0001067155 0.0005739984
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.014e-05 7.083e-05 9.422e-05 1.367e-04 1.461e-04 3.757e-03
Using spatial lag variables with row-standardized weights,
COVID19RATE.lag <- lag.listw(rswm_q_mexi_combined_filtered, mexi_combined_filtered$rate_of_covid)
COVID19RATE.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.677386
Append the spatially lag COVID19RATE.lag values onto mexi_combined_filtered SpatialPolygonDataFrame by using the code chunk below.
lag.list <- list(mexi_combined_filtered$NOMGEO, lag.listw(rswm_q_mexi_combined_filtered, mexi_combined_filtered$rate_of_covid))
lag.res <- as.data.frame(lag.list)
colnames(lag.res) <- c("NOMGEO", "lag COVID 19 Rate")
mexi_combined_filtered@data <- left_join(mexi_combined_filtered@data,lag.res)## CVEGEO CVE_ENT CVE_MUN NOMGEO Pop2010 Pop2020 new13 new14
## 1 09002 09 002 Azcapotzalco 414711 408441 18 27
## 2 09003 09 003 Coyoacán 620416 621952 15 19
## 3 09004 09 004 Cuajimalpa de Morelos 186391 199809 18 9
## 4 09005 09 005 Gustavo A. Madero 1185772 1176967 21 68
## 5 09006 09 006 Iztacalco 384326 393821 9 16
## 6 09007 09 007 Iztapalapa 1815786 1815551 32 59
## 7 09008 09 008 La Magdalena Contreras 239086 245147 9 12
## 8 09009 09 009 Milpa Alta 130582 139371 0 1
## 9 09010 09 010 Ã\201lvaro Obregón 727034 755537 18 34
## 10 09011 09 011 Tláhuac 360265 366586 2 8
## new15 new16 new17 new18 new19 new20 new21 new22 new23 new24 new25 new26
## 1 41 69 98 137 224 316 336 323 349 310 345 328
## 2 63 106 162 179 277 319 342 374 389 299 382 272
## 3 18 29 35 45 58 63 91 159 141 137 140 126
## 4 97 204 337 373 546 677 778 735 831 692 611 525
## 5 38 123 155 209 254 319 305 276 283 220 245 206
## 6 125 301 535 693 863 966 1076 873 884 818 747 665
## 7 18 38 86 51 74 133 144 142 126 122 181 179
## 8 22 36 78 70 94 119 195 163 183 223 170 137
## 9 56 124 142 197 238 322 385 387 383 458 435 367
## 10 19 49 96 135 203 303 352 233 303 237 279 277
## new27 new28 new29 new30 new31 new32 cumul13 cumul14 cumul15 cumul16 cumul17
## 1 316 301 300 332 271 24 19 46 87 156 254
## 2 382 313 490 427 396 30 25 44 107 213 375
## 3 168 159 171 220 114 7 66 75 93 122 157
## 4 573 512 625 663 505 59 28 96 193 397 734
## 5 165 172 186 233 196 62 11 27 65 188 343
## 6 762 670 744 728 620 135 41 100 225 526 1061
## 7 208 229 345 247 277 12 15 27 45 83 169
## 8 176 151 204 197 145 11 0 1 23 59 137
## 9 388 419 564 492 425 62 56 90 146 270 412
## 10 295 312 351 326 285 38 3 11 30 79 175
## cumul18 cumul19 cumul20 cumul21 cumul22 cumul23 cumul24 cumul25 cumul26
## 1 391 615 931 1267 1590 1939 2249 2594 2922
## 2 554 831 1150 1492 1866 2255 2554 2936 3208
## 3 202 260 323 414 573 714 851 991 1117
## 4 1107 1653 2330 3108 3843 4674 5366 5977 6502
## 5 552 806 1125 1430 1706 1989 2209 2454 2660
## 6 1754 2617 3583 4659 5532 6416 7234 7981 8646
## 7 220 294 427 571 713 839 961 1142 1321
## 8 207 301 420 615 778 961 1184 1354 1491
## 9 609 847 1169 1554 1941 2324 2782 3217 3584
## 10 310 513 816 1168 1401 1704 1941 2220 2497
## cumul27 cumul28 cumul29 cumul30 cumul31 cumul32 activ13 activ14 activ15
## 1 3238 3539 3839 4171 4442 4466 34 70 130
## 2 3590 3903 4393 4820 5216 5246 41 82 166
## 3 1285 1444 1615 1835 1949 1956 70 86 112
## 4 7075 7587 8212 8875 9380 9439 68 173 321
## 5 2825 2997 3183 3416 3612 3674 20 58 129
## 6 9408 10078 10822 11550 12170 12305 75 174 427
## 7 1529 1758 2103 2350 2627 2639 24 41 70
## 8 1667 1818 2022 2219 2364 2375 1 9 43
## 9 3972 4391 4955 5447 5872 5934 80 135 215
## 10 2792 3104 3455 3781 4066 4104 8 29 61
## activ16 activ17 activ18 activ19 activ20 activ21 activ22 activ23 activ24
## 1 199 332 516 776 1126 1456 1795 2143 2451
## 2 321 513 731 1054 1398 1751 2117 2496 2841
## 3 140 182 237 298 385 498 660 801 917
## 4 571 940 1438 2015 2755 3557 4275 5115 5756
## 5 286 482 722 1008 1328 1618 1893 2147 2394
## 6 854 1492 2315 3207 4207 5208 6095 6979 7723
## 7 113 206 280 386 507 653 786 925 1057
## 8 107 186 288 388 548 737 921 1129 1304
## 9 368 545 758 1053 1403 1807 2172 2651 3058
## 10 142 257 441 695 1023 1317 1579 1881 2111
## activ25 activ26 activ27 activ28 activ29 activ30 activ31 activ32 death13
## 1 2810 3117 3444 3746 4034 4340 4460 4466 0
## 2 3169 3496 3862 4227 4681 5075 5245 5246 0
## 3 1080 1181 1384 1540 1677 1885 1954 1956 0
## 4 6329 6842 7459 8022 8636 9182 9427 9439 2
## 5 2591 2766 2953 3137 3369 3553 3666 3674 1
## 6 8434 9100 9874 10528 11287 11950 12283 12305 2
## 7 1221 1408 1626 1877 2170 2446 2638 2639 0
## 8 1481 1626 1787 1964 2170 2322 2375 2375 0
## 9 3471 3812 4286 4744 5284 5705 5924 5934 2
## 10 2408 2681 3016 3325 3664 3940 4101 4104 0
## death14 death15 death16 death17 death18 death19 death20 death21 death22
## 1 2 4 6 15 14 38 33 44 47
## 2 1 7 8 19 11 26 30 36 34
## 3 2 1 3 5 1 4 11 8 9
## 4 5 26 19 37 74 82 98 117 114
## 5 1 3 7 26 20 30 46 49 37
## 6 6 15 33 76 115 127 174 176 127
## 7 1 4 3 0 9 4 16 12 5
## 8 0 0 1 3 5 7 4 3 7
## 9 2 11 11 23 29 31 50 40 60
## 10 1 1 2 5 16 14 14 22 17
## death23 death24 death25 death26 death27 death28 death29 death30 death31
## 1 44 43 47 38 27 19 31 20 19
## 2 41 34 37 33 31 25 23 20 16
## 3 11 18 6 5 5 9 10 7 8
## 4 121 148 112 77 72 56 57 52 33
## 5 29 35 26 22 19 25 23 17 12
## 6 120 110 103 63 67 65 52 29 29
## 7 19 10 10 7 12 8 5 5 6
## 8 4 6 7 11 4 2 2 0 2
## 9 56 61 41 32 39 26 28 20 24
## 10 20 17 24 12 11 8 10 3 8
## death32 wk13_32_cumul_count wk13_32_activ_count wk13_32_death_count
## 1 8 4465 41445 499
## 2 6 5236 48512 438
## 3 1 1908 17043 124
## 4 18 9432 92320 1320
## 5 3 3672 37794 431
## 6 17 12296 124517 1506
## 7 3 2633 21073 139
## 8 0 2375 21761 68
## 9 8 5896 53405 594
## 10 2 4103 36783 207
## rate_of_covid active_rate_of_covid death_rate_of_covid lag COVID 19 Rate
## 1 109.31811 1014.7120 12.217187 65.17350
## 2 84.18656 779.9959 7.042344 83.91900
## 3 95.49119 852.9646 6.205927 56.33468
## 4 80.13819 784.3890 11.215268 58.84120
## 5 93.24033 959.6746 10.944058 69.31115
## 6 67.72600 685.8359 8.295002 69.83156
## 7 107.40494 859.6067 5.670067 66.89785
## 8 170.40848 1561.3722 4.879064 54.46100
## 9 78.03721 706.8482 7.861958 78.95109
## 10 111.92462 1003.3935 5.646697 86.62845
Plotting both the COVID19Rate and spatial lag COVID19Rate for comparison using the code chunk below.
covid <- qtm(mexi_combined_filtered, "rate_of_covid")
lag_covid <- qtm(mexi_combined_filtered, "lag COVID 19 Rate")
tmap_arrange(covid, lag_covid, asp=1, ncol=2) # Measure of Global Spatial Autocorrelation
Perform Moran’s I statistics testing using moran.test() of spdep.
moran.test(mexi_combined_filtered$rate_of_covid, listw=rswm_q_mexi_combined_filtered, zero.policy = TRUE, na.action=na.omit)##
## Moran I test under randomisation
##
## data: mexi_combined_filtered$rate_of_covid
## weights: rswm_q_mexi_combined_filtered
##
## 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 Moran I statistic shows that the value of 0.495940195 is a positive spatial autocorrelation that of a sign of Clustering. The p value we got is lesser than the targeted 0.05 does we reject the null hypothesis where distribution of COVID 19 rate is randomly distributed.
The code chunk below performs permutation test for Moran’s I statistic by using moran.mc() of spdep.
set.seed(1234)
bperm= moran.mc(mexi_combined_filtered$rate_of_covid, listw=rswm_q_mexi_combined_filtered, nsim=999, zero.policy = TRUE, na.action=na.omit)
bperm##
## Monte-Carlo simulation of Moran I
##
## data: mexi_combined_filtered$rate_of_covid
## weights: rswm_q_mexi_combined_filtered
## number of simulations + 1: 1000
##
## statistic = 0.49594, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
## [1] -0.006244441
## [1] 0.002138681
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.133572 -0.039006 -0.009037 -0.006244 0.023391 0.193374
## Computing local Moran’s I
To compute local Moran’s I, the localmoran() function of spdep will be used. It computes Ii values, given a set of zi values and a listw object providing neighbour weighting information for the polygon associated with the zi values.
fips <- order(mexi_combined_filtered$NOMGEO)
localMI <- localmoran(mexi_combined_filtered$rate_of_covid, rswm_q_mexi_combined_filtered)
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-09
localmoran() function returns a matrix of values whose columns are:
Ii: the local Moran’s I statistics E.Ii: the expectation of local moran statistic under the randomisation hypothesis Var.Ii: the variance of local moran statistic under the randomisation hypothesis Z.Ii:the standard deviate of local moran statistic Pr(): the p-value of local moran statistic The code chunk below list the content of the local Moran matrix derived by using printCoefmat().
printCoefmat(data.frame(localMI[fips,], row.names=mexi_combined_filtered$NOMGEO[fips]), check.names=FALSE)## Ii E.Ii Var.Ii Z.Ii
## Ã\201lvaro Obregón 3.13141221 -0.00571429 0.13118547 8.66142143
## Acambay de RuÃz Castañeda 0.26390635 -0.00571429 0.23344785 0.55803082
## Acolman 0.01426135 -0.00571429 0.13118547 0.05515157
## Aculco 0.56632612 -0.00571429 0.23344785 1.18394565
## Almoloya de Alquisiras 0.30401140 -0.00571429 0.23344785 0.64103579
## Almoloya de Juárez 0.13581163 -0.00571429 0.13118547 0.39074473
## Almoloya del RÃo 0.67470812 -0.00571429 0.23344785 1.40826267
## Amacuzac 0.32183542 -0.00571429 0.23344785 0.67792598
## Amanalco -0.01484814 -0.00571429 0.13118547 -0.02521803
## Amatepec 0.65595393 -0.00571429 0.31298526 1.18271017
## Amecameca -0.05132161 -0.00571429 0.23344785 -0.09439296
## Apaxco 0.19887632 -0.00571429 0.47206007 0.29777420
## Atenco 0.05565255 -0.00571429 0.15391044 0.15642267
## Atizapán 0.74051905 -0.00571429 0.31298526 1.33386755
## Atizapán de Zaragoza 0.02399602 -0.00571429 0.15391044 0.07573090
## Atlacomulco -0.11524758 -0.00571429 0.18572540 -0.25416192
## Atlatlahucan 0.21459988 -0.00571429 0.13118547 0.60827443
## Atlautla 0.24540192 -0.00571429 0.18572540 0.58269205
## Axapusco 0.09806448 -0.00571429 0.23344785 0.21478976
## Axochiapan 0.05576287 -0.00571429 0.31298526 0.10988839
## Ayala 0.18529881 -0.00571429 0.10088550 0.60137970
## Ayapango 0.12507698 -0.00571429 0.15391044 0.33338396
## Azcapotzalco 3.71578995 -0.00571429 0.18572540 8.63540803
## Benito Juárez 2.04673506 -0.00571429 0.15391044 5.23164671
## Calimaya 0.02624586 -0.00571429 0.11414173 0.09459899
## Capulhuac -0.00302759 -0.00571429 0.31298526 0.00480237
## Chalco 0.92752879 -0.00571429 0.10088550 2.93819355
## Chapa de Mota 0.46452020 -0.00571429 0.23344785 0.97323907
## Chapultepec 0.23978916 -0.00571429 0.23344785 0.50811575
## Chiautla 0.12379509 -0.00571429 0.13118547 0.35756776
## Chicoloapan 0.05317621 -0.00571429 0.23344785 0.12188501
## Chiconcuac 0.09787837 -0.00571429 0.31298526 0.18516847
## Chimalhuacán 0.00927770 -0.00571429 0.23344785 0.03102874
## Coacalco de Berriozábal 0.04804869 -0.00571429 0.23344785 0.11127263
## Coatepec Harinas 0.45174810 -0.00571429 0.11414173 1.35404518
## Coatetelco 0.55220367 -0.00571429 0.15391044 1.42212017
## Coatlán del RÃo 0.30183344 -0.00571429 0.15391044 0.78393217
## Cocotitlán 0.43055795 -0.00571429 0.31298526 0.77982227
## Coyoacán 3.90730541 -0.00571429 0.18572540 9.07980204
## Coyotepec -0.04554039 -0.00571429 0.23344785 -0.08242763
## Cuajimalpa de Morelos 2.26893414 -0.00571429 0.23344785 4.70781452
## Cuauhtémoc 3.75433271 -0.00571429 0.15391044 9.58427430
## Cuautitlán -0.13585213 -0.00571429 0.11414173 -0.38519565
## Cuautitlán Izcalli -0.02486671 -0.00571429 0.15391044 -0.04881909
## Cuautla -0.06787939 -0.00571429 0.23344785 -0.12866242
## Cuernavaca 0.11106901 -0.00571429 0.13118547 0.32243180
## Donato Guerra 0.37624778 -0.00571429 0.23344785 0.79054264
## Ecatepec de Morelos 0.04146713 -0.00571429 0.08160370 0.16516423
## Ecatzingo 0.44738058 -0.00571429 0.23344785 0.93776541
## El Oro -0.09683642 -0.00571429 0.23344785 -0.18859446
## Emiliano Zapata 0.34151463 -0.00571429 0.15391044 0.88507861
## Gustavo A. Madero 1.89986950 -0.00571429 0.11414173 5.64034689
## Huehuetoca 0.06620238 -0.00571429 0.23344785 0.14884512
## Hueypoxtla 0.08465995 -0.00571429 0.31298526 0.16154098
## Huitzilac -0.64879100 -0.00571429 0.15391044 -1.63918793
## Huixquilucan 0.10760397 -0.00571429 0.18572540 0.26294459
## Isidro Fabela 0.21351661 -0.00571429 0.18572540 0.50870513
## Ixtapaluca -0.00068109 -0.00571429 0.15391044 0.01282951
## Ixtapan de la Sal 0.34631182 -0.00571429 0.18572540 0.81684417
## Ixtapan del Oro 0.42327324 -0.00571429 0.31298526 0.76680110
## Ixtlahuaca 0.24341734 -0.00571429 0.18572540 0.57808702
## Iztacalco 3.31053141 -0.00571429 0.18572540 7.69504290
## Iztapalapa 1.98037802 -0.00571429 0.11414173 5.87864445
## Jaltenco 0.00755830 -0.00571429 0.15391044 0.03383153
## Jantetelco 0.42730045 -0.00571429 0.23344785 0.89620578
## Jilotepec 0.31007892 -0.00571429 0.18572540 0.73276907
## Jilotzingo 0.01164928 -0.00571429 0.18572540 0.04029055
## Jiquipilco 0.36485488 -0.00571429 0.15391044 0.94457237
## Jiutepec 0.28604550 -0.00571429 0.23344785 0.60385198
## Jocotitlán 0.08604397 -0.00571429 0.13118547 0.25333913
## Jojutla 0.02376588 -0.00571429 0.23344785 0.06101477
## Jonacatepec de Leandro Valle 0.23448221 -0.00571429 0.23344785 0.49713202
## Joquicingo 0.09431862 -0.00571429 0.15391044 0.25498160
## Juchitepec -0.43747374 -0.00571429 0.10088550 -1.35933808
## La Magdalena Contreras 3.80785505 -0.00571429 0.31298526 6.81662978
## La Paz -0.02838077 -0.00571429 0.15391044 -0.05777635
## Lerma 0.01477010 -0.00571429 0.09028051 0.06817512
## Luvianos 0.10721642 -0.00571429 0.31298526 0.20185992
## Malinalco 0.27576759 -0.00571429 0.15391044 0.71749089
## Mazatepec 0.35469114 -0.00571429 0.15391044 0.91866524
## Melchor Ocampo 0.07343547 -0.00571429 0.31298526 0.14147758
## Metepec 0.18947611 -0.00571429 0.13118547 0.53890919
## Mexicaltzingo 0.44736022 -0.00571429 0.31298526 0.80985578
## Miacatlán 0.32815993 -0.00571429 0.13118547 0.92180705
## Miguel Hidalgo 2.88083576 -0.00571429 0.13118547 7.96959463
## Milpa Alta 4.54149654 -0.00571429 0.11414173 13.45931183
## Morelos 0.33742780 -0.00571429 0.15391044 0.87466137
## Naucalpan de Juárez 0.30194898 -0.00571429 0.11414173 0.91065402
## Nextlalpan -0.03645473 -0.00571429 0.11414173 -0.09098878
## Nezahualcóyotl 0.50115932 -0.00571429 0.11414173 1.50029770
## Nicolás Romero 0.12350691 -0.00571429 0.13118547 0.35677210
## Nopaltepec 0.43630779 -0.00571429 0.94928452 0.45367595
## Ocoyoacac -0.00101117 -0.00571429 0.11414173 0.01392076
## Ocuilan 0.16230965 -0.00571429 0.11414173 0.49733488
## Ocuituco 0.45881920 -0.00571429 0.23344785 0.96143979
## Otumba 0.00078201 -0.00571429 0.23344785 0.01344531
## Otzoloapan 0.30792007 -0.00571429 0.18572540 0.72775965
## Otzolotepec 0.01036922 -0.00571429 0.15391044 0.04099650
## Ozumba 0.11314263 -0.00571429 0.13118547 0.32815694
## Papalotla 0.09712461 -0.00571429 0.31298526 0.18382114
## Polotitlán 0.59015434 -0.00571429 0.47206007 0.86726512
## Puente de Ixtla -0.12795927 -0.00571429 0.15391044 -0.31159969
## Rayón 0.00363773 -0.00571429 0.23344785 0.01935576
## San Antonio la Isla 0.09152224 -0.00571429 0.15391044 0.24785370
## San Felipe del Progreso 0.25841641 -0.00571429 0.15391044 0.67326313
## San José del Rincón 0.45265169 -0.00571429 0.23344785 0.94867496
## San MartÃn de las Pirámides 0.00009449 -0.00571429 0.18572540 0.01347873
## San Mateo Atenco 0.09725230 -0.00571429 0.31298526 0.18404938
## San Simón de Guerrero 0.02798668 -0.00571429 0.31298526 0.06023936
## Santo Tomás 0.25920407 -0.00571429 0.31298526 0.47353284
## Soyaniquilpan de Juárez 0.22277735 -0.00571429 0.94928452 0.23451580
## Sultepec 0.51615315 -0.00571429 0.18572540 1.21094536
## Tecámac 0.00022609 -0.00571429 0.11414173 0.01758296
## Tejupilco -0.06464077 -0.00571429 0.13118547 -0.16269256
## Temamatla 0.20128425 -0.00571429 0.18572540 0.48032104
## Temascalapa -0.03381126 -0.00571429 0.23344785 -0.05815200
## Temascalcingo 0.05679423 -0.00571429 0.23344785 0.12937319
## Temascaltepec 0.23541759 -0.00571429 0.10088550 0.75917213
## Temixco 0.40424972 -0.00571429 0.18572540 0.95128373
## Temoac 0.63410025 -0.00571429 0.23344785 1.32421702
## Temoaya 0.18005254 -0.00571429 0.13118547 0.51289127
## Tenancingo -0.01925520 -0.00571429 0.15391044 -0.03451549
## Tenango del Aire -0.05415416 -0.00571429 0.18572540 -0.11240026
## Tenango del Valle 0.00271679 -0.00571429 0.10088550 0.02654414
## Teoloyucan -0.00619852 -0.00571429 0.23344785 -0.00100221
## Teotihuacán -0.09974196 -0.00571429 0.15391044 -0.23967441
## Tepalcingo 0.14441834 -0.00571429 0.23344785 0.31072782
## Tepetlaoxtoc 0.03383955 -0.00571429 0.13118547 0.10920582
## Tepetlixpa 0.27109306 -0.00571429 0.18572540 0.64230596
## Tepotzotlán 0.03150028 -0.00571429 0.13118547 0.10274724
## Tepoztlán -0.32562726 -0.00571429 0.13118547 -0.88326088
## Tequixquiac 0.12717342 -0.00571429 0.23344785 0.27503620
## Tetecala 0.40067870 -0.00571429 0.31298526 0.72641410
## Tetela del Volcán 0.47880979 -0.00571429 0.23344785 1.00281408
## Texcaltitlán 0.23391857 -0.00571429 0.15391044 0.61081869
## Texcalyacac -0.00817712 -0.00571429 0.15391044 -0.00627771
## Texcoco 0.03057381 -0.00571429 0.09028051 0.12077226
## Tezoyuca 0.02884594 -0.00571429 0.31298526 0.06177527
## Tianguistenco -0.13626139 -0.00571429 0.06301054 -0.52006860
## Timilpan 0.31583391 -0.00571429 0.15391044 0.81961905
## Tláhuac 6.24857236 -0.00571429 0.18572540 14.51249650
## Tlalmanalco 0.10702064 -0.00571429 0.13118547 0.31125448
## Tlalnepantla -0.24518797 -0.00571429 0.18572540 -0.55567664
## Tlalnepantla de Baz 0.25768033 -0.00571429 0.13118547 0.72721701
## Tlalpan 3.53735045 -0.00571429 0.10088550 11.15487509
## Tlaltizapán de Zapata 0.21236172 -0.00571429 0.13118547 0.60209502
## Tlaquiltenango 0.05996626 -0.00571429 0.15391044 0.16741823
## Tlatlaya 1.09905005 -0.00571429 0.94928452 1.13389136
## Tlayacapan -0.15746976 -0.00571429 0.18572540 -0.35213461
## Toluca -0.03825691 -0.00571429 0.10088550 -0.10245618
## Tonanitla 0.01782151 -0.00571429 0.18572540 0.05461266
## Tonatico 0.17498718 -0.00571429 0.47206007 0.26300441
## Totolapan 0.27301347 -0.00571429 0.18572540 0.64676210
## Tultepec -0.00134689 -0.00571429 0.18572540 0.01013414
## Tultitlán -0.01468854 -0.00571429 0.09028051 -0.02986768
## Valle de Bravo -0.02563639 -0.00571429 0.13118547 -0.05500377
## Valle de Chalco Solidaridad -0.53484420 -0.00571429 0.18572540 -1.22779726
## Venustiano Carranza 3.77273368 -0.00571429 0.23344785 7.82021168
## Villa de Allende 0.52745224 -0.00571429 0.23344785 1.10348881
## Villa del Carbón 0.43882908 -0.00571429 0.18572540 1.03152197
## Villa Guerrero 0.28384665 -0.00571429 0.18572540 0.67189951
## Villa Victoria 0.43743371 -0.00571429 0.18572540 1.02828412
## Xalatlaco -0.03817267 -0.00571429 0.23344785 -0.06717876
## Xochimilco 9.32722299 -0.00571429 0.18572540 21.65622193
## Xochitepec 0.40194830 -0.00571429 0.18572540 0.94594350
## Xonacatlán -0.00951967 -0.00571429 0.23344785 -0.00787596
## Xoxocotla 0.23390315 -0.00571429 0.15391044 0.61077938
## Yautepec 0.26680939 -0.00571429 0.11414173 0.80664418
## Yecapixtla 0.30399490 -0.00571429 0.10088550 0.97507879
## Zacatepec -0.18739238 -0.00571429 0.23344785 -0.37601712
## Zacazonapan 0.22850884 -0.00571429 0.18572540 0.54349321
## Zacualpan 0.65163149 -0.00571429 0.23344785 1.36050123
## Zacualpan de Amilpas 0.59318036 -0.00571429 0.23344785 1.23952559
## Zinacantepec 0.00901322 -0.00571429 0.15391044 0.03754009
## Zumpahuacán 0.19379577 -0.00571429 0.15391044 0.50854659
## Zumpango -0.01007570 -0.00571429 0.10088550 -0.01373133
## Pr.z...0.
## Ã\201lvaro Obregón 0.0000
## Acambay de RuÃz Castañeda 0.2884
## Acolman 0.4780
## Aculco 0.1182
## Almoloya de Alquisiras 0.2607
## Almoloya de Juárez 0.3480
## Almoloya del RÃo 0.0795
## Amacuzac 0.2489
## Amanalco 0.5101
## Amatepec 0.1185
## Amecameca 0.5376
## Apaxco 0.3829
## Atenco 0.4378
## Atizapán 0.0911
## Atizapán de Zaragoza 0.4698
## Atlacomulco 0.6003
## Atlatlahucan 0.2715
## Atlautla 0.2801
## Axapusco 0.4150
## Axochiapan 0.4562
## Ayala 0.2738
## Ayapango 0.3694
## Azcapotzalco 0.0000
## Benito Juárez 0.0000
## Calimaya 0.4623
## Capulhuac 0.4981
## Chalco 0.0017
## Chapa de Mota 0.1652
## Chapultepec 0.3057
## Chiautla 0.3603
## Chicoloapan 0.4515
## Chiconcuac 0.4265
## Chimalhuacán 0.4876
## Coacalco de Berriozábal 0.4557
## Coatepec Harinas 0.0879
## Coatetelco 0.0775
## Coatlán del RÃo 0.2165
## Cocotitlán 0.2177
## Coyoacán 0.0000
## Coyotepec 0.5328
## Cuajimalpa de Morelos 0.0000
## Cuauhtémoc 0.0000
## Cuautitlán 0.6500
## Cuautitlán Izcalli 0.5195
## Cuautla 0.5512
## Cuernavaca 0.3736
## Donato Guerra 0.2146
## Ecatepec de Morelos 0.4344
## Ecatzingo 0.1742
## El Oro 0.5748
## Emiliano Zapata 0.1881
## Gustavo A. Madero 0.0000
## Huehuetoca 0.4408
## Hueypoxtla 0.4358
## Huitzilac 0.9494
## Huixquilucan 0.3963
## Isidro Fabela 0.3055
## Ixtapaluca 0.4949
## Ixtapan de la Sal 0.2070
## Ixtapan del Oro 0.2216
## Ixtlahuaca 0.2816
## Iztacalco 0.0000
## Iztapalapa 0.0000
## Jaltenco 0.4865
## Jantetelco 0.1851
## Jilotepec 0.2318
## Jilotzingo 0.4839
## Jiquipilco 0.1724
## Jiutepec 0.2730
## Jocotitlán 0.4000
## Jojutla 0.4757
## Jonacatepec de Leandro Valle 0.3095
## Joquicingo 0.3994
## Juchitepec 0.9130
## La Magdalena Contreras 0.0000
## La Paz 0.5230
## Lerma 0.4728
## Luvianos 0.4200
## Malinalco 0.2365
## Mazatepec 0.1791
## Melchor Ocampo 0.4437
## Metepec 0.2950
## Mexicaltzingo 0.2090
## Miacatlán 0.1783
## Miguel Hidalgo 0.0000
## Milpa Alta 0.0000
## Morelos 0.1909
## Naucalpan de Juárez 0.1812
## Nextlalpan 0.5362
## Nezahualcóyotl 0.0668
## Nicolás Romero 0.3606
## Nopaltepec 0.3250
## Ocoyoacac 0.4944
## Ocuilan 0.3095
## Ocuituco 0.1682
## Otumba 0.4946
## Otzoloapan 0.2334
## Otzolotepec 0.4836
## Ozumba 0.3714
## Papalotla 0.4271
## Polotitlán 0.1929
## Puente de Ixtla 0.6223
## Rayón 0.4923
## San Antonio la Isla 0.4021
## San Felipe del Progreso 0.2504
## San José del Rincón 0.1714
## San MartÃn de las Pirámides 0.4946
## San Mateo Atenco 0.4270
## San Simón de Guerrero 0.4760
## Santo Tomás 0.3179
## Soyaniquilpan de Juárez 0.4073
## Sultepec 0.1130
## Tecámac 0.4930
## Tejupilco 0.5646
## Temamatla 0.3155
## Temascalapa 0.5232
## Temascalcingo 0.4485
## Temascaltepec 0.2239
## Temixco 0.1707
## Temoac 0.0927
## Temoaya 0.3040
## Tenancingo 0.5138
## Tenango del Aire 0.5447
## Tenango del Valle 0.4894
## Teoloyucan 0.5004
## Teotihuacán 0.5947
## Tepalcingo 0.3780
## Tepetlaoxtoc 0.4565
## Tepetlixpa 0.2603
## Tepotzotlán 0.4591
## Tepoztlán 0.8115
## Tequixquiac 0.3916
## Tetecala 0.2338
## Tetela del Volcán 0.1580
## Texcaltitlán 0.2707
## Texcalyacac 0.5025
## Texcoco 0.4519
## Tezoyuca 0.4754
## Tianguistenco 0.6985
## Timilpan 0.2062
## Tláhuac 0.0000
## Tlalmanalco 0.3778
## Tlalnepantla 0.7108
## Tlalnepantla de Baz 0.2335
## Tlalpan 0.0000
## Tlaltizapán de Zapata 0.2736
## Tlaquiltenango 0.4335
## Tlatlaya 0.1284
## Tlayacapan 0.6376
## Toluca 0.5408
## Tonanitla 0.4782
## Tonatico 0.3963
## Totolapan 0.2589
## Tultepec 0.4960
## Tultitlán 0.5119
## Valle de Bravo 0.5219
## Valle de Chalco Solidaridad 0.8902
## Venustiano Carranza 0.0000
## Villa de Allende 0.1349
## Villa del Carbón 0.1511
## Villa Guerrero 0.2508
## Villa Victoria 0.1519
## Xalatlaco 0.5268
## Xochimilco 0.0000
## Xochitepec 0.1721
## Xonacatlán 0.5031
## Xoxocotla 0.2707
## Yautepec 0.2099
## Yecapixtla 0.1648
## Zacatepec 0.6465
## Zacazonapan 0.2934
## Zacualpan 0.0868
## Zacualpan de Amilpas 0.1076
## Zinacantepec 0.4850
## Zumpahuacán 0.3055
## Zumpango 0.5055
Before mapping the local Moran’s I map, it is wise to append the local Moran’s I dataframe (i.e. localMI) onto hunan SpatialPolygonDataFrame. The code chunks below can be used to perform the task. The out SpatialPolygonDataFrame is called mexi_combined_filtered.localMI.
Using choropleth mapping functions of tmap package, to plot the local Moran’s I values by using the code chinks below.
tm_shape(mexi_combined_filtered.localMI) +
tm_fill(col = "Ii",
style = "pretty",
palette = "RdBu",
title = "local moran statistics") +
tm_borders(alpha = 0.5) ### Mapping local Moran’s I p-values The choropleth shows there is evidence for both positive and negative Ii values. However, it is useful to consider the p-values for each of these values, as consider above.
The code chunks below produce a choropleth map of Moran’s I p-values by using functions of tmap package.
tm_shape(mexi_combined_filtered.localMI) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5) Mapping both local Moran’s I values and p-values For effective interpretation, it is better to plot both the local Moran’s I values map and its correpence p-values map next to each other.
The code chunk below will be used to create such visualisation.
localMI.map <- tm_shape(mexi_combined_filtered.localMI) +
tm_fill(col = "Ii",
style = "pretty",
title = "local moran statistics") +
tm_borders(alpha = 0.5)
pvalue.map <- tm_shape(mexi_combined_filtered.localMI) +
tm_fill(col = "Pr.z...0.",
breaks=c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf),
palette="-Blues",
title = "local Moran's I p-values") +
tm_borders(alpha = 0.5)
tmap_arrange(localMI.map, pvalue.map, asp=1, ncol=2) ## LISA Cluster Map
The LISA Cluster Map shows the significant locations color coded by type of spatial autocorrelation. The first step before we can generate the LISA cluster map is to plot the Moran scatterplot.
Plotting Moran scatterplot The Moran scatterplot is an illustration of the relationship between the values of the chosen attribute at each location and the average value of the same attribute at neighboring locations.
nci <- moran.plot(mexi_combined_filtered$rate_of_covid, rswm_q_mexi_combined_filtered, labels=as.character(mexi_combined_filtered$NOMGEO), xlab="COVID 19 Rate (Pert 10,000)", ylab="Spatially Lag COVID 19 Rate") The pattern could be observed where the higher COVID 19 rate is being clustered together with the high COVID 19 rate area, while the lower rate are directly in the opposite.
Scaling the variable of interest First we will use scale() to centers and scales the variable. here centering is done by subtracting the mean (omitting NAs) the corresponding columns, and scaling is done by dividing the (centered) variable by their standard deviations.
The as.vector() added to the end is to make sure that the data type we get out of this is a vector, that map neatly into the dataframe.
Now, we are ready to plot the Moran scatterplot again by using the code chunk below.
nci2 <- moran.plot(mexi_combined_filtered$Z.rate_of_covid, rswm_q_mexi_combined_filtered, labels=as.character(mexi_combined_filtered$NOMGEO), xlab="z-rate_of_COVID19 (Per 10,000)", ylab="Spatially Lag z-rate_of_COVID19")The code chunks below show the steps to prepare a LISA cluster map. Next, we centers the variable of interest around its mean. This is follow by centering the local Moran’s around the mean. Next, we will set a statistical significance level for the local Moran. These four command lines define the high-high, low-low, low-high and high-low categories. Lastly, places non-significant Moran in the category 0.
quadrant <- vector(mode="numeric",length=nrow(localMI))
DV <- mexi_combined_filtered$rate_of_covid - mean(mexi_combined_filtered$rate_of_covid)
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] <- 0Now, we can build the LISA map by using the code chunks below.
mexi_combined_filtered.localMI$quadrant <- quadrant
colors <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusters <- c("insignificant", "low-low", "low-high", "high-low", "high-high")
tm_shape(mexi_combined_filtered.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) # Hot Spot and Cold Spot Area Analysis
Beside detecting cluster and outliers, localised spatial statistics can be also used to detect hot spot and/or cold spot areas.
The term ‘hot spot’ has been used generically across disciplines to describe a region or value that is higher relative to its surroundings (Lepers et al 2005, Aben et al 2012, Isobe et al 2015).
An alternative spatial statistics to detect spatial anomalies is the Getis and Ord’s G-statistics (Getis and Ord, 1972; Ord and Getis, 1995). It looks at neighbours within a defined proximity to identify where either high or low values clutser spatially. Here, statistically significant hot-spots are recognised as areas of high values where other areas within a neighbourhood range also share high values too.
The analysis consists of three steps:
First, we need to define a new set of neighbours. Whist the spatial autocorrelation considered units which shared borders, for Getis-Ord we are defining neighbours based on distance or also known as proximity.
There are two type of distance-based proximity matrix, they are:
But since fixed distance weight matrix is not so effective as seen or demonstrate from the Local Moran’s I result. Therefore, we will proceed to use adaptive approach.
Instead of using a fixed distance to deirved the proximity matrix, a fixed number of neighbours can be also used to derived the proximity matrix. This is called adaptive distance method.
The code chunk below will be used to derive an adaptive distance weight matrix.
coords <- coordinates(mexi_combined_filtered)
#knb <- knn2nb(knearneigh(coords, k=8, longlat = TRUE), row.names=row.names(mexi_combined_filtered$rate_of_covid))
knn6_lw <- nb2listw(knn6, style = 'B')
summary(knn6_lw)## Characteristics of weights list object:
## 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
## Link number distribution:
##
## 6
## 176
## 176 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 6 links
## 176 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 with 6 links
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 176 30976 1056 1874 26004
Next, we will plot the centroids of the counties and neighbours by using the code chunk below.
plot(mexi_combined_filtered, border="lightgrey")
plot(knn6, coords, pch = 19, cex = 0.6, add = TRUE, col = "red") In the choropleth map above, counties shaded in red are the hot spot areas and counties shaded in blue are the cold spot areas. The darkness of the colours representing the intensity of the Gi values.
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.
Gi statistics using adaptive distance
The code chunk below are used to compute the Gi values for Rate of COVID 19 by using an adaptive distance weight matrix (i.e knb_lw).
fips <- order(mexi_combined_filtered$NOMGEO)
gi.adaptive <- localG(mexi_combined_filtered$rate_of_covid, knn6_lw)
mexi_combined_filtered.gi <- cbind(mexi_combined_filtered, as.matrix(gi.adaptive))
names(mexi_combined_filtered.gi)[95] <- "gstat_adaptive"tm_shape(mexi_combined_filtered.gi) +
tm_fill(col = "gstat_adaptive",
style = "pretty",
palette="-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5) + tm_text("NOMGEO", size=0.2) As we can observed the result above, the hotspot are being highlighted in red and the cold spots have a color marked as blue. The color describe the density and the weight of the Gi values. Which allows us to better analyse the Moran’s I method. The area of Mexico States and City are being highlighted with red as it has a considerably high clustered population than the Morelos State area. But looking at Mexico States, the area are wide apart with the COVID spread around the corner while the critical zone were near the Mexico City. And even though Mexico City is small but the population are higher due to the area in measurement is smaller compared to Mexico States. Therefor it greatly affect the influence of spreading the COVID.
At last, we are able to identify the result from the outcome after analyze the Study Area with various techniques and methods like EDA, Thematic mapping, Weight Matrix and spatial analysis. With a higher intense COVID 19 range of 109.3 to 170.4 in Mexico City and have lesser municipalities as compared to Mexico State and Morelos State. It could be due to higher densely population and are more closely connected that increase the population to falls into COVID 19.