Procesos puntuales. Prevalencia de enfermedad de motoneurona y sus fenotipos en Comunidades Autónomas Españolas.
GRECS, Universidad de Girona, España
Maria Antònia Barceló y Marc Saez
17 de Octubre de 2019
La Enfermedad de Motoneurona (EMN, nomenclatura ORPHANET sobre enfermedades raras, ORPHA 98503) es una enfermedad neurodegenerativa caracterizada por una pérdida progresiva de las neuronas motoras superiores e inferiores que conduce a la atrofia muscular, parálisis y muerte del paciente, normalmente debido a un fallo respiratorio (Kiernan et al., 2011; Boylan, 2015; Ingre et al., 2015; Mancuso et al., 2015).
El término de Enfermedad de la Motoneurona (ENM) hace referencia a un grupo de enfermedades neurodegenerativas, frecuentemente esporádicas, que se caracterizan patológicamente por la presencia de agregados de proteínas en el citoplasma, en la inmensa mayoría de los casos de TDP-43. Bajo este término se engloban con frecuencia la Esclerosis Lateral Amiotrófica y enfermedades relacionadas, como la Esclerosis lateral Primaria y la Atrofia Muscular Progresiva. Esta subdivisión de ENM en las tres enfermedades anteriores se realiza según el grado de afectación clínica y neurofisiológica de la motoneurona superior e inferior (MNS y MNI, respectivamente). Aunque las tres comparten causas, anatomía patológica y mecanismos fisiopatológicos, se considera relevante esta distinción por el diferente pronóstico y por la variedad de diagnósticos diferenciales que se plantean en cada una de ellas.
La Esclerosis Lateral Amiotrófica (ELA, ORPHA 803) se caracteriza por una afectación combinada de primera y segunda motoneurona y una evolución rápida (media de supervivencia de 3 a 4 años).
La Esclerosis Lateral Primaria (ELP, ORPHA 35689) se caracteriza por una afectación clínica y neurofisiológica exclusiva de MNS durante al menos cuatro años, con curso lento y con frecuencia estacionario a largo plazo, pudiendo llegar a presentar supervivencias normales. El periodo de cuatro años se exige porque es frecuente que en ese periodo se desarrollen signos de MNI y el paciente pueda ser clasificado en ese momento como ELA. La ELP se presenta con frecuencia en forma de paraparesia espástica asimétrica por lo que el diagnóstico diferencial se amplía a otras causas de paraparesia progresiva (paraparesia espástica hereditaria, esclerosis múltiple primaria progresiva, paraparesia espástica tropical, etc.).
La Atrofia Muscular Progresiva (AMP, ORPHA 454706) se caracteriza por una afectación clínica exclusiva de MNI durante al menos cuatro años. Aunque algunos autores no utilizan el criterio de los cuatro años, éste parece exigible para igualarlo a la ELP y para diferenciar claramente la AMP de la ELA, ya que con frecuencia la ELA puede debutar sin que los signos de motoneurona superior sean evidentes.
La ELA es una enfermedad muy heterogénea (o son muchas enfermedades distintas que comparten en gran medida causas y mecanismos). Esto hace que lo que actualmente categorizamos como ELA pueda adquirir presentaciones muy diversas en la práctica clínica. Para tratar de recoger esta complejidad, se han ido describiendo distintos fenotipos de enfermedad que, por lo general, hacen referencia a las características clínicas presentes en el momento del diagnóstico. Al contrario que las categorías diagnósticas, los fenotipos no constituyen entidades propias con causas, patologías y cursos clínicos claramente diferenciados, sino que son subgrupos de pacientes agrupados por compartir determinadas características clínicas.
Se han identificado diferentes fenotipos y presentaciones de la ELA, definiéndose en función de dónde aparecen en primer lugar los síntomas de la enfermedad y/o las regiones afectas en el momento del diagnóstico (ELA espinal y ELA bulbar); en función de si la motoneurona más afectada es la superior, la inferior o ambas (ELA clásica, ELA MNI y ELA MNS); en función del grado de afectación cognitiva (ELA demencia frontotemporal, ELA deterioro cognitivo disejecutivo, ELA deterioro cognitivo no disejecutivo y ELA alteración leve del comportamiento); en función de la edad de inicio (ELA juvenil, en menores de 25 años y ELA de inicio en la juventud, en menores de 45 años) y en función de cuán rápidamente progresa la enfermedad (rápida, menos de 2 años de supervivencia desde el diagnóstico de la enfermedad, media, de 2 a 5 años de supervivencia desde el diagnóstico de la enfermedad, o lenta, más de 5 años de supervivencia desde el diagnóstico de la enfermedad).
Existen distintas formas de clasificar a los pacientes de ELA según la característica clínica en torno a la cuál los queramos agrupar. Estos son las más frecuentes:
- Según el grado de afectación de MNS y MNI, distinguiríamos tres fenotipos:
- ELA clásica: combinación de signos de MNS y MNI
- ELA MNI: mínimos signos clínicos o neurofisiológicos de MNS
- ELA MNS: mínimos signos clínicos o neurofisiológicos de MNI
- Según la región de inicio y/o las regiones afectas en el momento del diagnóstico:
- ELA bulbar. Se habla también de parálisis bulbar progresiva (PBP) o parálisis pseudobulbar progresiva para designar a aquellos pacientes con afectación aislada de la MNI o MNS respectivamente, en la región bulbar, en el momento del diagnóstico. Hay que tener en cuenta que la parálisis pseudobulbar progresiva no siempre va a evolucionar a una ELA ya que en algunos casos puede hacerlo a una ELP
- ELA espinal. Incluye a su vez otros subfenotipos:
- Flail leg: inicio asimétrico como pie caído con predominio de MNI y sin progresión a otro territorio durante al menos un año
- Pseudopolineurítica: inicio simétrico como pie caído bilateral con predominio de MNI y patrón distal proximal (patrón polineuropático)
- Flail arm: inicio proximal en MMSS, asimétrico y predominio MNI
- Hemipléjica (síndrome de Mills): inicio hemicorporal, predominio MNS
- Respiratoria/axial
- Según el grado de afectación cognitiva:
- ELA – demencia frontotemporal
- ELA – deterioro cognitivo disejecutivo
- ELA – deterioro cognitivo no disejecutivo
- ELA – alteración leve del comportamiento
- Según la edad de inicio:
- ELA juvenil: en el momento del diagnóstico de la enfermedad el paciente tiene menos de 25 años.
- ELA de inicio en la juventud: en el momento del diagnóstico de la enfermedad el paciente tiene menos de 45 años.
- Según la supervivencia:
- ELA supervivencia corta: tiempo entre diagnóstico y muerte menor a 2 años.
- ELA supervivencia media: tiempo entre diagnóstico y muerte entre 2 y 5 años.
- ELA supervivencia larga: tiempo entre diagnóstico y muerte mayor a 5 años.
- Según la progresión de la enfermedad:
- ELA rápidamente progresiva: retraso diagnóstico menor a 8 meses
- ELA media: retraso diagnóstico entre 8 y 18 meses
- ELA lentamente progresiva: retraso diagnóstico mayor a 18 meses
- Según la progresión de la enfermedad en función de la pérdida de puntos por mes en la Adaptación Española de la Escala revisada de Valoración Funcional de la Esclerosis Lateral Amiotrófica, ALSFRS-R (Salas et al., 2010):
- ELA rápida: pérdida >1-4 puntos
- ELA media:
- ELA lenta: pérdida <0-4 puntos
- Según la historia familiar:
- ELA familiar
- ELA esporádica
La principal importancia de estos diferentes fenotipos está en su pronóstico, debido a que algunos son menos incapacitativos y se extienden más lentamente que otros y, desafortunadamente, otros tienen un pronóstico peor en los problemas señalados (Chiò et al., 2013; Wolf et al., 2014). El fenotipo espinal representa el 42% del total de casos de ELA; el fenotipo bulbar, el 33,5% del total de casos; seguidos del fenotipo flail leg (8,5%), el fenotipo flail arm (6,5%), el piramidal (5%) y el respiratorio (4,5%) (Chiò et al., 2013; Wolf et al., 2014). Las mutaciones genéticas no explican esta enfermedad totalmente, ya que una misma mutación puede estar asociada a una gran variabilidad de fenotipos de ELA.
La ELA familiar heredada dominante (fALS), para la cuál la heredabilidad es principalmente dominante autosomal, sólo representa el 10% de todos los pacientes con ELA, mientras que la ELA esporádica (sALS), que no tiene una heredabilidad aparente, es la forma más común de la ELA. Aunque se ha hecho un gran progreso para entender la genética de fALS, las causas subyacentes de sALS permanecen desconocidas (Al-Chalabi et al., 2013; Schwartz et al., 2017). Hasta ahora los únicos factores de riesgo reconocidos para la ELA son la edad avanzada, ser hombre y tener una historia familiar de ELA (Belbasis et al., 2016).
Incidencia y prevalencia de la ELA
Marin et al., en un meta-análisis muy reciente, que incluyó 825 millones de personas-año de seguimiento y cubrió 45 áreas geográficas en 11 sub-continentes, estiman la incidencia cruda mundial global de la ELA en 1,75 por 100.000 personas/año de seguimiento (95% IC: 1,55-1,96) y una incidencia estandarizada de 1,68 (95% IC: 1,50-1,85) (Marin et al., 2017). Sin embargo, la distribución geográfica de la incidencia de la ELA es muy heterogénea. La incidencia de la ELA en Europa se ha estimado entre 2 y 4 casos por 100000 personas/año (Al-Chalabi et al., 2013). Este estimador es algo menor en el norte de Europa con una incidencia estandarizada igual a 1,89 por 100.000 personas-año (95% IC: 1,46-2,3) (Marin et al., 2017) y algo mayor en Europa occidental. Chiò et al., en una revisión sistemática de estudios observacionales con base poblacional, que incluía 37 artículos, estimaron la incidencia en Europa Occidental en el periodo comprendido entre 1995 y 2011, en 5,40 por 100.000 y 4,70 por 100.000 en 2012 (Chiò et al., 2013). En España, Santurtun et al. encontraron una incidencia anual en el rango de 1 a 3 casos por 100.000 habitantes en 2013 (Santurtún et al., 2016). Para Cataluña (en España), también para 2013, Pradas et al. estimaron una incidencia anual de 1,4 por 100.000 habitantes (Pradas et al., 2013). En Asia, la incidencia es mucho menor que en Europa. La incidencia estandarizada de la ELA es de 0,83 por 100000 personas-año (95% IC: 0,42-1,24) para Asia Oriental, y 0,73 por 100.000 personas-año (95% IC: 0,58-0,89) para Asia del sur (Marin et al., 2017). Sin embargo, en Japón y para el mismo periodo, la incidencia de la ELA se estimó en 2,2 por 100.000 personas-año (Doi et al., 2014).
Respecto a la prevalencia de la ELA se ha estimado entre 5 por 100.000 habitantes (Tesauro et al., 2017) y 5,4 por 100.000 habitantes (Belbasis et al., 2016). También en el caso de la prevalencia, existe una importante variabilidad geográfica. Las cifras más altas se encuentran en Japón, 9,9 por 100000 habitantes en 2013 (Doi et al., 2014), prácticamente el doble que en los Estados Unidos (3,9 por 100.000 habitantes en 2010-2011 - Doi et al., 2014-, 4,3 en 2013 - Mehta et al., 2016-) y Europa (4,06 por 100.000 habitantes en Europa Occidental - Chiò et al., 2013-), ambos con prevalencias similares. Pradas et al., estimaron una prevalencia de 5,4 por 100.000 habitantes para Cataluña, España (Pradas et al., 2013).
La heterogeneidad, tanto en la incidencia como en la prevalencia de la ELA, también se presenta en otras variables como el sexo y la edad. La incidencia es entre 1,3 y 1,5 veces mayor en hombres que en mujeres (Mehta et al, 2016; Marin et al., 2017; Pradas et al., 2013). La ELA es muy poco frecuente antes de los 40 años, con una edad media de aparición que varía entre los 58 y 63 años para la forma esporádica y entre los 40 y 60 años para la forma familiar (Belbasis et al., 2016), con poca variación entre fenotipos (Belbasis et al., 2016; Wolf et al., 2014). Después de la edad de aparición, la incidencia de la ELA aumenta alcanzando un máximo después de los 60 años. Chiò et al. encuentran que la incidencia de la ELA alcanza su máximo entre los 60 y 75 años de edad (Chiò et al., 2013). En el estudio Cima este máximo ocurre entre los 65 y 74 años de edad (Cima et al., 2009), y en varios estudios en Nueva Escocia, Canadá, entre los 70 y 79 años (Bonaparte et al., 2007).
La heterogeneidad en la incidencia y la prevalencia de la ELA puede ser una consecuencia del resultado de criterios de diagnóstico diferentes, diferente práctica clínica, formas diferentes de recoger los casos, así como el resultado de una interacción entre factores genéticos y no genéticos (Caller et al., 2012; Caller et al., 2013; Al-Chalabi et al., 2013; Trojsi et al., 2013; Modgil et al., 2014; Ingre et al., 2015; Oskarsson et al., 2105).
Existen en España un número muy reducido de estudios que estimen la prevalencia de la ELA (Pradas et al., 2013; Santurtún et al., 2016) y ninguno de ellos proporciona estimaciones de la misma según los diferentes fenotipos.
Nuestro principal objetivo es el de estimar la prevalencia de la enfermedad de la motoneurona, así como de sus diferentes fenotipos, en las Comunidades Autónomas de Cataluña, Valencia y Madrid en el periodo 2011-2019.
1 Datos
Indicaremos cuál es nuestro directorio de trabajo utilizando la función setwd().
setwd("~/Documents/Maria/Treballs Varis/ELA/Article Barcelona-Valencia")Importaremos los datos.
# Data
library(foreign)
ALSprevalence=read.spss("ELA_Valencia.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.
ALSprevalence es un objeto tipo data frame que contiene las siguientes variables:
- censustrac : una variable indicando la sección censal de cada una de las secciones censales que conforman cada uno de los municipios de Valencia
- 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
- GEO_x, GEO_y : las coordenadas geográficas del lugar de residencia de cada individuo
- UTM_x, UTM_y : las coordenadas UTM del lugar de residencia de cada individuo
- UTM_x_ct, UTM_y_ct : las coordenadas UTM del centroide de la sección censal donde reside el individuo
- id : identificador de cada uno de los individuos
- code : código individuo
- date_of_birth : edad del individuo
- province_of_birth : provincia de nacimiento del individuo
- gender : sexo del individuo
- age_diag, age_diagQ : edad del individuo, continua y en quintiles
- exitus : variable que indica si el individuo ha fallecido o no
- date_of_exitus : fecha de fallecimiento del individuo
- province_residence : provincia de residencia del individuo
- municipality_residence : municipio de residencia del individuo
- zip_code : código postal de residencia del individuo
- occupation : ocupación del individuo
- date_of_visit_clinical_data : fecha de la primera visita
- date_of_clinical_onset : fecha de inicio de los síntomas
- date_of_diagnosis_ALS : fecha de diagnóstico de ELA
- phenotype_at_diagnosis : fenotipo al diagnóstico
- site_of_onset : lugar de inicio de los síntomas
- definite_diagnosis : diagnóstico definitive (AMP, ELA, ELP, flail arm)
- phenotype_at_exitus : fenotipo al exitus
- genetics_studied : genética estudiada
- genotype : genotipo
- family_history_ALS : historia familiar de ELA
- degree_relative_familial_ALS : grado de parentesco
- codigo_familia : código de la familia. Este código es compartido para cada uno de los miembros de la familia.
- family_history_neurodegenerative : historia familiar de enfermedad neurodegenerativa
- family_history_neurodegenerative_who_what : ¿quién en la familia ha presentado alguna enfermedad neurodegenerativa? ¿de qué enfermedad neurodegenerativa se trataba?
- smoker : fumador (0 no, 1 sí)
- paquets_year : paquetes al año fumados
- tracheostomy : traqueotomía (0 no, 1 sí)
- date_placement_tracheostomy : fecha de colocación de la traqueotomía
- ethnicity : etnia
- ALSFRSR : Escala de valoración funcional para ELA revisada
- municipality_birth_code : código del municipio de nacimiento
- municipality_birth : municipio de nacimiento
- country_birth : país de nacimiento
- MND : una variable que indica si el individuo es un enfermo de motoneurona, (1) o no (0)
- cases_MND : una variable que indica el número de MND en la sección censal donde reside el individuo
- ALS : una variable que indica si el individuo es un enfermo de ELA, (1) o no (0)
- cases_ALS : una variable que indica el número de casos de ELA en la sección censal donde reside el individuo
- PMA : una variable que indica si el individuo es un enfermo de atrofia muscular progresiva, (1) o no (0)
- cases_PMA : una variable que indica el número de casos de PMA en la sección censal donde reside el individuo
- PLS : una variable que indica si el individuo es un enfermo de esclerosis lateral amiotrófica, (1) o no (0)
- cases_PLS : una variable que indica el número de casos de PLS en la sección censal donde reside el individuo
- flail_arm : una variable que indica si el individuo es un enfermo de flail arm, (1) o no (0)
- cases_flail_arm : una variable que indica el número de casos de flail arm en la sección censal donde reside el individuo
- ALS_spinal : una variable que indica si el individuo está diagnosticado de ELA espinal, (1) o no (0)
- cases_ALS_spinal : una variable que indica el número de casos de ELA espinal en la sección censal donde reside el individuo
- ALS_bulbar : una variable que indica si el individuo está diagnosticado de ELA bulbar, (1) o no (0)
- cases_ALS_bulbar : una variable que indica el número de casos de ELA bulbar en la sección censal donde reside el individuo
- ALS_dementia : una variable que indica si el individuo está diagnosticado de ELA que inicia con demencia, (1) o no (0)
- cases_ALS_dementia : una variable que indica el número de casos de ELA demencia en la sección censal donde reside el individuo
- ALS_respiratory : una variable que indica si el individuo está diagnosticado de ELA respiratoria, (1) o no (0)
- cases_ALS_respiratory : una variable que indica el número de casos de ELA respiratoria en la sección censal donde reside el individuo
- 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 (Duque et al., 2019).
- men, male_under_16, male_16to64, male_65ormore : variables que indican la población de hombres y 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 Valencia
- female, female_under_16, female_16to64, female_65ormore : variables que indican la población de mujeres y 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 Valencia
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)## Loading required package: spp<-shapefile("Valencia_ct.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 ENM 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 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 unidades funcionales de neurona motora presentes en cada uno de los hospitales de referencia de las Comunidades Autónomas que estudiamos 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 pacientes sobre la población cubierta por cada hospital, 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 cubierta por cada hospital.
Sin embargo, es muy probable que los individuos contactasen con las unidades (de las que provienen los datos), no sólo porque tuviesen alguno de los síntomas de inicio de la enfermedad, sino que también podrían existir factores no observados que influyesen en la utilización de la unidad y 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 pacientes hayan sido observados no es la misma para todos ellos. Así pues, 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 contactado con la unidad funcional del Hospital Politécnico y Universitario La Fé). En concreto, variables individuales, como el sexo y la edad; variables contextuales, como el índice de privación y la población estratificada por sexo y grupos de edad, de la sección censal donde reside el individuo.
En la segunda parte utilizamos un GLMM con vínculo Poisson truncada para modelizar el número de casos de ELA por sección censal. Cabe señalar que utilizamos este vínculo y no una Poisson, porque el número de secciones censales donde no existe ningún caso de ELA es muy elevado (86,13%). 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.first, id.second).
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 caso de la variable respuesta (en nuestro caso, ELA) 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 ELA 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}^{5}β_{2k}age\_diagQ_{ik}+\sum_{k=2}^{5}β_{3k}deprivationQ_{ik}\\ +β_4 log(men\_16to64_i)+β_5 log(men\_65ormore_i)+β_6 log(women\_16to64_i)+β_7 log(women\_65ormore_i)+offset(population\_over16_i)\)
donde el subíndice i indica la sección censal de los municipios de cada una de las Comunidades Autónomas; 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 la variable respuesta en cada una de las secciones censales, \(θ(x_i ), i=1,…,n\), se distribuye como una Poisson truncada.
\(Y_i |θ(x_i )\sim~Poisson\_truncated(θ(x_i ))\)
En este caso la función vínculo es la que sigue,
\(log(θ(x_i ))=β_0+η_{1i}+S(x_i)+β_1 gender_i+\sum_{k=2}^{5}β_{2k}age\_diagQ_{ik}+\sum_{k=2}^{5}β_{3k}deprivationQ_{ik}\\ +β_4 log(men\_16to64_i)+β_5 log(men\_65ormore_i)+β_6 log(women\_16to64_i)+β_7 log(women\_65ormore_i)\)
Siendo la definición de los parámetros, de las variables y de los efectos aleatorios, la misma que antes.
2.2 Preparación de los datos
Mediante la instrucción attach activaremos el data frame.
#Data Management
attach(ALSprevalence)## The following object is masked from package:datasets:
##
## womenComo 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, ALS y casos de ALS, respectivamente. También, definimos el identificador del efecto aleatorio para cada una de las partes del modelo, id.first y id.second.
#Data prepation
n<- dim(ALSprevalence)[1]
N<- n*2
Y = matrix(NA, N, 2)
Y[1:n, 1] = ALS
Y[1:n + n, 2] = cases_ALS
id.first <- c(id,id)
id.second <- c(id,id)Ahora haremos lo mismo para el resto de covariables (gender, age_diagQ, deprivationQ, men_16to64, men_65ormore, women_16to64, women_65ormore).
#Data prepation
#Construction of the matrices of the covariables
Gender= numeric(N)
Gender[1:n] = gender
Gender[1:n + n] = gender
Age_diagQ= numeric(N)
Age_diagQ[1:n] = age_diagQ
Age_diagQ[1:n + n] = age_diagQ
DeprivationQ = numeric(N)
DeprivationQ[1:n] = deprivationQ
DeprivationQ[1:n + n] = deprivationQ
Men_16to64=numeric(N)
Men_16to64[1:n] = men_16to64
Men_16to64[1:n + n] = men_16to64
Men_65ormore=numeric(N)
Men_65ormore[1:n] = men_65ormore
Men_65ormore[1:n + n] = men_65ormore
Women_16to64=numeric(N)
Women_16to64[1:n] = women_16to64
Women_16to64[1:n + n] = women_16to64
Women_65ormore=numeric(N)
Women_65ormore[1:n] = women_65ormore
Women_65ormore[1:n + n] = women_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 siguientes 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()mesh = inla.mesh.2d(cbind(UTM_x_ct,UTM_y_ct), max.edge=c(25000,25000), cutoff=0.3)
plot(mesh)
A continuación, añadiremos la localización de la residencia de cada sujeto.
plot(mesh)
points(cbind(UTM_x[ALS==1],UTM_y[ALS==1]), 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
# Spatial index
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 con ELA, 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_ALS<- data.frame(id.first,id.second,Field,Y,Gender,Age_diagQ,DeprivationQ,Men_16to64,Men_65ormore,Women_16to64,Women_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_diagQ) + factor(DeprivationQ) + log(Men_16to64)+log(Men_65ormore) + log(Women_16to64) + log(Women_65ormore) +
f(id.first, model="iid",
hyper= list( prec= list( prior="pc.prec", param=c(0.5, 0.01)))) + f(id.second,copy="id.first")+
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.first".
Estimamos ahora el modelo mediante la instrucción inla():
#Call INLA
result= inla(formula, data=data_ALS, family= c("binomial","zeroinflatedpoisson0"), 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))En 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 truncada (family= c("binomial","zeroinflatedpoisson0")); que utilice la base de datos data_ALS; 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\",
## \"zeroinflatedpoisson0\"), ", " data = data_ALS, 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 = 3.15, Running = 392, Post = 2.58, Total = 398
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) -8.053 0.778 -9.606 -8.044 -6.549 -8.027
## Gender 1.247 0.117 1.016 1.247 1.476 1.248
## factor(Age_diagQ)2 2.171 0.240 1.702 2.170 2.644 2.168
## factor(Age_diagQ)3 2.055 0.244 1.580 2.054 2.536 2.052
## factor(Age_diagQ)4 2.020 0.252 1.529 2.019 2.516 2.017
## factor(Age_diagQ)5 1.789 0.260 1.283 1.788 2.304 1.785
## factor(DeprivationQ)2 -0.038 0.158 -0.349 -0.038 0.271 -0.037
## factor(DeprivationQ)3 0.071 0.162 -0.247 0.071 0.389 0.072
## factor(DeprivationQ)4 -0.252 0.187 -0.624 -0.250 0.111 -0.246
## factor(DeprivationQ)5 -0.228 0.195 -0.617 -0.227 0.150 -0.223
## log(Men_16to64) -0.398 0.261 -0.909 -0.399 0.115 -0.400
## log(Men_65ormore) -0.038 0.103 -0.234 -0.041 0.172 -0.046
## log(Women_16to64) -0.208 0.271 -0.738 -0.209 0.324 -0.210
## log(Women_65ormore) -0.089 0.103 -0.285 -0.091 0.118 -0.095
## kld
## (Intercept) 0
## Gender 0
## factor(Age_diagQ)2 0
## factor(Age_diagQ)3 0
## factor(Age_diagQ)4 0
## factor(Age_diagQ)5 0
## factor(DeprivationQ)2 0
## factor(DeprivationQ)3 0
## factor(DeprivationQ)4 0
## factor(DeprivationQ)5 0
## log(Men_16to64) 0
## log(Men_65ormore) 0
## log(Women_16to64) 0
## log(Women_65ormore) 0
##
## Random effects:
## Name Model
## id.first IID model
## Field SPDE2 model
## id.second Copy
##
## Model hyperparameters:
## mean
## zero-probability parameter for zero-inflated poisson_0[2] 9.12e-01
## Precision for id.first 1.07e+04
## Range for Field 1.32e-01
## Stdev for Field 4.39e+01
## sd
## zero-probability parameter for zero-inflated poisson_0[2] 5.00e-03
## Precision for id.first 1.11e+05
## Range for Field 2.40e-01
## Stdev for Field 4.29e+01
## 0.025quant
## zero-probability parameter for zero-inflated poisson_0[2] 0.903
## Precision for id.first 132.496
## Range for Field 0.013
## Stdev for Field 3.434
## 0.5quant
## zero-probability parameter for zero-inflated poisson_0[2] 0.912
## Precision for id.first 1441.644
## Range for Field 0.067
## Stdev for Field 31.446
## 0.975quant
## zero-probability parameter for zero-inflated poisson_0[2] 9.21e-01
## Precision for id.first 6.94e+04
## Range for Field 6.61e-01
## Stdev for Field 1.58e+02
## mode
## zero-probability parameter for zero-inflated poisson_0[2] 0.913
## Precision for id.first 277.707
## Range for Field 0.028
## Stdev for Field 9.660
##
## Expected number of effective parameters(stdev): 18.26(4.33)
## Number of equivalent replicates : 384.54
##
## Deviance Information Criterion (DIC) ...............: 4023.14
## Deviance Information Criterion (DIC, saturated) ....: NaN
## Effective number of parameters .....................: 18.54
##
## Watanabe-Akaike information criterion (WAIC) ...: 4019.32
## Effective number of parameters .................: 14.17
##
## Marginal log-Likelihood: -2063.14
## Posterior marginals for the linear predictor and
## the fitted values are computed3 Representación en un mapa de la prevalencia de ELA
Por último, representaremos la prevalencia de ELA en el mapa de las secciones censales de la Comunidad Autónoma de Valencia, utilizando el seccionado del Censo de Población y Viviendas del 2011 (INE, 2019), así como la librería leaflet (Cheng, 2018). 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), expresados todos ellos por 100.000 habitantes.
# Mapping ALS prevalence
# Computation of the ALS prevalence and the 95% credibility interval
prev_ALS_mean= result$summary.fitted.values[(n+1):N,1]
prev_ALS_ll= result$summary.fitted.values[(n+1):N,3]
prev_ALS_ul= result$summary.fitted.values[(n+1):N,5]
summary(prev_ALS_mean*100000)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3449 1.2578 1.6004 8.8567 2.2125 423.8130summary(prev_ALS_ll*100000)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1529 0.7845 0.9944 5.4157 1.3440 188.5445summary(prev_ALS_ul*100000)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6633 1.8695 2.3908 13.5658 3.3795 875.1469Así pues, la estimación puntual de ELA por 100.000 habitantes en la Comunidad Autónoma de Valencia es de 8.8563 (Intervalo de credibilidad al 95% por 100.000 habitantes: 5.3848-13.623)
Representaremos en un mapa tipo raster la prevalencia. Para ello utilizaremos la librería raster (Hijmans, 2019).
# Mapping ALS prevalence
library(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 secciones censales de Valencia. 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 ALS.
#Rasterization of prevalence of ALS
r_ALS_mean=rasterize(p,rr,field=prev_ALS_mean*100000,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 100000 habitantes adultos (`prev_ALS_mean*100000).
Finalmente, representamos mediante leaflet el raster que hemos construido.
#Plot the raster
library(leaflet)
bins<-c(0,quantile(r_ALS_mean,probs=seq(0,1,0.25)),Inf)
pal<-colorBin(palette="YlOrRd",domain=prev_ALS_mean*100000,bins=bins,na.color="transparent")
leaflet(p)%>%addTiles()%>%
addRasterImage(r_ALS_mean,colors=pal,opacity=0.5)%>%
addLegend("bottomright",pal=pal,values=values(r_ALS_mean),title="Prevalence of ALS (per 100,000 inhabitants)")%>%
addScaleBar(position=c("bottomleft")) Nótese que en la representación hemos categorizado las prevalencias en cuartiles (bins<-quantile(r_ALS_mean,probs=seq(0,1,0.25)) y pal<-colorBin(…, bins=bins,…)).
Alternativamente, podemos utilizar la librería leaflet (Cheng, 2018) para representar la prevalencia (y su intervalo de credibilidad al 95%), en un mapa de secciones censales según el seccionado del Censo de Población y Viviendas del 2011 (INE, 2019).
Para agregar los polígonos de las secciones censales debemos utilizar un mapa con las coordenadas geográficas (longitud y latitud). Por esta razón, utilizando la librería sp (Pebesma y Bivand, 2005; Bivand et al., 2013) transformamos el mapa, originalmente con coordenadas UTM, a otro (que llamaremos p_src) con coordenadas geográficas (longitud y latitud).
# Add the polygons (of census tract)
# In order to add the polygons (addPolygons) you must use a map with the geographic coordinates (not UTM). For this reason, we transformed them before (now the map will be called p_src)
library(sp)
p_src<-spTransform(p,CRS("+init=epsg:4326"))
#summary(p_src)Añadiremos ahora a la base de datos del mapa p_src las variables que queremos representar: el nombre del municipio al que pertenece la sección censal (municipality_residence), el estimador de la prevalencia (para cada sección censal) (prev_ALS_mean) y sus intervalos de credibilidad al 95% (prev_ALS_ll y prev_ALS_ul) y la población a riesgo en la sección censal, la de 16 años o mayor (population_over16), además del código de la sección censal (censustract).
Previamente, sin embargo, debemos hacer algunos arreglos tanto en el mapa como en el data frame d, formado por las variables que queremos representar.
# Add data to map
d=data.frame(censustrac,municipality_residence,prev_ALS_mean,prev_ALS_ll,prev_ALS_ul,population_over16)En primer lugar, pondremos el código de la sección censal (CENSUSTRACT) como nombre de los polígonos del mapa. Utilizaremos éste como variable clave en la fusión de la data frame d y de la base de datos del mapa p_src@data.
# We must put the name to the polygons. This instruction does not put them.
# names(polygons(p_src))
#names(polygons(p_src))
for (i in 1:nrow(p_src@data)){
p_src@polygons[[i]]@ID=as.character(p_src@data$CENSUSTRAC[i])
}Para que el programa R reconozca la variable clave en ambas bases de datos (d y p_src@data), pondremos en minúscula los nombres de las variables de la base de datos del mapa.
#names(polygons(p_src))
names(p_src)<-casefold(names(p_src))
names(p_src)## [1] "objectid" "censustrac" "cumun" "cca" "cpro"
## [6] "cmun" "dist" "secc" "cudis" "obs"
## [11] "cnut0" "cnut1" "cnut2" "cnut3" "clau2"
## [16] "npro" "nca" "nmun" "shape_leng" "shape_area"
## [21] "shape_len"Utilizando la librería stringr (Wickham, 2019), eliminaremos los espacios en blanco de la variable alfanumérica (carácter), municipality_residence del data frame d.
library(stringr)
d$municipality_residence=str_trim(d$municipality_residence,side="both")Mediante la instrucción group_by de la librería dplyr (Wickham et al., 2019), agregaremos por sección censal (censustract), los valores de las variables municipality_residence (cogeremos el primer valor), prev_ALS_mean, prev_ALS_ll, prev_ALS_ul y population_over16 (en estas cuatro variables cogeremos la media de los valores).
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniond_group <- group_by(d,censustrac) %>%
summarize(municipality_residence = first(municipality_residence),
prev_ALS_mean = mean(prev_ALS_mean),
prev_ALS_ll = mean(prev_ALS_ll),
prev_ALS_ul = mean(prev_ALS_ul),
population_over16=mean(population_over16))Uniremos la base de datos del mapa (p_src@data) y el data frame agregado d_group, por la variable clave censustrac, mediante la instrucción left_join (la base de datos del mapa tiene más filas que el data frame agregado) de las librería tidyverse (Wickham, 2017).
dim(d_group)## [1] 3473 6dim(p_src@data)## [1] 3478 21library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.1.0 ✔ readr 1.1.1
## ✔ tibble 2.1.1 ✔ purrr 0.2.5
## ✔ tidyr 0.8.1 ✔ forcats 0.3.0## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()d_group2=left_join(p_src@data,d_group,by="censustrac")Arreglaremos la base de datos agregada para evitar duplicados.
dim(d_group2)## [1] 3478 26names(d_group2)## [1] "objectid" "censustrac"
## [3] "cumun" "cca"
## [5] "cpro" "cmun"
## [7] "dist" "secc"
## [9] "cudis" "obs"
## [11] "cnut0" "cnut1"
## [13] "cnut2" "cnut3"
## [15] "clau2" "npro"
## [17] "nca" "nmun"
## [19] "shape_leng" "shape_area"
## [21] "shape_len" "municipality_residence"
## [23] "prev_ALS_mean" "prev_ALS_ll"
## [25] "prev_ALS_ul" "population_over16"d_group2=d_group2[,-c(1,3:21)]
names(d_group2)## [1] "censustrac" "municipality_residence"
## [3] "prev_ALS_mean" "prev_ALS_ll"
## [5] "prev_ALS_ul" "population_over16"Finalmente, tras poner a las filas del data frame agregado el nombre de la sección censal (variable clave), fusionaremos éste con la base de datos del mapa, utilizando la instrucción SpatialPolygonsDataFrame de la librería sp (Pebesma y Bivand, 2005; Bivand et al., 2013).
library(sp)
rownames(d_group2)<-d_group2$censustrac
p_src<-SpatialPolygonsDataFrame(p_src,d_group2,match.ID=TRUE)
head(p_src@data)## censustrac municipality_residence prev_ALS_mean prev_ALS_ll
## 300101001 300101001 Adsubia 2.019237e-05 1.220356e-05
## 300201001 300201001 Agost 1.481559e-05 9.092329e-06
## 300201002 300201002 Agost 9.106325e-06 5.633249e-06
## 300201003 300201003 Agost 1.398717e-05 8.567335e-06
## 300301001 300301001 Agres 2.223084e-05 1.332861e-05
## 300401001 300401001 Aigües 1.928925e-05 1.217569e-05
## prev_ALS_ul population_over16
## 300101001 3.084567e-05 625
## 300201001 2.250987e-05 1345
## 300201002 1.367155e-05 1740
## 300201003 2.117250e-05 995
## 300301001 3.414076e-05 505
## 300401001 2.858629e-05 880Ahora representaremos la prevalencia de ELA en un mapa interactivo por secciones censales utilizando leaflet.
# Map per census tract without labels
library(leaflet)
bins<-c(0,quantile(prev_ALS_mean*100000,probs=seq(0,1,0.25)),Inf)
pal<-colorBin(palette="YlOrRd",domain=prev_ALS_mean*100000,bins=bins,na.color="transparent")
l=leaflet(p_src)%>%addTiles()
l%>%addPolygons(color="grey",weight=1,fillColor = ~pal(prev_ALS_mean*100000),fillOpacity=0.5)%>%
addLegend("bottomright",pal=pal,values=values(r_ALS_mean),title="Prevalence of ALS (per 100,000 inhabitants)")%>%
addScaleBar(position=c("bottomleft"))