Datos por areas

“To err is human, to forgive divine - but to include errors in your design is statistical” - Leslie Kish

0.1 Datos por areas

Se recogen tantos datos en las divisiones administrativas que existen técnicas especializadas para analizarlos.

Hasta ahora, los datos presentados han sido puntos, los cuales tienen la ventaja de ser la representación más precisa que se puede tener, no obstante, es más común encontrar los datos agregados por razones de confidencialidad. Desafortunadamente agregar datos hace que se pierda información y reduce la potencia de los análisis estadísticos.

Las divisiones regionales no están pensadas para beneficiar a los estadísticos. Responden a una mezcla de lìmites hostóricos, naturales como ríos y montañas y a razones políticas. Esto se traduce en que que en áreas pequeñas se puede concentrar mucha población y viceverza. Estas razones pueden producir mapas engañosos. Un remedio para esto es son los cartogramas.

UN cartograma es un mapa o diagrama que muestra datos de cantidad asociados a sus respectivas áreas, mediante la modificación de los tamaños de las unidades de enumeración. La información es aportada mediante la distorsión de las superficies reales, utilizando cada superficie de enumeración como un símbolo proporcional, el cual aumenta o disminuye en función de los valores correspondientes. Un ejemplo podría ser la representación de los países, donde su tamaño en el diagrama dependiera del número de habitantes.

library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

library(lattice)

# Create a new folder called "London Suicides" inside your working directory (my.dir)
# and save the data downloaded from
# https://sites.google.com/a/r-inla.org/stbook/datasets

# Import the shapefile 
my.dir<-"D:\\DriveW7\\2019\\Sabana\\Estadistica aplicada Salud Publica\\00. Datos\\Datos 2015\\LondonSuicides_data\\"

london.gen <- readShapePoly(paste(my.dir,"LDNSuicides.shp",sep=""))

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

# Load code for thinning map 
source(paste(my.dir,"naiveThinShp.R",sep=""))
london.gen.thin <- naive.thin.shp(london.gen,thin=5)

## [1] "WARNING: only use the thinned shapefile for plotting and always check against the original version!!"

# Load the data 

load(paste(my.dir,"LondonSuicides.RData",sep=""))
names <- sort(london.gen$NAME)
data <- data.frame(NAME=names, y=y, E=E, x1=x1, x2=x2)
Nareas <- length(data[,1])

# *** Code for Figure 1.7 
# Compute Standard Mortality Ratio
SMR <- y/E
SMR.cutoff <- c(0.6,0.8,1,1.2,1.4,1.8)
SMR.factor <- cut(SMR,breaks=SMR.cutoff,include.lowest=TRUE)
#Deprivation index
x.cutoff <- round(quantile(x1,seq(0,1,0.2)),1)
x1.factor <- cut(x1,breaks=x.cutoff,include.lowest=TRUE)
#Fragmentation index
x2.factor <- cut(x2,breaks=x.cutoff,include.lowest=TRUE)

data <- data.frame(NAME=names,SMR=SMR.factor,x1=x1.factor, x2=x2.factor)
boroughs <- london.gen.thin
data.boroughs <- attr(boroughs, "data")
attr(boroughs, "data") <- merge(data.boroughs, data, by="NAME")

#trellis.par.set(axis.line=list(col=NA))
# spplot(obj=boroughs, zcol="x1", col.regions=gray(4.5:0.5/5),main="",par.settings=list(fontsize=list(text=25)))
# spplot(obj=boroughs, zcol="x2", col.regions=gray(4.5:0.5/5),main="",par.settings=list(fontsize=list(text=25)))
spplot(obj=boroughs, zcol="SMR", col.regions=gray(4.5:0.5/5),main="",par.settings=list(fontsize=list(text=25)))

0.1.1 Cartograma

Las grandes áreas, como las ciudades o los países, a menudo se dividen en unidades administrativas más pequeñas, a menudo en zonas con una población aproximadamente igual. Pero el área de esas unidades puede variar considerablemente. Cuando se cartografían, las áreas grandes tienen más “peso” visual que las áreas pequeñas, aunque también hay mucha gente que vive en las áreas pequeñas.

Una técnica para corregir esto es el cartograma. Se trata de una distorsión controlada de las regiones, expandiendo algunas y contrayendo otras, de modo que el área de cada región es proporcional a una cantidad deseada, como la población. El cartograma también trata de mantener la geografía correcta tanto como sea posible, manteniendo las regiones en aproximadamente el mismo lugar relativo entre sí.

El paquete de cartogramas contiene funciones para crear cartogramas. Le da un marco de datos espaciales y el nombre de una columna, y obtiene un marco de datos similar pero con regiones distorsionadas para que el área de la región sea proporcional al valor de la columna de las regiones.

También usará el paquete rgeos para calcular las áreas de regiones individuales con la función gArea().

# Use the cartogram and rgeos packages
library(cartogram)
library(rgeos)

## rgeos version: 0.3-28, (SVN revision 572)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

# Make a scatterplot of electorate vs borough area
london_ref<-readRDS("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/Bibliografía/DataCamp Spatial Statistics in R/datasets/london_eu.rds", refhook = NULL)

names(london_ref)

##  [1] "NAME"             "TOTAL_POP"        "Electorate"      
##  [4] "Votes_Cast"       "Remain"           "Leave"           
##  [7] "Rejected_Ballots" "Pct_Remain"       "Pct_Leave"       
## [10] "Pct_Rejected"     "Assembly"

plot(london_ref$Electorate, gArea(london_ref, byid = TRUE))

# Make a cartogram, scaling the area to the electorate
carto_ref <- cartogram(london_ref, "Electorate")

## 
## Please use cartogram_cont() instead of cartogram().

## Mean size error for iteration 1: 1.5881743190908

## Mean size error for iteration 2: 1.32100446455657

## Mean size error for iteration 3: 1.18227723476121

## Mean size error for iteration 4: 1.10676057030171

## Mean size error for iteration 5: 1.0657703107641

## Mean size error for iteration 6: 1.04259259672006

## Mean size error for iteration 7: 1.02832326230708

## Mean size error for iteration 8: 1.01931941526112

## Mean size error for iteration 9: 1.01341424685212

## Mean size error for iteration 10: 1.00941370606418

## Mean size error for iteration 11: 1.00663364742297

## Mean size error for iteration 12: 1.00470553629914

## Mean size error for iteration 13: 1.00336434720465

## Mean size error for iteration 14: 1.00241457265516

## Mean size error for iteration 15: 1.00174179254187

plot(carto_ref)

# Check the linearity of the electorate-area plot
plot(carto_ref$Electorate, gArea(carto_ref, byid = TRUE))

# Make a fairer map of the Remain percentage
spplot(carto_ref, "Remain")

0.1.2 Ejericio

# Use the cartogram and rgeos packages
library(cartogram)
library(rgeos)
library(maptools)
library(lattice)

# Make a scatterplot of electorate vs borough area


desempeno <- readShapePoly(paste("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/shp/desempeno.shp",sep=""))

summary(desempeno$Desempen19)

desempeno$Desempen19<-as.numeric(as.character(desempeno$Desempen19))
desempeno$AREA<-as.numeric(as.character(desempeno$AREA))

# plot(desempeno$Desempen19, desempeno$AREA)

# Make a cartogram, scaling the area to the electorate
carto_ref <- cartogram(desempeno, "Desempen19")
plot(carto_ref)

# Check the linearity of the electorate-area plot
# plot(carto_ref$Electorate, gArea(carto_ref, byid = TRUE))
# 
# # Make a fairer map of the Remain percentage
# spplot(carto_ref, "Remain")

Desempen	mpio-dpto
Desempen1	Municipio
Desempen2	Departamento
Desempen3	Grupo dotaciones
Desempen4	Categoria de ruralidad
Desempen5	Educación - Cobertura media neta
Desempen6	Educación - SABER 11 Matematicas
Desempen7	Educación - SABER 11 Lenguaje
Desempen8	Educación - Cobertura Transición
Desempen9	Salud - Cobertura salud
Desempen10	Salud - Vacunación Pentavalente
Desempen11	Salud - Mortalidad Infantil
Desempen12	Servicios públicos - Cobertura eléctrica rural
Desempen13	Servicios públicos - Cobertura internet
Desempen14	Servicios públicos - Cobertura Acueducto
Desempen15	Servicios públicos - Cobertura Alcantarillado
Desempen16	Seguridad - Hurtos x 10000 hab
Desempen17	Seguridad - Homicidios x 10000 hab
Desempen18	Seguridad - Violencia Intrafamiliar x 10000 hab
Desempen19	Ingresos tributarios y no tributarios per cápita
Desempen20	Desempeno Municipal

0.1.3 Ejercicio

Realice el cartograma del PIB Departamental

# Use the cartogram and rgeos packages
library(cartogram)
library(rgeos)
library(maptools)
library(lattice)

# Make a scatterplot of electorate vs borough area


dpto <- readShapePoly("D:\\DriveW7\\2019\\Sabana\\Estadistica aplicada Salud Publica\\00. Datos\\shp\\deptPIB17a.shp",sep="")



dpto$pib1_PIB1<-as.numeric(as.character(dpto$pib1_PIB))

# plot(dpto$Desempen19, dpto$AREA)

# Make a cartogram, scaling the area to the electorate
carto_ref <- cartogram(dpto, "pib1_PIB1")
plot(carto_ref)

# Check the linearity of the electorate-area plot
# plot(carto_ref$Electorate, gArea(carto_ref, byid = TRUE))
# 
# # Make a fairer map of the Remain percentage
# spplot(carto_ref, "Remain")

0.1.3.1 matrices de contigüidad

0.1.3.1.1 Creación telesado (Partición)

x <- c(1,2,2,1,1)
y <- c(1,1,2,2,1)

vl1 <- cbind(x,y);vl2 <- cbind(x-1,y+1);vl3 <- cbind(x,y+1);vl4 <- cbind(x+1,y+1);vl5 <- cbind(x-1,y);vl6 <- cbind(x+1,y);vl7 <- cbind(x-1,y-1);vl8 <- cbind(x,y-1);vl9 <- cbind(x+1,y-1)

Sr1 = Polygon(vl1)
Sr2 = Polygon(vl2)
Sr3 = Polygon(vl3)
Sr4 = Polygon(vl4)
Sr5 = Polygon(vl5)
Sr6 = Polygon(vl6)
Sr7 = Polygon(vl7)
Sr8 = Polygon(vl8)
Sr9 = Polygon(vl9)

Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3), "s3")
Srs4 = Polygons(list(Sr4), "s4")
Srs5 = Polygons(list(Sr5), "s5")
Srs6 = Polygons(list(Sr6), "s6")
Srs7 = Polygons(list(Sr7), "s7")
Srs8 = Polygons(list(Sr8), "s8")
Srs9 = Polygons(list(Sr9), "s9")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3,Srs4,Srs5,Srs6,Srs7,Srs8,Srs9))
plot(SpP, col = 1:9, pbg="white")

library(spdep)

## Loading required package: Matrix

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`

## 
## Attaching package: 'spData'

## The following objects are masked _by_ '.GlobalEnv':
## 
##     x, y

(wr <- poly2nb(SpP, queen=FALSE))

## Neighbour list object:
## Number of regions: 9 
## Number of nonzero links: 24 
## Percentage nonzero weights: 29.62963 
## Average number of links: 2.666667

(wm <- nb2mat(wr, style='B', zero.policy = TRUE))

##    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## s1    0    0    1    0    1    1    0    1    0
## s2    0    0    1    0    1    0    0    0    0
## s3    1    1    0    1    0    0    0    0    0
## s4    0    0    1    0    0    1    0    0    0
## s5    1    1    0    0    0    0    1    0    0
## s6    1    0    0    1    0    0    0    0    1
## s7    0    0    0    0    1    0    0    1    0
## s8    1    0    0    0    0    0    1    0    1
## s9    0    0    0    0    0    1    0    1    0
## attr(,"call")
## nb2mat(neighbours = wr, style = "B", zero.policy = TRUE)

0.1.3.1.2 Graficamos los enlaces entre los polígonos:

plot(SpP, col='gray', border='blue')
xy <- coordinates(SpP)
plot(wr, xy, col='red', lwd=2, add=TRUE)

0.1.3.1.3 Creación lista de vecinos y grafica de Moran

wrl <- nb2listw(wr, style="B") 
# Número de vecinos por área
i <- rowSums(wm)
vals <- i/sum(i)
moran.plot(vals,wrl)

0.1.3.1.4 Otras formas de influencia entre polígonos:

0.1.3.1.4.1 Basado en la distancia:

wd10 <- dnearneigh(xy, 0, 10, longlat=TRUE)
wd25 <- dnearneigh(xy, 0, 25)
#, longlat=TRUE)

0.1.3.1.4.2 Vecino más cercano:

k3 <- knn2nb(knearneigh(xy, k=3, RANN=FALSE))
k6 <- knn2nb(knearneigh(xy, k=6, RANN=FALSE))

0.1.3.1.4.3 Rezago 2:

wr2 <- wr
for (i in 1:length(wr)) {
    lag1 <- wr[[i]]
    lag2 <- wr[lag1]
    lag2 <- sort(unique(unlist(lag2)))
    lag2 <- lag2[!(lag2 %in% c(wr[[i]], i))]
    wr2[[i]] <- lag2
}

0.1.3.1.4.4 Graficamos todas las influencias con la función plotit

plotit <- function(nb, lab='') {
  plot(SpP, col='gray', border='white')
  plot(nb, xy, add=TRUE, pch=20)
  text(2.5, 2.8, paste0('(', lab, ')'), cex=1.25)
}
par(mfrow=c(2, 3), mai=c(0,0,0,0))
plotit(wr, 'adjacency')
plotit(wr2, 'lag-2 adj.')
plotit(wd10, '10 km')
plotit(wd25, '25 km')
plotit(k3, 'k=3')
plotit(k6, 'k=6')

0.1.4 Prueba de autocorrelación espacial

0.1.4.1 I de moran

Es una medida de autocorrelación espacial. La autocorrelación espacial es más compleja que una dimensión de autocorrelación debido a que la correlación espacial es multi-dimensionales (es decir, 2 o 3 dimensiones del espacio) y multi-direccional.

La I de Moran se define:

\[{\displaystyle I={\frac {N}{\sum _{i}\sum _{j}w_{ij}}}{\frac {\sum _{i}\sum _{j}w_{ij}(X_{i}-{\bar {X}})(X_{j}-{\bar {X}})}{\sum _{i}(X_{i}-{\bar {X}})^{2}}}}\]

donde $N$ es el número de unidades espaciales indexados por $ i$ y $j$; $X$ es la variable de interés; ${\bar {X}}$ es la media de $B$; y ${\displaystyle w_{ij}}$ es un elemento de una matriz de pesos espaciales.

Es una versión del coeficiente de correlación de Pearson

$\displaystyle r={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{{\sqrt {\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}{\sqrt {\sum _{i=1}^{n}(y_{i}-{\bar {y}})^{2}}}}}$

AL ser una versión del coeficiente de correlación, se interpreta similarmente. Los valores negativos (positivos) indican negativo (positivo) de autocorrelación espacial. Los valores oscilan entre -1 (indicando dispersión perfecta) a 1 (correlación perfecta). Un valor de cero indica un patrón espacial aleatoria.

I de Moran es inversamente proporcional a C de Geary, pero no es idéntica. El I de Moran es una medida de autocorrelación espacial global, mientras que C de Geary es más sensible a la autocorrelación espacial local.

cual sería el valor del I de Moran del teselado.

https://www.researchgate.net/figure/Anselins-Moran-scatter-plot-interpretation-guide_fig1_294870442

El paquete spdep tiene funciones para medir la correlación espacial, también conocida como dependencia espacial. El cálculo de estas medidas requiere que primero se determine qué regiones son vecinas a través de la función poly2nb(), abreviatura de “polígonos a vecinos”. El resultado es un objeto de la clase nb. Luego se puede calcular la estadística de la prueba y ejecutar una prueba de significación sobre la hipótesis nula de que no hay correlación espacial. La prueba de significación puede ser hecha por Monte-Carlo o por modelos teóricos.

En este ejemplo se utilizará la estadística Moran “I” para probar la correlación espacial

\[ H_0: Aleatoriedad\ espacial\ completa \]

\[ H_a: autocorrelación \]

0.1.4.2 Prueba de Moran por Monte Carlo

# Use the spdep package
library(spdep)

# Make neighbor list
borough_nb <- poly2nb(london_ref)

# Get center points of each borough
borough_centers <- coordinates(london_ref)

# Show the connections
plot(london_ref); plot(borough_nb, borough_centers, add = TRUE)

# Map the total pop'n
spplot(london_ref, zcol = "TOTAL_POP")

# Run a Moran I test on total pop'n
moran.test(
  london_ref$TOTAL_POP, 
  nb2listw(borough_nb)
)

## 
##  Moran I test under randomisation
## 
## data:  london_ref$TOTAL_POP  
## weights: nb2listw(borough_nb)    
## 
## Moran I statistic standard deviate = 1.2124, p-value = 0.1127
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.11549264       -0.03225806        0.01485190

# Map % Remain
spplot(london_ref, zcol = "Pct_Remain")

# Run a Moran I MC test on % Remain
moran.mc(
  london_ref$Pct_Remain, 
  nb2listw(borough_nb), 
  nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  london_ref$Pct_Remain 
## weights: nb2listw(borough_nb)  
## number of simulations + 1: 1000 
## 
## statistic = 0.42841, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

0.1.4.3 Ejercicio

Calcule el I De Moran para el PIB departamental

library(spdep)
library(cartogram)
library(rgeos)
library(maptools)
library(lattice)

dpto <- readShapePoly("D:\\DriveW7\\2019\\Sabana\\Estadistica aplicada Salud Publica\\00. Datos\\shp\\deptPIB17a.shp")
dpto$pib1_PIB1<-as.numeric(as.character(dpto$pib1_PIB))

# 
# 
dpto_nb <- poly2nb(dpto)
# 
# # Get center points of each dpto
# dpto_centers <- coordinates(dpto)
# 
# # Show the connections
# plot(dpto); plot(dpto_nb, dpto_centers, add = TRUE)
# 
# dpto$pib1_PIB1s <- scale(dpto$pib1_PIB1)
# 
# # Map the pib
# spplot(dpto, zcol = "pib1_PIB1")
# 
# # Run a Moran I test on total pop'n
# moran.test(
#   dpto$pib1_PIB1, 
#   nb2listw(dpto_nb)
# )
# 
# # Map % Remain
# # spplot(dpto, zcol = "pib1_PIB1")
# 
# # Run a Moran I MC test on % Remain
# moran.mc(
# dpto$pib1_PIB1, 
# nb2listw(dpto_nb), 
#   nsim = 999)

mp <- moran.plot(dpto$pib1_PIB1, nb2listw(dpto_nb), labels = as.character(dpto$NOM_DPTO_x), 
    xlab = "PIB", ylab = "PIB rezagado")

0.1.5 LISA

El análisis espacial global es una estadística para resumir toda la zona de estudio. En otras palabras, el análisis global asume homogeneidad. Si este supuesto no se sostiene, el hecho de tener sólo una estadística, no tiene sentido, como sucede siempre en datos por áreas, las unidades de estudio no son homogeneas.

Pero si no hay autocorrelación global o ninguna agrupación, todavía podemos encontrar grupos a nivel local utilizando autocorrelación espacial local. El hecho de que de la Morán I es una suma de productos cruzados individuales es explotado por los “indicadores locales de asociación espacial” (LISA) para evaluar la agrupación de las unidades individuales mediante el cálculo de la I de Moran local para cada unidad espacial y la evaluación de la significación estadística para cada I.

dpto <- readShapePoly("deptPIB17a.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

dpto_nbq <- poly2nb(dpto)  #Queens neighbors
dpto_nbq_w <- nb2listw(dpto_nbq)  #row standardize

dpto$pib1_PIB1<-as.numeric(as.character(dpto$pib1_PIB))
dpto$pib1_PIB1s <- scale(dpto$pib1_PIB1)
dpto$pib1_PIB1s_lag <- lag.listw(dpto_nbq_w, dpto$pib1_PIB1s)

# Create a LISA map
locm <- localmoran(dpto$pib1_PIB1, dpto_nbq_w)  #Calculate the local Morann
# Organize the data
dpto$quad_sig <- NA
dpto@data[(dpto$pib1_PIB1s >= 0 & dpto$pib1_PIB1s_lag >= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 1
dpto@data[(dpto$pib1_PIB1s <= 0 & dpto$pib1_PIB1s_lag <= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 2
dpto@data[(dpto$pib1_PIB1s >= 0 & dpto$pib1_PIB1s_lag <= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 3
dpto@data[(dpto$pib1_PIB1s >= 0 & dpto$pib1_PIB1s_lag <= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 4
dpto@data[(dpto$pib1_PIB1s <= 0 & dpto$pib1_PIB1s_lag >= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 5

# Set up the map
breaks <- seq(1, 5, 1)
labels <- c("high-High", "low-Low", "High-Low", "Low-High", "Not Signif.")
np <- findInterval(dpto$quad_sig, breaks)
colors <- c("red", "blue", "lightpink", "skyblue2", "white")
plot(dpto, col = colors[np])
mtext("Local Moran's I", cex = 1.5, side = 3, line = 1)
legend("topleft", legend = labels, fill = colors, bty = "n")

0.2 Datos espaciales en salud

Es común encontrar datos de salud que pueden ser analizados por los investigadores, agregandole valor a estos datos por medio de representaciones cartográficas y herramientas analiticas. Es común tener en cuenta el tamaño de la población y no solo los datos, por lo cual se calculan las tasas.

\[\displaystyle R_i=\dfrac{N_i}{P_i}\]

Donde $i$ hace alusión a la zona sobre la que es calculada la tasa.

\[\displaystyle tasa =\dfrac{Número\ de\ casos}{Población\ en\ riesgo}\]

Usualmente se multiplica por potencias de 10, para ganar interpretabilidad y evitar los problemas en los decimales.

0.2.1 Tasa estandarizadas

\[SMR\ en\ la\ region\ i =\dfrac{R_i}{tasa\ general}\]

0.2.1.1 Casos esperados

\[E_i=R x P_i\]

Casos esperados = tasa general x población

\[SMR=\frac{N_i}{E_i}\] \[\displaystyle Tasa\ estandarizada=\frac{Número\ de\ casos}{Casos\ esperados}\]

Si la tasa es menor que uno, se aprecian menos casos de los esperados.

Usando modelos adecuados, es posible construir intervalos de confianza y es posible calcular zonas donde la SMR sea mayor a un umbral, por ejemplo el doble con algun nivel de confianza.

london<-readRDS("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/Bibliografía/DataCamp Spatial Statistics in R/datasets/london_2017_2.rds", refhook = NULL)
# Get a summary of the data set
summary(london)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##        min      max
## x 503574.2 561956.7
## y 155850.8 200933.6
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy
## +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894]
## Data attributes:
##      CODE               NAME              Flu_OBS          Vom_OBS     
##  Length:32          Length:32          Min.   :  0.00   Min.   : 0.00  
##  Class :character   Class :character   1st Qu.: 11.00   1st Qu.:14.00  
##  Mode  :character   Mode  :character   Median : 33.00   Median :40.00  
##                                        Mean   : 38.56   Mean   :37.28  
##                                        3rd Qu.: 61.00   3rd Qu.:57.50  
##                                        Max.   :112.00   Max.   :96.00  
##    Diarr_OBS        Gastro_OBS      TOTAL_POP        CCGcode         
##  Min.   :  0.00   Min.   :  0.0   Min.   :157711   Length:32         
##  1st Qu.: 22.50   1st Qu.: 48.0   1st Qu.:237717   Class :character  
##  Median : 62.00   Median :120.5   Median :272017   Mode  :character  
##  Mean   : 57.03   Mean   :113.7   Mean   :270780                     
##  3rd Qu.: 90.75   3rd Qu.:176.8   3rd Qu.:316911                     
##  Max.   :122.00   Max.   :251.0   Max.   :379691                     
##  CCG.geography.code   CCG.name         Asthma_Prevalence
##  Length:32          Length:32          Min.   :3.550    
##  Class :character   Class :character   1st Qu.:4.412    
##  Mode  :character   Mode  :character   Median :4.660    
##                                        Mean   :4.624    
##                                        3rd Qu.:4.925    
##                                        Max.   :5.470    
##  Obesity_Prevalence Cancer_Prevalence Diabetes_Prevalence     Income      
##  Min.   : 3.930     Min.   :0.870     Min.   :3.620       Min.   :0.0730  
##  1st Qu.: 5.855     1st Qu.:1.438     1st Qu.:5.265       1st Qu.:0.1308  
##  Median : 7.565     Median :1.605     Median :6.305       Median :0.1665  
##  Mean   : 7.585     Mean   :1.684     Mean   :6.245       Mean   :0.1655  
##  3rd Qu.: 8.810     3rd Qu.:1.903     3rd Qu.:7.067       3rd Qu.:0.1985  
##  Max.   :12.130     Max.   :2.540     Max.   :9.060       Max.   :0.2530  
##    Employment       Education      HealthDeprivation     Crime        
##  Min.   :0.0570   Min.   : 3.958   Min.   :-1.4100   Min.   :-0.1550  
##  1st Qu.:0.0905   1st Qu.:10.047   1st Qu.:-0.5055   1st Qu.: 0.3745  
##  Median :0.1095   Median :13.925   Median :-0.2050   Median : 0.5515  
##  Mean   :0.1092   Mean   :14.024   Mean   :-0.2044   Mean   : 0.5379  
##  3rd Qu.:0.1283   3rd Qu.:17.480   3rd Qu.: 0.2010   3rd Qu.: 0.7823  
##  Max.   :0.1560   Max.   :27.182   Max.   : 0.5430   Max.   : 1.0190  
##     Services      Environment          i        
##  Min.   :19.63   Min.   :13.37   Min.   : 0.00  
##  1st Qu.:24.43   1st Qu.:24.03   1st Qu.: 7.75  
##  Median :30.41   Median :28.20   Median :15.50  
##  Mean   :29.55   Mean   :31.38   Mean   :15.50  
##  3rd Qu.:34.74   3rd Qu.:40.15   3rd Qu.:23.25  
##  Max.   :41.89   Max.   :55.00   Max.   :31.00

# Map the OBServed number of flu reports
spplot(london, "Flu_OBS")

# Compute and print the overall incidence of flu
r <- sum(london$Flu_OBS) / sum(london$TOTAL_POP)
r

## [1] 0.0001424128

# Calculate the expected number for each borough
london$Flu_EXP <- london$TOTAL_POP * r

# Calculate the ratio of OBServed to EXPected
london$Flu_SMR <- london$Flu_OBS / london$Flu_EXP

# Map the SMR
spplot(london, "Flu_SMR")

0.2.2 Intervalos de confianza binomiales

Los SMR por encima de 1 representan altas tasas de enfermedad, pero ¿cuán alto debe ser un SMR antes de que pueda ser considerado estadísticamente significativo?

Dado un número de casos y una población, es posible calcular intervalos de confianza a algún nivel de la estimación de la proporción de casos por población utilizando las propiedades de la distribución binomial. El paquete epitools tiene una función binom.exact() que puede utilizar para calcular los intervalos de confianza para los datos de la gripe. Estos pueden ser escalados para que sean intervalos de confianza en el SMR dividiéndolos por la tasa global.

# For the binomial statistics function
# install.packages("epitools")
library(epitools)

# Get CI from binomial distribution
flu_ci <- binom.exact(london$Flu_OBS, london$TOTAL_POP)

# Add borough names
flu_ci$NAME <- london$NAME

# Calculate London rate, then compute SMR
r <- sum(london$Flu_OBS) / sum(london$TOTAL_POP)
flu_ci$SMR <- flu_ci$proportion / r

# Subset the high SMR data
flu_high <- flu_ci[flu_ci$SMR > 1, ]

# Plot estimates with CIs
library(ggplot2)
ggplot(flu_high, aes(x = NAME, y = proportion / r,
                     ymin = lower / r, ymax = upper / r)) +
  geom_pointrange() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

0.2.3 Probabilidades de superación

Las distribuciones y los intervalos de confianza pueden ser cosas difíciles de presentar a los no estadísticos. Una alternativa es presentar una probabilidad de que un valor esté por encima de un umbral. Por ejemplo, los equipos de salud pública podrían estar interesados en saber cuándo una Tasa estandarizada ha duplicado con creces, y como estadístico, usted puede calcular la probabilidad de que esto ocurra. Para que posteriormente, el equipo de salud pública podría decidir emitir una alerta cuando la probabilidad de que se duplique el tasa estandatizada sea superior a 0,95.

Las propiedades de la distribución binomial permiten calcularla para datos que son proporciones. Luego puede mapear estas probabilidades de superación para algún umbral, y usar una paleta de colores adecuada para resaltar probabilidades cercanas a 1.

# Probability of a binomial exceeding a multiple
binom.exceed <- function(observed, population, expected, e){
    1 - pbinom(e * expected, population, prob = observed / population)
}

# Compute P(rate > 2)
london$Flu_gt_2 <- binom.exceed(
            observed = london$Flu_OBS,
            population = london$TOTAL_POP,
            expected = london$Flu_EXP,
            e = 2)

# Use a 50-color palette that only starts changing at around 0.9
pal <- c(
  rep("#B0D0B0", 40),
  colorRampPalette(c("#B0D0B0", "orange"))(5), 
  colorRampPalette(c("orange", "red"))(5)
)

# Plot the P(rate > 2) map
spplot(london, "Flu_gt_2", col.regions = pal, at = seq(0, 1, len = 50))

0.2.3.1 Taller

Seleccione la causa de muerte empleada en el primer taller. Calcule la razon estandarizada de mortalidad

Calcule el I de Moran
Calcule los LISA
Realice la partición en cuartiles y compare
Calcule los intervalos de confianza Binomiales
Calcule las Probabilidades de superación
Interprete y Arriesgue una solución, si encuentra algún problema.

load("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/avpp13.Rdata")
# names(defun_avpp13)
# table(defun_avpp13$ops_un_digito)


# defun_avpp_neoplasia<-defun_avpp13[defun_avpp13$ops_un_digito==1,]
library(tidyverse)

## -- Attaching packages ------------------------------------------------------------------------------------------------------------ tidyverse 1.2.1 --

## v tibble  1.4.2     v purrr   0.2.5
## v tidyr   0.8.1     v dplyr   0.7.6
## v readr   1.1.1     v stringr 1.3.1
## v tibble  1.4.2     v forcats 0.3.0

## -- Conflicts --------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(readxl)

# defun_avpp13$ops_tres_digitos<-as.numeric(defun_avpp13$ops_desc_tres_digitos)

base<-defun_avpp13 %>% 
  #filter(ops_un_digito==5)%>%
  group_by(cod_mpio,ops_desc_un_digito)%>%
  #group_by(cod_mpio,ops_tres_digitos)%>%
  summarise(
    n=n())

base_rs<-reshape(data.frame(base), idvar = "cod_mpio", timevar = "ops_desc_un_digito", direction = "wide")

base_rs[is.na(base_rs)]<-0
pob<-read_xls("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/poblacion2013.xls",sheet = 1)

base<-merge(base_rs, pob, by="cod_mpio")

write.csv(x = base,"D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/667pob.csv")

cod_mpio
n.CAUSAS EXTERNAS
n.CIERTAS AFECCIONES ORIGINADA EN EL PERIODO PERINATAL
n.ENFERMEDADES DEL SISTEMA CIRCULATORIO
n.ENFERMEDADES TRANSMISIBLES
n.NEOPLASIAS (TUMORES)
n.SINTOMAS SIGNOS Y AFECCIONES MAL DEFINIDAS
n.TODAS LAS DEMAS CAUSAS
DPMP
MPIO
poblacion

mpio <- readShapePoly("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/shp/mpioopsundigito.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

## Field name: '667pob_fie' changed to: 'X667pob_fie'

## Field name: '667pob_n.C' changed to: 'X667pob_n.C'

## Field name: '667pob_n_1' changed to: 'X667pob_n_1'

## Field name: '667pob_n.E' changed to: 'X667pob_n.E'

## Field name: '667pob_n_2' changed to: 'X667pob_n_2'

## Field name: '667pob_n.N' changed to: 'X667pob_n.N'

## Field name: '667pob_n.S' changed to: 'X667pob_n.S'

## Field name: '667pob_n.T' changed to: 'X667pob_n.T'

## Field name: '667pob_DPM' changed to: 'X667pob_DPM'

## Field name: '667pob_MPI' changed to: 'X667pob_MPI'

## Field name: '667pob_pob' changed to: 'X667pob_pob'

“In God we trust, all else must bring data” - Art Dempster