Datos de área. Riesgo de ocurrencia de cáncer de colon y recto en el País Vasco.

GRECS, Universidad de Girona, España

En este tutorial, inspirado en el trabajo de Moraga (2019), utilizaremos métodos para representar y modelizar mapas de enfermedades con el objetivo de estimar el riesgo de ocurrencia de cáncer de colon y recto en el País Vasco, España, en el periodo 2007-2014. Utilizaremos una base de datos que contiene la población, los casos de cáncer de colon y recto, las variables de estratificación (sexo y grupos de edad) y como variable explicativa el índice de privación MEDEA 2011, para cada uno de los municipios del País Vasco. Los datos de la población, así como del sexo, de los grupos de edad y del índice de privación MEDEA 2011, se obtuvieron a partir del Censo Español de Población y Viviendas 2011 (https://www.ine.es/censos2011_datos/cen11_datos_microdatos.htm, País Vasco es la comunidad 16). Los casos de cáncer de colon y recto por municipio, sexo y grupos de edad se obtuvieron a partir del Registro de Cáncer del País Vasco (http://www.euskadi.eus/informacion/registros-de-cancer/web01-a3regepi/es/). Los datos de la cartografía se obtuvieron también a partir del Censo Español de Población y Viviendas 2011 (http://www.ine.es/censos2011_datos/cen11_datos_resultados_seccen.htm).

Mostraremos como calcular el número de casos observados y esperados, así como la Razón de Incidencia Estandarizada (RIE, SIR en inglés) en cada uno de los municipios del País Vasco. También como estimaremos el riesgo de ocurrencia de cáncer de colon y recto, controlando por factores de riesgo utilizando la aproximación de Laplace anidada integrada (INLA, por sus siglas en inglés) (Rue et al., 2009, 2017; Lindgren and Rue, 2015; Bakka et al., 2018). Por último, mostraremos como construir mapas interactivos del riesgo utilizando la librería leaflet de R (Cheng et al., 2018).

1 Datos y mapa

Indicaremos cuál es nuestro directorio de trabajo utilizando la función setwd().

setwd("~/Documents/Maria/Mis presentaciones/Areal data")

Cargaremos los datos.

# Data
load("~/Documents/Maria/Mis presentaciones/Areal data/PaisVascoCRC.RData")

Exploraremos los datos.

class(PaisVascoCRC)

## [1] "list"

names(PaisVascoCRC)

## [1] "geo"             "data"            "deprivationQ"    "explanatory"    
## [5] "spatial.polygon"

Vemos que PaisVascoCRC es un objeto tipo lista que contiene los siguientes elementos:

• geo : una tabla que contiene los identificadores de cada municipio, y las coordenadas UTM del centroide de cada una de ellos,
• data : una tabla que contiene además de los identificadores de cada municipio, el número de casos, la población y el estrato de sexo y grupo de edad,
• deprivationQ : una tabla de identificadores de cada municipio, junto con el quintil del índice de privación MEDEA 2011,
• explanatory: una tabla de identificadores de cada municipio, junto con las variables explicativas que hemos considerado factores de riesgo (sexo, grupo de edad, índice de privación MEDEA 2011 y los quintiles de este índice),
• spatial.polygon: un objeto tipo SpatialPolygons.

Si disponemos de la cartografía en un fichero vectorial tipo shape, deberemos convertirla en un objeto tipo polígono espacial. Debe tenerse en cuenta que en otras lenguas que no sean inglés existen caracteres (como la ñ, los acentos, los apóstrofos, etc.) que pueden originar problemas al importar el mapa. Por ello sugerimos que, además del nombre del área (COUNTY en nuestro caso), el fichero shape también contenga una variable identificadora numérica.

Nosotros disponemos de la cartografía en un fichero vectorial. Cargaremos la librería maptools (Bivand y Lewin-Koh, 2019) e importaremos el mapa como un objeto tipo polígono espacial.

library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

spatial.polygon=readShapePoly("municipios.shp", IDvar="COUNTY")

Nótese que hemos indicado que el identificador de cada municipio en el mapa está recogido en la variable COUNTY (IDvar="COUNTY").

Incluiremos este objeto en la lista mediante la siguiente instrucción.

PaisVascoCRC=list(geo=PaisVascoCRC$geo,data=PaisVascoCRC$data,deprivationQ=PaisVascoCRC$deprivationQ,explanatory=PaisVascoCRC$explanatory,spatial.polygon=spatial.polygon)

Visualizaremos las primeras filas de PaisVascoCRC$data, PaisVascoCRC$explanatory y PaisVascoCRC$deprivationQ.

PaisVascoCRC$data contiene el número de casos de cáncer de colon y recto (cases ) y la población de cada municipio (population), estratificado por sexo (hombres y mujeres) (gender: m y f ) y grupos de edad (age: <19, 20-64,>=65 ).

head(PaisVascoCRC$data)

##                                      county cases population gender   age
## 1 Abadiño                                       0        660      m   <19
## 2 Abadiño                                      12       2493      m 20-64
## 3 Abadiño                                      14        482      m  >=65
## 4 Abadiño                                       0        712      f   <19
## 5 Abadiño                                       4       2420      f 20-64
## 6 Abadiño                                       6        596      f  >=65

dim(PaisVascoCRC$data)

## [1] 1530    5

`PaisVascoCRC$explanatory contiene las variables de estratificación, sexo y grupos de edad, y como posible variable explicativa, el índice de privación MEDEA 2011 específico para cada municipio. Este índice de privación se presenta como variable contínua (Deprivation) y en quintiles (DeprivationQ).

head(PaisVascoCRC$explanatory)

##                                      county gender   age Deprivation
## 1 Abadiño                                        m   <19  -0.1457854
## 2 Abadiño                                        m 20-64  -0.1457854
## 3 Abadiño                                        m  >=65  -0.1457854
## 4 Abadiño                                        f   <19  -0.1457854
## 5 Abadiño                                        f 20-64  -0.1457854
## 6 Abadiño                                        f  >=65  -0.1457854
##        DeprivationQ
## 1 -0.1688 to 0.0031
## 2 -0.1688 to 0.0031
## 3 -0.1688 to 0.0031
## 4 -0.1688 to 0.0031
## 5 -0.1688 to 0.0031
## 6 -0.1688 to 0.0031

dim(PaisVascoCRC$explanatory)

## [1] 1530    5

PaisVascoCRC$deprivationQ contiene el índice de privación MEDEA 2011 específico para cada municipio. Este índice de privación se presenta en quintiles (DeprivationQ).

head(PaisVascoCRC$deprivationQ)

##                                       county       DeprivationQ
## 1  Abadiño                                    -0.1688 to 0.0031
## 2 Abaltzisketa                                -0.1688 to 0.0031
## 3  Abanto y Ciérvana-Abanto Zierbena           0.0031 to 0.3631
## 4 Aduna                                        0.0031 to 0.3631
## 5 Agurain/Salvatierra                          0.3631 to 1.2000
## 6 Aia                                        .0.5037 to -0.1688

dim(PaisVascoCRC$deprivationQ)

## [1] 272   2

El mapa de los municipios del País Vasco está contenido en el objeto tipo SpatialPolygons denominado PaisVascoCRC$spatial.polygon. Podemos representar el mapa como sigue:

map<-PaisVascoCRC$spatial.polygon
class(map)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

plot(map)

1.1 Transformación de las coordenadas

Para este tutorial utilizaremos la librería leaflet. Esta librería utiliza coordenadas geográficas (latitud/longitud) y no coordenadas UTM (Este/Norte). Transformaremos las coordenadas UTM a coordenadas geográficas, utilizando la función spTransform() de la librería sp (Pebesma and Bivand, 2005; Bivand et al., 2013).

Empezaremos creando un objeto de tipo SpatialPoints al que llamaremos sps utilizando la librería sp, especificando la proyección del País Vasco, es decir, WGS84, ETRS89/UTM zona 30 N. Después, transformaremos las coordenadas UTM de sps a coordenadas geográficas utilizando spTransform() donde en CRS escribiremos CRS("+proj=longlat +datum=WGS84").

#Transform coordinates
library(sp)
sps  <- SpatialPoints(PaisVascoCRC$geo[, c("UTM_x", "UTM_y")], proj4string = CRS("+proj=utm +zone=30N")) #create a SpatialPoints object. We need to specify that the projection of País Vasco is UTM zone 30N
head(sps)

## SpatialPoints:
##         UTM_x   UTM_y
## [1,] 530726.0 4774870
## [2,] 572849.8 4765502
## [3,] 492787.8 4795094
## [4,] 576654.6 4784869
## [5,] 550084.0 4743056
## [6,] 569058.1 4788158
## Coordinate Reference System (CRS) arguments: +proj=utm +zone=30N
## +ellps=WGS84

spst <- spTransform(sps, CRS("+proj=longlat +datum=WGS84")) #transform the UTM coordinates in the data to geographic coordinates (latitude/longitude) because leaflet need geographic coordinates
head(spst)

## SpatialPoints:
##         UTM_x    UTM_y
## [1,] -2.62226 43.12594
## [2,] -2.10567 43.03871
## [3,] -3.08893 43.30864
## [4,] -2.05629 43.21271
## [5,] -2.38714 42.83844
## [6,] -2.14939 43.24306
## Coordinate Reference System (CRS) arguments: +proj=longlat
## +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

head(coordinates(spst))

##         UTM_x    UTM_y
## [1,] -2.62226 43.12594
## [2,] -2.10567 43.03871
## [3,] -3.08893 43.30864
## [4,] -2.05629 43.21271
## [5,] -2.38714 42.83844
## [6,] -2.14939 43.24306

Finalmente, añadiremos las coordenadas geográficas al objeto geo de la lista PaisVascoCRC:

PaisVascoCRC$geo[, c("long", "lat")] <- coordinates(spst)  #add longitude and latitude variables to the data set PaisVascoCRC
head(PaisVascoCRC$geo)

##                                       county    UTM_x   UTM_y     long
## 1  Abadiño                                   530726.0 4774870 -2.62226
## 2 Abaltzisketa                               572849.8 4765502 -2.10567
## 3  Abanto y Ciérvana-Abanto Zierbena         492787.8 4795094 -3.08893
## 4 Aduna                                      576654.6 4784869 -2.05629
## 5 Agurain/Salvatierra                        550084.0 4743056 -2.38714
## 6 Aia                                        569058.1 4788158 -2.14939
##        lat
## 1 43.12594
## 2 43.03871
## 3 43.30864
## 4 43.21271
## 5 42.83844
## 6 43.24306

2 Preparación de los datos

Crearemos un data frame denominado d con las siguientes columnas:

• county: id de cada municipio,
• Y: número de casos observados en cada municipio,
• E: número de casos esperados en cada municipio,
• DeprivationQ: quintil del índice de privación de cada municipio,
• `SIR: SIR de cada municipio.

2.1 Casos observados

PaisVascoCRC$data contiene los casos en cada municipio estratificados por edad y sexo.

head(PaisVascoCRC$data)

##                                      county cases population gender   age
## 1 Abadiño                                       0        660      m   <19
## 2 Abadiño                                      12       2493      m 20-64
## 3 Abadiño                                      14        482      m  >=65
## 4 Abadiño                                       0        712      f   <19
## 5 Abadiño                                       4       2420      f 20-64
## 6 Abadiño                                       6        596      f  >=65

Obtendremos los casos en cada municipio, Y, agregando las filas de PaisVascoCRC$data por municipio y sumando el número observado de casos, mediante la librería dplyr (Wickham et al., 2019).

# Data preparation
# Observed cases
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

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

## # A tibble: 6 x 2
##   county                                           Y
##   <fct>                                        <dbl>
## 1 "Abadiño                                  "     36
## 2 "Abaltzisketa                              "     3
## 3 "Abanto y Ciérvana-Abanto Zierbena        "     73
## 4 "Aduna                                     "     1
## 5 "Agurain/Salvatierra                       "    36
## 6 "Aia                                       "    10

dim(d)

## [1] 255   2

2.1 Casos esperados

Calcularemos el número esperado de casos en cada municipio utilizando la estandarización indirecta. El número de casos esperado en cada municipio es el número total de casos que se esperaría si la población del municipio tuviese la misma estructura de edad y sexo que la población de referencia, en nuestro caso, el País Vasco. Calcularemos los esperados utlizando la función expected() de la librería SpatialEpi (Kim y Wakefield, 2018). Esta función tiene tres argumentos,

• population: un vector del número de habitantes en cada uno de los estratos de cada una de las áreas,
• cases: un vector con el número de casos observados en cada uno de los estratos de cada área,
• n.strata: número de estratos considerado.

Antes de hacer este cálculo, los vectores population y cases tienen que ser ordenados primero, en función del área y, en segundo lugar, los estratos, dentro de cada área, deben ser ordenados del mismo modo. Todos los estratos deben ser incluidos en los vectores, incluso aquéllos con 0 casos.

Utilizaremos para ello la función order() especificando el orden deseado, en nuestro caso, municipio, sexo y edad.

# Expected cases
library(SpatialEpi)
PaisVascoCRC$data <- PaisVascoCRC$data[order(PaisVascoCRC$data$county, PaisVascoCRC$data$gender, PaisVascoCRC$data$age), ]
head(PaisVascoCRC$data)

##                                      county cases population gender   age
## 1 Abadiño                                       0        660      m   <19
## 2 Abadiño                                      12       2493      m 20-64
## 3 Abadiño                                      14        482      m  >=65
## 4 Abadiño                                       0        712      f   <19
## 5 Abadiño                                       4       2420      f 20-64
## 6 Abadiño                                       6        596      f  >=65

A continuación, obtendremos los casos esperados E en cada municipio ejecutando la función expected(), donde population será igual a PaisVascoCRC$data$population y cases será igual a PaisVascoCRC$data$cases. Tenemos 2 sexos y 3 grupos de edad para cada municipio, por lo tanto, el número de estratos será igual a 2 x 3 = 6.

E<-expected(population=PaisVascoCRC$data$population,cases=PaisVascoCRC$data$cases,n.strata=6)
head(E)

## [1] 48.257586  2.314618 64.266157  2.779103 30.755510 13.570882

Ahora añadiremos el vector E a la data frame d que contiene los identificadores de los municipios (county ) y los casos observados (Y), asegurándonos que los elementos en E tienen el mismo orden que en d$county. Para hacer esto, usaremos match() construyendo el vector de posiciones que apareja d$county con unique(PaisVascoCRC$data$county) que se corresponden con los municipios en E. Después, ordenaremos `E utilizando ese vector.

d$E<-E[match(d$county,unique(PaisVascoCRC$data$county))]
head(d)

## # A tibble: 6 x 3
##   county                                           Y     E
##   <fct>                                        <dbl> <dbl>
## 1 "Abadiño                                  "     36 48.3 
## 2 "Abaltzisketa                              "     3  2.31
## 3 "Abanto y Ciérvana-Abanto Zierbena        "     73 64.3 
## 4 "Aduna                                     "     1  2.78
## 5 "Agurain/Salvatierra                       "    36 30.8 
## 6 "Aia                                       "    10 13.6

dim(d)

## [1] 255   3

Tal y como ya hemos dicho anteriormente, hemos calculado el número esperado de casos en cada municipio utilizando la estandarización indirecta mediante la función la función expected() de la librería SpatialEpi (Kim y Wakefield, 2018).

La estandarización indirecta se podría realizar manualmente de la siguiente manera.

Como ya sabemos, disponemos de información de 255 municipios. Consideremos tres grupos de edad, ≤19 años, 20-64 años y ≥65 años y dos grupos de sexo, hombres y mujeres. Nos centraremos en los dos primeros municipios, Abadiño y Abaltzisketa, pero el mismo proceso se repetiría para todos y cada uno de ellos.

Primero necesitamos calcular los riesgos de la enfermedad (riesgo de ocurrencia de cáncer de colon y recto) para cada grupo de edad y sexo en la población de referencia. Utilizaremos, como población de referencia, la población conjunta de las 255 áreas, es decir, la población del País Vasco.

Después multiplicaremos el riesgo de la enfermedad para cada grupo de edad y sexo, por la población del área en este grupo de edad y sexo y después sumaremos para todos los grupos de edad.

Recordemos que el número de casos esperado en cada municipio es el número total de casos que se esperaría si la población del municipio tuviese la misma estructura de edad y sexo que la población de referencia, en nuestro caso, el País Vasco. Como podemos ver, los casos esperados calculados manualmente para los municipios de Abadiño y Abaltzisketa coinciden con los calculados mediante la función expected() de la librería SpatialEpi (48,25759 vs 48,3 y 2,314619 vs 2,31).

Ahora añadiremos otra columna que corresponderá al quintil del índice de privación del municipio.

2.3 Quintil del índice de privación

Añadiremos a la d la variable deprivationQ que corresponde al quintil del índice de privación promedio del municipio (estaba recogido por sección censal antes). Añadiremos esta variable usando la función merge() en la que especificaremos county como la columna clave para la fusión.

# Deprivation Index
d<-merge(d,PaisVascoCRC$deprivationQ,by="county")
head(d)

##                                       county  Y         E
## 1  Abadiño                                   36 48.257586
## 2 Abaltzisketa                                3  2.314618
## 3  Abanto y Ciérvana-Abanto Zierbena         73 64.266157
## 4 Aduna                                       1  2.779103
## 5 Agurain/Salvatierra                        36 30.755510
## 6 Aia                                        10 13.570882
##         DeprivationQ
## 1  -0.1688 to 0.0031
## 2  -0.1688 to 0.0031
## 3   0.0031 to 0.3631
## 4   0.0031 to 0.3631
## 5   0.3631 to 1.2000
## 6 .0.5037 to -0.1688

2.4 Razón de incidencia estandarizada (RIE)

Los mapas de enfermedades (o de mortalidad) más simples representan los casos o muertes observadas en cada área geográfica. Sin embargo, cualquier interpretación de la estructura geográfica de los casos está limitada por la falta de información sobre la distribución espacial de la población que podría estar ‘a riesgo’ y, consecuentemente, ha generado esos casos observados. Por tanto, se prefiere la representación de tasas que permitan incorporar el efecto de la población a riesgo, en vez de representar los casos brutos. Sin embargo, el uso directo de las tasas brutas no permite la comparación entre distintas áreas, ya que las diferencias observadas entre ellas pueden deberse a factores de riesgo que no hayan sido tenidos en cuenta, como la edad. Una medida que tiene en cuenta la estructura de edad es la tasa estandarizada por edad. Hay dos métodos para la estandarización por edad, que se conocen como estandarización directa e indirecta. En la representación de los mapas de enfermedad se prefiere la utilización del método indirecto, que consiste en la comparación de los casos observados de la enfermedad en un área, con los casos esperados si los riesgos para cada grupo de edad fueran los mismos que en cierta población de referencia. La razón observada/esperada se llama razón de incidencia (o mortalidad) estandarizada (RIE, SIR por sus iniciales en inglés o RME, SMR en inglés), que no es más que un estimador del riesgo relativo del área, es decir, del riesgo de enfermad (o muerte) en relación con el grupo considerado de referencia (Barceló et al., 2008).

La RIE para cada área \(i=1,...,n\) (municipio en nuestro caso), se obtiene como el cociente entre el número de casos observados y el número de casos esperados en esa área, \(RIE_i=Y_i⁄E_i\) . El número de casos esperados representa el número total de casos de una determinada enfermedad que uno esperaría si el área tuviera el mismo riesgo de enfermedad que la población de referencia. Mediante estandarización indirecta, el número esperado de casos se calcula como: \(E_i=\sum_{j=1}^{m}r_j^{(s)}n_j^{(i)}\), donde \(r_j^{(s)}\) es el riesgo de enfermedad en el estrato \(j\) en la población de referencia (calculado como la división entre el número de casos de enfermedad en el estrato \(j\) en el conjunto de áreas \(i\) y la población en ese estrato en la población de referencia) y \(n_j^{(i)}\) es la población en el estrato \(j\) en el área \(i\). El estrato \(j\) suele referirse a un grupo de edad-sexo, pero también puede ser otra variable, como, por ejemplo, un grupo de población determinado.

En nuestro caso queremos calcular las RIE (por grupos de edad – tres grupos de edad - y sexo) de cáncer de colon y recto para cada área, es decir, para cada uno de los 255 municipios del País Vasco. Anteriormente, hemos calculado el número de casos esperados mediante estandarización indirecta (\(E_i\); \(i=1,⋯,255)\), tomando como población de referencia, la población conjunta de las 255 áreas, es decir, la población del País Vasco.

Por ejemplo, la RIE para Abadiño se calcularía dividiendo el número de casos observados (36 casos) por el número de casos esperados, \(RIE_1=36/48,25759= 0,7459966\), y para Abaltzisketa, cuyo número de casos observados es de 3, \(RIE_2=3/2,314619=1,29611\).

Si en un área determinada tenemos una RIE mayor que uno significa que en dicha área hay más casos observados que los que se esperaría en la población de referencia, es decir, el riesgo en dicha área es mayor que el esperado. Si la RIE es menor que uno significa que hay menos casos observados que esperados, es decir, el riesgo en dicha área es menor. Si la RIE es igual a uno significa que el número de casos observados en dicha área es el mismo que el número de casos que se esperaría en la población de referencia.

Siguiendo con nuestro ejemplo, el municipio de Abadiño tiene un 25,40% menos riesgo de ocurrencia de la enfermedad que el esperado y el municipio de Abaltzisketa, sin embargo, tiene un 29,61% más de riesgo.

Calcularemos el vector de RIEs como el cociente entre los casos observados y los casos esperados, y lo añadiremos a la data frame d.

# SIRs
d$SIR<-d$Y/d$E

2.5 Datos

La base de datos final es la siguiente:

head(d)

##                                       county  Y         E
## 1  Abadiño                                   36 48.257586
## 2 Abaltzisketa                                3  2.314618
## 3  Abanto y Ciérvana-Abanto Zierbena         73 64.266157
## 4 Aduna                                       1  2.779103
## 5 Agurain/Salvatierra                        36 30.755510
## 6 Aia                                        10 13.570882
##         DeprivationQ       SIR
## 1  -0.1688 to 0.0031 0.7459967
## 2  -0.1688 to 0.0031 1.2961099
## 3   0.0031 to 0.3631 1.1359011
## 4   0.0031 to 0.3631 0.3598284
## 5   0.3631 to 1.2000 1.1705220
## 6 .0.5037 to -0.1688 0.7368718

2.6 Añadir los datos al mapa

Hemos construido la data frame d que contiene el número de casos observados y esperados de cáncer de colon y recto, el quintil del índice de privación del municipio y el SIR para cada uno de los municipios.

Pretendemos añadir la data frame d al mapa, fusionándola con la base de datos del mapa. Utilizaremos, como variable clave, county. El problema es que la variable county contenida en d (d$county) no siempre coincide con la variable county contenida en la base de datos del mapa.

data.frame(d$county,map$COUNTY)[1:5,]

##                                     d.county
## 1  Abadiño                                  
## 2 Abaltzisketa                              
## 3  Abanto y Ciérvana-Abanto Zierbena        
## 4 Aduna                                     
## 5 Agurain/Salvatierra                       
##                          map.COUNTY
## 1                           Abadiño
## 2                      Abaltzisketa
## 3 Abanto y Ciérvana-Abanto Zierbena
## 4                             Aduna
## 5               Agurain/Salvatierra

Así, por ejemplo, nótese como el nombre del primer municipio en d$county, tiene un espacio en blanco el cual no aparece en map$COUNTY. Lo mismo ocurre con el nombre del tercer municipio, así como de muchos otros que no mostramos. Este hecho produciría un error que impediría fusionar las dos data frames.

Para hacerlos coincidir, eliminaremos de d$county todos los posibles espacios, tanto a la derecha como a la izquierda, utilizando la instrucción str_trim contenida en la librería stringr (Wickham, 2019).

# Add data to map

library(stringr)
d$county=str_trim(d$county,side="both")
data.frame(d$county,map$COUNTY)[1:5,]

##                            d.county                        map.COUNTY
## 1                           Abadiño                           Abadiño
## 2                      Abaltzisketa                      Abaltzisketa
## 3 Abanto y Ciérvana-Abanto Zierbena Abanto y Ciérvana-Abanto Zierbena
## 4                             Aduna                             Aduna
## 5               Agurain/Salvatierra               Agurain/Salvatierra

El mapa de los municipios del País Vasco es un objeto del tipo SpatialPolygons al que denominamos map. Utilizando este objeto y la data frame d crearemos un SpatialPolygonsDataFrame que nos permitirá representar mapas de las variables de d. Para hacerlo, llamaremos a las filas de la data frame d con los mismos nombres que d$county. A continuación, fusionaremos el SpatialPolygons map y la data frame d aparejando los polígonos ID del objeto SpatialPolygons con los nombres de las filas de d, es decir los nombres de county (match.ID = TRUE).

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

##                                                              county  Y
## Abadiño                                                     Abadiño 36
## Abaltzisketa                                           Abaltzisketa  3
## Abanto y Ciérvana-Abanto Zierbena Abanto y Ciérvana-Abanto Zierbena 73
## Aduna                                                         Aduna  1
## Agurain/Salvatierra                             Agurain/Salvatierra 36
## Aia                                                             Aia 10
##                                           E       DeprivationQ       SIR
## Abadiño                           48.257586  -0.1688 to 0.0031 0.7459967
## Abaltzisketa                       2.314618  -0.1688 to 0.0031 1.2961099
## Abanto y Ciérvana-Abanto Zierbena 64.266157   0.0031 to 0.3631 1.1359011
## Aduna                              2.779103   0.0031 to 0.3631 0.3598284
## Agurain/Salvatierra               30.755510   0.3631 to 1.2000 1.1705220
## Aia                               13.570882 .0.5037 to -0.1688 0.7368718

3 Representación gráfica de las RIE

Podemos representar en un mapa interactivo de tipo cloropleta los casos observados y esperados de cáncer de colon y recto, las RIE, así como el quintil del índice de privación de cada municipio. Construiremos el mapa utilizando la librería leaflet.

Para crear este mapa, primero, utilizaremos la instrucción leaflet(), añadiéndole como fondo un mapa de mosaicos procedente de OpenStreetMap (https://www.openstreetmap.org) mediante la instrucción addTiles(). Posteriormente, añadiremos los límites municipales del País Vasco con addPolygons() donde especificaremos el color del borde de las áreas (color) y su grosor (weight). Rellenaremos las áreas con los colores proporcionados por la función paleta de color generada por colorNumeric(), e indicaremos en fillOpacity un valor menor que 1 para poder visualizar el mapa de fondo. La función colorNumeric() crea una función de paleta de color que asigna a los valores de los datos diferentes colores de la paleta. Esta función tiene dos parámetros:

• palette: función de color con los valores que deben ser representados, y
• domain: posibles valores que pueden ser representados.

Finalmente, añadiremos una leyenda especificando la función de paleta de color (pal) y los valores utilizados para generar los colores de dicha función de paleta (values). Daremos a opacity el mismo valor que la transparencia de las áreas, y especificaremos un título y una posición para la leyenda.

# Mapping variables
library(leaflet)
l <- leaflet(map) %>% addTiles()
pal <- colorNumeric(palette = "YlOrRd", domain = map$SIR)
l %>% addPolygons(color = "grey", weight = 1, fillColor = ~pal(SIR), fillOpacity = 0.5) %>%
  addLegend(pal = pal, values = ~SIR, opacity = 0.5, title = "SIR", position = "bottomright")

Podemos mejorar el mapa resaltando los municipios cuando el ratón pase por encima, y mostrando la información sobre los casos observados y esperados, las SIR, y el quintil de índice de privación. Esto lo podemos hacer añadiendo los highlightOptions, label y labelOptions a addPolygons(). Señalaremos las áreas utilitzando un borde bastante grueso con (highlightOptions(weight = 4)).

Crearemos las etiquetas utilizando sintaxis HTML. Primero, crearemos el texto que se mostrará utilizando la función sprintf() que devuelve un vector alfanumérico que contiene una combinación de texto y valores de variables con un determinado formato y luego ejecutaremos la instrucción htmltools::HTML() que señala el texto. En labelOptions especificaremos el estilo (style), el tamaño del texto (textsize), y la dirección de las etiquetas (direction). Las posibles opciones para direction son left, right y auto y se refieren a la dirección de la etiqueta en relación al puntero del ratón. Hemos escogido auto ya que queremos que la dirección óptima se escoja dependiendo de la posición del puntero.

labels <- sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
Deprivation quintile: %s <br/> SIR: %s",
                  map$county, map$Y, round(map$E, 2), map$DeprivationQ, round(map$SIR, 2)) %>%
  lapply(htmltools::HTML)

l %>% addPolygons(color = "grey", weight = 1, fillColor = ~pal(SIR), 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 = ~SIR, opacity = 0.5, title = "SIR", position = "bottomright")

Examinaremos el mapa y veremos qué municipios tienen

• SIR = 1 que indica que los casos observados son los mismos que los esperados,
• SIR > 1 que indica que los casos observados son mayores que los esperados (es decir, son municipios con un riesgo mayor),
• SIR < 1 que indica que los casos observados son menores que los esperados (es decir, son municipios con un riesgo menor).

Este mapa muestra algunos municipios con un riesgo claramente mayor de padecer cáncer de colon y recto. Sin embargo, los SIR pueden ser erróneos y poco fiables en municipios con poca población. Esto es así porque la varianza de las SIR (y también de las SMR si de mortalidad se tratara) es proporcional al inverso de los valores esperados, los cuáles dependen de la estructura poblacional. Por lo tanto, poca población implica una gran varianza, lo que distorsiona la visualización del patrón geográfico si lo hubiese. Por contra, las aproximaciones basadas en modelos, o suavización, permiten incorporar covariables y tener en cuenta la dependencia espacial con lo que se mejorarían las estimaciones.

4 Suavización

En esta Sección especificaremos un modelo para los datos observados, y detallaremos los pasos necesarios para ajustar el modelo y obtener estimadores del riesgo utilizando la librería INLA (Rue et al., 2009, 2017, 2019; Lindgren and Rue, 2015; Bakka et al., 2018).

4.1 Modelo

Las RIE, aun cuando han sido ampliamente utilizadas, tienen algunas limitaciones. Dependen en gran medida del tamaño poblacional, puesto que la varianza de las RIE es inversamente proporcional a los valores esperados; así, áreas con poca población presentaran estimadores con gran variabilidad. En este sentido, las RIE extremas y, por tanto, dominantes en el aparente patrón geográfico, son las estimadas con menos precisión en las áreas que tienen pocos casos. Además, la variabilidad de los casos observados suele ser bastante mayor que la esperada, produciéndose lo que se denomina ‘extra variabilidad’. De hecho, cuando se dispone de datos espaciales es importante distinguir dos fuentes de extra variación. En primer lugar, la fuente más importante suele ser la denominada ‘dependencia espacial’ y es consecuencia de la correlación de la unidad espacial con unidades espaciales vecinas, generalmente las contiguas geográficamente. De este modo, las RIE de áreas contiguas, o cercanas, son más similares que las RIE de áreas distantes espacialmente. Parte de esta dependencia no es realmente una dependencia estructural, sino que se debe, principalmente, a la existencia de variables no controladas, es decir, no incluidas en el análisis. Respecto a la segunda fuente, debe asumirse la existencia de extra variación independiente e incorrelacionada espacialmente, denominada ‘heterogeneidad’, debida a variables no observadas sin estructura espacial que podrían influir en el riesgo (Barceló et al., 2008).

Con el fin de solucionar los problemas derivados de la utilización directa de las RIE se han propuesto varias alternativas para ‘suavizarlas’, es decir, reducir la extra variación. En concreto, para estimar los riesgos de enfermedad es preferible utilizar modelos (denominados en inglés ‘disease mapping models’) ya que éstos permiten incorporar variables explicativas y tomar prestada información de las áreas vecinas para mejorar los estimadores locales, suavizando los valores extremos consecuencia de tamaños de muestra pequeños (Barceló et al., 2008).

Un modelo comúnmente utilizado es el modelo de Besag-York-Mollié (BYM) (Besag et al., 1991). En nuestro caso, el modelo se especifica como sigue:

\(Y_i |θ_i \sim~Poisson(E_i×θ_i), i= 1,..., n\)

El número de casos observados (de cáncer de colon y recto) en el municipio \(i\), \(Y_i\), \(i=1,…,n\), se distribuye como una Poisson con media \(E_i×θ_i\), donde \(E_i\) es el número de casos esperados (de cáncer de colon y recto) en el municipio \(i\) y \(θ_i\) es el riesgo relativo para dicho municipio. Municipios con riesgos relativos \(θ_i>1\) y \(θ_i<1\) son municipios con mayor y menor riesgo que el que se esperaría en la población de referencia, respectivamente. Municipios con \(θ_i=1\) tienen el mismo riesgo que el que se esperaría en la población de referencia.

A su vez, el riesgo relativo se modelizaría como,

\(log(θ_i)=β_0+\sum_{k=2}^{5}β_{k}DeprivationQ_{ki}+u_i+v_i\)

• \(β_0\) es el término independiente
• \(β_k\) son los coeficientes para cada uno de los quintiles de privación (hemos escogido el primer quintil como categoría de referencia)
• \(u_i\) es un efecto espacial estructurado, \(u_i|u_{-i}\sim~N(\overline{u_{δ_i}},{\frac{1}{τ_vn_{δ_i}}})\) (un modelo condicionalmente autoregresivo (CAR))
• \(v_i\) es un efecto espacial no estructurado, también denominado heterogeneidad, \(v_i\sim~N(0,\frac{1}{τ_v})\)

En el modelo BYM se incluyen dos efectos aleatorios, uno no estructurado, que recoge la heterogeneidad, y otro estructurado, que recoge la dependencia espacial. El efecto espacial estructurado \(u_i\) es modelizado mediante la distribución condicional autoregresiva (CAR) que suaviza los datos de acuerdo con una determinada estructura de adyacencias dada por una matriz de vecindad que señala que dos áreas son vecinas si tienen algún límite común. La distribución CAR se expresa como: \(u_i|u_{-i}\sim~N(\overline{u_{δ_i}},{\frac{σ_u^2}{n_{δ_i}}})\) donde \(\overline{u_{δ_i}}=n_{δ_i}^{-1}\sum_{j∈δ_i}u_j\),y \(δ_i\) y \(n_{δ_i}\) representan respectivamente el conjunto de vecinos y el número de vecinos del área \(i\). La componente no estructurada \(v_i\) se modela mediante distribuciones independientes e idénticamente distribuidas como una normal de media cero y varianza igual a \(σ_v^2\).

4.2 Matriz de vecindad

Crearemos la matriz de vecindad necesaria para definir el efecto aleatorio espacial utilizando las funciones poly2nb() y nb2INLA() de la librería spdep (Bivand et al., 2018). En primer lugar, utilizaremos poly2nb() para crear una lista de vecinos constituida por los municipios adyacentes. Cada elemento de la lista nb representa un municipio y contiene los indicadores de sus vecinos (adyacentes).

# Neighbourhood matrix
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'))`

## Loading required package: sf

## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

library(INLA)

## Loading required package: parallel

## This is INLA_19.04.16 built 2019-04-16 07:32:56 UTC.
## See www.r-inla.org/contact-us for how to get help.
## To enable PARDISO sparse library; see inla.pardiso()

nb <- poly2nb(map)
head(nb)

## [[1]]
##  [1]  26  45  69  78  80  89 109 135 136 169 201
## 
## [[2]]
## [1]  16 105 170 195 237
## 
## [[3]]
## [1] 107 183 200 212 249
## 
## [[4]]
## [1]  21 232 252
## 
## [[5]]
## [1]  31  41  56 132 203 211
## 
## [[6]]
## [1]   7  42  95 114 197 229 243 248 252

Por ejemplo, nb[[2]] contiene los vecinos del municipio número 2 en el mapa (Abaltzisketa), que en nuestro caso serían los municipios números 16 (Amezketa), 105 (Gaintza), 170 (Mancomunidad de Amezketa y Ordizia), 195 (Orendain) y 237 (Zaldibia).

#To see the municipalities adjacent to municipality 2  (Abaltzisketa)
d$county[c(2,16,105,170,195,237)]

## [1] "Abaltzisketa"                      
## [2] "Amezketa"                          
## [3] "Gaintza"                           
## [4] "Mancomunidad de Amezketa y Ordizia"
## [5] "Orendain"                          
## [6] "Zaldibia"

Después, utilizaremos nb2INLA() para convertir esta lista en un fichero con la matriz de vecindad tal y como la requiere INLA. A continuación, leeremos este fichero utilizando la función inla.read.graph(), y la guardaremos en el objeto g que después usaremos para especificar el modelo con INLA.

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

4.3 Inferencia utilizando INLA

El modelo incluye dos efectos aleatorios, \(u_i\) que captura la variación residual espacial, y \(v_i\) que recoge la variación no estructurada. Necesitamos incluir dos vectores en la data frame que recojan los identificadores de estos efectos aleatorios. Los denominaremos \(re\_u\) para \(u_i\), y \(re\_v\) para \(v_i\). Ambos, \(re\_u\) y \(re\_v\) irán desde 1 hasta \(n\), siendo \(n\) el número de municipios. En nuestro ejemplo, \(n\)=255 y este puede ser obtenido mediante el número de filas de la data frame (nrow(map@data)).

# Inference using INLA
map$re_u <- 1:nrow(map@data)
map$re_v <- 1:nrow(map@data)

Especificaremos la fórmula del modelo incluyendo la variable respuesta en el lado izquierdo, y los efectos fijos y aleatorios en el lado derecho. Los efectos aleatorios se indican utilizando f() con parámetros igual al nombre de la variable indicadora del efecto y el modelo escogido. Para \(u_i\), utilizamos model = "besag" con matriz de vecindad dada por g (en el anexo, explicaremos otro modelo alternativo más utilizado actualmente). Para \(v_i\) escogeremos el modelo model="iid".

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

Nótese que al especificar el efecto aleatorio que captura la variabilidad espacial hemos pedido que nos lo ‘escale’ (scale.model=TRUE).

En el modelo de Besag, York y Mollié, BYM (Besag et al., 1991), el hiperparámetro del efecto aleatorio espacial (es decir, su varianza) se interpretaría como la medida de la dispersión de cada área respecto al riesgo de la enfermedad (de ocurrencia del cáncer de colon y recto en nuestro caso). El problema es que ni la variación estructurada espacialmente ni la variación no estructurada (recogidos por los efectos aleatorios espaciales y por la heterogeneidad, respectivamente) son variables observables, por lo que podrían no ser independientes (problema denominado de la no identificabilidad). Como consecuencia, parte de la dependencia espacial (variación estructurada) podría ser realmente heterogeneidad (variación no estructurada) y viceversa. Además, sus hiperparámetros (las varianzas de los efectos aleatorios) dependen de la estructura espacial del problema (la matriz de vecindad en nuestro caso) y, por tanto, no son interpretables de la misma forma cuando la estructura espacial no es la misma (Riebler et al., 2016). Por ejemplo, el hiperparámetro del efecto aleatorio espacial no sería el mismo si en vez de considerar municipios del País Vasco, trabajásemos con comarcas (agrupación de municipios), aun cuando lo único que ha variado es la estructura espacial. Por este motivo, Sørbye y Rue (2014) indican que los modelos deben escalarse a fin que la interpretación de las varianzas de los efectos aleatorios sea la misma con independencia de la estructura espacial utilizada.

Ajustaremos el modelo mediante la función inla() y utilizando sus priors por defecto (en el anexo explicamos los priors más utilizados actualmente). Especificaremos la formula, la familia, los datos y los casos esperados. Asimismo, indicaremos control.predictor igual a list(compute = TRUE) para calcular las medias a posteriori de los predictores.

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

4.4 Resultados

Podemos ver el objeto de resultados res utilizando summary().

# Results
summary(res)

## 
## Call:
##    c("inla(formula = formula, family = \"poisson\", data = map@data, 
##    ", " E = E, control.predictor = list(compute = TRUE))") 
## Time used:
##     Pre = 1.79, Running = 0.56, Post = 0.095, Total = 2.44 
## Fixed effects:
##                                          mean    sd 0.025quant 0.5quant
## (Intercept)                            -0.366 0.057     -0.483   -0.364
## factor(DeprivationQ).0.5037 to -0.1688  0.170 0.071      0.033    0.168
## factor(DeprivationQ)-0.1688 to 0.0031   0.232 0.072      0.094    0.230
## factor(DeprivationQ)0.0031 to 0.3631    0.292 0.071      0.154    0.291
## factor(DeprivationQ)0.3631 to 1.2000    0.279 0.072      0.143    0.277
##                                        0.975quant   mode kld
## (Intercept)                                -0.258 -0.361   0
## factor(DeprivationQ).0.5037 to -0.1688      0.314  0.165   0
## factor(DeprivationQ)-0.1688 to 0.0031       0.376  0.228   0
## factor(DeprivationQ)0.0031 to 0.3631        0.434  0.289   0
## factor(DeprivationQ)0.3631 to 1.2000        0.426  0.274   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 70.66 53.13      17.81    55.81     212.33 37.48
## Precision for re_v 37.10 15.90      14.22    34.49      75.88 29.36
## 
## Expected number of effective parameters(stdev): 91.86(10.87)
## Number of equivalent replicates : 2.78 
## 
## Marginal log-Likelihood:  -906.89 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

El estimador del término independiente es \(\widehat{β_0}\)=-0.3661 con un intervalo de credibilidad al 95% igual a (-0.4827, -0.2583), y el estimador del coeficiente para el quinto quintil de privación, por ejemplo, es \(\widehat{β_5}\)=0.2793 con un intervalo de credibilidad al 95% igual a (0.1433, 0.4255).

Sin embargo, el modelo es no lineal. De hecho, no es más que una regresión de Poisson con dos efectos aleatorios. Cuando el modelo no es lineal, la interpretación de los estimadores de los coeficientes sin transformar (y de sus intervalos de credibilidad), tiene poco sentido. Conviene transformarlos, tomando sus anti-logaritmos y obteniendo de este modo los riesgos relativos (RR) y los intervalos de credibilidad de los mismos. Utilizaremos el objeto summary.fixed contenido en el objeto de resultados res.

# Relative risks and 95% credible intervals
res$summary.fixed

##                                              mean         sd  0.025quant
## (Intercept)                            -0.3660527 0.05723350 -0.48270879
## factor(DeprivationQ).0.5037 to -0.1688  0.1696076 0.07141389  0.03312716
## factor(DeprivationQ)-0.1688 to 0.0031   0.2316577 0.07176184  0.09411558
## factor(DeprivationQ)0.0031 to 0.3631    0.2915312 0.07114125  0.15445874
## factor(DeprivationQ)0.3631 to 1.2000    0.2793040 0.07190022  0.14328632
##                                          0.5quant 0.975quant       mode
## (Intercept)                            -0.3644651 -0.2582806 -0.3612108
## factor(DeprivationQ).0.5037 to -0.1688  0.1682565  0.3138395  0.1654962
## factor(DeprivationQ)-0.1688 to 0.0031   0.2304535  0.3761569  0.2279973
## factor(DeprivationQ)0.0031 to 0.3631    0.2905583  0.4341684  0.2885787
## factor(DeprivationQ)0.3631 to 1.2000    0.2774877  0.4255228  0.2736823
##                                                 kld
## (Intercept)                            2.530811e-05
## factor(DeprivationQ).0.5037 to -0.1688 2.682822e-06
## factor(DeprivationQ)-0.1688 to 0.0031  5.225678e-07
## factor(DeprivationQ)0.0031 to 0.3631   2.590719e-07
## factor(DeprivationQ)0.3631 to 1.2000   5.766330e-06

Nos interesan los nombres de las filas del objeto (rownames(res$summary.fixed)), la columna mean (la primera columna), es decir la que contiene los estimadores de los coeficientes, de la que tomaremos antilogaritmos (exp(res$summary.fixed[,1])) y las columnas tercera y quinta las cuales contienen los límites inferior y superior de los intervalos de credibilidad al 95%, de los que también tomaremos sus anti-logaritmos (exp(res$summary.fixed[,3]) y exp(res$summary.fixed[,5])).

Coefficients=rownames(res$summary.fixed)
RR=exp(res$summary.fixed[,1])
Lower_RR=exp(res$summary.fixed[,3])
Upper_RR=exp(res$summary.fixed[,5])

data.frame(Coefficients=Coefficients,RR=RR, Lower_RR=Lower_RR,Upper_RR=Upper_RR)

##                             Coefficients        RR  Lower_RR  Upper_RR
## 1                            (Intercept) 0.6934663 0.6171095 0.7723785
## 2 factor(DeprivationQ).0.5037 to -0.1688 1.1848398 1.0336820 1.3686700
## 3  factor(DeprivationQ)-0.1688 to 0.0031 1.2606882 1.0986867 1.4566757
## 4   factor(DeprivationQ)0.0031 to 0.3631 1.3384754 1.1670261 1.5436789
## 5   factor(DeprivationQ)0.3631 to 1.2000 1.3222092 1.1540602 1.5303902

El riesgo relativo correspondiente al segundo quintil de privación, por ejemplo, es \(\widehat{RR_2}\)=1.1848 con un intervalo de credibilidad al 95% igual a (1.034, 1.369). Es decir, los municipios situados en el segundo quintil de privación tienen un 18,5% (intervalo de credibilidad al 95%, ICr95%=3,4% a 36.9%) más de riesgo de ocurrencia de cáncer de colon y recto que los municipios situados en el primer quintil. Los municipios situados en el quinto quintil tienen un 32,2% (ICr95%=15,4% a 53.0%) más de riesgo de ocurrencia de cáncer de colon y recto.

Dibujaremos ahora la distribución a posteriori del coeficiente asociado al quinto quintil del índice de privación, por ejemplo. Lo haremos calculando una suavización de la distribución marginal del coeficiente inla.smarginal() y luego la representaremos ggplot() de la librería ggplot2.

# Posterior distribution of the coefficient of the fifth quintile of the deprivation index

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

4.5 Añadir los resultados al mapa

Los estimadores del riesgo del cáncer de colon y recto, así como la incertidumbre para cada uno de los municipios vienen dados por la media a posteriori y por los intervalos de credibilidad al 95% de \(θ_i\), \(1,…,n\), que están en res$summary.fitted.values. La columna mean es la media a posteriori y 0.025quant y 0.975quant son los percentiles 2.5 y 97.5, respectivamente. Añadiremos estos datos a map para poder representarlos gráficamente. El estimador del riesgo relativo vendrá dado por mean y los límites inferior y superior de los intervalos de credibilidad al 95% de los riesgos, por 0.025quant y 0.975quant.

# Add the results to the map
head(res$summary.fitted.values)

##                           mean         sd 0.025quant  0.5quant 0.975quant
## fitted.Predictor.001 0.8179650 0.09858092  0.6325466 0.8152189   1.019490
## fitted.Predictor.002 0.9830483 0.19666314  0.6485814 0.9652710   1.420611
## fitted.Predictor.003 1.1119869 0.11290076  0.9023886 1.1078583   1.345150
## fitted.Predictor.004 0.9421554 0.19506603  0.6063501 0.9260254   1.371457
## fitted.Predictor.005 1.0071937 0.13875402  0.7620663 0.9977880   1.305542
## fitted.Predictor.006 0.7979653 0.12644992  0.5735121 0.7894453   1.071936
##                           mode
## fitted.Predictor.001 0.8102339
## fitted.Predictor.002 0.9335522
## fitted.Predictor.003 1.0996776
## fitted.Predictor.004 0.8966797
## fitted.Predictor.005 0.9791399
## fitted.Predictor.006 0.7743264

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

5 Representación riesgos suavizados de enfermedad en un mapa

Representaremos los riesgos suavizados de la enfermedad en un mapa interactivo utilizando leaflet. En el mapa, añadiremos etiquetas que aparecerán cuando el ratón pase por encima de los municipios, mostrando información de los casos observados y esperado, SIR, quintil del índice de privación, RR, y límites superior e inferior de los intervalos de credibilidad al 95%.

# Mapping disease risk
pal <- colorNumeric(palette = "YlOrRd", domain = map$RR)
labels <- sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
                  Deprivation quintile: %s <br/> SIR: %s <br/> RR: %s (%s, %s)",
                  map$county, map$Y, round(map$E, 2), map$DeprivationQ, round(map$SIR, 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")

En el mapa anterior, la escala no es la misma que en el mapa de los riesgos sin suavizar. Ello podría conducir a interpretaciones no del todo correctas. De hecho, si observamos el rango de los riesgos sin suavizar y de los riesgos suavizados, comprobaremos que son muy diferentes.

# Range of values of SMRs and RRs
range(map@data$SIR,na.rm=T)

## [1] 0.000000 2.898309

range(map@data$RR)

## [1] 0.3925484 1.2472133

Por lo tanto, deberíamos representarlos utilizando la misma escala.

# Map disease risk with the same scale as map of SIR

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

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")

Cuando representamos los riesgos suavizados en la misma escala que los riesgos sin suavizar, lo primero que se observa es que no se puede identificar un patrón geográfico. Sin embargo, si observamos más detalladamente podemos ver que el riesgo es un poco mayor en torno a la ría de Bilbao, así como en torno a Donostia y Vitoria. Quizás también en el oeste, en el área metropolitana de Bilbao, el riesgo podría ser mayor.

6 Representación probabilidades a posteriori en un mapa

Para ver más claras las conductas observadas anteriormente, así como para evaluar la existencia de aglomeraciones de exceso de casos (más conocidos por su nombre en inglés, clusters), podemos calcular las probabilidades a posteriori (PRP), denominadas también probabilidades de excedencia (Richardson et al., 2004). Las PRP se definen como la probabilidad de que los riesgos relativos sean mayores de 1.

En primer lugar, extraeremos la distribución marginal a posteriori de los efectos aleatorios que capturan la dependencia espacial (\(re\_v\)).

# PRPs
csi=res$marginals.random$re_v[1:nrow(d)]

A continuación, utilizando la instrucción inla.pmarginal, de INLA, obtendremos la probabilidad que RR sea mayor que uno (o equivalentemente que uno menos la probabilidad que sea mayor que cero), a la que denominaremos PRP (Blangiardo y Cameletti, 2015).

PRP=lapply(csi,function(x) {1-inla.pmarginal(0,x)})

El problema es que PRP se extrae como una lista,

head(PRP)

## $index.1
## [1] 0.2879799
## 
## $index.2
## [1] 0.5453543
## 
## $index.3
## [1] 0.6472886
## 
## $index.4
## [1] 0.4050367
## 
## $index.5
## [1] 0.809571
## 
## $index.6
## [1] 0.4511766

Por este motivo lo transformaremos como vector de PRP, añadiéndolo al mapa.

# Add the results to the map
map$PRP = as.numeric(PRP)

Una vez calculadas las PRP podremos localizar aquellas áreas con un riesgo elevado. Siguiendo a Richardson et al. (2004), éstas se corresponderán a aquellas con PRP mayor del 80%. De hecho, estas áreas tienen la máxima sensibilidad (probabilidad de detección superior al 80%) y la máxima especificidad (detección falsa inferior al 10%) (Richardson et al., 2004).

Representaremos las PRP en un mapa interactivo utilizando leaflet. Como antes, añadiremos etiquetas que aparecerán cuando el ratón pase por encima de los municipios, mostrando información de los casos observados y esperados, del quintil del índice de privación, de los RRs, y de las PRPs.

# Mapping PRPs
pal <- colorNumeric(palette = "GnBu", domain = map$PRP)
labels <- sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
                  Deprivation quintile: %s <br/> RR: %s <br/> PRP: %s",
                  map$county, map$Y, round(map$E, 2), map$DeprivationQ, round(map$RR, 2),
                  round(map$PRP, 3)) %>%
  lapply(htmltools::HTML)

leaflet(map) %>% addTiles() %>%
  addPolygons(color = "grey", weight = 1, fillColor = ~pal(PRP), 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 = ~PRP, opacity = 0.5, title = "PRR", position = "bottomright") 

Observamos como las áreas de riesgo más elevado (es decir, PRP mayor que 0.8) se corresponden con las indicadas anteriormente: ría de Bilbao, área metropolitana de Vitoria y Donostia/San Sebastián.

7 Anexo. El modelo BYM2

En el modelo clásico de Besag, York y Mollié, BYM (Besag et al., 1991), la variación estructurada espacialmente no es independiente de la variación no estructurada (problema denominado de la no identificabilidad). Como consecuencia parte de la dependencia espacial (variación estructurada) podría ser realmente heterogeneidad (variación no estructurada) y viceversa. Existen formulaciones alternativas al modelo de BYM, como los modelos de Leroux (Leroux et al., 2000) y de Dean (Dean et al., 2001), en las que se asegura que la variación espacial estructurada sea independiente de la no estructurada. Sin embargo, ninguno de los dos modelos escala la variación espacial. Como consecuencia, los hiperparámetros dependen de la estructura espacial del problema y no pueden ser interpretados correctamente (Riebler et al., 2016).

Por otra parte, hemos realizado las inferencias utilizando una aproximación Bayesiana. En este contexto, la elección de las distribuciones a priori de los hiperparámetros, conocidas como priors, puede tener un impacto considerable en los resultados. En los modelos de Leroux y de Dean se utilizan priors estándar que conducen a un sobreajuste, overfitting (Riebler et al., 2016). La principal consecuencia del overfitting (problema también conocido como multicolinealidad en el contexto de la regresión lineal múltiple), es que los estimadores de las varianzas son mayores que los reales y, por tanto, los intervalos de credibilidad serán mucho más amplios que lo esperado, lo que implica que no se rechazará la hipótesis nula (que los coeficientes sean igual a cero) más veces de lo que debería.

Simpson et al. (2017) propusieron una modificación del modelo BYM que soluciona estos problemas, por cuanto escala la variación estructurada espacialmente y utiliza priors que penalizan la complejidad (denominados PC priors). Estos priors son robustos, en el sentido que no tienen impacto en los resultados y, además, tienen una interpretación epidemiológica.

La especificación del modelo en INLA es la siguiente:

# BYM2
formula <- Y ~ factor(DeprivationQ) + f(re_u, model = "bym2", graph = g, scale.model = TRUE, hyper=list(
  phi=list(prior="pc", param=c(0.5,2/3)),
  prec=list(prior="pc.prec", param=c(0.75,0.05))))

Nótese que hemos especificado PC priors tanto para el hiperparámetro asociado a la dependencia espacial, \(∅\) (phi), como para el hiperparámetro asociado al efecto no estructurado \(τ_b\) (prec) (en realidad para la precisión, precisión=1/varianza, \(σ_b^2=1⁄τ_b\)).

El hiperparámetro \(∅\) se interpreta como el porcentaje de la varianza total que puede ser atribuida a la dependencia espacial (Riebler et al. 2016). En nuestro caso, hemos supuesto que \(Prob(∅<0,5)\)=2/3. Según Riebler et al. (2016) ésta es una elección conservadora, que supone que el efecto aleatorio no estructurado captura más variabilidad que el efecto aleatorio estructurado espacialmente.

En el modelo BYM2, \(τ_b\) es la precisión marginal y puede relacionarse con los riesgos relativos. En concreto, hemos especificado \(Prob(\frac{1}{\sqrt{τ_b}})\)>0,75=0,05). Hemos supuesto un 95% de probabilidad ((1-0,05)*100) de que los riesgos relativos sean menores que 2. El riesgo relativo medio igual a \(RR=\sqrt{-Ln(0,05))⁄0,75}=2\)). De nuevo, siguiendo a Riebler et al. (2016) es una elección conservadora.

Como antes, ajustaremos el modelo mediante la instrucción

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

El resultado puede obtenerse mediante,

summary(resBYM2)

## 
## Call:
##    c("inla(formula = formula, family = \"poisson\", data = map@data, 
##    ", " E = E, control.predictor = list(compute = TRUE))") 
## Time used:
##     Pre = 7.41, Running = 0.797, Post = 0.0986, Total = 8.31 
## Fixed effects:
##                                          mean    sd 0.025quant 0.5quant
## (Intercept)                            -0.383 0.058     -0.501   -0.381
## factor(DeprivationQ).0.5037 to -0.1688  0.180 0.074      0.037    0.179
## factor(DeprivationQ)-0.1688 to 0.0031   0.243 0.075      0.098    0.242
## factor(DeprivationQ)0.0031 to 0.3631    0.300 0.074      0.156    0.299
## factor(DeprivationQ)0.3631 to 1.2000    0.295 0.073      0.156    0.294
##                                        0.975quant   mode kld
## (Intercept)                                -0.272 -0.379   0
## factor(DeprivationQ).0.5037 to -0.1688      0.328  0.177   0
## factor(DeprivationQ)-0.1688 to 0.0031       0.392  0.240   0
## factor(DeprivationQ)0.0031 to 0.3631        0.448  0.298   0
## factor(DeprivationQ)0.3631 to 1.2000        0.443  0.291   0
## 
## Random effects:
##   Name     Model
##     re_u BYM2 model
## 
## Model hyperparameters:
##                      mean    sd 0.025quant 0.5quant 0.975quant   mode
## Precision for re_u 19.267 5.112     11.311    18.56     31.268 17.210
## Phi for re_u        0.422 0.181      0.116     0.41      0.788  0.361
## 
## Expected number of effective parameters(stdev): 98.42(10.00)
## Number of equivalent replicates : 2.59 
## 
## Marginal log-Likelihood:  -665.20 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Y los riesgos relativos de la variable explicativa junto con sus intervalos de credibilidad al 95%,

Coefficients=rownames(resBYM2$summary.fixed)
RR=exp(resBYM2$summary.fixed[,1])
Lower_RR=exp(resBYM2$summary.fixed[,3])
Upper_RR=exp(resBYM2$summary.fixed[,5])
data.frame(Coefficients=Coefficients,RR=RR, Lower_RR=Lower_RR,Upper_RR=Upper_RR)

##                             Coefficients        RR  Lower_RR  Upper_RR
## 1                            (Intercept) 0.6819472 0.6056506 0.7618299
## 2 factor(DeprivationQ).0.5037 to -0.1688 1.1973283 1.0380987 1.3887202
## 3  factor(DeprivationQ)-0.1688 to 0.0031 1.2745362 1.1029536 1.4805249
## 4   factor(DeprivationQ)0.0031 to 0.3631 1.3497745 1.1683854 1.5654311
## 5   factor(DeprivationQ)0.3631 to 1.2000 1.3436004 1.1685607 1.5571018

Los riesgos relativos de la variable explicativa (y sus intervalos de credibilidad al 95%) son muy parecidos en el BYM2 y en el BYM. Sin embargo, con el BYM2 podemos interpretar la variabilidad. En concreto, un 42,2% (intervalo de credibilidad al 95%, ICr95%: 11,6%-78,8%) puede ser atribuida a la dependencia espacial. La precisión marginal es igual a 19,267 (ICr95%: 11.311-31.268), la cual se corresponde a una desviación estándar marginal igual a 0,228 (ICr95%: 0,179-0,297).

8 Referencias

Bakka H, Rue H, Fuglstad A, Riebler A, Bolin D, Krainski ET, Simpson D, Lindgren F. Spatial modelling with R-INLA: A review. Invited extended review, arxiv: 1802.06350, 2018.

Barceló MA, Saez M, Cano-Serral G, Martínez-Beneito MA, Martínez JM, Borrell C, Ocaña-Riola R, Montoya I, Calvo M, López-Abente G, Rodríguez-Sanz M, Toro S, Alcalá JT, Saurina C, Sánchez-Villegas P, Figueiras A. Métodos para la suavización de indicadores de mortalidad: aplicación al análisis de desigualdades en mortalidad en ciudades del Estado español (Proyecto MEDEA). Gaceta Sanitaria 2008; 22(6):596-608.

Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 1991; 43(1):1–20.

Bivand RS, Pebesma E, Gomez-Rubio V. Applied spatial data analysis with R, Second edition. Nueva York: Springer, 2013.

Bivand, Roger S, Wong DWS. Comparing implementations of global and local indicators of spatial association. TEST 2018; 27(3):716-748.

Bivand RS, Lewin-Koh N. maptools: Tools for handling spatial objects. ‘Maptools’ Library. R package version 0.9-5. https://CRAN.R-project.org/package=maptools, 2019.

Blangiardo M, Cameletti M. Spatial and spatio-temporal Bayesian models with R-INLA. Chichester, UK: John Wiley & Sons, 2015.

Cheng J, Karambelkar B, Xie Y. leaflet: Create Interactive Web Maps with the JavaScript ‘Leaflet’ Library. R package version 2.0.2. https://CRAN.R-project.org/package=leaflet, 2018.

Dean C, Ugarte M, Militino A. Detecting interaction between random region and fixed age effects in disease mapping. Biometrics 2001; 57(1): 197–202.

Kim AY, Wakefield J. SpatialEpi: Methods and Data for Spatial Epidemiology. R package version 1.2.3. https://CRAN.R-project.org/package=SpatialEpi, 2018.

Leroux BG, Lei X, Breslow N. Estimation of disease rates in small areas: A new mixed model for spatial dependence. In Statistical Models in Epidemiology, the Environment, and Clinical Trials. Springer, 2000.

Lindgren F, Rue H. Bayesian spatial modelling with R-INLA. Journal of Statistical Software 2015; 63(19):1-25.

Moraga P. Geospatial health data: Modeling and visualization with R-INLA and Shiny. Chapman & Hall/CRC Biostatistical Series, 2019.

Pebesma EJ, Bivand RS. Classes and methods for spatial data in R. R News 5(2) 2005; https://cran.r-project.org/doc/Rnews/.

Riebler A, Sørbye SH , Simpson DP, Rue H. An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research 2016; 25(4):1145:1165 [Disponible en: https://arxiv.org/pdf/1601.01180.pdf, ultimo acceso el 16 de Abril de 2019].

Richardson S, Thomson A, Best N, Elliott P. Interpreting posterior relative risk estimates in disease-mapping studies. Environmental Health Perspectives 2004; 112(9):1016-1025 [Disponible en: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1247195/pdf/ehp0112-001016.pdf, ultimo acceso el 16 de Abril de 2019].

Rue H, Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society Series B, 2009; 71:319-392.

Rue H, Riebler A, Sørbye H, Illian JB, Simpson DP, Lindgren FK. Bayesian computing with INLA: A review. Annual Reviews of Statistics and its Applications 2017; 4(March):395-421.

Rue H, Lindgren F, Simpson D, Martino S, Krainski ET, Bakka H, Riebler A, Fuglstad GA. INLA: Functions with allow to perform full Bayesian analysis of latent Gaussian Models Using Integrated Nested Laplace Approximation. R package version INLA_19.05.19, 2019.

Simpson DP, Rue H, Martins TG, Riebler A, Sørbye SH. Penalising model component complexity: A principled, practical approach to constructing priors (with discussion). Statistical Science 2017; 32(1):1:46 [Disponible en: https://repository.kaust.edu.sa/bitstream/handle/10754/623413/euclid.ss.1491465621.pdf?sequence=1, último acceso el 16 de Abril de 2019].

Sørbye SH, Rue H. Scaling intrinsic Gaussian Markov random field priors in spatial modelling. Spatial Statistics 2014; 8:39–51.

Wickham H, François R, Henry R, Müller K. dplyr: A Grammar of Data Manipulation. R package version 0.8.1. https://CRAN.R-project.org/package=dplyr, 2019.

Wickham H. stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.4.0, https://CRAN.R-project.org/package=stringr, 2019.