En este tutorial usamos metodos espaciales para estimar el riesgo de cancer de pulmon en Pensilvania, Estados Unidos, en el anyo 2002. Usamos datos del paquete de R SpatialEpi con la poblacion, los casos de cancer de pulmon, y las proporciones de fumadores en cada uno de los condados de Pensilvania. Los datos de poblacion se han obtenido del censo del anyo 2000, y los casos de cancer y la proporcion de fumadores de la pagina web del Departamento de Salud de Pensilvania.

Mostramos como calcular el numero de casos observados y esperados, y las razones de mortalidad estandardizadas (Standardized Mortality Ratios (SMRs)) en cada uno de los condados de Pensilvania. Tambien mostramos como obtener estimaciones del riesgo relativo de la enfermedad y como cuantificar factores de riesgo usando INLA. Finalmente, mostramos como hacer mapas interactivos de las estimaciones del riesgo usando leaflet.

Una version extendida de este ejemplo aparece en el articulo Moraga, P. Small Area Disease Risk Estimation and Visualization Using R. The R Journal, 10(1):495-506, 2018, y el libro Moraga, P. Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny. Chapman & Hall/CRC, 2019.

1 Datos y mapa

Empezamos cargando el paquete SpatialEpi y adjuntado los datos.

library(SpatialEpi)
data(pennLC)

A continuacion, inspeccionamos los datos.

class(pennLC)
## [1] "list"
names(pennLC)
## [1] "geo"             "data"            "smoking"         "spatial.polygon"
?pennLC

Vemos que pennLC es una lista con los siguientes elementos:

  • geo: una tabla con los ids de los condados, y la longitud y latitud de los centroides de cada condado,
  • data: una tabla con los ids de los condados, numero de casos y numero de habitantes para cada uno de los stratos de la poblacion,
  • smoking: una tabla con los ids de los condados y las proporciones de fumadores,
  • spatial.polygon: un objeto SpatialPolygons con los condados de Pensilvania.

Inspeccionamos las primeras filas de pennLC$data y pennLC$smoking. pennLC$data contiene el numero de casos de cancer de pulmon y el numero de habitantes en cada condado, para cada raza (blanco y no blanco), genero (femenino y masculino) y edad (menos de 40, 40-59, 60-69 and 70+).

head(pennLC$data)

pennLC$smoking contiene la proporcion de fumadores en cada condado.

head(pennLC$smoking)

pennLC$spatial.polygon es un objeto de clase SpatialPolygons que contiene el mapa con los condados de Pensilvania. Podemos visualizar el mapa de la siguiente manera:

map <- pennLC$spatial.polygon
plot(map)

2 Preparacion de los datos

Ahora creamos un data frame llamado d con columnas que contengan los ids de los condados, el numero de casos observados y esperados, las proporciones de fumadores y las SMRs. En concreto, d tendra las siguientes columnas:

  • county: id de cada condado,
  • Y: numero observado de casos en cada condado,
  • E: numero esperado de casos en cada condado,
  • smoking: proporcion de fumadores en cada condado,
  • SMR: SMR de cada condado.

2.1 Casos observados

pennLC$data contiene los casos en cada condado estratificados por raza, genero y edad. Podemos obtener el numero de casos en cada condado, Y, agregando las filas de pennLC$data por condado y sumando el numero de casos observados usando el paquete dplyr.

library(dplyr)
d <- group_by(pennLC$data, county) %>% summarize(Y = sum(cases))
head(d)

En lugar de usar dplyr, tambien podemos calcular el numero de casos en cada condado usando la funcion aggregate(). aggregate() acepta los siguientes argumentos:

  • x: data frame con los datos,
  • by: lista con los elementos de agrupacion (cada elemento necesita tener la misma longitud que las variables en el data frame en x),
  • FUN: funcion para calcular los estadisticos resumen que se aplican a cada uno de los subconjuntos de datos.

Podemos usar aggregate() con x igual al vector del numero de casos, by igual a una lista con los condados, y FUN igual a sum. A continuacion, ponemos nombres county e Y a las columnas del data frame obtenido.

d <- aggregate(
  x = pennLC$data$cases,
  by = list(county = pennLC$data$county),
  FUN = sum
)
names(d) <- c("county", "Y")

2.2 Casos esperados

Ahora calculamos el numero de casos esperados a cada condado. Los casos esperados representan el numero de casos de la enfermedad que se esperaria obtener si la poblacion de cada condado se comportara de la misma manera que la poblacion de Pensilvania. El numero de casos esperados \(E_i\) en cada condado \(i\) se calcula usando estandardizacion indirecta como \[E_i=\sum_{j=1}^m r_j^{(s)} n_j^{(i)},\] donde \(r_j^{(s)}\) es la tasa (numero de casos divido por la poblacion) en el estrato \(j\) de la poblacion de Pensilvania, y \(n_j^{(i)}\) es la poblacion en el estrato \(j\) del condado \(i\).

Podemos calcular el numero de casos esperados usando la funcion expected() del paquete SpatialEpi. Esta funcion tiene tres argumentos:

  • population: vector con la poblacion en cada estrato en cada area,
  • cases: vector con el numero de casos es cada estrato en cada area,
  • n.strata: numero de estratos considerados.

Los vectores population and cases tienen que ordenarse primero por area, y luego, dentro de cada area, el numero de casos de cada estrato tienen que listarse en el mismo orden. Todos los estratos necesitan incluirse en los vectores, incluyendo estratos con 0 casos.

Para obtener el numero de casos esperados primero ordenamos los datos usando la funcion order() donde especificamos el orden como condado, raza, genero y finalmente edad.

pennLC$data <- pennLC$data[order(pennLC$data$county, pennLC$data$race, pennLC$data$gender, pennLC$data$age), ]

A continuacion, obtenemos el numero de casos esperados E en cada condado llamando a la funcion expected() donde asignamos population igual a pennLC$data$population y cases igual a pennLC$data$cases. Hay 2 razas, 2 generos y 4 grupos de edad en cada condado, asi que el numero de estratos es 2 x 2 x 4 = 16.

E <- expected(population = pennLC$data$population, cases = pennLC$data$cases, n.strata = 16)

Ahora anyadimos el vector E al data frame d que pasa a tener los ids de los condados (county), el numero de casos observados (Y), y el numero de casos esperados (E).

d$E <- E
head(d)

2.3 Proporciones de fumadores

Ahora anyadimos a d la variable smoking que representa la proporcion de fumadores en cada condado. Anyadimos esta variable usando la funcion merge() donde especificamos county como la columna de union.

d <- merge(d, pennLC$smoking, by = "county")

2.4 SMRs

Para cada condado \(i\), la SMR es el numero de casos observados dividido el numero de casos esperados, \[\mbox{SMR}_i=Y_i/E_i.\] Calculamos el vector de SMRs como el numero de casos observados dividido el numero de casos esperados, y lo anyadimos al data frame d.

d$SMR <- d$Y/d$E

2.5 Datos

El data frame d construido es el siguiente:

head(d)

2.6 Anyadir datos al mapa

Hemos construidos un data frame d con el numero de casos observados y esperados de la enfermedad, las proporciones de fumadores y la SMR para cada uno de los condados. El mapa de los condados de Pensilvania esta en el objeto map de clase SpatialPolygons. Usando este objeto y el data frame d podemos crear un SpatialPolygonsDataFrame que nos permitira crear mapas de las variables en d. Para hacer eso, primero asignamos el nombre de las filas del data frame d igual a d$county. A continuacion, unimos el SpatialPolygons map y el data frame d uniendo los valores de los slots de los ID de los Polygons SpatialPolygons con el nombre de las filas del data frame (match.ID = TRUE).

library(sp)
rownames(d) <- d$county
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE)
head(map@data)

3 Mapa de las SMR

Podemos visualizar el numero de casos observados y esperados de la enfermedad, las SMRs, y las proporciones de fumadores en un mapa interactivo. Construimos el mapa con el paquete leaflet. Primero usamos la funcion leaflet() y anyadimos al mapa las tiles OpenStreetMap que vienen por defecto usando addTiles(). A continuacion, anyadimos los condados de Pensilvania con addPolygons() donde especificamos el color de los bordes de las areas (color) y el grosor del borde (weight). Rellenamos las areas con los colores de la funcion paleta de colores generada con colorNumeric(), y asignamos fillOpacity igual a un numero menor que 1 para poder ver el mapa del background. Usamos colorNumeric() para crear una funcion paleta de colores que mapea los valores de los datos a los colores de la paleta. Creamos la funcion paleta de colores usando los siguientes dos parametros:

  • palette: funcion color con la que se mostraran los valores, y
  • domain: valores que se pueden mostrar.

Finalmente, anyadimos la leyenda especificando la funcion paleta de color (pal) y los valores usados para generar con la funcion paleta (values). Asignamos opacity el mismo valor que la opacidad de las areas, y especificamos un titulo y una posicion para la leyenda.

library(leaflet)
l <- leaflet(map) %>% addTiles()

pal <- colorNumeric(palette = "YlOrRd", domain = map$SMR)

l %>% addPolygons(color = "grey", weight = 1, fillColor = ~pal(SMR), fillOpacity = 0.5) %>%
  addLegend(pal = pal, values = ~SMR, opacity = 0.5, title = "SMR", position = "bottomright")

Podemos mejorar este mapa resaltando los condados cuando pasamos el raton por encima, y mostrando informacion de los casos observados y esperados, SMRs, y proporciones de fumadores. Para hacer esto, anyadimos los argumentos highlightOptions, label y labelOptions en addPolygons(). Escogemos resaltar las areas usando un grosor mayor (highlightOptions(weight = 4)). Creamos etiquetas de texto con los valores que queremos mostrar usando sintaxis HTML. Primero, creamos el texto que mostraremos usando la funcion sprintf() la cual devuelve un vector que contiene una combinacion de texto y valores de la variables, y despues aplicamos htmltools::HTML() para marcar el texto como HTML. En labelOptions especificamos style, textsize, and direction. Los valores posibles para direction son left, right y auto y especifican la direccion en que la etiqueta se mostrara en relacion al marker. Escogemos auto para que la direccion optima se escoja dependiendo de la posicion del marker.

labels <- sprintf("<strong>%s</strong><br/>Observed: %s <br/>Expected: %s <br/>Smokers proportion: %s <br/>SMR: %s",
  map$county, map$Y,  round(map$E, 2), map$smoking, round(map$SMR, 2)) %>%
  lapply(htmltools::HTML)

l %>% addPolygons(color = "grey", weight = 1, fillColor = ~pal(SMR), fillOpacity = 0.5,
    highlightOptions = highlightOptions(weight = 4),
    label = labels,
    labelOptions = labelOptions(style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px", direction = "auto")) %>%
    addLegend(pal = pal, values = ~SMR, opacity = 0.5, title = "SMR", position = "bottomright")

Ahora podemos examinar el mapa y ver que areas tienen

  • SMR \(=\) 1 indicando mismos casos observados que esperados,
  • SMR \(>\) 1 indicando mas casos observados que esperados,
  • SMR \(<\) 1 indicando menos casos observados que esperados.

El mapa muestra las SMRs y da una idea del riesgo del cancer de pulmon en Pensilvania. Sin embargo, las SMRs pueden ser enganyosas y poco fiables en los condados con poblaciones pequenyas. A continuacion, veremos como podemos usar INLA para ajustar un modelo que permite incorporar covariables y utilizar informacion de condados vecinos para mejorar las estas estimaciones y producir estimaciones de riesgo de la enfermedad suavizados.

4 Modelizacion

En esta seccion especificamos el modelo para estimar el riesgo de cancer de pulmon en Pensilvania y detallamos los pasos requeridos para ajustar el modelo y visualizar los resultados usando INLA.

4.1 Modelo

Sean \(Y_i\) y \(E_i\) el numero de casos de cancer observados y esperados, respectivamente, y sea \(\theta_i\) el riesgo relativo del condado \(i=1,\ldots,n\). El modelo se especifica como sigue:

\[Y_i|\theta_i \sim Poisson(E_i \times \theta_i),\ i=1,\ldots,n,\]

\[\log(\theta_i) = \beta_0 + \beta_1 \times smoking_i + u_i + v_i.\]

  • \(\beta_0\) intercept,
  • \(\beta_1\) coeficiente de la covariable proporcion de fumadores,
  • \(u_i\) efecto espacial estructurado, \(u_i|\mathbf{u_{-i}} \sim N(\bar{u}_{\delta_i}, \frac{1}{\tau_u n_{\delta_i}})\) (Conditionally autoregressive model (CAR)),
  • \(v_i\) efecto no estructurado, \(v_i \sim N(0, 1/\tau_v)\).

4.2 Matriz de vecindad

Creamos la matriz de vecindad que se necesita para definir el efecto aleatorio espacial usando las funciones poly2nb() y nb2INLA() del paquete spdep. Primero, usamos poly2nb() para crear una lista de vecinos donde los vecinos de una determinada area son las areas con las que comparte bordes. Cada elemento de la lista nb representa un area y contiene los indices de sus vecinos. Por ejemplo, nb[[2]] contiene los vecinos del area 2.

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

A continuacion, usamos la funcion nb2INLA() para transformar la lista nb en un archivo llamado map.adj con la representacion de la matriz vecindad que se requiere para usar INLA. Esta funcion guarda el archivo map.adj en el directorio de trabajo que puede obtenerse con getwd(). A continuacion, leemos el archivo map.adj usando la funcion inla.read.graph() de INLA, y lo asignamos al objeto g que se usara mas tarde para especificar el modelo espacial con INLA.

nb2INLA("map.adj", nb)
g <- inla.read.graph(filename = "map.adj")

4.3 Inferencia usando INLA

El modelo incluye dos efectos aleatorios, concretamente, \(u_i\) para modelizar variabilidad residual espacial, y \(v_i\) para modelizar ruido no estructurado. Necesitamos incluir dos vectores en el data frame para representar los indices de los efectos aleatorios. Llamamos re_u al vector de los indices para \(u_i\), y re_v al vector de los indices para \(v_i\). Asignamos re_u y re_v igual a \(1,\ldots,n\), donde \(n\) es el numero de condados. En nuestro ejemplo, \(n\)=67 y puede ser obtenido con el numero de filas en los datos (nrow(map@data)).

map$re_u <- 1:nrow(map@data)
map$re_v <- 1:nrow(map@data)

Especificamos la formula del modelo incluyendo la respuesta a la izquierda, y los efectos fijos y aleatorios a la derecha. Los efectos aleatorios se asignan usando la funcion f() con argumentos igual al nombre de la variable y el modelo escogido. Para \(u_i\), usamos model = "besag" con matriz de vecindad dada por g. Para \(v_i\) escogemos model = "iid".

formula <- Y ~ smoking + f(re_u, model = "besag", graph = g, scale.model = TRUE) + f(re_v, model = "iid")

Finalmente, ajustamos el modelo ejecutando la funcion inla() y usando las priors por defecto de INLA. Especificamos la formula, familia, datos, y el numero de casos esperados, y asignamos control.predictor igual a list(compute = TRUE) para obtener las distribuciones a posteriori de las predicciones.

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

4.4 Resultados

Inspeccionamos el objecto de los resultados res usando summary().

summary(res)
## 
## Call:
## c("inla(formula = formula, family = \"poisson\", data = map@data, ",  "    E = E, control.predictor = list(compute = TRUE))")
## 
## Time used:
##  Pre-processing    Running inla Post-processing           Total 
##           0.953           1.056           0.369           2.378 
## 
## Fixed effects:
##                mean     sd 0.025quant 0.5quant 0.975quant    mode kld
## (Intercept) -0.3233 0.1502    -0.6207  -0.3231    -0.0277 -0.3228   0
## smoking      1.1556 0.6244    -0.0808   1.1568     2.3841  1.1600   0
## 
## Random effects:
## Name   Model
##  re_u   Besags ICAR model 
## re_v   IID model 
## 
## Model hyperparameters:
##                        mean       sd 0.025quant 0.5quant 0.975quant
## Precision for re_u   230.78   124.81      77.34   202.17     551.13
## Precision for re_v 18148.62 18199.95    1164.16 12725.30   66367.47
##                       mode
## Precision for re_u  156.91
## Precision for re_v 3144.96
## 
## Expected number of effective parameters(std dev): 18.68(4.433)
## Number of equivalent replicates : 3.587 
## 
## Marginal log-Likelihood:  -289.45 
## Posterior marginals for linear predictor and fitted values computed

Vemos que el intercept \(\hat \beta_0=\) -0.323 con un intervalo de credibilidad del 95% igual a (-0.621, -0.028), y el coeficiente de proporcion de fumadores es \(\hat \beta_1=\) 1.156 con un intervalo de credibilidad al 95% igual a (-0.081, 2.384). Hacemos un plot de la distribucion a posteriori del coeficiente de fumadores calculando un suavizado de la distribucion marginal del coeficiente con inla.smarginal() y haciendo un plot con la funcion ggplot() del paquete ggplot2.

library(ggplot2)
marginal <- inla.smarginal(res$marginals.fixed$smoking)
marginal <- data.frame(marginal)
ggplot(marginal, aes(x = x, y = y)) + geom_line() + labs(x = expression(beta[1]), y = "Density") +
  geom_vline(xintercept = 0, col = "blue") + theme_bw()

4.5 Anyadir resultados al mapa

Las estimaciones del riesgo de cancer de pulmon y la incertidumbre para cada uno de los condados vienen dadas por las medias de las distribuciones a posteriori y los intervalos de credibilidad al 95% de \(\theta_i\), \(i=1,\ldots,n\) que estan en res$summary.fitted.values. La columna mean es la media a posteriori y 0.025quant and 0.975quant son los percentiles 2.5 and 97.5, respectivamente.

head(res$summary.fitted.values)

Anyadimos estos datos a map para poder hacer mapas de estas variables. Asignamos mean a la estimacion del riesgo relativo (RR), y 0.025quant and 0.975quant a los limites inferior y superior de los intervalos de credibilidad al 95% de los riesgos relativos.

map$RR <- res$summary.fitted.values[, "mean"]
map$LL <- res$summary.fitted.values[, "0.025quant"]
map$UL <- res$summary.fitted.values[, "0.975quant"]

5 Mapas del riesgo relativo

Mostramos los riesgos relativos de cancer de pulmon en un mapa interactivo usando leaflet. En el mapa, anyadimos etiquetas de texto que aparecen cuando el raton pasa por encima de los condados mostrando informacion sobre el numero de casos observados y esperados, SMRs, proporcion de fumadores, RRs, y limites inferiores y superiores de intervalos de credibilidad al 95%.

Vemos que los condados con mayor riesgo de cancer de pulmon estan localizados en el oeste y el sudeste de Pensilvania, y los condados con riesgo mas bajo estan en el centro. Los intervalos de credibilidad al 95% indican la incertidumbre de las estimaciones del riesgo.

pal <- colorNumeric(palette = "YlOrRd", domain = map$RR)

labels <- sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
                  Smokers proportion: %s <br/>SMR: %s <br/>RR: %s (%s, %s)",
                  map$county, map$Y,  round(map$E, 2),  map$smoking, round(map$SMR, 2),
                  round(map$RR, 2), round(map$LL, 2), round(map$UL, 2)) %>%
  lapply(htmltools::HTML)

leaflet(map) %>% addTiles() %>%
    addPolygons(color = "grey", weight = 1, fillColor = ~pal(RR),  fillOpacity = 0.5,
    highlightOptions = highlightOptions(weight = 4),
    label = labels,
    labelOptions = labelOptions(style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px", direction = "auto")) %>%
    addLegend(pal = pal, values = ~RR, opacity = 0.5, title = "RR", position = "bottomright")

6 Ejercicio

Modelizacion espacial del riesgo de cancer de labio en hombres en Escocia, 1975-1980.

Datos:

map contiene un mapa con los ids de los condados (areaid), numero de casos observados (Y), numero de casos esperados (E), y una covariable (covariate) para cada uno de los condados de Escocia. La covariable representa la proporcion de la poblacion que se dedica a agricultura, pesca, o ingenieria forestal (AFF).

  • Calcula las razones de incidencia estandardizada (Standardized Incidence Ratios (SIRs)) en cada uno de los condados. Crea un mapa leaflet que muestre los SIRs.
  • Ajusta un modelo espacial para obtener los riesgos relativos en cada uno de los condados. El predictor linear del modelo debe incluir un intercept, una covariable, un efecto aleatorio espacial, y un efecto aleatorio no estructurado.
  • Cual es el efecto de la covariable en el riesgo de la enfermedad?
  • Crea una tabla con columnas id de area, numero de casos observados, numero de casos esperados, covariable, SIR, RR, and limites inferiores y superiores de intervalos de credibilidad al 95%.
  • Crea un mapa leaflet map que muestre los valores de. RR.
  • Crea un dashboard interactivo con el paquete flexdashboard que incluya un titulo, un mapa leaflet, y una tabla con los resultados.
library(SpatialEpi)
data(scotland)
d <- scotland$data
map <- scotland$spatial.polygon
# Map is in projection OSGB 1936/British National Grid which has EPSG code 27700
proj4string(map) <- "+proj=tmerc +lat_0=49 +lon_0=-2
+k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36
+units=km +no_defs"
# Transform to latitude/longitude
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84 +no_defs"))

library(sp)
rownames(d) <- d$county.names
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE)

map$county <- map$county.names
d$county <- d$county.names
d$Y <- d$cases
d$E <- d$expected

order <- match(map$county, d$county)
map@data <- d[order, c("county", "Y", "E", "AFF")]
head(map@data)

7 Referencias

Moraga, P. Small Area Disease Risk Estimation and Visualization Using R. The R Journal, 10(1):495-506, 2018 https://journal.r-project.org/archive/2018/RJ-2018-036/index.html

Moraga, P. (2019). Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny. Chapman & Hall/CRC Biostatistics Series, http://www.paulamoraga.com/book-geospatial


Creative Commons License
Este trabajo realizado por se publica bajo la licencia Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License.