Procesos puntuales. Prevalencia de enfermos crónicos avanzados en las Áreas Básicas de Salud gestionadas por el Instituto de Asistencia Sanitaria Girona.
GRECS, Universidad de Girona, España
Maria Antònia Barceló y Marc Saez
27 de Septiembre de 2019
Desde hace más de diez años se han ido implementando diversas estrategias coordinadas de atención a personas con enfermedades crónicas en España. La principal razón de este protagonismo de la cronicidad en la agenda global de las políticas en salud, es que supone que ésta es la principal causa de morbilidad y mortalidad (Ministerio de Sanidad, Servicios Sociales e Igualdad, 2012). Por otra parte, la estratificación se ha convertido en el eje más importante sobre el que gravitan esas estrategias. Estratificar consiste en segmentar las personas con enfermedades crónicas en función del uso que realizan de los servicios asistenciales y de su coste. En este sentido, se identifican diferentes niveles de necesidades y se ponen en marcha intervenciones en cada grupo con el fin principal de manejar mejor la presencia de múltiples enfermedades crónicas y de cronicidad avanzada, disminuir los ingresos hospitalarios y las demandas de urgencia, así como realizar un uso más adecuado de los medicamentos.
La Estrategia para el Abordaje de la Cronicidad en el Sistema Nacional de Salud (Blay y Limón, 2017) propone la priorización de pacientes en una pirámide tipo Kaiser (Feachem et al., 2002; Nuño, 2007):
- Nivel 0: Personas sanas con o sin factores de riesgo.
- Nivel 1: Presencia de 1 o 2 enfermedades crónicas. Pacientes de bajo riesgo, con condiciones todavía en estados incipientes.
- Nivel 2: Pacientes crónicos complejos (PCC). Presencia de múltiples enfermedades, 3 o 4 enfermedades crónicas. Son pacientes de alto riesgo, pero de menor complejidad.
- Nivel 3: Enfermos crónicos avanzados (MACA). Pacientes con múltiples enfermedades crónicas y mayor carga de fragilidad. Presencia de 5 o más enfermedades crónicas.
La estratificación se realiza fundamentalmente, bien a partir de la presencia de criterios previamente definidos, o bien a partir de modelos predictivos. Entre ellos, destacamos los basados en el índice PROFUND (Díez-Manglano et al., 2015), aplicado principalmente en Andalucía; en el cuestionario CARS (Avilés et al., 2014), utilizado en Valencia; en el ambulatory care groups, ACG-case mix (Weiner et al., 1991; Starfield et al., 1991; Bolíbar et al., 1996), aplicado en el País Vasco; y en los grupos de morbilidad ajustada, GMA (Monterde et al., 2016) aplicado en Cataluña. El Sistema Nacional de Salud tomó este último como referencia en su proyecto de estratificación de la población (Blay y Limón, 2017) y ha sido validado en Cataluña (Monterde et al., 2016; Estupiñán et al., 2019).
Por otra parte, hay mucha evidencia de que existen desigualdades en la salud (Bouchard et al., 2015) y, como ya se apuntaba en el Informe Ladonde (Lalonde, 1974) o en el informe Black (Black, 1980), fue el informe Acheson (Acheson, 1998) el que concluyó firmemente que las desigualdades en la salud tienen una explicación socioeconómica (Marmot, 2005; Kunst, 2007; Mackenbach et al., 2008; Bouchard et al., 2015). Una proporción no insignificante de ellas está causada por problemas medioambientales. Estos factores están generalmente, aunque no únicamente, vinculados a condiciones socioeconómicas (Bowen, 2001; Evans y Kantrowitz, 2002; Brulle y Pellow, 2006; Laurent et al., 2007; Deguen y Zmirou-Navier, 2010; Hajat et al., 2015).
Sin embargo, la mayoría de modelos predictivos de estratificación no tienen en cuenta las desigualdades socioeconómicas (de hecho, sólo el del País Vasco) y ninguno de ellos considera la desigualdad medioambiental. Si no se tienen en cuenta las desigualdades socioeconómicas y medioambientales, algunas personas podrían estar mal estratificadas. Esto podría tener importantes implicaciones respecto a: la estrategia principal de atención a la cronicidad y a los cuidados de salud, por parte del profesional; la eficiencia de los servicios de atención hospitalaria y primaria, así como la evaluación de modelos de continuidad asistencial y de atención socio-sanitaria, por parte de la persona con enfermedades crónicas; y, por último, la asignación de los recursos disponibles, la mejora de procesos y la reducción del gasto, por parte del sistema sanitario.
Nuestro objetivo en este tutorial consiste en estimar la prevalencia de la enfermedad crónica avanzada (MACA) en las Áreas Básicas de Salud (ABS) que estuvieron gestionadas por el Instituto de Asistencia Sanitaria (IAS), controlando por desigualdades socioeconómicas y medioambientales. El IAS gestionó estas ABS (Anglès, Cassà de la Selva y Breda-Hostalrich), las cuales se corresponden con la comarca de la Selva Interior, Girona, España, hasta el año 2012 inclusive. A partir de 2013 estas áreas han pasado a ser gestionadas por el Instituto Catalán de la Salud (ICS), tras su fusión con el IAS.
Para cumplir nuestro objetivo utilizamos un diseño transversal obtenido a partir de una cohorte retrospectiva de base poblacional (detalles adicionales pueden verse en Barceló et al., 2016) 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 prevalencia de los MACA utilizando la librería leaflet de R (Cheng et al., 2018).
La base de datos está constituida por 21.061 sujetos adultos (mayores de 14 años de edad), asignados a alguna de las ABS gestionadas por el IAS y que tuvieron algún contacto con los servicios de salud entre el 1 de enero de 2005 y el 31 de diciembre de 2012, ambos inclusive.
En este tutorial analizamos conjuntamente dos variables respuesta: una a nivel individual que recoge si el individuo fue o no MACA (es decir, tenía un diagnóstico de cinco o más enfermedades crónicas al final de su seguimiento – entre 2005 y 2012-) y otra a nivel agregado de sección censal, que se corresponde con el número de MACAs en la sección censal en que residía el individuo en cuestión durante su seguimiento.
Como covariables consideramos, en primer lugar, variables individuales como el sexo, la edad que tenía el individuo el 1 de enero de 2005, así como un indicador de niveles elevados de contaminación atmosférica en el domicilio de cada individuo. Este indicador fue construido estimando los niveles de contaminación atmosférica a los que estaba expuesto el individuo en su domicilio, a partir de la información obtenida en las 136 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-2012 (http://dtes.gencat.cat/icqa/start.do?lang=en). En concreto, utilizamos la media anual de las medias diarias de partículas (PM10 y PM2.5), dióxido de nitrógeno (NO2), dióxido de azufre (SO2), benceno (C6H6), ozono (O3), arsénico (As) y níquel (Ni). Categorizamos en quintiles las estimaciones en el domicilio de cada sujeto y construimos un indicador que toma el valor 1 si los niveles de todos los contaminantes estaban situados en los quintiles cuarto o quinto y el valor 0 en otro caso.
Además, utilizamos variables contextuales a nivel de la sección censal donde residía el individuo, tales como el índice de privación basado en el utilizado en el proyecto IneqCities, construido combinando siete indicadores socioeconómicos a nivel de la sección censal en la que el individuo estaba domiciliado (Hoffman et al., 2014) y ajustado con el índice de privación 2011 de la Sociedad Española de Epidemiología, IP2011 (Duque et al., 2019); la tasa de paro; y la estructura de población por sexo y grupos de edad de esa sección censal. Estos indicadores fueron obtenidos a partir del Censo Español de Población y Viviendas del 2011 (https://www.ine.es/censos2011_datos/cen11_datos_microdatos.htm, Cataluña es la comunidad 9).
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).
1 Datos
Indicaremos cuál es nuestro directorio de trabajo utilizando la función setwd().
setwd("~/Documents/Maria/Mis presentaciones/Procesos puntuales/Procesos puntuales_prevalencia crónicas")Importaremos los datos.
# Data
library(foreign)
SelvaInteriorMACA=read.spss("SelvaInteriorMACA.sav",to.data.frame=TRUE)Con la opción to.data.frame=TRUE estamos indicando que importe los datos en un formato data frame.
Exploraremos los datos.
class(SelvaInteriorMACA)## [1] "data.frame"names(SelvaInteriorMACA)## [1] "nhc" "maca" "cases_maca"
## [4] "pcc" "cases_pcc" "UTM_x_ct"
## [7] "UTM_y_ct" "UTM_y" "UTM_x"
## [10] "gender" "age" "age_cat"
## [13] "deprivation" "deprivationQ" "unemployment_rate"
## [16] "unemployment_rateQ" "polluted" "male_under_16"
## [19] "male_16to64" "male_65ormore" "female_under16"
## [22] "female_16to64" "female_65ormore" "population_under16"
## [25] "population_over16" "censustrac" "cpro"
## [28] "cmun" "dist" "secc"
## [31] "municipality" "name_municipality"SelvaInteriorMACA es un objeto tipo data frane que contiene las siguientes variables:
- nhc : identificador de cada uno de los individuos
- maca : una variable que indica si el individuo es un enfermo crónico avanzado, MACA, (1) o no (0)
- cases_maca : una variable que indica el número de MACAs en la sección censal donde reside el individuo
- pcc : una variable que indica si el individuo es un enfermo crónico complejo, PCC, (1) o no (0)
- cases_pcc : una variable que indica el número de PCCs en la sección censal donde reside el individuo
- UTM_x_ct, UTM_y_ct : las coordenadas UTM del centroide de la sección censal donde reside el individuo
- UTM_x, UTM_y : las coordenadas UTM del lugar de residencia de cada individuo
- gender : sexo del individuo (male, hombre; female, mujer)
- age, age_cat : edad del individuo, continua y en grupos de edad de cinco años, respectivamente
- deprivation, deprivationQ : índice de privación IneqCities 2011 ajustado con el índice IP 2011 correspondiente a la sección censal de residencia del individuo, continua y en quintiles respectivamente
- unemployment_rate, unemployment_rate_Q : tasa de paro correspondiente a la sección censal de residencia del individuo, continua y en quintiles respectivamente
- polluted : una variable que indica si el lugar de residencia del individuo está muy contaminado (1) o no (0).
- male_under_16, male_16to64, male_65ormore : variables que indican la población de hombres estratificada en menores de 16 años, de 16 a 64 años y 65 años o más de cada una de las secciones censales de los municipios de las ABS gestionadas por el IAS
- female_under_16, female_16to64, female_65ormore : variables que indican la población de mujeres estratificada en menores de 16 años, de 16 a 64 años y 65 años o más de cada una de las secciones censales de los municipios de las ABS gestionadas por el IAS
- population_under_16, population_over_16 : variables que indican la población total estratificada en menores de 16 años y 16 años o más de cada una de las secciones censales de los municipios de las ABS gestionadas por el IAS
- censustrac : una variable indicando la sección censal de cada uno de los municipios de las ABS gestionadas por el IAS
- cpro, cmun, dist, secc : variables indicadoras de la provincia, el municipio, el distrito y la sección censal del lugar de residencia de cada uno de los individuos
- municipality : una variable indicando el municipio de residencia de cada uno de los individuos
- name_municipality : una variable indicando el nombre del municipio de residencia de cada uno de los individuos
Disponemos también de la cartografía en un fichero vectorial. Cargaremos la librería raster (Hijmans, 2019) e importaremos el mapa como un objeto tipo polígono espacial.
#From shapefile to raster
library(raster)## Warning: package 'raster' was built under R version 3.5.2## Loading required package: spp<-shapefile("ABS i municipis.shp")Visualizaremos el área objeto de estudio.
plot(p)
2 Modelización
En esta Sección especificaremos un modelo para los datos observados, y detallaremos los pasos necesarios para ajustar el modelo y obtener estimadores de la prevalencia de MACA utilizando INLA (Rue et al., 2009; Rue et al., 2017; Rue et al., 2019; Lindgren y Rue, 2015; Bakka et al., 2018) e INLA-SPDE (Lindgren et al., 2011; Krainski et al., 2019).
2.1 Modelo
En los diseños observacionales, como el que utilizamos en este tutorial, es muy probable que ocurra un sesgo de selección. Algunos individuos tienen más probabilidad de ser observados y, por lo tanto, de estar presentes en la muestra que otros. En este sentido, los individuos que contactaron con alguna de las ABS que estaban gestionadas por el IAS están sobrerrepresentados. Si la selección fuese exógena, es decir, si la probabilidad de que se observase a un individuo fuese idéntica para todos los individuos, bastaría con ponderar la muestra de tal manera que se otorgase menor peso a los individuos efectivamente observados. La ponderación utilizada sería igual al inverso del porcentaje de los usuarios sobre la población asignada a cada ABS, estratificando por género y grupo de edad. En epidemiología esta manera de ponderar se denomina estandarización (o ajuste) por edad y sexo (Cricelli et al., 2003), considerando como población estándar la población asignada a cada ABS.
Sin embargo, es muy probable que los individuos contactasen con los servicios de atención primaria (de los que provienen los datos), bien porque tenían un problema de salud o bien porque creen que lo tenían, o podrían existir factores no observados que influyesen en el uso de los servicios de atención primaria que estarían correlacionados con los factores no observables que afectasen a la variable respuesta (Heckman, 1979). En todo caso, la probabilidad de que estos individuos hayan sido observados no es la misma para todos ellos. En este caso, la ponderación (estandarización) por edad y sexo no corregiría el sesgo de selección (Heckman, 1979).
En este caso (denominado selección endógena), se debe utilizar un modelo en dos partes (Heckman, 1976). En la primera parte se estima la probabilidad de que un individuo sea observado. Estas probabilidades se utilizan como ponderaciones en la segunda parte del modelo con el objetivo de corregir la no aleatoriedad (es decir, el sesgo de selección).
Estas partes se pueden estimar de forma separada, como la aproximación utilizada muy recientemente para estimar la prevalencia y la incidencia del Plasmodium falciparum, causante de la malaria (Bhatt et al., 2015; Battle et al., 2019; Weiss et al., 2019). El problema de esta aproximación (de hecho, de todas las aproximaciones en varias etapas) es que el error, inherente a la estimación, cometido en la primera parte se arrastra a la segunda parte. Si ese error es aleatorio, los estimadores (de la prevalencia) son insesgados pero ineficientes. Es decir, los intervalos de confianza serán muy amplios y será difícil de encontrar asociaciones estadísticamente significativas. Sin embargo, si el error no fuese aleatorio, por ejemplo, cuando no se especifique bien la primera parte del modelo, los estimadores serán sesgados y sus varianzas estarán mal calculadas.
Preferimos, por tanto, estimar las dos partes de forma conjunta, como la aproximación propuesta por Saez et al. (2009), en la que se especifica un modelo en dos partes denominado Hurdle.
En la primera parte estimamos la probabilidad de que un individuo sea observado utilizando un modelo lineal generalizado mixto (GLMM) con vínculo binomial. Incluimos como variables asociadas a esta probabilidad aquellas que podrían explicar que un individuo fuese observado (es decir, hubiese utilizado los servicios de atención primaria). En concreto, variables individuales, como el sexo y la edad; variables contextuales, como el índice de privación, la tasa de paro y la población estratificada por sexo y grupos de edad, de la sección censal donde reside el individuo; y un indicador de niveles elevados de contaminación atmosférica en la residencia del individuo.
En la segunda parte utilizamos un GLMM con vínculo Poisson para modelizar el número de MACAs por sección censal. Cabe señalar que utilizamos este vínculo y no una Poisson truncada como en Saez et al. (2009), porque el número de secciones censales donde no existe ningún MACA es muy reducido. Incluimos las mismas variables explicativas que en la primera parte del modelo.
En ambas partes, incluimos también efectos aleatorios que recogen confusores no observados. Uno de los efectos (field) recoge la dependencia espacial (y es común para las dos partes) y el otro recoge la heterogeneidad individual y es específico para cada una de las partes (id.primera, id.segunda).
Especificaremos el modelo para la primera parte como sigue:
Condicional al verdadero riesgo en la localización x_i, la probabilidad de que ocurra un MACA en esa localización, \(P(x_i ), i=1,…,n\), se distribuye como una binomial.
\(Y_i |P(x_i )\sim~Binomial(n.trials_i,P(x_i ))\)
\(n.trials_i\) es la población a riesgo de ser MACA en la localización \(x_i\).
La función vínculo es la siguiente,
\(log(\frac {P(x_i)}{1-P(x_i)})=β_0+η_{1i}+S(x_i)+β_1 Gender_i+\sum_{k=2}^{15}β_{2k}Age\_cat_{ik}+\sum_{k=2}^{5}β_{3k}DeprivationQ_{ik}\\ +\sum_{k=2}^{5}β_{4k}Unemployment\_rateQ_{ik}+β_5 Polluted_i+β_6 log(Male\_16to64_i)+β_7 log(Male\_65ormore_i)\\ +β_8 log(Female\_16to64_i)+β_9 log(Female\_65ormore_i)\)
donde el subíndice i indica la sección censal donde se encuentra la residencia del individuo; y \(β\) son los coeficientes de las variables explicativas (\(e^β\) es el riesgo relativo asociado a cada una de ellas). \(η\) y \(S\) son efectos aleatorios. \(η\) recoge la heterogeneidad individual no estructurada espacialmente. Es decir, recoge aquellos confusores no observables asociados a cada individuo que no varían en el tiempo. \(S\) es un efecto aleatorio estructurado espacialmente que se distribuye normalmente con media cero y una función de covarianza Mátern:
\(Cov(S(x_i),S(x_i'))=\frac{σ^2}{2^{(ν-1)} Γ(ν)}(k||x_i-x_i')^νK_ν(k||x_i-x_i'||)\)
donde \(Κ_ν\) es la función de Bessel modificada de segundo tipo y orden \(ν>0\). \(ν\) es un parámetro de suavización, \(σ^2\) es la varianza y \(κ>0\) se relaciona con el rango (\(ρ=\sqrt{8 ν}⁄κ\)), la distancia a la cuál la correlación espacial es próxima a 0.1.
Por lo que se refiere a la segunda parte, la especificación de la misma será como sigue,
Condicional a verdadero riesgo en la localización \(x_i\), la esperanza matemática de que ocurran casos de MACA en la sección censal donde vive cada individuo, \(θ(x_i ), i=1,…,n\), se distribuye como una Poisson.
\(Y_i |θ(x_i )\sim~Poisson(θ(x_i ))\)
En este caso la función vínculo es la que sigue,
\(log(θ(x_i ))=β_0+η_{2i}+S(x_i)+β_1 Gender_i+\sum_{k=2}^{15}β_{2k}Age\_cat_{ik}+\sum_{k=2}^{5}β_{3k}DeprivationQ_{ik}\\ +\sum_{k=2}^{5}β_{4k}Unemployment\_rateQ_{ik}+β_5 Polluted_i+β_6 log(Male\_16to64_i)\\ +β_7 log(Male\_65ormore_i)+β_8 log(Female\_16to64_i)+β_9 log(Female\_65ormore_i)\)
Siendo la definición de los parámetros, de las variables y de los efectos aleatorios, la misma que antes.
Obsérvese que en las variables que se refieren a la población de cada una de las secciones censales se les ha aplicado la transformación logarítmica por dos motivos. En primer lugar, para estabilizar las varianzas, puesto que el rango de estas variables puede ser muy elevado. En segundo lugar, para facilitar la interpretación de los coeficientes. En este caso, los estimadores de los parámetros correspondientes se interpretan como elasticidades.
2.2 Preparación de los datos
Mediante la instrucción attach activaremos el data frame.
#Data Management
attach(SelvaInteriorMACA)Como estimamos un modelo en dos partes, necesitamos preparar los datos en forma matricial con dos columnas, utilizándose cada una de ellas en cada una de las dos partes del modelo, respectivamente.
Empezamos definiendo las dos variables respuesta que constituirán las dos columnas de la matriz de variables respuesta, MACA y casos de MACA, respectivamente. También, definimos el identificador del efecto aleatorio para cada una de las partes del modelo, id.primera y id.segunda.
#Data prepation
n<- dim(SelvaInteriorMACA)[1]
N<- n*2
Y = matrix(NA, N, 2)
Y[1:n, 1] = maca
Y[1:n + n, 2] = cases_maca
id.primera <- c(nhc, nhc)
id.segunda <- c(nhc,nhc)Ahora haremos lo mismo para el resto de covariables (Gender, Age_cat, DeprivationQ, Unemployment_rateQ, Polluted, Male_16to64, Male_65ormore, Female_16to64, Female_65ormore).
#Data prepation
#Construction of the matrices of the covariables
Gender= numeric(N)
Gender[1:n] = gender
Gender[1:n + n] = gender
Age_cat= numeric(N)
Age_cat[1:n] = age_cat
Age_cat[1:n + n] = age_cat
DeprivationQ = numeric(N)
DeprivationQ[1:n] = deprivationQ
DeprivationQ[1:n + n] = deprivationQ
Unemployment_rateQ = numeric(N)
Unemployment_rateQ[1:n] = unemployment_rateQ
Unemployment_rateQ[1:n + n] = unemployment_rateQ
Polluted = numeric(N)
Polluted[1:n] = polluted
Polluted[1:n + n] = polluted
Male_16to64=numeric(N)
Male_16to64[1:n] = male_16to64
Male_16to64[1:n + n] = male_16to64
Male_65ormore=numeric(N)
Male_65ormore[1:n] = male_65ormore
Male_65ormore[1:n + n] = male_65ormore
Female_16to64=numeric(N)
Female_16to64[1:n] = female_16to64
Female_16to64[1:n + n] = female_16to64
Female_65ormore=numeric(N)
Female_65ormore[1:n] = female_65ormore
Female_65ormore[1:n + n] = female_65ormoreDebemos señalar que en este tutorial suponemos que el efecto de las covariables observadas y de los efectos aleatorios (con excepción del asociado a cada uno de los individuos) es el mismo para cada una de las variables respuesta. Es por este motivo que, en cada una de las matrices, la segunda parte contiene la misma variable que la primera parte.
2.3 Construcción de la rejilla
Con el objetivo de recoger la dependencia espacial, necesitamos triangular la región de estudio 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, 2017, 2019; Lindgren and Rue, 2015 y 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.
# Modelling
#Spatial Dependence
# Mesh
library(INLA)## Loading required package: Matrix## 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(UTM_x_ct, UTM_y_ct)
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 la residencia de cada sujeto.
# Modelling
#Spatial Dependence
# Mesh
library(INLA)
coo <- cbind(UTM_x_ct, UTM_y_ct)
mesh <- inla.mesh.2d(loc = coo, max.edge = c(25000,25000), cutoff = 0.3)
plot(mesh)
points(cbind(UTM_x,UTM_y), col = "red")
2.4 Construcción del modelo SPDE en la rejilla
Ahora, construiremos el 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.
Una vez construidos la rejilla y el modelo SPDE, construiremos la matriz de efectos aleatorios espaciales (Field).
# Build the SPDE model on the mesh
Field=numeric(N)
Field[1:n]=mesh$idx$loc
Field[1:n+n]=mesh$idx$locEn cada una de las partes especificaremos un modelo diferente. En la primera de ellas, utilizamos un modelo binomial y es por este motivo por lo que construiremos la matriz de denominadores (n.trials).
#Construction of the matrix of the denominators
n.trials = rep(NA, N)
n.trials[1:n]= population_over16
n.trials[1:n+n]= 1Sólo utilizamos el denominador, que representa la población “a riesgo” de la que se observan los individuos MACA, para la primera de las partes del modelo. Para la segunda parte, el denominador será la unidad.
En la librería INLA,los datos faltantes, si los hay, se imputan de una forma óptima utilizando el mismo modelo especificado por el investigador. Cuando los modelos no son lineales, como en nuestro caso, debe incluirse la opción link=1 para que cuando se estime el modelo se utilicen las mismas distribuciones no lineales especificadas por el investigador. Si no se incluye esta opción, la opción por defecto es asumir una distribución Gaussiana. Como especificamos modelos no lineales, construimos la matriz correspondiente.
#Construction of the link matrix
link = rep(NA, N)
link[1:n] = 1
link[1:n+n] = 1Finalmente, construimos un data frame que contenga todos los elementos (matrices y vectores) construidos anteriormente.
#Construction of the data frame
data_MACA<- data.frame(id.primera,id.segunda,Field,Y,Gender,Age_cat,DeprivationQ,Unemployment_rateQ,Polluted,Male_16to64,Male_65ormore,Female_16to64,Female_65ormore,link,n.trials)2.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, la cuál utilizamos en este tutorial, se utiliza para evaluar la asociación de una variable respuesta con diversas variables explicativas, controlando por variables confusoras. La segunda, 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.
Especificamos el modelo:
#Model formula
formula = Y~ Gender+factor(Age_cat) + factor(DeprivationQ) + factor(Unemployment_rateQ) + Polluted+ log(Male_16to64)+log(Male_65ormore) + log(Female_16to64) + log(Female_65ormore) +
f(id.primera, model="iid",
hyper= list( prec= list( prior="pc.prec", param=c(0.5, 0.01)))) + f(id.segunda,copy="id.primera")+
f(Field, model=spde) Las dos partes del modelo se conectan a través de los efectos aleatorios que recogen la heterogeneidad individual, utilizando la opción copy="id.primera".
Estimamos ahora el modelo mediante la instrucción inla():
#Call INLA
result= inla(formula, data=data_MACA, family= c("binomial","poisson"), Ntrials=n.trials, control.compute = list(dic=T,waic=T), control.predictor=list(link=link,compute=TRUE), control.fixed = list(expand.factor.strategy='inla', prec = 0.01, prec.intercept = 0.01))## Warning in log(Male_65ormore): NaNs produced
## Warning in log(Male_65ormore): NaNs producedsummary(result)##
## Call:
## c("inla(formula = formula, family = c(\"binomial\", \"poisson\"),
## data = data_MACA, ", " Ntrials = n.trials, control.compute =
## list(dic = T, waic = T), ", " control.predictor = list(link =
## link, compute = TRUE), control.fixed = list(expand.factor.strategy
## = \"inla\", ", " prec = 0.01, prec.intercept = 0.01))")
## Time used:
## Pre = 4.82, Running = 3638, Post = 3.4, Total = 3646
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -1.265 0.078 -1.417 -1.265 -1.113
## Gender 0.059 0.006 0.048 0.059 0.070
## factor(Age_cat)2 0.383 0.024 0.335 0.383 0.430
## factor(Age_cat)3 0.392 0.022 0.349 0.392 0.435
## factor(Age_cat)4 0.399 0.020 0.359 0.399 0.439
## factor(Age_cat)5 0.405 0.020 0.366 0.405 0.444
## factor(Age_cat)6 0.424 0.020 0.386 0.424 0.463
## factor(Age_cat)7 0.417 0.020 0.378 0.417 0.455
## factor(Age_cat)8 0.398 0.020 0.359 0.398 0.436
## factor(Age_cat)9 0.396 0.020 0.357 0.396 0.435
## factor(Age_cat)10 0.393 0.020 0.354 0.393 0.432
## factor(Age_cat)11 0.382 0.020 0.343 0.382 0.422
## factor(Age_cat)12 0.386 0.020 0.346 0.386 0.427
## factor(Age_cat)13 0.375 0.021 0.334 0.375 0.417
## factor(Age_cat)14 0.349 0.023 0.303 0.349 0.395
## factor(Age_cat)15 0.353 0.028 0.298 0.353 0.408
## factor(DeprivationQ)2 -1.078 0.057 -1.189 -1.078 -0.966
## factor(DeprivationQ)3 0.679 0.057 0.568 0.679 0.790
## factor(DeprivationQ)4 0.467 0.058 0.355 0.467 0.581
## factor(DeprivationQ)5 0.312 0.058 0.199 0.312 0.426
## factor(Unemployment_rateQ)2 1.121 0.035 1.053 1.120 1.189
## factor(Unemployment_rateQ)3 1.127 0.035 1.059 1.127 1.196
## factor(Unemployment_rateQ)4 1.223 0.037 1.152 1.223 1.295
## factor(Unemployment_rateQ)5 1.735 0.037 1.663 1.735 1.807
## Polluted 0.428 0.009 0.411 0.428 0.446
## log(Male_16to64) -1.823 0.038 -1.898 -1.823 -1.748
## log(Male_65ormore) -0.264 0.021 -0.305 -0.264 -0.223
## log(Female_16to64) 0.894 0.034 0.828 0.894 0.960
## log(Female_65ormore) 0.251 0.014 0.224 0.251 0.278
## mode kld
## (Intercept) -1.264 0
## Gender 0.059 0
## factor(Age_cat)2 0.383 0
## factor(Age_cat)3 0.392 0
## factor(Age_cat)4 0.399 0
## factor(Age_cat)5 0.405 0
## factor(Age_cat)6 0.424 0
## factor(Age_cat)7 0.417 0
## factor(Age_cat)8 0.398 0
## factor(Age_cat)9 0.396 0
## factor(Age_cat)10 0.393 0
## factor(Age_cat)11 0.382 0
## factor(Age_cat)12 0.386 0
## factor(Age_cat)13 0.375 0
## factor(Age_cat)14 0.349 0
## factor(Age_cat)15 0.353 0
## factor(DeprivationQ)2 -1.079 0
## factor(DeprivationQ)3 0.678 0
## factor(DeprivationQ)4 0.467 0
## factor(DeprivationQ)5 0.311 0
## factor(Unemployment_rateQ)2 1.120 0
## factor(Unemployment_rateQ)3 1.127 0
## factor(Unemployment_rateQ)4 1.223 0
## factor(Unemployment_rateQ)5 1.735 0
## Polluted 0.428 0
## log(Male_16to64) -1.823 0
## log(Male_65ormore) -0.264 0
## log(Female_16to64) 0.894 0
## log(Female_65ormore) 0.251 0
##
## Random effects:
## Name Model
## id.primera IID model
## Field SPDE2 model
## id.segunda Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for id.primera 68.97 1.42 66.70 68.79 72.17
## Range for Field 111.82 16.08 81.69 111.50 144.60
## Stdev for Field 40.04 14.94 13.55 39.79 66.78
## mode
## Precision for id.primera 68.13
## Range for Field 111.59
## Stdev for Field 36.27
##
## Expected number of effective parameters(stdev): 6740.47(25.03)
## Number of equivalent replicates : 6.25
##
## Deviance Information Criterion (DIC) ...............: 2980608.81
## Deviance Information Criterion (DIC, saturated) ....: 195978.91
## Effective number of parameters .....................: 6693.15
##
## Watanabe-Akaike information criterion (WAIC) ...: 3161230.98
## Effective number of parameters .................: 162331.10
##
## Marginal log-Likelihood: -1491405.00
## Posterior marginals for the linear predictor and
## the fitted values are computedEn esta instrucción, hemos indicado que se estimen de forma simultánea dos modelos: una regresión logística y una regresión de Poisson (family= c("binomial","poisson")); que utilice la base de datos data_MACA; que calcule el DIC y el WAIC como medidas de ajuste del modelo; que utilice la aproximación gaussiana en la distribución a posteriori de los parámetros; y que el modelo impute los missings utilizando su propia especificación (expand.factor.strategy='inla'). Nótese que hemos añadido la opción link comentada anteriormente.
El resultado del modelo se puede visualizar utilizando la siguiente instrucción:
summary(result)##
## Call:
## c("inla(formula = formula, family = c(\"binomial\", \"poisson\"),
## data = data_MACA, ", " Ntrials = n.trials, control.compute =
## list(dic = T, waic = T), ", " control.predictor = list(link =
## link, compute = TRUE), control.fixed = list(expand.factor.strategy
## = \"inla\", ", " prec = 0.01, prec.intercept = 0.01))")
## Time used:
## Pre = 4.82, Running = 3638, Post = 3.4, Total = 3646
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant
## (Intercept) -1.265 0.078 -1.417 -1.265 -1.113
## Gender 0.059 0.006 0.048 0.059 0.070
## factor(Age_cat)2 0.383 0.024 0.335 0.383 0.430
## factor(Age_cat)3 0.392 0.022 0.349 0.392 0.435
## factor(Age_cat)4 0.399 0.020 0.359 0.399 0.439
## factor(Age_cat)5 0.405 0.020 0.366 0.405 0.444
## factor(Age_cat)6 0.424 0.020 0.386 0.424 0.463
## factor(Age_cat)7 0.417 0.020 0.378 0.417 0.455
## factor(Age_cat)8 0.398 0.020 0.359 0.398 0.436
## factor(Age_cat)9 0.396 0.020 0.357 0.396 0.435
## factor(Age_cat)10 0.393 0.020 0.354 0.393 0.432
## factor(Age_cat)11 0.382 0.020 0.343 0.382 0.422
## factor(Age_cat)12 0.386 0.020 0.346 0.386 0.427
## factor(Age_cat)13 0.375 0.021 0.334 0.375 0.417
## factor(Age_cat)14 0.349 0.023 0.303 0.349 0.395
## factor(Age_cat)15 0.353 0.028 0.298 0.353 0.408
## factor(DeprivationQ)2 -1.078 0.057 -1.189 -1.078 -0.966
## factor(DeprivationQ)3 0.679 0.057 0.568 0.679 0.790
## factor(DeprivationQ)4 0.467 0.058 0.355 0.467 0.581
## factor(DeprivationQ)5 0.312 0.058 0.199 0.312 0.426
## factor(Unemployment_rateQ)2 1.121 0.035 1.053 1.120 1.189
## factor(Unemployment_rateQ)3 1.127 0.035 1.059 1.127 1.196
## factor(Unemployment_rateQ)4 1.223 0.037 1.152 1.223 1.295
## factor(Unemployment_rateQ)5 1.735 0.037 1.663 1.735 1.807
## Polluted 0.428 0.009 0.411 0.428 0.446
## log(Male_16to64) -1.823 0.038 -1.898 -1.823 -1.748
## log(Male_65ormore) -0.264 0.021 -0.305 -0.264 -0.223
## log(Female_16to64) 0.894 0.034 0.828 0.894 0.960
## log(Female_65ormore) 0.251 0.014 0.224 0.251 0.278
## mode kld
## (Intercept) -1.264 0
## Gender 0.059 0
## factor(Age_cat)2 0.383 0
## factor(Age_cat)3 0.392 0
## factor(Age_cat)4 0.399 0
## factor(Age_cat)5 0.405 0
## factor(Age_cat)6 0.424 0
## factor(Age_cat)7 0.417 0
## factor(Age_cat)8 0.398 0
## factor(Age_cat)9 0.396 0
## factor(Age_cat)10 0.393 0
## factor(Age_cat)11 0.382 0
## factor(Age_cat)12 0.386 0
## factor(Age_cat)13 0.375 0
## factor(Age_cat)14 0.349 0
## factor(Age_cat)15 0.353 0
## factor(DeprivationQ)2 -1.079 0
## factor(DeprivationQ)3 0.678 0
## factor(DeprivationQ)4 0.467 0
## factor(DeprivationQ)5 0.311 0
## factor(Unemployment_rateQ)2 1.120 0
## factor(Unemployment_rateQ)3 1.127 0
## factor(Unemployment_rateQ)4 1.223 0
## factor(Unemployment_rateQ)5 1.735 0
## Polluted 0.428 0
## log(Male_16to64) -1.823 0
## log(Male_65ormore) -0.264 0
## log(Female_16to64) 0.894 0
## log(Female_65ormore) 0.251 0
##
## Random effects:
## Name Model
## id.primera IID model
## Field SPDE2 model
## id.segunda Copy
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant
## Precision for id.primera 68.97 1.42 66.70 68.79 72.17
## Range for Field 111.82 16.08 81.69 111.50 144.60
## Stdev for Field 40.04 14.94 13.55 39.79 66.78
## mode
## Precision for id.primera 68.13
## Range for Field 111.59
## Stdev for Field 36.27
##
## Expected number of effective parameters(stdev): 6740.47(25.03)
## Number of equivalent replicates : 6.25
##
## Deviance Information Criterion (DIC) ...............: 2980608.81
## Deviance Information Criterion (DIC, saturated) ....: 195978.91
## Effective number of parameters .....................: 6693.15
##
## Watanabe-Akaike information criterion (WAIC) ...: 3161230.98
## Effective number of parameters .................: 162331.10
##
## Marginal log-Likelihood: -1491405.00
## Posterior marginals for the linear predictor and
## the fitted values are computedNosotros proponemos visualizarlo de una forma más simplificada:
#Simplified representation of results
matriz=matrix(NA, length(names(result$marginals.fixed)),1)
for(i in 1: length(names(result$marginals.fixed))){
matriz[i,]=ifelse(result$summary.fixed[i,1]>0,1-inla.pmarginal(0,result$marginals.fixed[[i]]), inla.pmarginal(0,result$marginals.fixed[[i]]))}
prob.coef.dif.zero=matriz[,1]
data.frame(RR=exp(result$summary.fixed[1:25,1]), exp(result$summary.fixed[1:25,c(3,5)]), prob.coef.dif.zero[1:25])## RR X0.025quant X0.975quant
## (Intercept) 0.2823358 0.242376 0.3285582
## Gender 1.0610076 1.049198 1.0729418
## factor(Age_cat)2 1.4661353 1.398273 1.5372622
## factor(Age_cat)3 1.4803569 1.418091 1.5454030
## factor(Age_cat)4 1.4904009 1.431979 1.5513009
## factor(Age_cat)5 1.4996475 1.442024 1.5596808
## factor(Age_cat)6 1.5285776 1.470356 1.5892192
## factor(Age_cat)7 1.5169871 1.459527 1.5768273
## factor(Age_cat)8 1.4886992 1.432506 1.5472131
## factor(Age_cat)9 1.4860849 1.429632 1.5448797
## factor(Age_cat)10 1.4816972 1.425097 1.5406543
## factor(Age_cat)11 1.4658188 1.408483 1.5255822
## factor(Age_cat)12 1.4715557 1.413617 1.5319595
## factor(Age_cat)13 1.4556254 1.396659 1.5171529
## factor(Age_cat)14 1.4179077 1.354270 1.4845323
## factor(Age_cat)15 1.4235228 1.347016 1.5041913
## factor(DeprivationQ)2 0.3401968 0.304519 0.3804594
## factor(DeprivationQ)3 1.9714418 1.765261 2.2043153
## factor(DeprivationQ)4 1.5959019 1.426454 1.7875359
## factor(DeprivationQ)5 1.3660893 1.219916 1.5314915
## factor(Unemployment_rateQ)2 3.0666101 2.865446 3.2834770
## factor(Unemployment_rateQ)3 3.0873737 2.883635 3.3070853
## factor(Unemployment_rateQ)4 3.3970737 3.162960 3.6501309
## factor(Unemployment_rateQ)5 5.6680357 5.275108 6.0928740
## Polluted 1.5349178 1.508657 1.5616031
## prob.coef.dif.zero.1.25.
## (Intercept) 1
## Gender 1
## factor(Age_cat)2 1
## factor(Age_cat)3 1
## factor(Age_cat)4 1
## factor(Age_cat)5 1
## factor(Age_cat)6 1
## factor(Age_cat)7 1
## factor(Age_cat)8 1
## factor(Age_cat)9 1
## factor(Age_cat)10 1
## factor(Age_cat)11 1
## factor(Age_cat)12 1
## factor(Age_cat)13 1
## factor(Age_cat)14 1
## factor(Age_cat)15 1
## factor(DeprivationQ)2 1
## factor(DeprivationQ)3 1
## factor(DeprivationQ)4 1
## factor(DeprivationQ)5 1
## factor(Unemployment_rateQ)2 1
## factor(Unemployment_rateQ)3 1
## factor(Unemployment_rateQ)4 1
## factor(Unemployment_rateQ)5 1
## Polluted 1data.frame(elasticity=result$summary.fixed[26:29,1], result$summary.fixed[26:29,c(3,5)], prob.coef.dif.zero[26:29])## elasticity X0.025quant X0.975quant
## log(Male_16to64) -1.8229973 -1.8981476 -1.7479055
## log(Male_65ormore) -0.2636486 -0.3045878 -0.2226603
## log(Female_16to64) 0.8940343 0.8275386 0.9604843
## log(Female_65ormore) 0.2508534 0.2240730 0.2776283
## prob.coef.dif.zero.26.29.
## log(Male_16to64) 1
## log(Male_65ormore) 1
## log(Female_16to64) 1
## log(Female_65ormore) 1Como vemos, con la instrucción mostramos los riesgos relativos (RR ) y las elasticidades (elasticity ), sus intervalos de credibilidad al 95% (X0.025quant;X0.975quant ), así como la probabilidad de que el coeficiente estimado sea diferente de cero (prob.coef.dif.zero ). Cuando esta probabilidad sea mayor que 0,95, el riesgo relativo y la elasticidad serán significativos al 95% y cuando sean mayores a 0,90, al 90%.
Ser mujer implica un riesgo 6,10% mayor de ser enfermo crónico avanzado (MACA ) que ser hombre. Este riesgo aumenta, respecto a los individuos de 20 a 24 años de edad, hasta los 45 años de edad en términos relativos. Es a partir de aquí que va disminuyendo, también en términos relativos, excepto para los individuos de 70 a 74 años de edad y para aquéllos con edades superiores a 85 años. Respecto a la interpretación del índice de privación, en los quintiles tercero y superiores hay un mayor riesgo de ser MACA respecto al primer quintil, el grupo de población más favorecido. Nótese como en el segundo quintil, el riesgo es incluso menor que en el primero. Por lo que ser refiere a la tasa de desempleo, existe un riesgo creciente respecto al primero, suave hasta el quinto quintil en el que este riesgo se dispara. Cuando el individuo reside en una zona altamente contaminada, el riesgo de ser MACA es un 53,49% mayor que en otra zona. Por último, las elasticidades son negativas para todos los grupos de edad en el caso de los hombres y positivas para todos los grupos de edad en el caso de las mujeres. Las elasticidades se interpretarían, por ejemplo, como sigue: al aumentar en un 1%, por encima de la media, la población de mujeres de 16 a 64 años de edad en la sección censal donde reside el individuo, el número de MACAs en esa sección censal aumenta un 0,89%.
Por tanto, podríamos concluir que el riesgo de ser MACA es mayor en mujeres que viven en zonas deprimidas y altamente contaminadas.
3 Representación en un mapa de la prevalencia de MACA
Por último, representaremos la prevalencia de MACA en el mapa de las ABS gestionadas por el IAS, la Selva Interior, Girona, utilizando la librería leaflet. Para ello, extraeremos la estimación de la media de la prevalencia, así como su intervalo de credibilidad al 95%, del objeto resultado del modelo (result).
#Mapping MACA prevalence
prev_MACA_mean= result$summary.fitted.values[(n+1):N,1]
prev_MACA_ll= result$summary.fitted.values[(n+1):N,3]
prev_MACA_ul= result$summary.fitted.values[(n+1):N,5]Representaremos en un mapa tipo raster la prevalencia. Para ello utilizaremos la librería raster.
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 las ABS gestionadas por el IAS. 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’ la prevalencia de MACA.
#Rasterization of prevalence of MACA
r_MACA_mean=rasterize(p,rr,field=prev_MACA_mean*1000,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 prevalencia (fun=mean).
Nótese que representamos la prevalencia por 1000 habitantes adultos (`prev_MACA_mean*1000).
Finalmente, representamos mediante leaflet el raster que hemos construido.
#Plot the raster
library(leaflet)
bins<-quantile(r_MACA_mean,probs=seq(0,1,0.2))
pal<-colorBin(palette="YlOrRd",domain=prev_MACA_mean*1000,bins=bins,na.color="transparent")
leaflet()%>%addTiles()%>%
addRasterImage(r_MACA_mean,colors=pal,opacity=0.5)%>%
addLegend("topright",pal=pal,values=values(r_MACA_mean),title="Prevalence of MACA")%>%
addScaleBar(position=c("bottomleft"))## Warning in colors(.): Some values were outside the color scale and will be
## treated as NANótese que en la representación hemos categorizado las prevalencias en quintiles (bins<-quantile(r_MACA_mean,probs=seq(0,1,0.2)) y pal<-colorBin(…, bins=bins,…)).
4 Referencias
Acheson D. Independent Inquiry into Inequalities in Health Report. London, England: The Stationary Office, 1998 [Disponible en: https://www.gov.uk/government/uploads/system/uploads/attachment_data/file/265503/ih.pdf, último acceso el 30 de Julio de 2019].
Avilés MJ, Cuevas MD, Zafra E (coords). Estrategia para la atención a pacientes crónicos en la Comunitat Valenciana. Valencia: Generalitat, Conselleria de Sanitat, 2014 [Disponible en: http://publicaciones.san.gva.es/publicaciones/documentos/V.2792-2014.pdf, último acceso el 3 de Agosto de 2019].
Bhatt S, Weiss DJ, Cameron E, Bisanzio D, Mappin B, Dalrympie U, Battle K, Moyes CL, Henry A, Eckhoff PA, Wenger EA, Briët O, Penny MA, Smith TA, Bennett A, Yukich J, Eisele TP, Griffin JT, Fergus CA, Lynch M, Lindgren F, Cohen JM, Murray CLJ, Smith DL, Hay SI, Cibulskis RE, Gething PW. The effect of malaria control on Plasmodium falciparum in Africa between 2000 and 2015. Nature 2015; 526(7572):207-211.
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, Coll-Negre M, Coll-de-Tuero G, Saez M. Effects of the financial crisis on psychotropic drug consumption in a cohort from a semi-urban region in Catalonia, Spain. PLoS One 2016; 11(2):e0148594.
Battle KE, Lucas TCD, Nguyen M, Howes RE, Nandi AK, Twohig KA, Pfeffer DA, Cameron E, Rao PC, Casey D, Gibson HS, Rozier JA, Dalrymple U, Keddie SH, Collins EL, Harris JR, Guerra CA, Thorn MP, Bisanzio D, Fullman N, Huynh CK, Kulikoff X, Kutz MJ, Lopez AD, Mokdad AH, Naghavi M, Nguyen G, Shackelford KA, Vos T, Wang H, Lim SS, Murray CJL, Price RN, Baird JK, Smith DL, Bhatt S, Weiss DJ, Hay SI, Gething PW. Mapping the global endemicity and clinical burden of Plasmodium vivax, 2000-17: a spatial and temporal modelling study. Lancet. 2019; 394(10195):332-343.
Black D. Health inequalities. Report of a Research Working Group. London, England: Department of Health and Social Security, 1980.
Blay C, Limón E (coords). Bases para un modelo catalán de atención a las personas con necesidades complejas. Conceptualización e introducción a los elementos operativos, 2017 [Disponible en: http://salutweb.gencat.cat/web/.content/_ambits-actuacio/Linies-dactuacio/Estrategies-de-salut/Cronicitat/Documentacio-cronicitat/arxius/bases_modelo_personas_complejidad_v_6.pdf, último acceso el 30 de Julio de 2019].
Bolíbar B, Prados A, Gervás J, Juncosa S, Carrillo E. Sistemas de clasificación en grupos iso-consumo (case-mix) en atención ambulatoria. Perspectivas para nuestra atención primaria. Atención Primaria 1996; 17(1):74-83.
Bouchard L, Albertini M, Batista R, de Montigny J. Research on health inequalities: A bibliometric analysis (1966-2014). Social Science & Medicine 2015; 141:100-108.
Bowen W. An analytical review of environmental justice research: what do we really know? Environmental Management 2002; 29(1):3-15.
Brulle RJ, Pellow DN. Environmental justice: human health and environmental inequalities. Annual Review of Public Health 2006; 27:103-124.
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.
Cricelli C, Mazzaglia G, Samani F, Marchi M, Sabatini A, Nardi R, Ventrigila G, Caputti AP. Prevalence estimates for chronic diseases in Italy: exploring the differences between self-report and primary care databases. Journal of Public Health 2003; 25(3):254–257.
Deguen S, Zmirou-Navier D. Social inequalities resulting from health risks related to ambient air quality – A European review. European Journal of Public Health 2010; 28(1):27-35.
Díez-Manglano J, Cabrerizo-García JL, García-Arilla E, Jimeno A, Calvo E, Martínez-Álvarez RM, Bejarano E, Caudevilla A. External validation of the PROFOUND index in polypathological patients from internal medicine and acuate geriatrics departments in Aragón. Intern Emerg Med. 2015; 10(8):915-926.
Duque I, Domínguez‐Berjón MF, Cebrecos A, Prieto‐Salceda MD, Esnaola S, Calvo Sánchez M, Marí‐Dell’Olmo M. Índice de privación en España por sección censal en 2011: hacia un seguimiento de la desigualdad en áreas pequeñas durante la crisis económica, 2019 [Disponible en: https://www.seepidemiologia.es/documents/dummy/ManualdeusodelIP2011.pdf, último acceso el 30 de Julio de 2019].
Estupiñán-Ramírez M, Tristancho-Ajamil R, Company-Sancho MC, Sánchez-Janáriz H. Comparación de modelos predictivos para la selección de pacientes de alta complejidad. Gac Sanit. 2019; 33(1):60-65.
Evans GW, Kantrowitz E. Socioeconomic status and health: the potential role of environmental risk exposure. Annual Review of Public Health 2002; 23:303-331.
Feachem RG, Sekhri NK, White KL. Getting more for their dollar: a comparison of the NHS with California’s Kaiser Permanente. British Medical Journal 2002;324:135-41.
Hajat A, Hsia C, O’Neill MS. Socioeconomic disparities and air pollution exposure: A global review. Current Environmental Health Reports 2015; 2(4):440-450.
Heckman J. The common structure of statistical models of truncation, sample selection, and limited dependent variables and a simple estimator for such models. Annals of Economic and Social Measurement 1976; 5:475–492.
Heckman J. Sample selection bias as a specification error. Econometrica 1979; 47(1):153–162.
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.
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.
Kunst AE. Describing socioeconomic inequalities in health in European countries: an overview of recent studies. Revue d’Epidemiologie et de Sante Publique 2007;55:3–11.
Lalonde M. A new perspective on the health of Canadians. A working document. Ottawa: Government of Canada, 1974.
Laurent O, Bard D, Filleul L, Segala C. Effect of socioeconomic status on the relationship between atmospheric pollution and mortality. Journal of Epidemiology and Community Health 2007; 61(8):665-675.
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.
Mackenbach JP, Stirbu I, Roskam AJ, Schaap MM, Menvielle G, Leinsalu M, Kunst AE. Socioeconomic inequalities in health in 22 European countries. New England Journal of Medicine 2008; 358(23):2468-2481.
Marmot M. Social determinants of health inequalities. Lancet 2005; 365(9464):1099-1104.
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].
Ministerio de Sanidad, Servicios Sociales e Igualdad. Estrategia para el abordaje de la cronicidad en el Sistema Nacional de Salud. Madrid: Ministerio de Sanidad, Servicios Sociales e Igualdad, 2012 [Disponible en: https://www.mscbs.gob.es/organizacion/sns/planCalidadSNS/pdf/ESTRATEGIA_ABORDAJE_CRONICIDAD.pdf, último acceso el 30 de Julio de 2019].
Monterde D, Vela E, Clèries M y grupo colaborativo GMA. Los grupos de morbilidad ajustados: nuevo agrupador de morbilidad poblacional de utilidad en el ámbito de la atención primaria. Aten Primaria. 2016; 48(10):674-682.
Nuño R. Buenas prácticas en gestión sanitaria: el caso Kaiser Permanente. Rev Adm San. 2007; 5(2):283-292
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.
Saez M, Barceló MA, Coll-de-Tuero G. A selection-bias free method to estimate the prevalence of hypertension from an administrative primary health care database in the Girona Health Region, Spain. Comput Methods Programs Biomed. 2009; 93Ç(3):228-240.
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 1 de Agosto de 2019].
Starfield BH, Weiner JP, Mumford LM, Steinwachs DM. Ambulatory Care Groups: a categorization of diagnoses for research and management. Health Ser Res. 1991; 26(1):53-74.
Weiner JP, Starfield BH, Steinwachs DM, Mumford LM. Development and application of a population-oriented measure of ambulatory care case-mix. Med Care. 1991; 29:452-472.
Weiss DJ, Lucas TCD, Nguyen M, Nandi AK, Bisanzio D, Battle KE, Cameron E, Twohig KA, Pfeffer DA, Rozier JA, Gibson HS, Rao PC, Casey D, Bertozzi-Villa A, Collins EL, Dalrymple U, Gray N, Harris JR, Howes RE, Kang SY, Keddie SH, May D, Rumisha S, Thorn MP, Barber R, Fullman N, Huynh CK, Kulikoff X, Kutz MJ, Lopez AD, Mokdad AH, Naghavi M, Nguyen G, Shackelford KA, Vos T, Wang H, Smith DL, Lim SS, Murray CJL, Bhatt S, Hay SI, Gething PW. Mapping the global prevalence, incidence, and mortality of Plasmodium falciparum, 2000-17: a spatial and temporal modelling study. Lancet 2019; 394(10195):322-331.