En epidemiología ambiental a menudo se necesitan estimar los niveles de factores de riesgo medioambientales, tales como los contaminantes atmosféricos, a los que está expuesto cada uno de los individuos. Habitualmente, se dispone de medidas de los niveles de exposición observados en las estaciones captadoras de contaminación atmosférica. El problema es que las localizaciones de estas estaciones casi nunca coinciden con las localizaciones de los individuos. Éstas casi siempre suelen referirse al lugar de residencia de cada uno de ellos. En todo caso, se necesita predecir los niveles de los contaminantes atmosféricos en esas localizaciones a partir de los datos observados en otras.

Nuestro objetivo en este tutorial consistirá en predecir los niveles de los contaminantes atmosféricos, en concreto los de dióxido de nitrógeno (NO2), a los que están expuestos los pacientes de enfermedad de neurona motora diagnosticados en la Unidad Funcional de la Motoneurona (UFMNA) del Hospital Universitario de Bellvitge, Cataluña, en el período 2011-2019. En total, 613 pacientes cumplieron los criterios diagnósticos de ‘El Escorial-Airlie’ (Brooks et al., 2000).

La información de los niveles de contaminación atmosférica (media mensual de medias diarias) se obtuvieron a partir de las observaciones en 144 estaciones captadoras de la Red Catalana para el Control y Prevención de la Contaminación (XVPCA) localizadas en toda Cataluña y correspondientes al periodo 2011-2018 (http://dtes.gencat.cat/icqa/start.do?lang=en).

Los datos de la cartografía se obtuvieron a partir del Censo Español de Población y Viviendas 2011 (http://www.ine.es/censos2011_datos/cen11_datos_resultados_seccen.htm).

Para cumplir nuestros objetivos utilizamos una cohorte retrospectiva de base poblacional (Povedano et al., 2018) y aplicamos métodos estadísticos Bayesianos bajo la aproximación de Laplace anidada integrada (INLA, por sus siglas en inglés) (Rue et al., 2009; Rue et al., 2017; Lindgren and Rue, 2015; Bakka et al., 2018; Martino y Riebler, 2019) y de las ecuaciones diferenciales parciales estocásticas (INLA, SPDE) (Lindgren et al., 2011; Krainski et al., 2019). Además, representaremos mapas interactivos de la exposición a NO2 utilizando la librería leaflet de R (Cheng et al., 2018).

1 Datos

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

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

Importaremos los datos desde ficheros en formato SPSS.

# Data
library(foreign)
CataloniaALS_data=read.spss("CataloniaALS_data_2011-2019.sav",to.data.frame=T)
Catalonia_airpollutants=read.spss("Catalonia_airpollutants_2011-2018.sav",to.data.frame=T)

Con la opción to.data.frame=TRUE estamos indicando que importe los datos en un formato data frame.

class(CataloniaALS_data)
## [1] "data.frame"
class(Catalonia_airpollutants)
## [1] "data.frame"
names(CataloniaALS_data)
## [1] "id"        "UTM_x"     "UTM_y"     "year_diag"
names(Catalonia_airpollutants)
##  [1] "id_station"    "year"          "month"         "PM10"         
##  [5] "PM2.5"         "NO2"           "NO"            "SO2"          
##  [9] "CO"            "C6H6"          "O3"            "As"           
## [13] "Ni"            "Cd"            "Pb"            "municipality" 
## [17] "id_ABS"        "UTM_x_station" "UTM_y_station"

Vemos que CataloniaALS_data y Catalonia_airpollutants son objetos tipo data frame.

CataloniaALS_data contiene las siguientes variables:

  • id : identificador de cada uno de los individuos
  • UTM_x,UTM_y : las coordenadas UTM del lugar de residencia de cada individuo
  • year_diag : año de diagnóstico de la enfermedad de neurona motora

Catalonia_airpollutants contiene las siguientes variables:

  • id_station : identificador de cada una de las estaciones captadoras de contaminación atmosférica
  • year : una variable que indica el año de observación de los niveles de contaminación
  • month : una variable que indica el mes de observación de los niveles de contaminación
  • PM10 : nivel de contaminación atmosférica observado en la estación por partículas de diámetro menor que 10 \(\mu\)m
  • PM2.5 : nivel de contaminación atmosférica observado en la estación por partículas de diámetro menor que 2,5 \(\mu\)m
  • NO2 : nivel de contaminación atmosférica observado en la estación por dióxido de nitrógeno
  • NO : nivel de contaminación atmosférica observado en la estación por monóxido de nitrógeno
  • SO2 : nivel de contaminación atmosférica observado en la estación por dióxido de azufre
  • CO : nivel de contaminación atmosférica observado en la estación por monóxido de carbono
  • C6H6 : nivel de contaminación atmosférica observado en la estación por benceno
  • O3 : nivel de contaminación atmosférica observado en la estación por ozono
  • As : nivel de contaminación atmosférica observado en la estación por arsénico
  • Ni : nivel de contaminación atmosférica observado en la estación por níquel
  • Cd : nivel de contaminación atmosférica observado en la estación por cadmio
  • Pb : nivel de contaminación atmosférica observado en la estación por plomo
  • municipality : una variable indicando el código del municipio donde está situada la estación captadora
  • id_ABS : identificador de cada una de las ABS donde está situada la estación captadora
  • UTM_x_station, UTM_y_station : las coordenadas UTM del lugar donde se encuentra situada la estación captadora

2 Preparación de los datos

Para estimar la contaminación por NO2 a la que ha estado expuesto el individuo consideraremos la exposición desde el 2011 (primer año para el cuál disponemos de datos observados para los diferentes contaminantes en Cataluña) hasta el año en el cuál el individuo fue diagnosticado de ENM.

Por este motivo, utilizaremos un proceso secuencial para cada uno de los años de diagnóstico, empezando por el año 2011 (Year=2011). Este mismo procedimiento lo realizaremos, posteriormente, para cada uno de los años de diagnóstico.

#Data preparation
Year=2011
ALS_data=CataloniaALS_data[CataloniaALS_data$year_diag==Year,]
airpollutants=Catalonia_airpollutants[Catalonia_airpollutants$year<=Year,]

3 Construcción de la rejilla

Necesitamos triangular el territorio de Cataluña mediante una rejilla con el objetivo de discretizar el campo aleatorio. Para construir la rejilla utilizaremos la función inla.mesh.2d() de la librería INLA (Rue et al., 2009; Rue et al., 201; Rue et al., 2019; Lindgren y Rue, 2015; Bakka et al., 2018), con los siguiente parámetros:

  • loc : coordenadas de la localización de los sujetos que se utilizarán como los nodos iniciales de la rejilla.
  • max.edge : valores que indican la longitud máxima de los bordes de cada uno de los triángulos. Este parámetro es un vector bidimensional. El primer componente indica la longitud máxima para la rejilla interior y el segundo, para la extensión. La unidad de escala de los dos componentes es la misma que la de las coordenadas (en nuestro caso, metros).
  • cutoff : mínima distancia permitida entre nodos. Si la distancia entre dos nodos es menor que la indicada, estos dos nodos se reemplazan por uno sólo. Es muy útil en el caso de nodos que provengan de datos agrupados, ya que evita la construcción de muchos triángulos pequeños alrededor de los mismos.
#Build Mesh
library(INLA)
## Loading required package: Matrix
## Loading required package: sp
## 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()
coo <- cbind(airpollutants$UTM_x_station,airpollutants$UTM_y_station)
mesh <- inla.mesh.2d(loc = coo, max.edge = c(25000,25000), cutoff = 0.3)
plot(mesh)

A continuación, añadiremos la localización de las estaciones captadoras.

#Build Mesh
library(INLA)
coo <- cbind(airpollutants$UTM_x_station,airpollutants$UTM_y_station)
mesh <- inla.mesh.2d(loc = coo, max.edge = c(25000,25000), cutoff = 0.3)
plot(mesh)
points(coo, col = "red")

4 Construcción del modelo SPDE en la rejilla

Utilizaremos la función inla.spde2.pcmatern de la librería INLA,

# Build the SPDE model on the mesh
spde <- inla.spde2.pcmatern(mesh = mesh,constr=TRUE,prior.range=c(0.01,0.01),prior.sigma=c(100,0.1))

En este tutorial, utilizamos priors que penalizan la complejidad (denominados PC priors). Estos priors son robustos, en el sentido que no tienen un impacto en los resultados y, además, tienen una interpretación epidemiológica (Simpson et al., 2017). Hemos especificado PC priors para el rango y la desviación típica marginal del campo aleatorio.

Para el rango especificamos que \(Prob(ρ<ρ_0 )=p_ρ\), donde \(ρ\) denota el rango; \(ρ_o\) es el primer parámetro de la función prior.range (en nuestro caso, 0.01); y \(p_ρ\) es el segundo parámetro de la función (en nuestro caso, también sería 0.01). Por lo que se refiere a la desviación típica marginal, especificamos \(Prob(σ>σ_0 )=p_σ\), donde \(σ\) denota la desviación típica marginal; \(σ_0\) es el primer parámetro de la función prior.sigma (en nuestro caso, 100) y \(p_σ\) es el segundo parámetro (en nuestro caso, 0.01). Si en \(p_ρ\) o en \(p_σ\) especificamos NA en lugar de un valor concreto, entonces, el rango será un parámetro fijo (no aleatorio), igual a \(ρ_0\) y la desviación típica marginal será fija e igual a \(σ_0\), respectivamente.

5 Inferencia utilizando INLA

Hay dos formas de especificar el modelo para el que hacer inferencias a partir de la aproximación SPDE. La primera se utiliza para evaluar la asociación de una variable respuesta con diversas variables explicativas, controlando por variables confusoras. La segunda, la cuál utilizamos en este tutorial, se utiliza para predecir valores de la variable respuesta en localizaciones (en el espacio y/o en el tiempo) en las que ésta no se observa.

En este tutorial, realizamos conjuntamente la estimación de los niveles de NO2 en las localizaciones de las estaciones captadoras y la predicción de los mismos en las localizaciones de los individuos. Esta forma conjunta de estimar y predecir al mismo tiempo evita arrastrar los errores de la estimación a la predicción (Krainsky et al., 2019).

Alternativamente, tal y como hace Moraga (2019), la estimación y la predicción podrían realizarse por separado. El problema de esta aproximación es una consecuencia del traslado del error de estimación a la predicción. Si ese error fuese aleatorio, los estimadores serían insesgados pero ineficientes. Es decir, los intervalos de confianza serían muy amplios. Sin embargo, si el error no fuese aleatorio, por ejemplo, cuando no se especifica bien el modelo y/o la dependencia espacial, los estimadores serían sesgados y sus varianzas estarían mal calculadas.

Especificaremos las matrices de proyección:

# Projectors matrices
# Projector matrix of the estimation of the observed air pollutants
A=inla.spde.make.A(mesh,coo)
# Projector matrix of the prediction of the air pollutants in the ENM locations
poo=cbind(ALS_data$UTM_x, ALS_data$UTM_y)
projgrid <- inla.mesh.projector(mesh, loc=poo)
prdcoo <- projgrid$lattice$loc
Aprd <- projgrid$proj$A

En primer lugar, hemos especificado la matriz de proyección para la estimación del NO2 observado (A). A continuación, utilizando la malla de los contaminantes observados, proyectamos las localizaciones de los individuos con ENM. Con esta última proyección especificamos la matriz de proyección para la predicción de los niveles de NO2 en las localizaciones de los pacientes con ENM (Aprd).

Utilizaremos ahora la función inla.stack() para organizar los datos, los efectos y las matrices de proyección para los contaminantes observados y sus predicciones.

Esta función tiene los siguientes argumentos:

  • data: lista de los vectores de datos
  • A : lista de las matrices de proyección
  • tag: etiqueta para identificar el stack
  • effects : lista de los efecto fijos y aleatorios

El argumento data contiene, en el stack de estimación, la variable respuesta (en nuestro caso, NO2) y en el stack de predicción, un vector de NA de dimensión igual al número de localizaciones de los pacientes con ENM, que es el lo que queremos predecir. A es una lista que contiene dos elementos, la matriz de proyección (A para la estimación y Aprd para la predicción) y un 1 para indicar que los efectos fijos se corresponden uno a uno con la variable respuesta. Los efectos aleatorios, indicados en el subargumento list(i=1:spde$n.spde), contienen los efectos espaciales modelizados por un campo aleatorio Gaussiano (Gaussian Random Field). Por último, los efectos fijos, indicados dentro del subargumento data.frame, contienen un término independiente (Intercept=rep(1,length(airpollutants$NO2))) para el stack de la estimación y (Intercept=rep(1,dim(Aprd)[1])) para el stack de la predicción. En los dos casos son vectores que contienen unos. En el primer caso, la dimensión de este vector es igual a la longitud de la variable respuesta (NO2) y en el segundo caso, al número de filas de la matriz de proyección para la predicción. El stack para la estimación contiene, también como efectos fijos, las covariables. El stack para la predicción no contiene nunca covariables. En nuestro caso, consideramos como covariables todos los otros contaminantes que se han observado en la misma estación captadora y el mes de observación del contaminante (por si hubiese efectos estacionales).

#Stack data for the estimation and prediction
#Stack for estimation
stk.dat=inla.stack(data=list(y=airpollutants$NO2),A=list(A,1),tag='dat',effects=list(list(i=1:spde$n.spde),data.frame(Intercept=rep(1,length(airpollutants$NO2)),PM10=airpollutants$PM10,PM2.5=airpollutants$PM2.5,NO=airpollutants$NO,SO2=airpollutants$SO2,CO=airpollutants$CO,C6H6=airpollutants$C6H6,O3=airpollutants$O3,As=airpollutants$As,Ni=airpollutants$Ni,Cd=airpollutants$Cd,Pb=airpollutants$Pb,month=airpollutants$month)))
#Stack for prediction
stk.prd<-inla.stack(data=list(y=NA),A=list(Aprd,1),effects=list(i=1:spde$n.spde,data.frame(Intercept=rep(1,dim(Aprd)[1]))),tag='prd') 
#Joining of the two stacks
stk.all<-inla.stack(stk.dat,stk.prd)

Nótese que hemos juntado los dos stacks.

A continuación, especificaremos la fórmula del modelo.

#Model formula
formula = y ~ 0 +Intercept+PM10+PM2.5+NO+SO2+CO+C6H6+O3+As+Ni+Cd+Pb+f(month, model= "rw1",scale.model=T,hyper=list(theta=list(prior="pc.prec",param=c(0.5,0.01))))+f(i,model=spde)

En la parte izquierda indicamos la variable respuesta (y) y en la parte derecha especificamos los efectos fijos y aleatorios. Nótese que excluimos el término independiente por defecto (indicado mediante el 0) ya que el stack conjunto ya lo contiene (Intercept). Los efectos aleatorios controlan la dependencia espacial (f(i,model=spde)) y una posible relación no lineal entre el mes en el cual se observa el contaminante y la variable respuesta (f(month,…)). Cuando el modelo es intrínseco, como el camino aleatorio especificado para el mes en nuestro modelo, debe escalarse la variación estructurada temporalmente (scale.model=T). Si ésta no se escalase, los hiperparámetros dependerían de la estructura temporal del problema y no podrían ser interpretados correctamente (Riebler et al., 2016). Nótese además que utilizamos priors que penalizan la complejidad (PC priors) (Simpson et al., 2017). Estos priors son robustos, en el sentido que no tienen impacto en los resultados y, además, tienen una interpretación epidemiológica.

Procederemos a ajustar el modelo mediante la instrucción inla().

#Call inla()
result <-inla(formula,family='gaussian',data=inla.stack.data(stk.all),control.predictor=list(A=inla.stack.A(stk.all),compute=TRUE,link=1),control.results=list(return.marginals.random=FALSE,return.marginals.predictor=FALSE))

Especificamos la fórmula, la familia (Gaussiana ya que la variable respuesta es una variable contínua), los datos y otras opciones. Los datos están en el stack conjunto. Este no tiene una estructura de data frame y por esto debemos utilizar la instrucción inla.stack.data. En control.predictor utilizamos las opciones compute=TRUE para calcular las medias posteriores de las predicciones; A=inla.stack.A(stk.all) para predecir las medias en la proyección conjunta; y link=1 para que cuando se estime el modelo, se utilicen las mismas distribuciones especificadas por el investigador para imputar de forma óptima los datos faltantes si los hubiese. En control.results indicamos que no se guarden las distribuciones marginales para que así, si la malla es compleja, podemos aumentar la velocidad de computación. Debe tenerse en cuenta que si nuestro objetivo fuese el de estimar asociaciones esta última opción no debería incluirse. Por último, extraeremos del objeto result la media y la desviación típica a posteriori de las predicciones.

#Posterior mean and standard deviation of the predictions
id.prd <- inla.stack.index(stk.all,'prd')$data
mean.prd <- result$summary.fitted.values$mean[id.prd]
sd.prd <- result$summary.fitted.values$sd[id.prd]

Los indicadores del stack conjunto que corresponden a las predicciones están etiquetados con tag='pred'. Estas predicciones las hemos obtenido utilizando la instrucción inla.stack.index. Como vemos, a continuación hemos obtenido la media y la desviación típica indicando a qué valores corresponden en el objeto summary.fitted.values del objeto result.

Guardaremos el id de cada individuo, el año de diagnóstico y las coordenadas de su domicilio, así como la media y la desviación típica de las predicciones de NO2 en las mismas coordenadas, en el objeto predictions.

#Getting the predictions
predictions=cbind(ALS_data$id,ALS_data$year_diag,poo,mean.prd,sd.prd)

Una vez obtenidas las predicciones para aquellos individuos diagnosticados en el año 2011, deberá repetirse todo el proceso para el resto de individuos y con sus correspondientes años de diagnóstico. A fin de simplificar el proceso, en el anexo mostramos una función que realiza este proceso iterativo.

6 Representación en un mapa de los niveles predichos de NO2

Por último, representaremos los niveles predichos de NO2 (medias y desviaciones típicas) en el mapa de Cataluña, utilizando la librería leaflet.

Disponemos de la cartografía de Cataluña en un fichero vectorial (Catalonia_ct.shp). Cargaremos la librería raster (Hijmans, 2019) e importaremos el mapa como un objeto tipo polígono espacial.

#Mapping of the predictions of NO2
#From shapefile to raster
library(raster)
## Warning: package 'raster' was built under R version 3.5.2
p<-shapefile("Catalonia_ct.shp")

Visualizaremos el área objeto de estudio.

plot(p)

Representaremos en un mapa tipo raster las predicciones.

ext<-extent(p)
rr<-raster(ext,res=250)

Mediante la instrucción extent indicamos la extensión espacial del raster que representa las coordenadas “X,Y” de las esquinas norte, sur, este y oeste del área que queremos representar, en nuestro caso el mapa vectorial de Cataluña. La resolución del raster, opción res, es la longitud, en metros, de cada uno de los píxeles del raster. En nuestro caso, hemos indicado que represente píxeles de 250x250 metros. A más resolución (menor longitud del píxel), más memoria se consume. Por lo que se debe tener cautela a la hora de escoger dicha resolución.

Ahora ‘rasterizaremos’ las predicciones.

#Rasterization of the predictions
#mean
r_mean=rasterize(p,rr,field=predictions[,5],fun=mean)
#standard deviations
r_sd=rasterize(p,rr,field=predictions[,6],fun=mean)

Rasterizar significa asignar a cada píxel un valor de la variable (indicada como field en la instrucción). Nosotros le hemos indicado que en cada píxel calcule la media de la estimación de la media y de la desviación estándard de las predicciones (fun=mean).

Finalmente, representamos mediante leaflet el raster que hemos construido.

#Plot the raster (mean)
library(leaflet)
bins<-c(0,quantile(r_mean,probs=seq(0,1,0.2)),Inf)
pal<-colorBin(palette="YlOrRd",domain=predictions[,5],bins=bins,na.color="transparent")
leaflet()%>%addTiles()%>%addRasterImage(r_mean,colors=pal,opacity=0.5)%>%addLegend("bottomright",pal=pal,values=values(r_mean),title="Mean of the predictions of NO2")%>%addScaleBar(position=c("bottomleft"))
#Plot the raster (sd)
library(leaflet)
bins<-c(0,quantile(r_sd,probs=seq(0,1,0.2)),Inf)
pal<-colorBin(palette="GnBu",domain=predictions[,6],bins=bins,na.color="transparent")
leaflet()%>%addTiles()%>%addRasterImage(r_sd,colors=pal,opacity=0.5)%>%addLegend("bottomright",pal=pal,values=values(r_sd),title="Standard deviation of the predictions of NO2")%>%addScaleBar(position=c("bottomleft"))

Nótese que en la representación hemos categorizado en quintiles y hemos incluido el 0 como límite inferior y el infinito como límite superior (bins<-c(0,quantile(r_mean,probs=seq(0,1,0.2)),Inf) y pal<-colorBin(…, bins=bins,…)).

7 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.

Brooks BR, Miller RG, Swash M, Munsat TL. El Escorial revisited:revised criteria for the diagnosis of amyotrophic lateral sclerosis. Amyotroph Lateral Scler Other Motor Neuron Disord 2000; 1: 293–299.

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.

Krainski ET, Gómez-Rubio V, Bakka H, Lenzi A, Castro-Camilo D, Simpson D, Lindgren F, Rue H. Advanced spatial modelling with stochastic partial differential equations using R and INLA. Boca Raton, Florida, USA: CRC Press, Taylor & Francis Group, 2019.

Lindgren F, Rue H, Lindström J. An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach (with discussion). Journal of the Royal Statistical Society, Series B. 2011; 73(4):423-498.

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

Martino S, Riebler A. Integrated Nested Laplace Approximations (INLA). 2019 [Disponible en: https://arxiv.org/abs/1907.01248, último acceso el 30 de Julio de 2019].

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

Povedano M, Saez M, Martínez-Matos JA, Barceló MA. Spatial assessment of the association between long-term exposure to environmental factors and the occurrence of amytrophic lateral sclerosis in Catalonia, Spain: A population-based nested case-control study. Neuroepidemiology 2018; 51(1-2):33-49.

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].

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.

Barceló MA, Varga D, Tobías A, Díaz J, Linares C, Saez M. Long term effects of traffic noise on mortality in the city of Barcelona, 2004-2007. Environmental Research 2016; 147:193-206.

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

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.

Diggle PJ, Moraga P, Rowlingson B, Taylor BM. Spatial and spatio-temporal log-Gaussian Cox processes: Extending the geostatistical paradigm. Statistical Science 2013; 28(4):542-563.

Hijmans RJ. raster: Geographic Data Analysis and Modeling. R package version 2.9-5. https://CRAN.R-project.org/package=raster, 2019.

Hoffmann R, Borsboom G, Saez M, Mari Dell’Olmo M, Burström B, Corman D, Costa D, Deboosere P, Domínguez-Berjón MF, Dzúrová D, Gandarillas A, Gotsens M, Kovács K, Mackenbach J, Martikainenn P, Maynou L, Morrison J, Palència L, Pérez G, Pikhart H, Rodríguez-Sanz M, Santana P, Saurina C, Tarkiainen L, Borrell C. Social differences in avoidable mortality between small areas of 15 European cities: an ecological study. International Journal of Health Geographics 2014; 13:8.

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.

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

Wickham H. tidyverse: Easily Install and Lad the ‘Tidyverse’. R package version 1.2.1. https://CRAN.R-project.org/package=tidyverse, 2017.

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.

Anexo. Función para estimar y predecir los niveles de contaminación atmosférica por NO2 de forma iterativa.

#Function to estimate and predict NO2 air pollution levels iteratively
for(j in 2012:2019){
#Data preparation
Year=j
ALS_data=CataloniaALS_data[CataloniaALS_data$year_diag==Year,]
airpollutants=Catalonia_airpollutants[Catalonia_airpollutants$year<=Year,]
#Build Mesh
coo<- cbind(airpollutants$UTM_x_station,airpollutants$UTM_y_station)
mesh<- inla.mesh.2d(loc =coo,max.edge=c(25000,25000),cutoff=0.3)
# Build the SPDE model on the mesh
spde<-inla.spde2.pcmatern(mesh = mesh,constr=TRUE,prior.range=c(0.01,0.01),prior.sigma=c(100,0.1))
# Projectors matrices
# Projector matrix of the estimation of the observed air pollutants
A=inla.spde.make.A(mesh,coo)
# Projector matrix of the prediction of the air pollutants in the ENM locations
poo=cbind(ALS_data$UTM_x, ALS_data$UTM_y)
projgrid<-inla.mesh.projector(mesh, loc=poo)
prdcoo<-projgrid$lattice$loc
Aprd<-projgrid$proj$A
#Stack data for the estimation and prediction
#Stack for estimation
stk.dat=inla.stack(data=list(y=airpollutants$NO2),A=list(A,1), tag='dat',effects=list(list(i=1:spde$n.spde),data.frame(Intercept=rep(1,length(airpollutants$NO2)),PM10=airpollutants$PM10,PM2.5=airpollutants$PM2.5,NO=airpollutants$NO,SO2=airpollutants$SO2,CO=airpollutants$CO,C6H6=airpollutants$C6H6,O3=airpollutants$O3,As=airpollutants$As,Ni=airpollutants$Ni,Cd=airpollutants$Cd,Pb=airpollutants$Pb,month=airpollutants$month)))
#Stack for prediction
stk.prd<-inla.stack(data=list(y=NA), A=list(Aprd,1), effects=list(i=1:spde$n.spde,data.frame(Intercept=rep(1,dim(Aprd)[1])),tag='prd') 
#Joining of the two stacks
stk.all<-inla.stack(stk.dat, stk.prd)
#Model formula
formula=y~0+Intercept+PM10+PM2.5+NO+SO2+CO+C6H6+O3+As+Ni+Cd+Pb+f(month, model= "rw1", scale.model=T,hyper=list(theta=list(prior="pc.prec",param=c(0.5,0.01))))+f(i,model=spde)
#Call inla()
result<-inla(formula, family='gaussian', data=inla.stack.data(stk.all),
control.predictor=list(A=inla.stack.A(stk.all),compute=TRUE, link=1),
control.results=list(return.marginals.random=FALSE,return.marginals.predictor=FALSE))
#Posterior mean and standard deviation of the predictions
id.prd<-inla.stack.index(stk.all,'prd')$data
mean.prd<-result$summary.fitted.values$mean[id.prd]
sd.prd<-result$summary.fitted.values$sd[id.prd]
#Getting the predictions
predictions1=cbind(ALS_data$id,ALS_data$year_diag,poo,mean.prd,sd.prd)
predictions=rbind(predictions,predictions1)}