Processing math: 100%
  • Advertencias
  • Contexto
  • Objetivos
  • INICIO
  • 07 - análisis espacial
    • instalar paquetes
    • librerias a usar
    • importar base de datos
    • _preparar base de datos
      • casos observados
      • variables exposición
      • agregar polígonos
    • _calcular casos esperados
    • _calcular SMR
      • crear un gráfico con ggplot2
    • _modelo de regresión espacial
      • extraer matriz de vecindad
      • realizar inferencia
      • visualizar resultados
  • RETROALIMENTACIÓN

Advertencias

Este documento es una continuación de Análisis de datos en Vigilancia Epidemiológica I: tiempo, espacio, persona y curva epidémica. Retorna a dicho documento para recordar las advertencias.

Este documento está vinculado a las siguientes diapositivas

Material original de “Small Area Disease Risk Estimation and Visualization Using R” Paula Moraga, The R Journal (2018) 10:1, pages 495-506 actualizado por la misma autora aquí. Adaptación disponible aquí

Contexto

En este tutorial usamos metodos espaciales para estimar el riesgo de cancer de pulmon en Pensilvania, Estados Unidos, en el año 2002.

Usamos datos del paquete de R SpatialEpi con la poblacion, los casos de cancer de pulmón, y las proporciones de fumadores en cada uno de los condados de Pensilvania.

Los datos de poblacion se han obtenido del censo del año 2000, y los casos de cancer y la proporción de fumadores de la pagina web del Departamento de Salud de Pensilvania.

Objetivos

  • Calcular por cada Condado/Distrito dentro de un Estado/Departamento el: (i )número de casos observados y esperados usando el paquete SpatialEpi, y (ii) las razones de mortalidad estandardizadas (Standardized Mortality Ratios o SMR).

  • Obtener estimaciones del riesgo relativo de la enfermedad y cuantificar factores de riesgo usando el paquete INLA.

  • Crear un mapa con las estimaciones de riesgo usando ggplot2.

INICIO

07 - análisis espacial

instalar paquetes

  • Aquí la lista de los paquetes requeridos:
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("sf")) install.packages("sf")
if(!require("SpatialEpi")) install.packages("SpatialEpi")
if(!require("spdep")) install.packages("spdep")
if(!require("leaflet")) install.packages("leaflet")

  • El paquete INLA no está en CRAN porque utiliza algunas librarías externas que dificultan la construcción de los binarios.

  • Por lo tanto, necesitamos instalarlo añadiendo la URL del repositorio INLA:

if(!require("INLA")) install.packages("INLA", repos = "https://inla.r-inla-download.org/R/stable", dep = TRUE)

librerias a usar

library(tidyverse)
library(sf)

importar base de datos

  • Extraer los datos pennLC que vienen en el paquete SpatialEpi
library(SpatialEpi)

# llamar a base de datos dentro de paquete
data(pennLC)

P: ¿Qué tipo de objeto es pennLc? ¿Qué elementos contiene?

# ¿qué tipo de objeto es?


# identificar sus atributos 
attributes(pennLC)

_preparar base de datos

casos observados

  • Usa group_by() %>% summarise() para crear una base de datos resumen con la suma sum de casos (Y) y suma total de población (pop) por condado/distrito
pennLC$data %>% 
  as_tibble() %>% 
  group_by(county) %>% 
  summarise(Y=sum(cases),
            pop=sum(population)) %>% 
  ungroup()

variables exposición

  • Unir con left_join() la proporción de fumadores smoking por condado/distrito
pennLC$smoking

agregar polígonos

  • Unir con left_join() los polígonos de cada condado/distrito

  • Utilizar el número de orden de cada condado/distrito para unir la información

(en el mejor de los casos, deberían disponer de los ubigeos)

pennLC$spatial.polygon %>% 
              st_as_sf() %>% 
              as_tibble() %>% 
              rownames_to_column()
  • Asignar el resultado final al objeto d
#outcome at county level + covariates (race,sex,age_strata)
d <- pennLC$data %>% 
  as_tibble() %>% 
  group_by(county) %>% 
  summarise(Y=sum(cases),
            pop=sum(population)) %>% 
  ungroup() %>% 
  
  #add exposure at county level
  left_join(pennLC$smoking) %>% 
  rownames_to_column() %>% 
  
  #add polygon to tibble
  left_join(pennLC$spatial.polygon %>% 
              st_as_sf() %>% 
              as_tibble() %>% 
              rownames_to_column())
  • Visualizar casos observados
d %>% 
  st_as_sf() %>% 
  st_transform(crs = 27700) %>% 
  ggplot(aes(fill=Y)) +
  geom_sf() +
  scale_fill_gradient(low = "white",high = "red") +
  guides(fill = guide_colorbar(reverse=F)) +
  theme_bw() +
  labs(title = "Observed cases",
       subtitle = "Political Map: Polygon boundaries for each county")

_calcular casos esperados

  • Ojo: El cálculo tomará en cuenta todos los estratos

  • Asigna la base de datos al objeto e

# Expected cases (using stratas!) ---------------------------------
e <- pennLC$data %>% 
  as_tibble() %>% 
  select(county,cases,population,race,gender,age)
  • Usar count() para identificar la cantidad de combinación de estratos por variable categórica
e %>% count(county) #67

e %>% count(race,gender,age) #16
  • Calcular los casos esperados usando la función expected()
population <- e$population
cases <- e$cases
n.strata <- 16 # = 2 races * 2 genders * 4 age bands
E <- expected(population, cases, n.strata)

P: Explorar E, ¿qué longitud tiene? ¿a qué nivel de agregación corresponden estos resultados?

E

_calcular SMR

  • Agregar con mutate() el número de casos esperados (E) a una lísta de condados/distritos
e %>% 
  select(county) %>% 
  distinct() %>% 
  mutate(E=E)
  • Unir con left_join() la base resumen d con la cantidad de casos esperados E
d %>% 
  left_join(
    e %>% 
      select(county) %>% 
      distinct() %>% 
      mutate(E=E)
  )
  • Usar mutate() para crear la variable SMR
d %>% 
  left_join(
    e %>% 
      select(county) %>% 
      distinct() %>% 
      mutate(E=E)
  ) %>% 
  mutate(SMR=Y/E)

  • Usar la función st_as_sf() para crear un objeto tipo sf

  • Usar la función as('Spatial') para transformar objeto sf a SpatialPolygonDataFrame

  • Asignar resultado final al objeto map

map <- d %>% 
  left_join(
    e %>% 
      select(county) %>% 
      distinct() %>% 
      mutate(E=E)
  ) %>% 
  mutate(SMR=Y/E) %>% 
  # add data to map 
  st_as_sf() %>% 
  # transform to SpatialPolygonDataFrame
  as('Spatial')
## Joining, by = "county"

P: ¿Qué información clave caracteriza a un objeto espacial?

crear un gráfico con ggplot2

  • Usar la plantilla de ggplot2 para crear un mapa
# Mapping SMR ---------------------------------
map %>% 
  st_as_sf() %>% 
  st_transform(crs = 27700) %>% 
  ggplot(aes(fill=SMR)) +
  geom_sf() +
  scale_fill_gradient2(high = "#e34a33",low = "#2b8cbe",midpoint = 1) +
  guides(fill = guide_colorbar(reverse=F)) +
  theme_bw() +
  labs(title = "Standardize Mortality Ratio",
       subtitle = "Political Map: Polygon boundaries for each county")

# alternative - cartogram
map_carto <- map %>% 
  st_as_sf() %>% 
  st_transform(crs = 27700) %>% 
  cartogram::cartogram_cont("pop", itermax=5)

map_carto %>% 
  st_as_sf() %>% 
  ggplot(aes(fill=SMR)) +
  geom_sf() +
  scale_fill_gradient2(high = "#e34a33",low = "#2b8cbe",midpoint = 1) +
  guides(fill = guide_colorbar(reverse=F)) +
  theme_bw() +
  labs(title = "Standardize Mortality Ratio",
       subtitle = "Cartogram: polygon area proportional to the population size")

P: ¿Qué problema observa con esta visualización?

_modelo de regresión espacial

extraer matriz de vecindad

  • Usar función poly2nb() de paquete spdep para extraer polígonos vecinos

  • nb2INLA() de paquete INLA para generar archivo local a ser usado por inla()

# Neighbourhood matrix ---------------------------------
library(spdep)
library(INLA)
nb <- poly2nb(map)
head(nb)
## [[1]]
## [1] 21 28 67
## 
## [[2]]
## [1]  3  4 10 63 65
## 
## [[3]]
## [1]  2 10 16 32 33 65
## 
## [[4]]
## [1]  2 10 37 63
## 
## [[5]]
## [1]  7 11 29 31 56
## 
## [[6]]
## [1] 15 36 38 39 46 54
nb2INLA("map.adj", nb)
g <- inla.read.graph(filename = "map.adj")

realizar inferencia

  • indexar efectos aleatorios

  • crear fórmula

  • ejecutar la función inla

# Inference using INLA ---------------------------------

#index both random effects
map$re_u <- 1:nrow(map@data) #spatial residual variation
map$re_v <- 1:nrow(map@data) #modeling unstructure noise

#create formula
#iid: independent and identically distributed
formula <- 
  Y ~ # outcome
  smoking + # exposure
  f(re_u, model = "besag", graph = g) + #random effect - spatial
  f(re_v, model = "iid") #random effect - noise

res <- inla(formula, 
            family = "poisson", 
            data = map@data, 
            E = E,
            control.predictor = list(compute = TRUE))

visualizar resultados

  • Inspeccionar la salida resumen en crudo del modelo
summary(res)
## 
## Call:
##    c("inla(formula = formula, family = \"poisson\", data = map@data, ", " 
##    E = E, control.predictor = list(compute = TRUE))") 
## Time used:
##     Pre = 1.98, Running = 2.64, Post = 1.16, Total = 5.78 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.323 0.150     -0.621   -0.323     -0.028 -0.323   0
## smoking      1.156 0.624     -0.080    1.157      2.384  1.160   0
## 
## Random effects:
##   Name     Model
##     re_u Besags ICAR model
##    re_v IID model
## 
## Model hyperparameters:
##                        mean       sd 0.025quant 0.5quant 0.975quant    mode
## Precision for re_u    92.67    50.19      31.00    81.16     221.44   62.96
## Precision for re_v 18097.70 18216.37    1165.24 12665.38   66267.31 3144.29
## 
## Expected number of effective parameters(stdev): 18.66(4.43)
## Number of equivalent replicates : 3.59 
## 
## Marginal log-Likelihood:  -320.54 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

  • Extraer la distribución posterior del coeficiente β1

  • Visualizar con ggplot2

marginal <- inla.smarginal(res$marginals.fixed$smoking) %>% 
  data.frame()

ggplot(marginal, aes(x = x, y = y)) + 
  geom_line() +
  geom_vline(xintercept = 0, col = "blue") + 
  labs(x = expression(beta[1]), 
       y = "Density") +
  theme_bw()

  • Extraer predictores lineales del desenlace del modelo

  • Visualizar con ggplot2

head(res$summary.fitted.values)
##                          mean         sd 0.025quant  0.5quant 0.975quant
## fitted.Predictor.01 0.8796107 0.05846419  0.7639557 0.8798268  0.9947103
## fitted.Predictor.02 1.0598951 0.02765373  1.0069360 1.0594517  1.1153543
## fitted.Predictor.03 0.9627956 0.05193450  0.8552488 0.9644781  1.0611516
## fitted.Predictor.04 1.0268012 0.05129911  0.9266835 1.0264757  1.1289543
## fitted.Predictor.05 0.9078911 0.05495643  0.7984740 0.9082307  1.0161501
## fitted.Predictor.06 0.9953488 0.04033642  0.9184137 0.9944626  1.0774326
##                          mode
## fitted.Predictor.01 0.8807642
## fitted.Predictor.02 1.0585487
## fitted.Predictor.03 0.9679679
## fitted.Predictor.04 1.0260002
## fitted.Predictor.05 0.9094453
## fitted.Predictor.06 0.9928306
map$RR <- res$summary.fitted.values[, "mean"]
map$LL <- res$summary.fitted.values[, "0.025quant"]
map$UL <- res$summary.fitted.values[, "0.975quant"]

map %>% 
  st_as_sf() %>% 
  st_transform(crs = 27700) %>% 
  ggplot(aes(fill=RR)) +
  geom_sf() +
  scale_fill_gradient2(high = "#e34a33",low = "#2b8cbe",midpoint = 1) +
  guides(fill = guide_colorbar(reverse=F)) +
  labs(fill = "linear\npredictor") +
  theme_bw()

map_carto <- map %>% 
  st_as_sf() %>% 
  st_transform(crs = 27700) %>% 
  cartogram::cartogram_cont("pop", itermax=5)

map_carto %>% 
  st_as_sf() %>% 
  ggplot(aes(fill=RR)) +
  geom_sf() +
  scale_fill_gradient2(high = "#e34a33",low = "#2b8cbe",midpoint = 1) +
  guides(fill = guide_colorbar(reverse=F)) +
  labs(fill = "linear\npredictor") +
  theme_bw()

RETROALIMENTACIÓN

08 - cierre

¡Gracias por su atención!

Andree Valle Campos

@avallecam click aquí

código desarrollado