Geoestadistic analysis of COVID cases distribution in Buenos Aires by age / Distribución de casos de COVID por edad Ciudad de Buenos Aires
Data: cases report since March 2020 until October 2021 from the Department of public health/ Data de datos publicos del ministerio de salud Marzo 2020- Octubre 2021
Ines Patop
25/10/2021
1 English
1.1 Introduction
I plan to analyze the age distribution of those tested for COVID in the Autonomous City of Buenos Aires (CABA) accumulated as of October 25, 2021.
1.2 Packages to use
library(sf)
library(tmap)
library(tidyverse)
library(spdep)1.3 Data loading and exploratory analysis
1.3.0.1 Geographic data - Administrative units
As was done in class, I will use the National Institute of Census and Statistics data and filter the departments of CABA and create a new variable that keeps only the commune number and another with the number of inhabitants in number format.
deptos <- st_read("pxdptodatosok.shp")## Reading layer `pxdptodatosok' from data source
## `/Users/inespatop/Documents/resume/clase2/pxdptodatosok.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 527 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.02985 ymin: -90 xmax: -25.02314 ymax: -21.74506
## Geodetic CRS: WGS 84depcaba <- deptos %>%
filter(codpcia == "02") %>%
mutate(departamento_id = substr(link,3, 5), .keep ="unused") %>%
mutate(pob = as.numeric(personas), .keep ="unused") %>%
select(departamento_id, departamen, pob)
depcabaI make a plot just to see the map
qtm(depcaba, fill = "pob",text = "departamento_id")
1.3.0.2 Loading data from COVID tests in the City of Buenos Aires
Since the data is very large, I load an already filtered dataset from the dataset only for CABA
##covid <- read_csv("./Covid19Cases.csv")
##glimpse(covid)
##covid_caba<-covid[covid$residence_province_id == "02",]
##rm(covid)
##save(covid_caba,"covid_caba.RData")
load("covid_caba.RData")
head(covid_caba)names(covid_caba)## [1] "id_evento_caso" "sexo"
## [3] "edad" "edad_años_meses"
## [5] "residencia_pais_nombre" "residencia_provincia_nombre"
## [7] "residencia_departamento_nombre" "carga_provincia_nombre"
## [9] "fecha_inicio_sintomas" "fecha_apertura"
## [11] "sepi_apertura" "fecha_internacion"
## [13] "cuidado_intensivo" "fecha_cui_intensivo"
## [15] "fallecido" "fecha_fallecimiento"
## [17] "asistencia_respiratoria_mecanica" "carga_provincia_id"
## [19] "origen_financiamiento" "clasificacion"
## [21] "clasificacion_resumen" "residencia_provincia_id"
## [23] "fecha_diagnostico" "residencia_departamento_id"
## [25] "ultima_actualizacion"##3 Exploring data and generating summary measures
I now make a summary measure of the number of cases by age, separating by confirmed and discarded cases.
covid_age <- covid_caba %>%
group_by( clasificacion_resumen) %>%
select(residencia_departamento_nombre, residencia_departamento_id, clasificacion_resumen, edad ) %>%
group_by(residencia_departamento_id,clasificacion_resumen) %>%
summarise(Edad_media = mean(edad,na.rm = T) )%>%
filter (residencia_departamento_id != "0")%>%
filter (clasificacion_resumen != "Sospechoso")## `summarise()` has grouped output by 'residencia_departamento_id'. You can
## override using the `.groups` argument.covid_ageI look at the distribution of the data both before and after generating the summary measure
#Total
ggplot(covid_caba%>% filter (residencia_departamento_id != "0")%>% filter (clasificacion_resumen != "Sospechoso"), aes(y=edad,x=residencia_departamento_id, fill=clasificacion_resumen))+geom_boxplot()+labs(title = "Distribución de edad acumulada por comuna")## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
#medida resumen
ggplot(covid_age, aes(y=Edad_media,x=clasificacion_resumen, fill=clasificacion_resumen))+geom_boxplot()+labs(title = "Distribución de edad promerdio acumulada")
1.3.0.3 Union with geographic data and graph of average age tested
Union
coviddepcaba <- left_join(depcaba, covid_age, by = c("departamento_id"="residencia_departamento_id"))Mean age in confirmed cases
qtm(coviddepcaba%>% filter (clasificacion_resumen == "Confirmado"), fill = "Edad_media")
Mean age in discarded cases
qtm(coviddepcaba%>% filter (clasificacion_resumen == "Descartado"), fill = "Edad_media")
At first glance, it can already be seen that in the southern zone there is more testing of young people and a greater number of positive tests on young people.
Now I make a graph by quantiles
tm_shape(coviddepcaba) +
tm_fill("Edad_media", style = "quantile", n = 5, palette = "Reds") +
tm_borders(alpha = 0.1) +
tm_layout(main.title = "Casos confirmados al 22/10/2021", main.title.size = 0.7 ,
legend.position = c("right", "bottom"), legend.title.size = 0.8)
## Neighbor Analysis
1.3.0.4 List of neighbors (of class nb) from polygons with poly2nb
For the rest of the work I will use only the data from positive covid tests.
coviddepcaba_pos<- coviddepcaba%>% filter (clasificacion_resumen == "Confirmado")I calculate the neighbors using poly2nb as it was done in class-
w <- poly2nb(coviddepcaba_pos, row.names = "departamento_id")I explore the neighbors
summary(w)## Neighbour list object:
## Number of regions: 15
## Number of nonzero links: 64
## Percentage nonzero weights: 28.44444
## Average number of links: 4.266667
## Link number distribution:
##
## 3 4 5 6 7
## 6 4 2 1 2
## 6 least connected regions:
## 1 8 9 10 12 13 with 3 links
## 2 most connected regions:
## 5 7 with 7 linksI map the neighbors
plot(st_geometry(coviddepcaba_pos), border="grey")
plot(w, coordinates(as(coviddepcaba_pos, "Spatial")), add=TRUE, col="blue")
title("Vecinos calculados con poly2nb")
View neighboring weight arrays
In this case the binary weights were done using style='B' as an option in particular this generates a matrix of 0 and 1 in which each neighbor has the same weight.
wmatb <- nb2mat(w, style='B')
wmatb## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0
## 2 1 0 1 0 1 0 0 0 0 0 0 0 0 1
## 3 1 1 0 1 1 0 0 0 0 0 0 0 0 0
## 4 1 0 1 0 1 0 1 1 0 0 0 0 0 0
## 5 0 1 1 1 0 1 1 0 0 0 0 0 0 1
## 6 0 0 0 0 1 0 1 0 0 0 1 0 0 0
## 7 0 0 0 1 1 1 0 1 1 1 1 0 0 0
## 8 0 0 0 1 0 0 1 0 1 0 0 0 0 0
## 9 0 0 0 0 0 0 1 1 0 1 0 0 0 0
## 10 0 0 0 0 0 0 1 0 1 0 1 0 0 0
## 11 0 0 0 0 0 1 1 0 0 1 0 1 0 0
## 12 0 0 0 0 0 0 0 0 0 0 1 0 1 0
## 13 0 0 0 0 0 0 0 0 0 0 0 1 0 1
## 14 0 1 0 0 1 0 0 0 0 0 0 0 1 0
## 15 0 0 0 0 1 1 0 0 0 0 1 1 1 1
## [,15]
## 1 0
## 2 0
## 3 0
## 4 0
## 5 1
## 6 1
## 7 0
## 8 0
## 9 0
## 10 0
## 11 1
## 12 1
## 13 1
## 14 1
## 15 0
## attr(,"call")
## nb2mat(neighbours = w, style = "B")For standardized weights per row (style=‘W’) instead each neighbor has a proportional weight. For example 3 neighbors will have a weight of 0.33
wmatW <- nb2mat(w, style='W')
wmatW## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## 1 0.00 0.3333333 0.3333333 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.25 0.0000000 0.2500000 0.0000000 0.2500000 0.0000000 0.0000000 0.0000000
## 3 0.25 0.2500000 0.0000000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000
## 4 0.20 0.0000000 0.2000000 0.0000000 0.2000000 0.0000000 0.2000000 0.2000000
## 5 0.00 0.1428571 0.1428571 0.1428571 0.0000000 0.1428571 0.1428571 0.0000000
## 6 0.00 0.0000000 0.0000000 0.0000000 0.2500000 0.0000000 0.2500000 0.0000000
## 7 0.00 0.0000000 0.0000000 0.1428571 0.1428571 0.1428571 0.0000000 0.1428571
## 8 0.00 0.0000000 0.0000000 0.3333333 0.0000000 0.0000000 0.3333333 0.0000000
## 9 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3333333 0.3333333
## 10 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000
## 11 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.2000000 0.2000000 0.0000000
## 12 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 13 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 14 0.00 0.2500000 0.0000000 0.0000000 0.2500000 0.0000000 0.0000000 0.0000000
## 15 0.00 0.0000000 0.0000000 0.0000000 0.1666667 0.1666667 0.0000000 0.0000000
## [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.2500000 0.0000000
## 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1428571 0.1428571
## 6 0.0000000 0.0000000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000
## 7 0.1428571 0.1428571 0.1428571 0.0000000 0.0000000 0.0000000 0.0000000
## 8 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 9 0.0000000 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 10 0.3333333 0.0000000 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000
## 11 0.0000000 0.2000000 0.0000000 0.2000000 0.0000000 0.0000000 0.2000000
## 12 0.0000000 0.0000000 0.3333333 0.0000000 0.3333333 0.0000000 0.3333333
## 13 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000 0.3333333 0.3333333
## 14 0.0000000 0.0000000 0.0000000 0.0000000 0.2500000 0.0000000 0.2500000
## 15 0.0000000 0.0000000 0.1666667 0.1666667 0.1666667 0.1666667 0.0000000
## attr(,"call")
## nb2mat(neighbours = w, style = "W")Let’s do the same now in list format
lwb <- nb2listw(w, style='B')
lwb## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 15
## Number of nonzero links: 64
## Percentage nonzero weights: 28.44444
## Average number of links: 4.266667
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 15 225 64 128 1208head(lwb)## $style
## [1] "B"
##
## $neighbours
## Neighbour list object:
## Number of regions: 15
## Number of nonzero links: 64
## Percentage nonzero weights: 28.44444
## Average number of links: 4.266667
##
## $weights
## $weights[[1]]
## [1] 1 1 1
##
## $weights[[2]]
## [1] 1 1 1 1
##
## $weights[[3]]
## [1] 1 1 1 1
##
## $weights[[4]]
## [1] 1 1 1 1 1
##
## $weights[[5]]
## [1] 1 1 1 1 1 1 1
##
## $weights[[6]]
## [1] 1 1 1 1
##
## $weights[[7]]
## [1] 1 1 1 1 1 1 1
##
## $weights[[8]]
## [1] 1 1 1
##
## $weights[[9]]
## [1] 1 1 1
##
## $weights[[10]]
## [1] 1 1 1
##
## $weights[[11]]
## [1] 1 1 1 1 1
##
## $weights[[12]]
## [1] 1 1 1
##
## $weights[[13]]
## [1] 1 1 1
##
## $weights[[14]]
## [1] 1 1 1 1
##
## $weights[[15]]
## [1] 1 1 1 1 1 1
##
## attr(,"mode")
## [1] "binary"
## attr(,"B")
## [1] TRUE1.4 Exploration of Spatial Association/Correlation Indices
1.4.0.1 Moran I Global
# todos los test
moran.test(coviddepcaba_pos$Edad_media, lwb)##
## Moran I test under randomisation
##
## data: coviddepcaba_pos$Edad_media
## weights: lwb
##
## Moran I statistic standard deviate = 2.1678, p-value = 0.01509
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.20243360 -0.07142857 0.01596030We see that the Global Moran’s i is significant (p value < 0.05) and since it is greater than zero even though it is slight (0.2), we have a grouped distribution of the data. This correlates with what we were seeing on the charts.
1.4.0.1.0.1 Moran Correlogram
Until the third neighbor we see the i of Morán
Icorr3 <-sp.correlogram(neighbours=w,var=coviddepcaba_pos$Edad_media,order=3,method="I",zero.policy=TRUE)
plot(Icorr3)
We see here that with the first neighbor there is a higher correlation than with the second and so
1.4.0.1.0.2 Moran scatterplot
mp <- moran.plot(coviddepcaba_pos$Edad_media, lwb, zero.policy=TRUE, labels=as.character(coviddepcaba_pos$"departamento_id"))
mpIn this case we see that there is an important or influential value, which is clearly the low age values, which are in the southern part of the city. However when we look at the rest of the values we see a lot of data, in fact many values in the anticorrelation zone. As expected, the slope is very low (i de moran).
Just as a reminder: on the x-axis we have the average age of covid-positive covid tests in each game. On the y-axis we see the spatial lag of x which is the average of values of neighbors to x. The slope of the regression line is Moran I.
1.4.0.2 Moran local Ii
We calculate the local Moran index. For this we use the localmoran function
The domain of Moran’s i is -1 to 1 and we now hope to see if there are areas in which there is autocorrelation. I then look for the values with p value less than 0.05
propc_lm_lwb <- localmoran(coviddepcaba_pos$Edad_media, lwb)
propc_lm_lwb_df<-as.data.frame(propc_lm_lwb)
propc_lm_lwb_dfpropc_lm_lwb_df[propc_lm_lwb_df$`Pr(z != E(Ii))`<0.05,]Ii: Moran i local E.Ii: local i Moran expect Var.Ii: variance of the local i Moran statistic Z.Ii: standard deviation of the local i Moran Pr(): p-value of the local i Moran
We see that for commune 4 and 8, the neighbors have similar values and for commune 15 dissimilar values. If we see this again on the map, we see that commune 4 and 8 have similar ages and are in the south of the city, while commune 15 is surrounded by lower age values than its own.
qtm(coviddepcaba_pos, fill = "Edad_media",text = "departamento_id")
1.4.0.3 LISA clusters and exploration of Moran II with maps
I add the output of the local Moran to the data by department (community) to be able to map them
moran_map_Edad_media <- cbind(coviddepcaba_pos, propc_lm_lwb)
head(moran_map_Edad_media)I copy the p-values into the variable $p to make it easier to call
moran_map_Edad_media$p <- moran_map_Edad_media$Pr.z....E.Ii..
head(moran_map_Edad_media)Quadrant for the LISA clusters (with positive autocorrelation) (Moran’s scatterplor hotspots on the map) It must be done with standardized or scaled values (value - mean)
I calculate the standardized values of the values and the averages of the neighbors
z_Edad_media <- moran_map_Edad_media$Edad_media - mean(moran_map_Edad_media$Edad_media) # Z propcasos
lag_Edad_media <- lag.listw(lwb, moran_map_Edad_media$Edad_media)#lag de Edad_media
z_lag_Edad_media <-lag_Edad_media - mean(lag_Edad_media) # Z de lag de Edad_mediaI create a coloring vector
quad <- vector(mode="numeric",length=nrow(moran_map_Edad_media))
quad[z_Edad_media > 0 & z_lag_Edad_media > 0] <- 1 #alto-alto HH rojo
quad[z_Edad_media > 0 & z_lag_Edad_media < 0] <- 2 #alto-bajo HL rosado
quad[z_Edad_media < 0 & z_lag_Edad_media > 0] <- 3 #bajo-alto LH celeste
quad[z_Edad_media < 0 & z_lag_Edad_media < 0] <- 4 #bajo-bajo LL azul
quad[moran_map_Edad_media$p > 0.05] <- 5 #No Signif blancoLet’s see the result
quad## [1] 5 5 5 3 5 5 5 4 5 5 5 5 5 5 3moran_map_Edad_media$quad <- quad
moran_map_Edad_mediaLISA <- c("red", "lightpink","skyblue2","blue", "white")I make the map of LISA with the numbers of the communes
tm_shape(moran_map_Edad_media) +
tm_fill(col="quad", palette = LISA, style = "fixed", breaks = c(1,2,3,4,5,6),
labels = c("HH", "HL", "LH","LL", "no signif"),title="Agrupamientos LISA") +
tm_legend(outside=TRUE)+
tm_borders(col = "Grey")+
tm_text(text = "departamento_id")
1.5 Conclusions
In the southern part of the city (Comunas 8 and 4) the people who are tested for covid are, on average, younger than in the rest of the city. In Comuna 15 in particular we see that it is a commune of low values surrounded by higher values.
What we saw on the map “by eye” we could also see when making the local Moran’s i calculations. With values of I of 4.7, 4.8 and -1.3 for communes 4, 8 and 15 respectively. However, for commune 4 we see that when doing the LISA we find that “Low surrounded by High” values appear as values. This makes sense when we see that commune 4 is elongated and has many more neighbors with higher average age values than it.
2 Spanish
2.1 Introducción
Me planteo analizar la distribución de edad de las testeadas para COVID en la Ciudad Autonoma de Buenos Aires (CABA) acumulados al día 25 de Octubre de 2021.
2.2 Paquetes a usar
library(sf)
library(tmap)
library(tidyverse)
library(spdep)2.3 Carga de datos y análisis exploratorio
2.3.0.1 Datos geográficos - Unidades administrativas
Tal como se hizo en clase utilizaré los datos del INDEC y filtro los departamentos de la CABA y creo una variable nueva que se queda solo con el numero de comuna y otra con el numero de habitantes en formato número.
deptos <- st_read("pxdptodatosok.shp")## Reading layer `pxdptodatosok' from data source
## `/Users/inespatop/Documents/resume/clase2/pxdptodatosok.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 527 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.02985 ymin: -90 xmax: -25.02314 ymax: -21.74506
## Geodetic CRS: WGS 84depcaba <- deptos %>%
filter(codpcia == "02") %>%
mutate(departamento_id = substr(link,3, 5), .keep ="unused") %>%
mutate(pob = as.numeric(personas), .keep ="unused") %>%
select(departamento_id, departamen, pob)
depcabaHago un plot simplemente para ver el mapa
qtm(depcaba, fill = "pob",text = "departamento_id")
2.3.0.2 Carga de datos de testeos de COVID en la Ciudad de Buenos Aires
Dado que los datos son muy grandes cargo un dataset ya filtrado del dataset solo para CABA
##covid <- read_csv("./Covid19Casos.csv")
##glimpse(covid)
##covid_caba<-covid[covid$residencia_provincia_id == "02",]
##rm(covid)
##save(covid_caba,"covid_caba.RData")
load("covid_caba.RData")
head(covid_caba)##3 Exploración de datos y generación de medidas resumen
Hago ahora una medida resumen de catidad de casos por edad separando por casos confirmados y descartados.
covid_age <- covid_caba %>%
group_by( clasificacion_resumen) %>%
select(residencia_departamento_nombre, residencia_departamento_id, clasificacion_resumen, edad ) %>%
group_by(residencia_departamento_id,clasificacion_resumen) %>%
summarise(Edad_media = mean(edad,na.rm = T) )%>%
filter (residencia_departamento_id != "0")%>%
filter (clasificacion_resumen != "Sospechoso")## `summarise()` has grouped output by 'residencia_departamento_id'. You can
## override using the `.groups` argument.covid_ageMiro la distribucion de los datos tanto antes como luego de generar la medida resumen
##Total
ggplot(covid_caba%>% filter (residencia_departamento_id != "0")%>% filter (clasificacion_resumen != "Sospechoso"), aes(y=edad,x=residencia_departamento_id, fill=clasificacion_resumen))+geom_boxplot()+labs(title = "Distribución de edad acumulada por comuna")## Warning: Removed 24 rows containing non-finite values (stat_boxplot).
##medida resumen
ggplot(covid_age, aes(y=Edad_media,x=clasificacion_resumen, fill=clasificacion_resumen))+geom_boxplot()+labs(title = "Distribución de edad promerdio acumulada")
2.3.0.3 Unión con datos geog´raficos y grafico de edad promedio testeada
Unión
coviddepcaba <- left_join(depcaba, covid_age, by = c("departamento_id"="residencia_departamento_id"))Edad media en casos confirmados
qtm(coviddepcaba%>% filter (clasificacion_resumen == "Confirmado"), fill = "Edad_media") Edad media en casos descartados
qtm(coviddepcaba%>% filter (clasificacion_resumen == "Descartado"), fill = "Edad_media") A simple vista ya se ve que en la zona sur hay mayor testeo de personas jóvenes y mayor cantidad de test positivos a personas jóvenes.
Ahora hago un gráfico por cunatiles
tm_shape(coviddepcaba) +
tm_fill("Edad_media", style = "quantile", n = 5, palette = "Reds") +
tm_borders(alpha = 0.1) +
tm_layout(main.title = "Casos confirmados al 22/10/2021", main.title.size = 0.7 ,
legend.position = c("right", "bottom"), legend.title.size = 0.8)
2.4 Análisis de Vecinos
2.4.0.1 Lista de vecinos (de clase nb) a partir de los poligonos con poly2nb
Para el resto del trabajo usaré solo los datos de testeos de covid positivos.
coviddepcaba_pos<- coviddepcaba%>% filter (clasificacion_resumen == "Confirmado")Calulo los vecinos usando poly2nb tal y como se hizo en clase-
w <- poly2nb(coviddepcaba_pos, row.names = "departamento_id")Exploro los vecinos
summary(w)## Neighbour list object:
## Number of regions: 15
## Number of nonzero links: 64
## Percentage nonzero weights: 28.44444
## Average number of links: 4.266667
## Link number distribution:
##
## 3 4 5 6 7
## 6 4 2 1 2
## 6 least connected regions:
## 1 8 9 10 12 13 with 3 links
## 2 most connected regions:
## 5 7 with 7 linksMapeo los vecinos
plot(st_geometry(coviddepcaba_pos), border="grey")
plot(w, coordinates(as(coviddepcaba_pos, "Spatial")), add=TRUE, col="blue")
title("Vecinos calculados con poly2nb")
Ver las matrices de pesos vecinos
En este caso los pesos binarios se realizaron usando style='B' como opción en particualr esto genera una matriz de 0 y 1 en la cual cada vecino tiene el mismo peso.
wmatb <- nb2mat(w, style='B')
wmatb## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0
## 2 1 0 1 0 1 0 0 0 0 0 0 0 0 1
## 3 1 1 0 1 1 0 0 0 0 0 0 0 0 0
## 4 1 0 1 0 1 0 1 1 0 0 0 0 0 0
## 5 0 1 1 1 0 1 1 0 0 0 0 0 0 1
## 6 0 0 0 0 1 0 1 0 0 0 1 0 0 0
## 7 0 0 0 1 1 1 0 1 1 1 1 0 0 0
## 8 0 0 0 1 0 0 1 0 1 0 0 0 0 0
## 9 0 0 0 0 0 0 1 1 0 1 0 0 0 0
## 10 0 0 0 0 0 0 1 0 1 0 1 0 0 0
## 11 0 0 0 0 0 1 1 0 0 1 0 1 0 0
## 12 0 0 0 0 0 0 0 0 0 0 1 0 1 0
## 13 0 0 0 0 0 0 0 0 0 0 0 1 0 1
## 14 0 1 0 0 1 0 0 0 0 0 0 0 1 0
## 15 0 0 0 0 1 1 0 0 0 0 1 1 1 1
## [,15]
## 1 0
## 2 0
## 3 0
## 4 0
## 5 1
## 6 1
## 7 0
## 8 0
## 9 0
## 10 0
## 11 1
## 12 1
## 13 1
## 14 1
## 15 0
## attr(,"call")
## nb2mat(neighbours = w, style = "B")Para pesos estandarizados por fila (style=‘W’) en cambio cada vecino tiene un peso proporcional. Por ejemplo 3 vecinos tendran un peso de 0.33
wmatW <- nb2mat(w, style='W')
wmatW## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## 1 0.00 0.3333333 0.3333333 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.25 0.0000000 0.2500000 0.0000000 0.2500000 0.0000000 0.0000000 0.0000000
## 3 0.25 0.2500000 0.0000000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000
## 4 0.20 0.0000000 0.2000000 0.0000000 0.2000000 0.0000000 0.2000000 0.2000000
## 5 0.00 0.1428571 0.1428571 0.1428571 0.0000000 0.1428571 0.1428571 0.0000000
## 6 0.00 0.0000000 0.0000000 0.0000000 0.2500000 0.0000000 0.2500000 0.0000000
## 7 0.00 0.0000000 0.0000000 0.1428571 0.1428571 0.1428571 0.0000000 0.1428571
## 8 0.00 0.0000000 0.0000000 0.3333333 0.0000000 0.0000000 0.3333333 0.0000000
## 9 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3333333 0.3333333
## 10 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000
## 11 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.2000000 0.2000000 0.0000000
## 12 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 13 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 14 0.00 0.2500000 0.0000000 0.0000000 0.2500000 0.0000000 0.0000000 0.0000000
## 15 0.00 0.0000000 0.0000000 0.0000000 0.1666667 0.1666667 0.0000000 0.0000000
## [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.2500000 0.0000000
## 3 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1428571 0.1428571
## 6 0.0000000 0.0000000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000
## 7 0.1428571 0.1428571 0.1428571 0.0000000 0.0000000 0.0000000 0.0000000
## 8 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 9 0.0000000 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 10 0.3333333 0.0000000 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000
## 11 0.0000000 0.2000000 0.0000000 0.2000000 0.0000000 0.0000000 0.2000000
## 12 0.0000000 0.0000000 0.3333333 0.0000000 0.3333333 0.0000000 0.3333333
## 13 0.0000000 0.0000000 0.0000000 0.3333333 0.0000000 0.3333333 0.3333333
## 14 0.0000000 0.0000000 0.0000000 0.0000000 0.2500000 0.0000000 0.2500000
## 15 0.0000000 0.0000000 0.1666667 0.1666667 0.1666667 0.1666667 0.0000000
## attr(,"call")
## nb2mat(neighbours = w, style = "W")Vamos lo mismo ahora en formato lista
lwb <- nb2listw(w, style='B')
lwb## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 15
## Number of nonzero links: 64
## Percentage nonzero weights: 28.44444
## Average number of links: 4.266667
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 15 225 64 128 1208head(lwb)## $style
## [1] "B"
##
## $neighbours
## Neighbour list object:
## Number of regions: 15
## Number of nonzero links: 64
## Percentage nonzero weights: 28.44444
## Average number of links: 4.266667
##
## $weights
## $weights[[1]]
## [1] 1 1 1
##
## $weights[[2]]
## [1] 1 1 1 1
##
## $weights[[3]]
## [1] 1 1 1 1
##
## $weights[[4]]
## [1] 1 1 1 1 1
##
## $weights[[5]]
## [1] 1 1 1 1 1 1 1
##
## $weights[[6]]
## [1] 1 1 1 1
##
## $weights[[7]]
## [1] 1 1 1 1 1 1 1
##
## $weights[[8]]
## [1] 1 1 1
##
## $weights[[9]]
## [1] 1 1 1
##
## $weights[[10]]
## [1] 1 1 1
##
## $weights[[11]]
## [1] 1 1 1 1 1
##
## $weights[[12]]
## [1] 1 1 1
##
## $weights[[13]]
## [1] 1 1 1
##
## $weights[[14]]
## [1] 1 1 1 1
##
## $weights[[15]]
## [1] 1 1 1 1 1 1
##
## attr(,"mode")
## [1] "binary"
## attr(,"B")
## [1] TRUE2.5 Exploración de Indices de asociación/correlación espacial
2.5.0.1 Moran I Global
## todos los test
moran.test(coviddepcaba_pos$Edad_media, lwb)##
## Moran I test under randomisation
##
## data: coviddepcaba_pos$Edad_media
## weights: lwb
##
## Moran I statistic standard deviate = 2.1678, p-value = 0.01509
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.20243360 -0.07142857 0.01596030Vemos que el i de Moran Global da significativo (pvalor < 0.05) y como es mayor a cero aunque sea leve (0.2), tenemos una distribución de los datos agrupada. Esto se correlaciona con lo que veíamos en los gráficos.
2.5.0.1.0.1 Correlograma de Moran
Hasta el tercer vecino vemos el i de Morán
Icorr3 <-sp.correlogram(neighbours=w,var=coviddepcaba_pos$Edad_media,order=3,method="I",zero.policy=TRUE)
plot(Icorr3)
Vemos aqui que con el primer vecino hay una correlación mayor que con el segundo y asi
2.5.0.1.0.2 Diagrama de dispersión de Morán
mp <- moran.plot(coviddepcaba_pos$Edad_media, lwb, zero.policy=TRUE, labels=as.character(coviddepcaba_pos$"departamento_id"))
mpEn este caso vemos que hay un valor importante, o influyente, que es claramente los valores de valor de edad bajos, que vimos antes están en la zona sur de la ciudad. Sin embargo cuando vemos el resto de los valores vemos muchos datos de hecho muchos valores en la zona de anticorrelación. Como esperamos la pendiente es muy baja (i de moran).
Solo a forma de recordar: en el eje x tenemos la edad promedio de los test de covid positivo de covid en cada partido. En el eje y vemos el lag espacial de x que es el promedio de valores de vecinos a x. La pendiente de recta de regresion es Moran I.
2.5.0.2 Moran local Ii
Calculamos el indice de Moran Local. Para esto usamos la función localmoran
El dominio del i de Moran es -1 a 1 y esperamos ver ahora si hay zonas en las cuales haya autocorrelación. Busco entonces los valores con p valor menos a 0.05
propc_lm_lwb <- localmoran(coviddepcaba_pos$Edad_media, lwb)
propc_lm_lwb_df<-as.data.frame(propc_lm_lwb)
propc_lm_lwb_dfpropc_lm_lwb_df[propc_lm_lwb_df$`Pr(z != E(Ii))`<0.05,]Ii: Moran local E.Ii: experanza del Moran local Var.Ii: varianza del local Moran statistic Z.Ii: desvío estándar del Moran local Pr(): p-value del Moran local
Vemos que para la comuna 4 y 8, los vecinos tienen valores similares y para la comuna 15 valores disímiles. Si vemos esto nuevamente en el mapa vemos que la comuna 4 y 8 tienen edades similares y están en el sur de la ciudad mientras que la 15 esta rodeada de valores de edad más bajos que los propios.
qtm(coviddepcaba_pos, fill = "Edad_media",text = "departamento_id")
2.5.0.3 LISA clusters y exploracion de Moran Ii con mapas
Agrego la salida del Moran local a los datos por depto (comuna) para poder mapearlos
moran_map_Edad_media <- cbind(coviddepcaba_pos, propc_lm_lwb)
head(moran_map_Edad_media)Copio los p-valores en la variable $p para que sea mas facil de llamar
moran_map_Edad_media$p <- moran_map_Edad_media$Pr.z....E.Ii..
head(moran_map_Edad_media)Cuadrante para los LISA clusters (con autocorrelación positiva) (hotspots del scatterplor de Moran en el mapa) Hay que hacerlo con valores estadarizados o escalados (valor - media)
Calculo los valores estandarizados de los valores y los promedios de los vecinos
z_Edad_media <- moran_map_Edad_media$Edad_media - mean(moran_map_Edad_media$Edad_media) ## Z propcasos
lag_Edad_media <- lag.listw(lwb, moran_map_Edad_media$Edad_media)##lag de Edad_media
z_lag_Edad_media <-lag_Edad_media - mean(lag_Edad_media) ## Z de lag de Edad_mediaCreo un vector de coloración
quad <- vector(mode="numeric",length=nrow(moran_map_Edad_media))
quad[z_Edad_media > 0 & z_lag_Edad_media > 0] <- 1 ##alto-alto HH rojo
quad[z_Edad_media > 0 & z_lag_Edad_media < 0] <- 2 ##alto-bajo HL rosado
quad[z_Edad_media < 0 & z_lag_Edad_media > 0] <- 3 ##bajo-alto LH celeste
quad[z_Edad_media < 0 & z_lag_Edad_media < 0] <- 4 ##bajo-bajo LL azul
quad[moran_map_Edad_media$p > 0.05] <- 5 ##No Signif blancoVemos el resutlado
quad## [1] 5 5 5 3 5 5 5 4 5 5 5 5 5 5 3moran_map_Edad_media$quad <- quad
moran_map_Edad_mediaLISA <- c("red", "lightpink","skyblue2","blue", "white")Hago el mapa de LISA con los números de las comunas
tm_shape(moran_map_Edad_media) +
tm_fill(col="quad", palette = LISA, style = "fixed", breaks = c(1,2,3,4,5,6),
labels = c("HH", "HL", "LH","LL", "no signif"),title="Agrupamientos LISA") +
tm_legend(outside=TRUE)+
tm_borders(col = "Grey")+
tm_text(text = "departamento_id")
2.6 Conclusiones
En la zona sur de la ciudad (Comunas 8 y 4) las persona que se testean para covid tienen en promedio una edad menos que en el resto de la ciduad. En la Comuna 15 en particular vemos que es una comuna de valores bajos rodeada de valores más altos.
Esto que vimos en el mapa “a ojo” también pudimos verlo al realizar los cáculos de i de Moran local. Con valores de I de 4.7, 4.8 y -1.3 para las comunas 4, 8 y 15 respectivamente. Sin embargo, para la comuna 4 vemos que al hacer el LISA encontramos que aprace como valores “Bajos rodeados de Altos”. Esto tiene sentido cuando vemos que la comuna 4 es alargado y tiene muchos más vecinos con valores de edad promedio más altos que ella.