Repositorio Institucional de Documentos

M´aster en Modelizaci´on Matem´atica, Estad´ıstica y Computaci´on Trabajo fin de m´aster An´alisis espacial de riesgos de mortalidad por c´ancer de pr´ostata con INLA Autora: Paula Navarro Esteban Directora: Mar´ıa Dolores Ugarte Mart´ınez (Universidad P´ublica de Navarra) 3 de septiembre de 2012 Agradecimientos En primer lugar, me gustar´ıa expresar mi m´as sincero agradecimiento a mi directora Mar´ıa Dolores Ugarte Mart´ınez por darme la oportunidad de trabajar con ella y por la gran cantidad de horas que me ha dedicado. Me gustar´ıa tambi´en darle las gracias a mis padres, a mi hermana y a mis amigos, por estar a mi lado apoy´andome cuando m´as lo necesitaba. iii ´ Indice general Pr´ologo 1 Introducci´on 5 1. Informaci´on general sobre el c´ancer de pr´ostata 7 1.1. ¿ Qu´e es el c´ancer de pr´ostata? ............................ 7 1.2. ¿C´omo se diagnostica? ................................. 9 1.3. Factores de riesgo .................................... 10 2. Creaci´on de la base de datos 13 2.1. Definiciones b´asicas ................................... 13 2.1.1. Salud p´ublica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.2. Registro de c´ancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2. Instalaci´on y carga de librer´ıas ............................ 14 2.3. Lectura de los datos .................................. 16 2.4. C´alculo de valores esperados .............................. 17 2.4.1. Tasas de referencia o tasa global de mortalidad en Espa˜na . . . . . . . . . 19 2.4.2. Esperados.................................... 23 2.4.3. Observados ................................... 24 3. Medidas cl´asicas de estimaci´on de riesgos: raz´on de mortalidad estandarizada 27 v Paula Navarro Esteban 3.1. C´alculo de la SMR del c´ancer de pr´ostata en Espa˜na . . . . . . . . . . . . . . . . 29 3.1.1. C´alculo de la SMR del primer periodo respecto del primer periodo . . . . 29 3.1.2. C´alculo de la SMR del segundo periodo respecto del segundo periodo . . . 30 3.1.3. C´alculo de la SMR del segundo periodo respecto del primer periodo . . . 33 4. Aplicaci´on del modelo BYM con INLA 41 4.1. Descripci´on del modelo BYM ............................. 42 4.1.1. Estad´ıstica bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.1.2. ModelodeBYM ................................ 44 4.2. Lectura de datos y cartograf´ıa ............................. 47 4.2.1. Importancia de los mapas . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2.2. Cartograf´ıa ................................... 48 4.2.3. Construcci´on de la matriz de vecindad . . . . . . . . . . . . . . . . . . . . 51 4.3. Suavizaci´on de Riesgos ................................. 53 4.3.1. INLA ...................................... 53 4.3.2. An´alisis del primer periodo (estandarizaci´on interna con tasas espec´ıficas de referencia del primer periodo) . . . . . . . . . . . . . . . . . . . . . . . 58 4.3.3. An´alisis del segundo periodo (estandarizaci´on interna con tasas espec´ıficas de referencia del segundo periodo) . . . . . . . . . . . . . . . . . . . . . . 64 4.3.4. An´alisis del segundo periodo (estandarizaci´on interna con tasas espec´ıficas de referencia del primer periodo) . . . . . . . . . . . . . . . . . . . . . . . 69 5. Conclusiones 81 vi Indice de tablas 2.1. N´umero de muertes y poblaci´on total seg´un grupos de edad para cada periodo . 22 3.1. SMR e intervalos de confianza al 95 % para el periodo 1975-1991. Parte 1. . . . . 31 3.2. SMR e intervalos de confianza al 95 % para el periodo 1975-1991. Parte 2. . . . . 32 3.3. SMR e intervalos de confianza al 95 % para el periodo 1992-2008. Parte 1. . . . . 34 3.4. SMR e intervalos de confianza al 95 % para el periodo 1992-2008. Parte 2. . . . . 35 3.5. SMR e intervalos de confianza al 95 % para el periodo 1992-2008 respecto al periodo 1975-1991. Parte 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.6. SMR e intervalos de confianza al 95 % para el periodo 1992-2008 respecto al periodo 1975-1991. Parte 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1. Poblaci´on total seg´un provincias y periodos. Parte 1. . . . . . . . . . . . . . . . . 73 4.2. Poblaci´on total seg´un provincias y periodos. Parte 2. . . . . . . . . . . . . . . . . 74 4.3. Casos de mortalidad seg´un provincias y periodos. Parte 1. . . . . . . . . . . . . . 75 4.4. Casos de mortalidad seg´un provincias y periodos. Parte 2. . . . . . . . . . . . . . 76 vii ´ Indice de figuras 1.1. Anatom´ıa del sistema reproductor masculino . . . . . . . . . . . . . . . . . . . . 8 1.2. Pr´ostata normal y HPB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.1. Raz´on de mortalidad estandarizada para el periodo 1975-1991 . . . . . . . . . . . 33 3.2. Raz´on de mortalidad estandarizada para el periodo 1992-2008 . . . . . . . . . . . 36 3.3. Raz´on de mortalidad estandarizada para el periodo 1992-2008 respecto al periodo 1975-1991 ........................................ 39 4.1. ModeloBYM...................................... 46 4.2. Mapa de Espa˜na sin modificar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.3. Mapa de Espa˜na con las islas Canarias trasladadas . . . . . . . . . . . . . . . . . 50 4.4. Riesgos suavizados para el periodo 1975-1991 . . . . . . . . . . . . . . . . . . . . 61 4.5. Probabilidad de los riesgos estimados del periodo 1975-1991 sean superiores a 100 63 4.6. Funci´on de densidad de la RMEs del periodo 1975-1991 . . . . . . . . . . . . . . 64 4.7. Riesgos suavizados para el periodo 1992-2008 . . . . . . . . . . . . . . . . . . . . 66 4.8. Probabilidad de que los riesgos estimados del periodo 1992-2008 sean superiores a100........................................... 67 4.9. Funci´on de densidad de la RMEs del periodo 1992-2008 . . . . . . . . . . . . . . 68 4.10. Riesgos suavizados para el periodo 1992-2008 seg´un las tasas del periodo 1975-1991 71 4.11. Probabilidad de que los riesgos estimados del periodo 1992-2008 seg´un tasas del periodo 1975-1991 sean superiores a 100 . . . . . . . . . . . . . . . . . . . . . . . 78 ix Cap´ıtulo 1 Informaci´on general sobre el c´ancer de pr´ostata Antes de profundizar en detalles matem´aticos veamos qu´e es el c´ancer de pr´ostata y c´omo se detecta. 1.1. ¿ Qu´e es el c´ancer de pr´ostata? El c´ancer de pr´ostata es una enfermedad en la cual se forman c´elulas malignas en los tejidos de la pr´ostata. ´ Esta es una gl´andula que pertenece al sistema reproductor masculino localizada justo por debajo de la vejiga y por delante del recto, tal y como se muestra en la Figura (1.1). Su tama˜no es como el de una nuez y rodea una parte de la uretra. La gl´andula prost´atica produce un fluido que forma parte del semen. El c´ancer de pr´ostata es el tercer tumor m´as frecuente en varones espa˜noles y supone la tercera causa de muerte por c´ancer en Espa˜na [4]. Por ejemplo, en el a˜no 2005 exist´ıa una tasa cruda de 25.89 casos por 100000 habitantes, en dicho a˜no fallecieron 5500 personas por este c´ancer [5]. Este c´ancer constituye aproximadamente el 11 % de todas las neoplasias (c´anceres) en los varones de Europa, y es el responsable del 9 % de las muertes por c´ancer entre los hombres dentro de la Uni´on Europea [6]. Su incidencia (n´umero de casos nuevos de c´ancer que se diagnostican cada a˜no) es muy variable en las distintas zonas del mundo pudiendo tener esto relaci´on con distintos factores como son: exposici´on ambiental, dieta, estilo de vida y factores gen´eticos [7]. La probabilidad de padecer un c´ancer de pr´ostata se incrementa con la edad y en el 90 % de 7 Paula Navarro Esteban Figura 1.1: Anatom´ıa del sistema reproductor masculino los casos se produce en hombres de m´as de 65 a˜nos. Esto es debido a que a medida que los hombres envejecen, la pr´ostata se puede agrandar y bloquear la vejiga o la uretra tal y como se muestra en la Figura (1.2). Esto puede ocasionar dificultades para orinar o interferir con la funci´on sexual. La afecci´on se llama hiperplasia prost´atica benigna (HPB). Figura 1.2: Pr´ostata normal y HPB En la raza negra existe un riesgo mayor de padecerlo y tambi´en se incrementa de forma importante en personas con parientes de primer grado afectados pudiendo llegar a ser siete veces superior en personas con dos o tres familiares de primer grado afectados [8]. La incidencia del c´ancer de pr´ostata est´a aumentando en pa´ıses desarrollados entre otras causas por el aumento de la esperanza de vida y principalmente por el desarrollo de nuevas t´ecnicas 8 An´alisis espacial de riesgos de mortalidad diagn´osticas, ver apartado 1.2. Desde el punto de vista anatomo-patol´ogico la mayor parte de los c´anceres de pr´ostata son adenocarcinomas, es decir, aparecen en las c´elulas glandulares que revisten ´organos internos siendo frecuentemente de localizaci´on multifocal y de predominio en la zona perif´erica de la gl´andula prost´atica [9]. Algunos de los s´ıntomas ocasionados por el c´ancer de la pr´ostata pueden ser: - Disminuci´on del calibre o interrupci´on del flujo urinario. - Aumento de la frecuencia de la micci´on (especialmente por la noche). - Dificultad para orinar. - Dolor o ardor durante la micci´on. - Presencia de sangre en la orina o en el semen. - Dolor constante en la espalda, las caderas o la pelvis. - Eyaculaci´on dolorosa. 1.2. ¿C´omo se diagnostica? Seg´un [4] se pueden utilizar las siguientes pruebas y procedimientos: -Examen digital del recto (EDR): examen del recto. El m´edico o el enfermero inserta un dedo protegido por un guante lubricado en el recto y palpa la pr´ostata a trav´es de la pared del recto en busca de bultos o ´areas anormales. -Prueba del ant´ıgeno prost´atico espec´ıfico (APE): prueba de laboratorio que mide las concentraciones del APE en la sangre. El APE es una sustancia elaborada por la pr´ostata que se puede encontrar en una mayor cantidad en la sangre de los hombres que tienen c´ancer de pr´ostata. La concentraci´on de APE tambi´en puede ser elevada en los hombres que sufren una infecci´on o una inflamaci´on de la pr´ostata, o que tienen HPB (pr´ostata agrandada, pero no cancerosa). Esta prueba junto con la anterior son las m´as utilizadas. -Ecograf´ıa transrectal: procedimiento por el cual se inserta en el recto una sonda que tiene aproximadamente el tama˜no de un dedo para examinar la pr´ostata. La sonda se utiliza para hacer rebotar ultrasonidos (ondas de sonido de alta energ´ıa) en los tejidos internos de la pr´ostata y crear ecos. Los ecos forman una imagen de los tejidos corporales que se llama sonograma. La ecograf´ıa transrectal se puede usar durante una biopsia. 9 Paula Navarro Esteban -Biopsia: extracci´on de c´elulas o tejidos realizada por un pat´ologo para observarlos bajo un microscopio. El pat´ologo examina la muestra en busca de c´elulas cancerosas. Hay dos tipos de biopsia utilizados para diagnosticar el c´ancer de pr´ostata. -Biopsia transrectal: se extrae tejido de la pr´ostata mediante la inserci´on de una aguja fina a trav´es del recto hasta la pr´ostata. Este procedimiento se suele realizar mediante ecograf´ıa transrectal para ayudar a guiar la aguja. -Biopsia transperineal: se extrae una muestra de tejido prost´atico mediante la inserci´on de una aguja fina a trav´es de la piel entre el escroto y el recto hasta la pr´ostata. 1.3. Factores de riesgo Un factor de riesgo es toda circunstancia o situaci´on que afecta la posibilidad de tener una enfermedad, en nuestro caso el c´ancer de pr´ostata. Hay factores de riesgo que se pueden modificar, como beber alcohol, y otros, como la edad, que no se pueden cambiar. El tenerlos no es una condici´on necesaria ni suficiente para padecer la enfermedad, es decir, hay personas con uno o m´as factores de riesgo que nunca padecen c´ancer, mientras que otras que lo padecen puede que hayan tenido pocos factores de riesgo conocidos o ninguno [10]. Todav´ıa no se entiende completamente las causas del c´ancer de pr´ostata, pero se han encontrado varios factores que pueden cambiar el riesgo de padecer esta enfermedad. Son los siguientes: -Edad. La probabilidad de tener c´ancer de pr´ostata aumenta r´apidamente despu´es de los 50 a˜nos, se conocen muy pocos casos antes de los 40. Casi dos de cada tres casos de c´ancer de pr´ostata se detectan en hombres mayores de 65 a˜nos. -Raza. Los hombres de raza negra es el grupo ´etnico que con m´as frecuencia padece c´ancer de pr´ostata respecto a otras razas. Adem´as ´estos tienen una mayor probabilidad de ser diagnosticados en una etapa avanzada, y tienen m´as del doble de probabilidad de morir de c´ancer de pr´ostata en comparaci´on con los hombres blancos. Al contrario, los hombres asi´atico-americanos y los hispanos o latinos tienen una menor frecuencia de padecer c´ancer de pr´ostata respecto a los hombres blancos. El porqu´e de estas diferencias no est´a claro todav´ıa. -Nacionalidad. El c´ancer de pr´ostata es m´as com´un en Norteam´erica y en la regi´on noroeste de Europa, Australia, y en las islas del Caribe. Es menos com´un en Asia, Centroam´erica y Sudam´erica. Estas diferencias tampoco est´an claras, es probable que se deba a los distintos estilos de vida o a que en los pa´ıses desarrollados se usan m´as pruebas de detecci´on. -Antecedentes familiares. Se sabe que se duplica el riesgo de que un hombre padezca el c´ancer si su padre o su hermano lo padecen tambi´en (siendo mayor el riesgo para los hombres que 10 An´alisis espacial de riesgos de mortalidad tienen un hermano con la enfermedad que para aquellos con un padre afectado por este c´ancer). Asimismo, el riesgo es mucho mayor en el caso de los hombres que tienen varios familiares afectados, especialmente si tales familiares eran j´ovenes en el momento en el que se les diagnostic´o el c´ancer. -Otras causas. Para factores como la alimentaci´on, la obesidad, el tabaquismo, la inflamaci´on de la pr´ostata, las infecciones de transmisi´on sexual o la vasectom´ıa, se han realizado estudios pero no se han encontrado conclusiones claras. 11 Cap´ıtulo 2 Creaci´on de la base de datos En este cap´ıtulo crearemos las bases de datos a partir del archivo “prostata.txt” (cedido por el Centro Nacional de Epidemiolog´ıa), que recoge los casos de mortalidad por c´ancer de pr´ostata en cada provincia de Espa˜na desde el a˜no 1975 hasta 2008. Para nuestro caso vamos a considerar dos periodos de estudio, 1975-1991 y 1992-2008. En primer lugar calcularemos los valores o muertes esperadas (mediante estandarizaci´on indirecta) seg´un la mortalidad por c´ancer de pr´ostata en hombres, describiendo la sintaxis del software R. Tambi´en calcularemos las tasas de mortalidad de referencia o tasa global de mortalidad en Espa˜na, y finalmente obtendremos la Raz´on de Mortalidad Estandarizada (RMEi) o “Standardized Mortality Ratio” (SMRi), en ingl´es, que definiremos matem´aticamente m´as adelante. Para entender los c´alculos mencionados vamos a definir primeramente una serie de conceptos. 2.1. Definiciones b´asicas 2.1.1. Salud p´ublica La salud p´ublica es la disciplina encargada de la protecci´on de la salud a nivel poblacional. De ella se encargan especialistas como m´edicos, bi´ologos, estad´ısticos, veterinarios, etc, aunque su desarrollo depende de los gobiernos. Su objetivo es mejorar las condiciones de salud de las comunidades mediante la promoci´on de estilos de vida saludables, las campa˜nas de concienciaci´on, la educaci´on y la investigaci´on. Los organismos de la salud p´ublica deben evaluar las necesidades de salud de la poblaci´on, los riesgos para la salud y analizar las causas de dichos riesgos. De acuerdo a lo detectado, 13 Paula Navarro Esteban deben establecer las prioridades y desarrollar programas y planes que permitan responder a las necesidades de la poblaci´on. La salud p´ublica tambi´en debe gestionar de manera ´optima los recursos para asegurar que sus servicios lleguen a la mayor cantidad de gente posible, ya que parte de unos principios comunitarios y no personales. 2.1.2. Registro de c´ancer Una de las mayores necesidades de la salud p´ublica es conocer, lo mejor posible, las caracter´ısticas de salud y enfermedad de la poblaci´on. Esta informaci´on es un instrumento fundamental para la planificaci´on de las actuaciones en el ´ambito de la salud en general y de la salud p´ublica en particular. Para ello se ha creado el sistema de registro de mortalidad, que es una fuente de datos exhaustiva, representativa tanto a nivel nacional como internacional, y que ha sido recopilada durante un largo periodo de tiempo. Todas las muertes ocurridas en una zona o pa´ıs quedan recogidas en este sistema, junto con informaci´on de la persona, el tiempo, el lugar y la causa de la defunci´on. Gracias a este registro se pueden construir indicadores de mortalidad que pueden ayudar a detectar problemas de salud importantes as´ı como desigualdades entre regiones. No obstante, el objetivo m´as efectivo de los programas de salud p´ublica es prevenir la causa que da origen a todos los dem´as trastornos o afecciones que conducen a la muerte. Con esa base se dise˜n´o el Modelo Internacional de Certificado M´edico de Causa de Defunci´on, en el que est´a basado el dise˜no del Bolet´ın Estad´ıstico de Defunci´on (B.E.D.). Fue en 1967 cuando la Asamblea Mundial de la Salud estableci´o que las causas de defunci´on (definidas como “todas aquellas enfermedades, estados morbosos o lesiones que produjeron la muerte o contribuyeron a ella y las circunstancias del accidente o de la violencia que produjo dichas lesiones”) fueran registradas en el certificado m´edico de causas de defunci´on. M´as tarde se implant´o la Clasificaci´on Internacional de Enfermedades (CIE), por la cual las causas de muerte eran codificadas de una forma homog´enea, una vez inscritas. De este modo se obten´ıa una fuente de datos universal. 2.2. Instalaci´on y carga de librer´ıas Debemos definir un directorio de trabajo (en este caso una carpeta llamada “tfm”en la unidad D:/) donde se guardar´an las bases de datos que creemos con ayuda del software R. > setwd("D:/Paula/master/TFM/tfm") 14 An´alisis espacial de riesgos de mortalidad Tras su instalaci´on en R, cargamos las librer´ıas necesarias con las siguientes instrucciones. (Nota: se ha escrito a continuaci´on de cada librer´ıa una breve descripci´on.) > library(plyr) # plyr es un conjunto de herramientas cuya finalidad es resolver > # un problema. Consiste en dividir un problema grande en varios > # m´as peque~nos y manejables, para luego volverlos a unir y poder > # asi resolver el problema principal. > > library(reshape) # reshape facilita la reestructuraci´on y la agregaci´on de > # datos usando solamente dos funciones: melt y cast. > > library(INLA) # INLA contiene funciones que permiten trabajar con modelos > # aditivos completamente bayesianos usando Integrated Nested > # Laplace Approximation. > > library(maptools) # maptools es un conjunto de herramientas que manipulan y > # leen datos geogr´aficos, en particular archivos ESRI. > > library(spdep) # spdep es una colecci´on de funciones que sirven para crear > # matrices espaciales a partir de pol´ıgonos contiguos y de > # funciones que estiman intervalos y errores de modelos > # espaciales simult´aneos autorregresivos (SAR) y de modelos > # espaciales de regresi´on (CAR), entre otros. > > library(RColorBrewer) # RColorBrewer proporciona paletas para colorear mapas > # de acuerdo a una variable. > > library(Hmisc) # Hmisc contiene numerosas funciones ´utiles para el an´alisis de > # datos, gr´aficos de alto nivel, funciones para el c´alculo del > # tama~no de una muestra, importar conjuntos de datos,etc. > > library(gplots) # gplots consta de varias herramientas para la representaci´on > # gr´afica de los datos. > > library(pixmap) # pixmap contiene funciones para importar, exportar, > # representar gr´aficamente y otras manipulaciones de im´agenes > # pixeladas. > > library(mvtnorm) # mvtnorm calcula probabilidades de normales multivariantes > # y de t de student, as´ı como cuantiles y densidades. 15 Paula Navarro Esteban 1975-1991 poblaci´on casos edad1 23783116 3 edad2 26650153 0 edad3 28289365 3 edad4 27812155 10 edad5 25733004 14 edad6 23405963 11 edad7 21409388 12 edad8 19591498 23 edad9 18897184 46 edad10 18474795 170 edad11 18123387 527 edad12 16715265 1317 edad13 14107813 3276 edad14 11457767 6392 edad15 8747924 10793 edad16 5989069 13892 edad17 3335444 12847 edad18 1713562 9326 1992-2008 poblaci´on casos edad1 17585302 0 edad2 17808776 0 edad3 19972259 0 edad4 23486777 3 edad5 27140860 5 edad6 29658756 3 edad7 29397065 8 edad8 27443972 11 edad9 24997207 47 edad10 22559250 154 edad11 20002762 568 edad12 18190253 1459 edad13 16964722 3555 edad14 15528149 7560 edad15 13273003 13217 edad16 9479828 18670 edad17 5670309 21230 edad18 3641927 25464 Tabla 2.1: N´umero de muertes y poblaci´on total seg´un grupos de edad para cada periodo > pob.h.p1 <- pob.h[pob.h$periodo==1,] #separamos el primer periodo > pob.h.p2 <- pob.h[pob.h$periodo==2,] #separamos el segundo periodo > #calculamos la poblaci´on total por grupos de edad > hombres.ce.p1<-data.frame(poblacion=apply(pob.h.p1[,3:20],2,sum),hombres.ce.p1) > hombres.ce.p2<-data.frame(poblacion=apply(pob.h.p2[,3:20],2,sum),hombres.ce.p2) Por lo tanto para el periodo 1 y 2 respectivamente obtenemos la Tabla 2.1. A la vista de la Tabla 2.1 podemos comprobar que para cada periodo se verifica que aproximadamente el 90 % de los casos de c´ancer de pr´ostata se produce en hombres mayores de 65 a˜nos, es decir, a partir del grupo edad14, tal y como hab´ıamos mencionado en el cap´ıtulo 1. Concretamente para el primer periodo el n´umero de casos de c´ancer de pr´ostata en hombres mayores de 65 a˜nos representa un 90.77 % del total de casos en ese periodo y para el segundo periodo los hombres mayores de 65 a˜nos que padecen c´ancer representan un 93.68 % del total de este mismo periodo. Calculamos las tasas espec´ıficas por edad del pa´ıs dividiendo los muertos por c´ancer en cierto grupo de edad por la poblaci´on en el grupo de edad. 22 An´alisis espacial de riesgos de mortalidad > tasa.h.ce.p1 <- hombres.ce.p1[2]/hombres.ce.p1$poblacion > tasa.h.ce.p2 <- hombres.ce.p2[2]/hombres.ce.p2$poblacion 2.4.2. Esperados El n´umero de muertes esperadas estar´an contenidas en un vector de dimensi´on 18. Para calcular esta cantidad por estandarizaci´on indirecta debemos multiplicar la poblaci´on de cada provincia por edad por las tasas de referencia en ese grupo de edad. Esto equivale a multiplicar la matriz de poblaciones, que contiene tantas filas como provincias y tantas columnas como grupos de edad por el vector de tasas calculado en el apartado anterior, esto es: Eij =X j nijRj donde Rij es la tasa espec´ıfica por edad y nij es la poblaci´on en la provincia idel grupo de edad j. Hay que tener en cuenta que habr´a 4 vectores de este tipo: los esperados para el periodo 1 tomando las tasas de referencia del periodo 1, los esperados para el periodo 1 tomando las tasas de referencia del periodo 2, los esperados para el periodo 2 tomando las tasas de referencia del periodo 1 y los esperados para el periodo 2 tomando las tasas de referencia del periodo 2. En la sintaxis que se muestra de R, los esperados tomando como referencia el periodo 1 est´an guardados en el objeto “esperados.p1” y los esperados tomando como referencia el periodo 2 en el objeto “esperados.p2”. Se puede acceder al periodo en que se esperan a trav´es de la variable periodo. Pasamos a ver c´omo se calculan los esperados. Vamos a ordenar las poblaciones por periodo y provincia para que nos quede as´ı en los esperados. > pob.h <- pob.h[order(pob.h$periodo, pob.h$SECCION),] Multiplicamos las tasas espec´ıficas por las poblaciones, con esto calculamos el denominador de las RMEs, es decir, los esperados. > #Esperados tomando como referencia las tasas del primer periodo > esperados.p1 <- data.matrix(pob.h[,3:20])%*%data.matrix(tasa.h.ce.p1) > esperados.p1 <- data.frame(periodo=pob.h$periodo,SECCION=pob.h$SECCION, + esperados.p1) > #Esperados tomando como referencia las tasas del segundo periodo > esperados.p2 <- data.matrix(pob.h[,3:20])%*%data.matrix(tasa.h.ce.p2) > esperados.p2 <- data.frame(periodo=pob.h$periodo,SECCION=pob.h$SECCION, + esperados.p2) 23 Paula Navarro Esteban Si pidi´eramos los valores extremos de los esperados obtendr´ıamos: Tomando como referencia el primer periodo, el m´ınimo de los esperados es: 301.9957 Tomando como referencia el primer periodo, el m´aximo de los esperados es: 10804.71 Tomando como referencia el segundo periodo, el m´ınimo de los esperados es: 294.8144 Tomando como referencia el segundo periodo, el m´aximo de los esperados es: 10421.61 Luego los esperados son muy variables, con rangos 301.9957-10804.71, tomando como referencia el primer periodo, y 294.8144-10421.61, tomando como referencia el segundo periodo. Esta variabilidad se ver´a reflejada cuando calculemos las SMRs. Obtendremos SMRs altas en las provincias escasamente pobladas o en provincias donde el c´ancer de pr´ostata sea una enfermedad muy rara y haya muy pocos casos. Por eso tenemos que usar m´etodos de suavizado de los riesgos de mortalidad que usen la informaci´on de todas las provincias y as´ı las estimaciones ser´an m´as fiables y los estimadores m´as estables. 2.4.3. Observados Los observados estar´an contenidos en una vector de igual estructura que los esperados pero en cada elemento deberemos tener los muertos observados por secci´on (provincia). Vamos a calcular los observados en cada uno de los periodos. Creamos la variable “secci´on” (an´alogamente a lo anterior) y la incluimos en la base de datos de “mort”. > seccion <- mort$prov > mort.h <- data.frame(mort[9], SECCION=seccion, mort[1:7],mort[10:27]) > #Notar que aqui no estamos cogiendo el a~no Sumamos las muertes por secci´on (provincia) y periodo. > mort.agr<-aggregate(data.frame(casos=mort.h$casos), + list(periodo=mort.h$periodo,SECCION=mort.h$SECCION), + FUN=sum) 24 An´alisis espacial de riesgos de mortalidad Las ordenamos y las renombramos por “observados”. > observados <- mort.agr[order(mort.agr$periodo,mort.agr$SECCION),] Guardamos todos los objetos que hemos ido creando. > hombres_ce <- list(observados=observados, esperados.p1=esperados.p1, + esperados.p2=esperados.p2) > save(hombres_ce,file="espana_interna.RData") 25 Cap´ıtulo 3 Medidas cl´asicas de estimaci´on de riesgos: raz´on de mortalidad estandarizada La aproximaci´on cl´asica asume que para el estudio de datos de mortalidad, el n´umero de muertes en cada provincia yi, sigue una distribuci´on de Poisson. Estas distribuciones son independientes y vienen definidas como: yi∼ P(λi=eiri),para i= 1,...,50 (3.1) donde eies el n´umero esperado de muertes en la provincia i-´esima (que es conocido) y ries el riesgo relativo asociado a la provincia i-´esima (que es desconocido). Para estimar el par´ametro rique aparece en la distribuci´on de probabilidad yi, usaremos el m´etodo de estimaci´on por m´axima verosimilitud que elige como estimador a aqu´el valor del par´ametro que tiene la propiedad de maximizar la probabilidad de la muestra aleatoria observada [13]. Definici´on 3.0.3. Sea Xuna variable aleatoria discreta con funci´on de probabilidad denotada por p(X;θ), es decir, que depende del par´ametro θy de la m.a.s. de tama˜no n X = (x1, . . . , xn). La funci´on de probabilidad conjunta viene dada por p(X;θ) = n Y i=1 p(xi;θ). Cuando θes conocido, esta funci´on determina la probabilidad de aparici´on de cada muestra, sin embargo en un problema de estimaci´on se conoce el valor de la muestra pero no el del par´ametro. Luego se puede ver esta funci´on de probabilidad como una funci´on del par´ametro desconocido. A esta nueva funci´on se le llama funci´on de verosimilitud y se denota por L(θ|X), donde L(θ|X) = p(X;θ). 27 Paula Navarro Esteban Es decir, que la verosimilitud es una funci´on que describe la dependencia de un/os par´ametro/s en los valores de la muestra. Al valor del par´ametro que maximiza esta funci´on se le llama estimador de m´axima verosimilitud o estimador MV y lo denotaremos por ˆ θ. En nuestro caso el par´ametro desconocido es riy denotaremos a su estimador MV como ˆri. Como los casos observados, yi, siguen una distribuci´on de Poisson de par´ametro λi=eiri, su funci´on de probabilidad ser´a: P(Yi=yi) = e−λiλyi i yi!,para i= 1,...,50. Adem´as, suponemos que los yi’s son independientes y por tanto la funci´on de probabilidad conjunta ser´a el producto de las funciones de probabilidad marginales, es decir, L(ri|yi) = 50 Y i=1 e−λiλyi i yi!= 50 Y i=1 e−eiri(eiri)yi yi!,(3.2) donde hemos tomado n= 50, ya que es el n´umero de provincias que vamos a considerar. Como el logaritmo es una funci´on mon´otona creciente, maximizar Les lo mismo que maximizar log(L), y esto ´ultimo es m´as sencillo. Por lo tanto, aplicando las propiedades del logaritmo a (3.2) obtenemos: log L(ri|yi) = 50 X i=1 (−eiri+yilog eiri−log yi!) . Ahora bien, los posibles candidatos a estimadores MV son los valores que resuelven ∂log L ∂ri = 0. Es decir, ∂log L ∂ri =−ri+yiei eiri =−ri+yi ri = 0. Por tanto, yi=riei⇔ˆri=yi ei . Es importante notar que el ˆricalculado es candidato a estimador MV, ya que el que la primera derivada sea igual a 0 es s´olo una condici´on necesaria pero no suficiente para obtener un m´aximo. Por lo tanto vamos a calcular la segunda derivada: ∂2log L ∂r2 iri=ˆri =−e2 i y2 i <0. Por lo tanto ˆri=yi eies el estimador MV de rique se denomina SMRi. 28 An´alisis espacial de riesgos de mortalidad Veamos cu´al es la varianza de este estimador MV. Para ello tomamos varianzas en la igualdad anterior usando que eies constante y las propiedades de la varianza : V ar(ˆri) = V ar yi ei=V ar(yi) e2 i . Como las yisiguen una distribuci´on de Poisson, su varianza es igual a la media, en este caso a riei, luego: V ar(ˆri) = riei e2 i =ri ei . La V ar(ˆri) depende de ri, luego un estimador de la varianza es: d V ar(ˆri) = ˆri ei =yi e2 i .(3.3) La varianza de este estimador ser´a grande si eies peque˜no, por tanto para enfermedades muy raras o ´areas peque˜nas la SMR podr´ıa ser muy inestable. Por ello se consideran modelos m´as sofisticados que suavizan los riesgos en el sentido de que los hacen m´as estables, es decir, menos variables. 3.1. C´alculo de la SMR del c´ancer de pr´ostata en Espa˜na A continuaci´on vamos a determinar la SMR del c´ancer de pr´ostata en Espa˜na. Como ya hemos comentado en el cap´ıtulo dos, esta raz´on es el cociente entre las muertes observadas y las esperadas en cada provincia. Frecuentemente las SMR suelen aparecer multiplicadas por 100, y as´ı lo vamos a considerar a continuaci´on. Para su c´alculo y representaci´on vamos a tener en cuenta los esperados seg´un cada periodo de tiempo. La explicaci´on de c´omo se realiza el mapa de Espa˜na en Rse muestra en el cap´ıtulo 4. 3.1.1. C´alculo de la SMR del primer periodo respecto del primer periodo Seleccionamos los datos del primer periodo y calculamos la SMR correspondiente. > observados1<- observados[observados$periodo==1,] > esperados1<- esperados.p1[esperados.p1$periodo==1,] > SMRperiodo1<-observados1$casos/esperados1$total*100 Calculamos los intervalos de confianza del riesgo relativo al 95 % (que es el m´as habitual). IC95 =ˆri−1.96qd V ar(ˆri),ˆri+ 1.96qd V ar(ˆri) Para calcular la d V ar(ˆri) usamos la expresi´on (3.3). En Rse calcular´ıa como sigue: 29 Paula Navarro Esteban > var.periodo1 <-observados1$casos/(esperados1$total)^2 #varianza estimada > l.i.periodo1 <-SMRperiodo1-1.96*(var.periodo1)^(1/2) > #extremo inferior del intervalo de confianza > u.i.periodo1 <-SMRperiodo1+1.96*(var.periodo1)^(1/2) > #extremo superior del intervalo de confianza Recogemos la SMR y los intervalos de confianza en la Tablas (3.1) y (3.2), en las cuales podemos obervar que no hay ning´un intervalo de confianza que contenga al 100, no obstante hay algunos que est´an muy pr´oximos a 100, como es el caso de la provincia de La Coru˜na y de Segovia. En estas tablas se han escrito con rojo las provincias que tienen un exceso de riesgo mayor que el global de Espa˜na, es decir, son provincias que tienen una SMR>100 y cuyo intervalo de confianza no contiene al 100. Las provincias que presentan un defecto de riesgo menor que el global de Espa˜na, es decir, provincias con un SMR<100 y cuyo intervalo de confianza no contiene al 100, se han escrito con verde. En la Figura (3.1) se muestran las SMR calculadas para el primer periodo tomando como valores esperados los valores del primer periodo. 3.1.2. C´alculo de la SMR del segundo periodo respecto del segundo periodo Seleccionamos los datos del segundo periodo y calculamos la SMR correspondiente. > observados2<- observados[observados$periodo==2,] > esperados2<- esperados.p2[esperados.p2$periodo==2,] > SMRperiodo2<-observados2$casos/esperados2$total*100 Calculamos los intervalos de confianza del riesgo relativo al 95 %, para ello calculamos la d V ar(ˆri) usando la expresi´on (3.3). > var.periodo2 <-observados2$casos/(esperados2$total)^2 #varianza estimada > l.i.periodo2 <-SMRperiodo2-1.96*(var.periodo2)^(1/2) > #extremo inferior del intervalo de confianza > u.i.periodo2 <-SMRperiodo2+1.96*(var.periodo2)^(1/2) > #extremo superior del intervalo de confianza Recogemos la SMR y los intervalos de confianza para el periodo 1992-2008 en las Tablas (3.3) y (3.4), en los cuales podemos obervar qu´e intervalos contienen al 100 y cu´ales no. An´alogamente al periodo 1, se han escrito con rojo las provincias que tienen un exceso de riesgo mayor que el global de Espa˜na, es decir, son provincias que tienen una SMR>100 y cuyo intervalo de confianza 30 An´alisis espacial de riesgos de mortalidad periodo 1975-1991 Provincia SMR Extremo inferior del IC Extremo superior del IC 1´ Alava 92.02732 91.92318 92.13146 2Albacete 110.44717 110.36488 110.52947 3Alicante 108.21751 108.16860 108.26641 4Almer´ıa 94.07112 93.99299 94.14924 5´ Avila 77.47897 77.39882 77.55912 6Badajoz 106.66730 106.60720 106.72740 7Baleares 106.66730 115.30053 115.42398 8Barcelona 106.98268 106.95751 107.00784 9Burgos 82.65539 82.58920 82.72159 10 C´aceres 89.90307 89.83803 89.96811 11 C´adiz 111.22778 111.16298 111.29257 12 Castell´on 114.90356 114.83234 114.97477 13 Ciudad Real 102.64547 102.57750 102.71345 14 C´ordoba 84.74963 84.69627 84.80299 15 La Coru˜na 100.08295 100.03636 100.12953 16 Cuenca 78.47246 78.39963 78.54529 17 Gerona 112.44738 112.37658 112.51818 18 Granada 86.55596 86.50169 86.61024 19 Guadalajara 81.88871 81.80049 81.97693 20 Guip´uzcoa 87.32163 87.25971 87.38355 21 Huelva 109.73550 109.65302 109.81798 22 Huesca 101.55544 101.47279 101.63809 23 Ja´en 79.34486 79.29164 79.39807 24 Le´on 96.61224 96.55350 96.67099 25 L´erida 99.03057 98.96212 99.09903 Tabla 3.1: SMR e intervalos de confianza al 95 % para el periodo 1975-1991. Parte 1. 31 Paula Navarro Esteban periodo 1992-2008 respecto al periodo 1975-1991 Provincia SMR Extremo inferior del IC Extremo superior del IC 26 La Rioja 109.33051 109.25744 109.40358 27 Lugo 101.75593 101.70501 101.80685 28 Madrid 90.44334 90.42491 90.46177 29 M´alaga 88.82532 88.78755 88.86308 30 Murcia 96.89576 96.85559 96.93593 31 Navarra 100.72587 100.67422 100.77752 32 Orense 109.80886 109.75372 109.86399 33 Princip. de Asturias 107.60661 107.56990 107.64332 34 Palencia 98.59275 98.51080 98.67470 35 Las Palmas 125.42646 125.36618 125.48674 36 Pontevedra 115.54346 115.49632 115.59060 37 Salamanca 104.13039 104.07296 104.18782 38 Sta Cruz de Tenerife 100.51233 100.46191 100.56274 39 Cantabria 103.45516 103.40184 103.50849 40 Segovia 98.54571 98.46158 98.629855 41 Sevilla 94.39442 94.35961 94.42924 42 Soria 85.33804 85.24760 85.42849 43 Tarragona 91.56986 91.52312 91.61659 44 Teruel 97.07215 96.99295 97.15136 45 Toledo 93.93642 93.88803 93.98480 46 Valencia 100.16652 100.13838 100.19466 47 Valladolid 107.23896 107.18068 107.29723 48 Vizcaya 98.40796 98.37004 98.44589 49 Zamora 100.23226 100.16519 100.29933 50 Zaragoza 105.76047 105.71949 105.80146 Tabla 3.6: SMR e intervalos de confianza al 95 % para el periodo 1992-2008 respecto al periodo 1975-1991. Parte 2. 38 An´alisis espacial de riesgos de mortalidad En la Figura (3.3) se muestran las SMR calculadas para el segundo periodo tomando como valores esperados los valores del primer periodo. RMS para el periodo 1992−2008 respecto al periodo 1975−1991 [ 70, 80) [ 80, 90) [ 90,100) [100,110) [110,120) [120,130] Figura 3.3: Raz´on de mortalidad estandarizada para el periodo 1992-2008 respecto al periodo 1975-1991 39 Cap´ıtulo 4 Aplicaci´on del modelo BYM con INLA Los modelos m´as utilizados son los modelos jer´arquicos bayesianos, concretamente el propuesto por Besag, York y Mollie (BYM) [12] el cual supone la principal referencia en m´etodos de suavizaci´on de riesgos e incluye dos efectos aleatorios: uno recoge la dependencia espacial (se sabe que regiones vecinas presentan riesgos parecidos debido a que tienen estilos de vida similares, as´ı como mismo tipo de alimentaci´on o igual contaminaci´on atmosf´erica, por ejemplo) y el otro recoge la dependencia no estructurada. Dicho modelo es estimado por el m´etodo INLA. La estimaci´on obtenida mediante el modelo BYM corresponde a una suavizaci´on de los riesgos de mortalidad que ser´a m´as pronunciada en aquellas ´areas con menor poblaci´on. Este modelo establece una ponderaci´on entre dos fuentes de variaci´on para obtener una estimaci´on del riesgo relativo de cada ´area. La primera de estas fuentes de variaci´on, de estructura espacial, comparte informaci´on entre unidades vecinas para modelizar la dependencia geogr´afica de los riesgos. Consideraremos que un ´area es vecina de otra si comparte frontera. Este efecto hace posible que regiones vecinas tengan estimaciones similares y de esta forma el riesgo var´ıe geogr´aficamente de forma suave. El segundo de los efectos, espacialmente heterog´eneo, toma valores independientes en todas las unidades geogr´aficas lo que permite que localizaciones vecinas presenten riesgos diferentes. En el caso de que el ´area peque˜na tenga una poblaci´on de gran tama˜no tendr´a mayor peso la informaci´on proporcionada por este ´area; en cambio, si presenta una poblaci´on de reducido tama˜no tendr´a mayor peso la informaci´on del res- to de ´areas (o ´areas vecinas). Mediante la ponderaci´on de ambos tipos de informaci´on el modelo BYM minimiza el problema relativo a la estabilidad de los riesgos estimados. A la estimaci´on del riesgo obtenida a partir del modelo BYM se la denominar´a de aqu´ı en adelante Raz´on de Mortalidad Estandarizada Suavizada (RMEs) [32]. 41 Paula Navarro Esteban En este cap´ıtulo vamos a exponer las ideas b´asicas de la modelizaci´on bayesiana antes de aplicar el modelo de BYM. 4.1. Descripci´on del modelo BYM 4.1.1. Estad´ıstica bayesiana Fundamentalmente la estad´ıstica bayesiana permite que los par´ametros de cualquier funci´on de probabilidad tengan distribuciones, es decir, sean variables aleatorias. A estas distribuciones se les llama distribuciones a priori. Debido a esto, los par´ametros de las distribuciones a priori mencionadas pueden ser tambi´en variables aleatorias, de aqu´ı que se establezca una jerarqu´ıa. Estos modelos jer´arquicos son la base de la inferencia bayesiana. Combinando la verosimilitud con las distribuciones a priori de los par´ametros, se forma una distribuci´on llamada distribuci´on a posteriori que describe el comportamiento de los par´ametros despu´es de haber visto los datos. Resumiendo, la distribuci´on a priori es la idea que tenemos de la variaci´on de los datos. Los datos se modelizan con ayuda de la funci´on de verosimilitud y a trav´es de la distribuci´on a posteriori actualizamos el conocimiento que tenemos de la variaci´on de los datos. Adem´as esta distribuci´on a posteriori puede convertirse en la distribuci´on a priori de los par´ametros antes de experimentar con los siguientes datos. Como ventajas de la estad´ıstica bayesiana podemos decir que permite la formulaci´on de modelos m´as complejos y flexibles que la estad´ıstica frecuentista y que el software bayesiano ha alcanzando en la actualidad un grado de madurez que permite la utilizaci´on cotidiana de esta metodolog´ıa. No obstante, la estad´ıstica bayesiana es criticada por los estad´ısticos frecuentistas por conciderarla subjetiva ya que dicen que la estad´ıstica bayesiana se encuentra en el observador, mientras que la estad´ıstica cl´asica o frecuentista es un concepto objetivo, que se encuentra en la naturaleza. De este modo, en estad´ıstica cl´asica s´olo se toma como fuente de informaci´on las muestras finitas obtenidas. En el caso bayesiano, sin embargo, adem´as de la muestra tambi´en juega un papel fundamental la informaci´on previa de los fen´omenos que se tratan de modelizar [16]. Veamos las diferencias entre ambas: - En estad´ıstica frecuentista la probabilidad de un suceso es la proporci´on de veces que se dar´ıa dicho suceso si se repitiera indefinidamente un experimento. La probabilidad es objetiva y es igual para todos los observadores, mientras que para un bayesiano la probabilidad de un suceso es una medida de la credibilidad de los posibles valores que puede tomar dicho 42 An´alisis espacial de riesgos de mortalidad suceso. La probabilidad es subjetiva y depende del observador [17]. - En estad´ıstica cl´asica los par´ametros de un modelo estad´ıstico se suponen constantes que habremos de determinar. Sin embargo en estad´ıstica bayesiana los par´ametros de un modelo estad´ıstico se suponen variables aleatorias cuyas distribuci´on de probabilidad habremos de determinar. - La inferencia de la estad´ıstica frecuentista se basa en la funci´on de verosimilitud, la cual mide la credibilidad de los valores de los par´ametros dados los datos que hemos observado: Verosim(Par´ametros)=P(datos|Par´ametros). La principal herramienta de la estad´ıstica bayesiana es el teorema de Bayes que actualiza la informaci´on, a priori, de los par´ametros con la informaci´on proporcionada por los datos. π(θ|datos) = π(datos|θ)π(θ) π(datos)∝π(datos|θ)π(θ), donde π(datos|θ) es la funci´on de verosimilitud, que ya hemos definido anteriormente, π(θ) es la distribuci´on a priori y π(θ|datos) es la distribuci´on a posteriori que vamos a definir a continuaci´on. Definici´on 4.1.1. Sean ylos datos observados y θel/los par´ametro/s del modelo. Sea π(θ)la distribuci´on a priori y π(y|θ)la verosimilitud. Se define la distribuci´on a posteriori como π(θ|y) = π(y|θ)π(θ)/C, donde C=Zπ(y|θ)π(θ)dθ (4.1) Notaci´on: π(θ|y)∝π(y|θ)π(θ). Observaciones. - En la inferencia bayesiana, las distribuciones a posteriori aparecen continuamente. Por ejemplo, los cuantiles o momentos se pueden escribir en funci´on de la media a posteriori de funciones de θ, siendo la media a posteriori de una funci´on f(θ): E[f(θ|y)] = Rf(θ)π(θ)π(y|θ)dθ Rπ(θ)π(y|θ)dθ . No obstante, el c´alculo de estas integrales ha sido uno de los principales problemas pr´acticos en la inferencia bayesiana, sobre todo en dimensiones altas. Como dicho c´alculo es imposible realizarlo de forma exacta, se aproxima num´ericamente mediante aproximaciones de Laplace o integraci´on de Monte Carlo, como veremos m´as adelante [21]. - La verosimilitud informa sobre los par´ametros a trav´es de los datos, mientras que las distribuciones a priori informan a trav´es de las asunciones, por eso cuando haya una gran cantidad de datos, la verosimilitud contribuir´a en mayor medida a la estimaci´on del riesgo. Sin embargo cuando la muestra de datos sea peque˜na, la/s distribuci´on/es a priori dominar´an el an´alisis. 43 Paula Navarro Esteban 4.1.2. Modelo de BYM Vamos a escribir el modelo de BYM bas´andonos en [12] y [15]. Sean xi= log(ri) el logaritmo del riesgo relativo, que es desconocido, en la provincia i-´esima para i= 1,...,50 e yiel n´umero de casos observados de muertes durante el periodo del estudio. Cuando la enfermedad no es contagiosa y rara, se suele suponer que los yi’s son variables de Poisson independientes con medias eiri, donde eies el n´umero de casos esperados en la provincia i-´esima suponiendo riesgo constante, es decir, basado s´olo en la tasa de incidencia global y en la poblaci´on de riesgo en la provincia i-´esima ajustada por edad. Definici´on 4.1.2. El modelo de BYM se define en dos etapas: yi|ri∼P(eiri) log(ri) = t+ui+vi,para i= 1,...,50,(4.2) donde yies el n´umero de casos observados de muertes durante el periodo del estudio, ries el riesgo relativo en la provincia i-´esima, eies el n´umero de casos esperados en la provincia i-´esima, tes un t´ermino est´andar asociado a las covariables medidas que son conocidas o sospechadas para ser relevantes en la enfermedad que normalmente tiene la forma de un modelo lineal t=Aθ con al menos Adesconocido, uiyvison t´erminos adicionales que se pueden interpretar como sustitutos para variables desconocidas o no observadas. En particular, uimuestra la estructura espacial en la que dos ´areas que son vecinas tendr´ıan valores m´as “parecidos” que cualquier par arbitrario de ´areas que tom´asemos, mientras que virepresenta las variabilidad no estructurada, es decir, como un “ruido” que recoge los efectos aleatorios no espaciales (llamado heterogenidad). La segunda etapa de este modelo tambi´en se puede expresar como: ri=et+ui+vi,para i= 1,...,50. Por simplicidad y sin p´erdida de informaci´on podemos ignorar ty denotaremos a uipor uy avipor v. Hasta el momento no tenemos mas informaci´on, as´ı que podemos suponer que uyvson independientes y que v∼ N(0, σ2 v) siendo σ2 vsu varianza que es desconocida. Para use elige una funci´on de densidad de la siguiente familia: p(u)∝exp  −X i<j ωijφ(ui−uj) con u∈R, donde ωij son los pesos no negativos tales que ωij = 0 excepto si iyjest´an en ´areas vecinas y φes una funci´on par de zcreciente con |z|. 44 An´alisis espacial de riesgos de mortalidad Por lo tanto la densidad condicional de uies: pi(ui| · · · )∝exp  −X j∈∂i ωijφ(ui−uj) ,(4.3) donde ∂idenota a las ´areas vecinas de i. Por simplicidad tomamos ωij = 1 y φ(z) = z2/(2σ2 u) con kconstante positiva desconocida. As´ı (4.3) se transforma en: p(u|σ2 u)∝1 (σ2 u)n/2exp  −1 2σ2 uX i∼j (ui−uj)2 , donde i∼jindica que iyjson vecinos. Notar que esta aproximaci´on especifica la dependencia espacial como un promedio del efecto de sus ´areas vecinas limitando la vecindad a sus ´areas contiguas. Esto recibe el nombre de autorregresi´on intr´ınseca gaussiana y una de sus ventajas es que los momentos condicionales est´an definidos como funciones simples de los “valores vecinos” y el n´umero de vecinos. E(ui|. . .) = ¯ui=X i∼j ui/ni, V ar(ui|. . .) = σ2 u/ni siendo nila cardinalidad de ∂i, es decir, el n´umero de vecinos de i. Como resultado final se obtienen para cada provincia, la media a posteriori de la distribuci´on del riesgo as´ı como la probabilidad de que ´esta sea mayor de 100. Para calcular las distribuciones a posteriori no se usa la definici´on previa sino que, en la pr´actica, se aproximan con los m´etodos de Monte Carlo. En nuestro caso, la densidad a posteriori de u,v,σ2 uyσ2 ven la cual basamos las inferencias es: P(u, v, σ2 u, σ2 v|y)∝ n Y i=1 (exp(−ciexi)(ciexi)yi/yi!) ×σ2 u −n/2exp  −1 2σ2 uX i∼j (ui−uj)2  ×σ2 v −n/2exp −1 2σ2 v n X i=1 v2 i!×prior(σ2 u, σ2 v) donde prior(σ2 u, σ2 v) es la densidad a priori de los dos hiperpar´ametros. Para evitar singularidades y lograr que la distribuci´on sea correcta realizamos la siguiente modificaci´on: prior(σ2 u, σ2 v)∝exp−/2σ2 uexp−/2σ2 v,con σ2 u, σ2 v>0 45 Paula Navarro Esteban donde es una constante peque˜na positiva que tomamos =0.001. Las densidades condicionales de uiyvino se ven modificadas, mientras que las de σ2 uyσ2 v permanecen en la familia de distribuciones gamma-inversas. Es decir, σ2 v∼Gamma−1(0.001,0.001); σ2 u∼Gamma−1(0.001,0.001). Un resumen de este modelo se recoge en la Figura (4.1), en la cual Eirepresenta los casos observados, en vez de ei. Figura 4.1: Modelo BYM En la pr´actica una vez obtenidas las estimaciones a posteriori de los par´ametros, se estiman los riesgos de cada provincia a trav´es de la siguiente expresi´on: RME = (exp(ui+vi)) ∗100 De esta forma se excluye al par´ametro tde la f´ormula (4.2) de la estimaci´on, por lo que la RME compara el riesgo de cada provincia con el nivel medio de Espa˜na. En esta secci´on hablaremos tambi´en de la matriz de vecindades ya que el modelo BYM suaviza 46 An´alisis espacial de riesgos de mortalidad los riesgos prestando informaci´on entre regiones vecinas en lugar de hacer uso de todas las unidades, independientemente de su ubicaci´on geogr´afica. 4.2. Lectura de datos y cartograf´ıa En el siguiente paso importaremos la cartograf´ıa necesaria para calcular la matriz de vecindades y representar gr´aficamente los mapas. Veamos primero la importancia de los mapas en epidemiolog´ıa. 4.2.1. Importancia de los mapas Tradicionalmente, los mapas o atlas se han utilizado para representar la geograf´ıa de un pa´ıs o de una regi´on. Pueden ser pol´ıticos, f´ısicos, econ´omicos,. . . En cambio, los mapas o atlas que vamos a usar, ilustrar´an las diferencias de mortalidad que existen entre las distintas provincias espa˜nolas con respecto a Espa˜na. Este trabajo se encamina en el ´ambito de la representaci´on cartogr´afica de enfermedades o “disease mapping” (en ingl´es) [30]. En este campo se suelen aplicar modelos a nivel de ´area, donde un ´area puede ser una secci´on censal, una provincia (como es nuestro caso), un pa´ıs, etc. Generalmente se asume que el n´umero de eventos (casos de mortalidad, por ejemplo) localizados en cada una de estas ´areas sigue una distribuci´on de Poisson con media igual al riesgo de mortalidad (ri) por el n´umero de casos esperados (ei). Actualmente el riesgo se estima mediante modelos jer´arquicos bayesianos, siendo el modelo de BYM [12] uno de los m´as utilizados para el estudio de la distribuci´on geogr´afica de riesgos en ´areas peque˜nas. Una vez realizada la estimaci´on, los valores del riesgo estimado se representan en un mapa, que permite identificar las zonas geogr´aficas con mayor riesgo de mortalidad [22]. Los primeros mapas sanitarios que se realizaron fueron a finales del siglo XVIII y principios del XIX, cuando la fiebre amarilla era uno de los mayores problemas en Salud P´ublica de Estados Unidos. Para conocer la causa de esta enfermedad se represent´o el mapa de Nueva York junto con los casos detectados de esta enfermedad. Se realizaron m´as mapas como ´este en los a˜nos sucesivos y gracias a ellos se detect´o el patr´on geogr´afico de la fiebre amarilla donde se mostraba que esta enfermedad era m´as com´un en las zonas en donde las temperaturas eran altas, ten´ıan una humedad elevada y hab´ıa una gran concentraci´on de insectos en el aire [22]. Desde entonces se ha seguido realizando este tipo de mapas con resultados exitosos, como el descubrimiento del origen del c´olera en Londres a mital del siglo XIX, cuando se observ´o (gracias a los mapas que dibuj´o el doctor Snow) que la falta de drenaje en un determinado distrito era la causa de esta epidemia. 47 Paula Navarro Esteban converjan (bajo ciertas condiciones de regularidad) a la distribuci´on de inter´es π(θ|y). De este modo, estos m´etodos permiten muestrear la distribuci´on a posteriori, aunque ´esta sea desconocida, gracias a la construcci´on de una cadena de Markov, mediante simulaci´on Monte Carlo, cuya distribuci´on estacionaria sea, precisamente π(θ|y). “MCMC es, esencialmente, integraci´on Monte Carlo, haciendo correr por largo tiempo una inteligentemente construida cadena de Markov .”[21] Precisando a´un m´as, vamos a introducir MCMC como un m´etodo para evaluar expresiones del siguiente tipo: E[f(X)] = Rf(x)π(x)dx Rπ(x)dx , donde Xun vector de kvariables aleatorias con distribuci´on π(.), en nuestro caso Xes el vector de los par´ametros del modelo y π(.) es su distribuci´on a posteriori. Para ello, vamos a describir las dos partes que componen este m´etodo: integraci´on Monte Carlo y Cadenas de Markov [21]. Cadena de Markov Una cadena de Markov es una sucesi´on de variables aleatorias, {X1, X2, . . . , Xt, . . .}tal que ∀t≥0, P(Xt+1|Xt, Xt−1, . . . , X1) = P(Xt+1|Xt). Es decir, dado Xt, el siguiente estado Xt+1 no depende de la historia de la cadena {X1, . . . , Xt−1}. Integraci´on Monte Carlo Es un m´etodo que mediante modelos matem´aticos intenta reproducir el comportamiento aleatorio de sistemas reales no din´amicos. Veamos c´omo lo logra. La integraci´on Monte Carlo eval´ua E[f(X)] extrayendo muestras {Xt, t = 1, . . . , n}de π(.) y entonces aproxima E[f(X)] ≈1 n n X t=1 f(Xt). Es decir, que la media poblacional de f(X) se estima a trav´es de la media muestral. Si las muestras {Xt}son independientes, por la ley de los grandes n´umeros podemos asegurar que las aproximaciones pueden ser tan precisas como queramos, basta con aumentar el tama˜no de la muestra, n. Notar que aqu´ı el valor de nlo decide el analista. 54 An´alisis espacial de riesgos de mortalidad En general, no es factible extraer muestras {Xt}independientes, sin embargo no es necesario que las muestras sean independientes. Las {Xt}pueden ser generadas por cualquier proceso, el cual extraiga muestras a trav´es de un conjunto de puntos de π(.), donde ´esta no se anule, en las proporciones adecuadas. Una forma de hacer esto es a trav´es de las cadenas de Markov, teniendo π(.) como su distribuci´on estacionaria. Esto es “Markov Chain Monte Carlo”. Observaci´on. Ley de los grandes n´umeros. Sea (Ω, F, P) un espacio de probabilidad. Sean X1, X2, . . . , Xnvariables aleatorias independientes e id´enticamente distribuidas definidas en dicho espacio de probabilidad y con media µ,E[Xi] = µ. Entonces se tiene: Pω∈Ω : l´ım n→+∞ X1+· · · Xn n=µ= 1. Gaussian Markov random field (GMRF) AGaussian Markov random field (GMRF),x={xi|i∈ν}, es un vector aleatorio normal |ν|- dimensional que satisface propiedades de independencia condicional de Markov. Para simplificar podemos asumir que ν={1, . . . , n}. Las propiedades de independencia condicional se pueden representar usando un grafo no dirigido G= (ν, ) con νnodos y arcos. Diremos que xiexjson condicionalmente independientes si y s´olo si i, j /∈, es decir, si no existe ning´un arco entre los nodos xiexj. Entonces diremos que xes un GMRF con respecto a G. Cuando {i, j} ∈ diremos que iyjson vecinos [24]. Los GMRFs son tambi´en conocidos en ingl´es como “conditional auto-regressions (CARs)” siguiendo el trabajo de Besag [23]. Esta familia de modelos cubre un buen n´umero de modelos jer´arquicos bayesianos, incluyendo la mayor´ıa de los usados en estad´ıstica espacial. Adem´as, las propiedades de Markov pueden ser usadas para modelizar la dependencia local. Modelo de variable latente Un modelo de variable latente es un modelo estad´ıstico que relaciona un conjunto de variables (las llamadas variables manifiestas) con un conjunto de varibles latentes (que son las variables que no son emp´ıricamente observables). INLA Recientemente se ha implementado el m´etodo conocido en ingl´es como Integrated Nested Laplace Approximation(INLA) para estimar modelos del tipo“Latent Gaussian Markov Random Fields”, en ingl´es [25]. Entre estos modelos que puede estimar INLA se encuentra el modelo BYM, 55 Paula Navarro Esteban para el cual proporciona una buena aproximaci´on a las distribuciones posteriores marginales de los par´ametros en el modelo. Entendemos por distribuciones posteriores marginales, denot´andolas por π(θi|y), a aquellas distribuciones a posteriori que dependen de alguno de los par´ametros del modelo. En nuestro caso, como queremos estimar el riesgo relativo en diferentes provincias, no nos interesa calcular la distribuci´on a posteriori conjunta sino que es suficiente con calcular las marginales. De esta forma es m´as f´acil hacer las aproximaciones o c´alculos necesarios que si tom´asemos la conjunta. Adem´as de proporcionar una aproximaci´on de las marginales posteriores de los par´ametros de un modelo, INLA puede calcular tambi´en criterios para la selecci´on del modelo, como el DIC (Deviance Information Criterion, en ingl´es), que es el an´alogo al AIC de la estad´ıstica frecuentista, o el histograma PIT (Probability Integral Transform, en ingl´es), que sirve para detectar valores at´ıpicos y para comparar y validar los modelos. A continuaci´on vamos a explicar qu´e aproximaciones usa INLA. En primer lugar, supongamos que hemos observado nvariables yi, i = 1, . . . , n con una distribuci´on (en nuestro caso Poisson) con media λique est´a relacionada con un predictor lineal ηia tav´es de una funci´on. A su vez, ηise modeliza como sigue: ηi=α+ nf X j=1 f(j)(uji) + nβ X k=1 βkzki +i, donde f(j)representa una funci´on no lineal o efecto aleatorio (estamos suponiendo que hay nf)en un conjunto de covariables u,βkson los coeficientes para los efectos lineales de un vector de covariables zyison los t´erminos no estructurados. Los efectos latentes x={{ηi}, α, {βk}, . . .} se suponen normales con media 0 y matriz de precisi´on Q(θ1), donde θ1es vector de hiperpar´ametros. Por lo tanto, las observaciones tendr´an una verosimilitud que depender´a de los efectos latentes xy de los par´ametros θ2. De aqu´ı que las obervaciones yiser´an supuestas independientes dados xyθ2. Este tipo de modelos GMRF son los que puede estimar INLA. En particular, para los modelos espaciales (como es nuestro caso), las f(j)(uji) representan los efectos aleatorios en una localizaci´on espacial i. Se supone que el cuerpo latente tiene propiedades de independencia condicionales, lo que convierte a xen un GMRF. Recordar que el objetivo de INLA es calcular las distribuciones a posteriori marginales. ´ Estas se pueden escribir como: π(xi|y)∝Zπ(xi|θ, y)π(θ|y)dθ 56 An´alisis espacial de riesgos de mortalidad y π(θj|y)∝Zπ(θ|y)dθ−j, donde θ−jdenota θmenos la componente j. Las aproximaciones se har´an de las integrales de los miembros de la derecha de las anteriores expresiones. Dichas aproximaciones ser´an posibles si la dimensi´on de θes peque˜na, pero esto es lo que suele suceder en la pr´actica. Notar que tambi´en se necesita una aproximaci´on de π(θ|y). De π(x, θ, y) = π(x|θ, y)π(θ|y)π(y) se sigue que la marginal posterior π(θ|y) de los hiperpar´ametros se puede obtener usando la aproximaci´on de Laplace: ˜π(θ|y)∝π(x, θ, y) ˜πG(x|θ, y)x=x∗(θ) , donde ˜πG(x|θ, y) es la aproximaci´on gaussiana (es decir, aproximamos la distribuci´on de una variable no normal a una gaussiana uniendo la moda y la curvatura de la moda [24]) a la condicional de xyx∗(θ) es la moda de la condicional dado un valor de θ. Por lo tanto, las marginales que nos interesan se pueden calcular usando la siguiente suma finita evaluada en los puntos θk: ˜π(xi|y) = X k ˜π(xi|θk, y)טπ(θk|y)×∆k, donde ∆krepresenta los pesos para cada vector de los valores de θken la red, ˜π(xi|θ, y) y ˜π(θ|y) representan las aproximaciones de π(xi|θ, y) y π(θ|y), respectivamente [14]. Los puntos θkse pueden elegir de dos maneras: usando la estrategia “GRID” o bien usando “CCD” (Central Composite Design). La primera se recomienda si la dimensi´on del vector de los hiperpar´ametros, h, es media (h= 6/12), y lo que hace es para definir una red de puntos que cubre el ´area donde est´a la mayor parte de la masa de ˜π(θk|y). La segunda estrategia, CCD, tiene su origen en el campo del dise˜no y consiste en colocar una peque˜na cantidad de “puntos” en un espacio m-dimensional para estimar la curvatura de ˜π(θk|y). Para m´as informaci´on ver [25]. Dicho art´ıculo recomienda usar la estrategia CCD para problemas cuyo vector de hiperpar´ametros tiene una dimensi´on hgrande, siendo h > 12, aproximadamente. Para calcular la aproximaci´on ˜π(xi|θk, y) existen tres formas posibles: aproximaciones gaussianas, totales de Laplace y simplificadas de Laplace. Cada una de ellas presenta caracter´ısticas diferentes, as´ı como tiempos de c´alculo y precisi´on. Aunque la aproximaci´on de Gauss es la m´as r´apida computacionalmente, puede generar errores en la localizaci´on de la media posterior y/o errores debido a la falta de asimetr´ıa [24]. Por otra parte, la aproximaci´on de Laplace es la m´as precisa, pero es muy costosa computacionalmente. Por ´ultimo, la aproximaci´on simplificada de Laplace es la m´as r´apida para calcular y por lo general es lo suficientemente precisa. Para concluir, INLA es una alternativa a los m´etodos MCMC (Markov Chain Monte Carlo) computacionalmente m´as r´apida que proporciona resultados satisfactorios en un corto periodo 57 Paula Navarro Esteban de tiempo computacional [14]. Esto se debe principalmente a que hace uso de aproximaciones de Laplace para obtener las distribuciones marginales a posteriori y utiliza computaci´on en paralelo. Veamos ahora la aplicaci´on de INLA en nuestro caso particular. En primer lugar, realizaremos un an´alisis descriptivo donde estudiaremos cada periodo de estudio (1975-1991 y 1992-2008) de forma independiente, es decir, consideraremos como si tuvi´eramos dos estudios de dise˜no transversal distintos, uno para cada periodo. 4.3.2. An´alisis del primer periodo (estandarizaci´on interna con tasas espec´ıficas de referencia del primer periodo) Vamos a obtener las RME suavizadas para el periodo 1975-1991 y estandarizaremos seg´un las tasas espec´ıficas por edad del mismo periodo. Como las tasas espec´ıficas se han calculado teniendo en cuenta toda Espa˜na (estandarizaci´on interna), las RMEs tendr´an como referencia la mortalidad de todo el pa´ıs. Creamos la matriz de datos necesaria para ejecutar el modelo. Seleccionamos los datos del primer periodo (1975-1991): > obs <- observados[observados$periodo==1,] > esp <- esperados.p1[esperados.p1$periodo==1,] Para simplificar la notaci´on, asignamos los valores observados y esperados a los objetos O y E, respectivamente. Adem´as creamos el objeto “region” que necesitaremos para el an´alisis: > O <- obs$casos > E <- esp$total > region <- 1:length(x.nb) Finalmente, creamos un objeto “data.frame” con todos los datos necesarios, para ejecutar el primer modelo: > Datos <- data.frame(region=region, region.struct=region, O=O, E=E) Definimos la f´ormula del modelo (formula1) y ejecutamos la funci´on “inla” que estima el modelo BYM, comentado anteriormente. As´ı, obtendremos los riesgos suavizados para provincia. Tambi´en obtendremos la probabilidad de que dichos riesgos sea mayor que 100 (probabilidad de exceso de riesgo) para cada provincia. Todos los resultados de este modelo se guardar´an en el objeto “resultado.p1”. 58 An´alisis espacial de riesgos de mortalidad > formula1 <- O ~ f(region.struct,model="besag",graph=esp_prov_nb, + param=c(1,0.001))+f(region,model="iid", param=c(1,0.01)) > resultado.p1<-inla(formula1, family="poisson", data=Datos, E=E, + control.compute=list(dic=T, cpo=TRUE), + control.predictor=list(compute=TRUE,cdf=c(log(1))), + control.inla=list(strategy="laplace",int.strategy="grid")) En graph hemos escrito la matriz de vecindad. El efecto de suavizado es especificado en “f()”, siendo model=“iid” para un efecto aleatorio i.i.d. y “besag” para un IGMRF (Intrinsic Gaussian Random Field). Luego estamos ante un modelo BYM ya que tenemos la suma de un“iid”m´as un “besag”. Adem´as con hyper=c(1,0.001) indicamos los par´ametros a y b de la distribuci´on a priori (LogGamma) del hiperpar´ametro (Log-Precision).Por defecto: a=1, b=0.001. Como dic=TRUE ycpo=TRUE, se calcular´an el DIC y las medidas de validaci´on CPO y PIT, respectivamente. Con compute=TRUE se calculan las marginales del predictor lineal log(λi) y con cdf=c(log(1)) se calcula para el predictor lineal la“cumulative distribution function”(cdf), es decir, se calculan las probabilidades Prob(log(λi)< x(p)). En este caso x(p) = log(1). Con strategy=laplace se utiliza la estrategia para obtener las aproximaciones de las distribuciones marginales a posteriori, como es “laplace” tiene mayor precisi´on aunque tambi´en m´as coste computacional que si hubi´eramos usado “gaussian” o “simplified.laplace”. Por ´ultimo con “int.strategy” indicamos la estrategia de integraci´on utilizada, puede ser “ccd”, “grid” (m´as precisa pero con m´as coste computacional) y “eb” (empirical bayes). Si quisi´eramos ver todas las componentes que podemos obtener con inla basta con escribir: > names(resultado.p1) [1] "names.fixed" "summary.fixed" [3] "marginals.fixed" "summary.lincomb" [5] "marginals.lincomb" "size.lincomb" [7] "summary.lincomb.derived" "marginals.lincomb.derived" [9] "size.lincomb.derived" "mlik" [11] "cpo" "model.random" [13] "summary.random" "marginals.random" [15] "size.random" "summary.linear.predictor" [17] "marginals.linear.predictor" "summary.fitted.values" [19] "marginals.fitted.values" "size.linear.predictor" [21] "summary.hyperpar" "marginals.hyperpar" [23] "internal.summary.hyperpar" "internal.marginals.hyperpar" [25] "si" "total.offset" [27] "model.spde2.blc" "summary.spde2.blc" 59 Paula Navarro Esteban [29] "marginals.spde2.blc" "size.spde2.blc" [31] "misc" "dic" [33] "mode" "neffp" [35] "joint.hyper" "nhyper" [37] "version" "Q" [39] "graph" "cpu.used" [41] ".args" "call" [43] "model.matrix" La ventaja de este modelo es que permite obtener las distribuciones a posteriori marginales de la suma del efecto espacial y el heterog´eneo. Los riesgos suavizados de cada ´area las obtenemos mediante la siguiente sintaxis: > RMEs.p1 <- resultado.p1$summary.fitted.values$mean*100 Es decir, que las RMEs se muestran con su valor multiplicado por 100. Por consiguiente, aquellas provincias con valores de RMEs superiores a 100 indican un exceso de riesgo de mortalidad respecto a la mortalidad de la poblaci´on espa˜nola, mientras que valores inferiores evidencian un menor riesgo. Para obtener la probabilidad de que los riesgos suavizados sea superior a 100 para cada ´area ejecutaremos la siguiente sintaxis: > PRP.p1 <- 1-resultado.p1$summary.linear.predictor[,"0 cdf"] Representaremos los riesgos suavizados mediante mapas de cortes fijos donde los colores verdes representan las ´areas con mayor defecto de mortalidad, mientras que las ´areas con mayor exceso de mortalidad son representadas en tonos marrones. Podemos observar que tras la suavizaci´on, las zonas despobladas no presentan comportamientos an´omalos y que las RME ahora no reflejan la distribuci´on de la poblaci´on en la regi´on de estudio sino el riesgo estudiado. En la Figura (4.4) se muestran los riesgos suavizados de mortalidad por c´ancer de pr´ostata para Espa˜na (concretamente, la media de las distribuciones). Como ya hemos mencionado anteriormente observamos las provincias cuyos riesgos estimados est´en por encima de 100, lo que indicar´a que tienen un exceso de riesgo de mortalidad respecto a la mortalidad de la poblaci´on espa˜nola, y las provincias que presenten riesgos menores de 100, lo que se traducir´a como un menor riesgo respecto a la mortalidad de la poblaci´on espa˜nola. En este caso se observa que se dan m´as casos de mortalidad de c´ancer de pr´ostata, comparando con la media de Espa˜na, en 60 An´alisis espacial de riesgos de mortalidad Riesgos suavizados del periodo 1975−1991 [ 80, 90) [ 90,100) [100,110) [110,120) [120,130) [130,140] Figura 4.4: Riesgos suavizados para el periodo 1975-1991 las Palmas, seguido por Valladolid, Baleares, Santa Cruz de Tenerife y en algunas zonas de la costa como Castell´on, Gerona, Vizcaya, C´adiz y M´alaga. Las provincias de Pontevedra, Asturias, Cantabria, Palencia, La Rioja, Navarra, Huesca, Zaragoza, Barcelona, Tarragona, Valencia, Alicante, Albacete, Ciudad Real, Badajoz y Huelva tienen un riesgo muy poco por encima de 100 (entre 100 y 110), luego tienen un mayor riesgo que en Espa˜na en global. En cambio puede observarse un patr´on dominado por las provincias ubicadas en el interior (´ Avila, Burgos, Soria, Madrid, Guadalajara, Cuenca, Toledo, Teruel, C´ordoba, Ja´en, etc) donde los riesgos son menores que en Espa˜na en su totalidad, ya que recordar que el valor 100 indica que el riesgo de la provincia es como el de Espa˜na. Si comparamos la Figura (4.4) con la Figura (3.1), realizada en el cap´ıtulo 3, vemos que pr´acticamente los mapas que representan la SMR y los riesgos suavizados son muy parecidos, aunque se encuentran diferencias en las provincias de Soria (42), ´ Avila (5), Cuenca (16) y Ja´en (23) que tienen una SMR entre 70 y 80, mientras que sus riesgos suavizados var´ıan entre 80 y 61 Paula Navarro Esteban 90. Cantabria tiene una SMR entre 110 y 120, mientras que sus riesgos suavizados var´ıan entre 100 y 110.Esto se debe a que estas provincias no son muy pobladas o tienen pocos casos de mortalidad por c´ancer. En efecto, si calculamos con Rel n´umero de observados y esperados para estas provincias en el primer periodo obtenemos: ´ Avila: n º de observados: 359 , n º de esperados es: 463.3516 Poblaci´on total de ´ Avila en el a~no 1991: 87608 Cuenca: n º de observados: 446 , n º de esperados es: 568.3522 Poblaci´on total de Cuenca en el a~no 1991: 102195 Soria: n º de observados: 235 , n º de esperados es: 301.9957 Poblaci´on total de Soria en el a~no 1991: 47149 Ja´en: n º de observados: 854 , n º de esperados es: 1076.314 Poblaci´on total de Ja´en en el a~no 1991: 315168 Con el objeto de cuantificar la evidencia estad´ıstica que proporcionan las estimaciones del riesgo en cada provincia crearemos mapas de la probabilidad de exceso de riesgo (RMEs>100). Para representar estas probabilidades en el mapa consideraremos la siguiente segmentaci´on: (0,1; 0,2; 0,8 y 0,9). Utilizaremos tonalidades verdes para los riesgos relativos con baja probabilidad de ser superior a 100, rojas para los riesgos relativos con alta probabilidad (´areas con exceso de riesgo, es decir, provincias donde el riesgo es mayor que el riesgo global de Espa˜na ) y el rango intermedio lo representaremos en amarillo. Consideraremos que una probabilidad es alta cuando sea mayor que 0.8. El mapa representado en la Figura (4.5) representa las probabilidades de las RMEs en el periodo 1975-1991. Veamos en ´el, qu´e probabilidades tienen las provincias que en el mapa representado en la Figura (4.4) presentaban un riesgo estimado por encima de 100 (que eran entre otras: las Palmas, Santa Cruz de Tenerife, Valladolid, Baleares, Castell´on, Gerona, Vizcaya, C´adiz, M´alaga, Pontevedra, Cantabria, Asturias, Palencia, Zaragoza, Barcelona, Valencia, Alicante, Albacete, Huelva y Badajoz). Observamos que todas estas provincias tienen una probabidad superior a 0.9, lo que significa que con una alta probabilidad el riesgo de mortalidad por c´ancer de pr´ostata de estas provincias est´a por encima del riesgo global de Espa˜na. Adem´as con alta probabilidad, entre 0.8 y 0.9, tenemos tambi´en Navarra. 62 An´alisis espacial de riesgos de mortalidad Probabilidad de que los riesgos estimados sean superiores a 100 Periodo 1975−1991 [0.0,0.1) [0.1,0.2) [0.2,0.8) [0.8,0.9) [0.9,1.0] Figura 4.5: Probabilidad de los riesgos estimados del periodo 1975-1991 sean superiores a 100 Para observar la forma de la distribuci´on de la RMEs representaremos su funci´on de densidad en la Figura (4.6), ya que bajo la filosof´ıa Bayesiana no se obtiene una estimaci´on puntual del riesgo sino toda la distribuci´on de probabilidad. En su interior se muestran los colores correspondientes a los grupos de valores utilizados en la representaci´on geogr´afica de las RMEs. La funci´on de densidad permite apreciar con detalle el rango de valores de la RMEs, el cual no es visible en el mapa. La forma de la funci´on de densidad de las RMEs es como la de una distribuci´on normal de media en torno a 115 pero con dos m´aximos, uno local y otro absoluto. Tambi´en vemos que la cola derecha de esta funci´on es m´as larga que la de la izquierda. 63 Paula Navarro Esteban Los riesgos suavizados de cada ´area las obtenemos mediante la siguiente sintaxis: > RMEs.p2.resp.p1 <- resultado.p2.resp.p1$summary.fitted.values$mean*100 Los riesgos suavizados del segundo periodo (1992-2008) teniendo en cuenta las tasas del primer periodo se representan en la Figura (4.10). En ella se representan los mapas de cortes fijos donde los colores verdes representan las ´areas con mayor defecto de mortalidad, mientras que los tonos marrones representan las ´areas con mayor exceso. Recordemos que los valores de los riesgos suavizados est´an en referencia a la media de la mortalidad en todo el pa´ıs para el primer periodo (1975-1991). > Cortes.fijos <- cut2(RMEs.p1, g=7, onlycuts=T) > Cortes.fijos[1] <- min(RMEs.p1,RMEs.p2.resp.p1) > Cortes.fijos[length(Cortes.fijos)] <- max(RMEs.p1,RMEs.p2.resp.p1) 70 An´alisis espacial de riesgos de mortalidad Riesgos suavizados del periodo 1992−2008 según tasas del periodo 1975−1991 [ 70, 80) [ 80, 90) [ 90,100) [100,110) [110,120) [120,130] Figura 4.10: Riesgos suavizados para el periodo 1992-2008 seg´un las tasas del periodo 1975-1991 Notar que en la Figura (4.10) donde se han representado los riesgos suavizados del segundo periodo (1992-2008) seg´un las tasas del primer periodo, se han fijado los cortes seg´un los cortes fijos de los riesgos suavizados del primer periodo, ver Figura (4.4). De este modo se puede comparar la tendencia en la distribuci´on de los riesgos estimados entre los dos periodos analizados. Se utiliza la misma gama de colores divergente que la utilizada en la representaci´on geogr´afica de los riesgos suavizados de las figuras anteriores. En esta figura observamos que Las Palmas, Huesca, Castell´on y Pontevedra presentan el mayor exceso de riesgo de mortalidad (concretamente mayor de 110), mientras que las provincias de Lugo, Orense, La Coru˜na, Asturias, Cantabria, Navarra, La Rioja, Zaragoza, L´erida, Gerona, Baleares, Valencia, Albacete, Valladolid, Zamora, Salamanca y Santa Cruz de Tenerife tienen tambi´en un exceso de riesgo pero en menor medidad (concretamente entre 100 y 110). Si comparamos la Figura (4.10) con la Figura (3.3) observamos que son exactamente iguales. Esto puede deberse al aumento de la poblaci´on o al descenso del n´umero de casos de mortalidad 71 Paula Navarro Esteban por c´ancer de pr´ostata. Veamos, por ejemplo, como ha variado la poblaci´on de cada provincia. Para ello mostraremos como ejemplo la poblaci´on en el a˜no 1991 y en el a˜no 2008 en las Tablas 4.1 y 4.2. Estas tablas se han realizado usando la siguiente orden en R: > poblaciontotal<-aggregate(mort$pob,by=list(provincia=mort$provincia, + agno=mort$agno),FUN=sum,na.rm=T) > casostotal<-aggregate(mort$casos,by=list(provincia=mort$provincia, + agno=mort$agno),FUN=sum,na.rm=T) En estas tablas aparecenen azul las provincias que han aumentado su poblaci´on, como puede verse la mayor´ıa de ellas la ha incrementado. A continuaci´on se han realizado las tablas 4.3 y 4.4, donde se muestran los casos de mortalidad por c´ancer en el a˜no 1991 y en el a˜no 2008 por provincias. Se han se˜nalado con azul las provincias cuya mortalidad por c´ancer ha disminuido. Observamos que el descenso se ha producido en un n´umero muy peque˜no de provincias. No obstante, en general las provincias que en el a˜no 2008 tienen m´as casos de mortalidad que en el a˜no 1991 han aumentado su poblaci´on de forma m´as significativa que el n´umero de casos de muertes. Cabe destacar tambi´en, que las provincias de ´ Avila, Le´on y Lugo presentan en el a˜no 2008 menos casos de mortalidad que en 1991, asi como n´umero de habitantes. En cambio, Asturias tiene menor n´umero de habitantes (ha pasado de 527875 a 506429 habitantes) y mayor n´umero de casos (en 1991 ten´ıa 133 y en 2008 ten´ıa 196). Eso mismo les ha pasado a las provincias de Orense, Palencia, Salamanca, Soria, Vizcaya y Zamora. Sin embrago, a Barcelona, Navarra y Valladolid les ha ocurrido lo contrario, han aumentado su poblaci´on pero han logrado reducir el n´umero de muertes por c´ancer. 72 An´alisis espacial de riesgos de mortalidad provincia pob. a˜no 1991 pob. a˜no 2008 ´ Alava 135846 152290 Albacete 170574 197900 Alicante 635766 934825 Almer´ıa 227068 348428 Asturias 527875 506429 ´ Avila 87608 85107 Badajoz 321359 333706 Barcelona 2267372 2619950 Burgos 177000 184488 C´aceres 203664 201717 C´adiz 537052 597355 Cantabria 257997 280573 Castell´on 221145 293266 Ciudad Real 232734 256289 C´ordoba 369279 385648 Coru˜na 528875 538275 Cuenca 102195 108062 Girona 253635 364375 Granada 388961 443619 Guadalajara 74240 118748 Guip´uzcoa 332616 339503 Huelva 219204 248844 Huesca 104774 112492 Islas Baleares 351142 530773 Ja´en 315168 324254 Tabla 4.1: Poblaci´on total seg´un provincias y periodos. Parte 1. 73 Paula Navarro Esteban provincia pob. a˜no 1991 pob. a˜no 2008 Las Palmas 385546 535774 Le´on 257647 235557 Lleida 176165 215834 Lugo 188022 168060 Madrid 2389806 3031261 M´alaga 571728 764375 Murcia 517273 725827 Navarra 258524 304285 Ourense 170094 157489 Palencia 92003 84218 Pontevedra 430940 454793 Rioja 130613 157673 Salamanca 174739 169100 Segovia 73647 81527 Sevilla 797993 903410 Soria 47149 47142 Sta. Cruz de Tenerife 358889 495613 Tarragona 270012 396502 Teruel 72134 74997 Toledo 245004 327341 Valencia 1034610 1237337 Valladolid 243236 255239 Vizcaya 564772 551484 Zamora 105214 96242 Zaragoza 410467 463864 Tabla 4.2: Poblaci´on total seg´un provincias y periodos. Parte 2. 74 An´alisis espacial de riesgos de mortalidad provincia casos a˜no 1991 casos a˜no 2008 ´ Alava 21 39 Albacete 39 65 Alicante 137 181 Almer´ıa 40 52 Asturias 133 196 ´ Avila 33 29 Badajoz 75 92 Barcelona 583 524 Burgos 46 69 C´aceres 52 63 C´adiz 81 114 Cantabria 69 69 Castell´on 70 74 Ciudad Real 59 71 C´ordoba 68 89 Coru˜na 129 161 Cuenca 28 45 Girona 75 88 Granada 70 75 Guadalajara 32 36 Guip´uzcoa 80 112 Huelva 36 63 Huesca 47 54 Islas Baleares 102 115 Ja´en 51 95 Tabla 4.3: Casos de mortalidad seg´un provincias y periodos. Parte 1. 75 Paula Navarro Esteban provincia casos a˜no 1991 casos a˜no 2008 Las Palmas 71 113 Le´on 98 95 Lleida 56 62 Lugo 75 71 Madrid 432 544 M´alaga 97 147 Murcia 107 172 Navarra 84 83 Ourense 79 91 Palencia 27 32 Pontevedra 109 130 Rioja 43 53 Salamanca 51 56 Segovia 21 28 Sevilla 111 201 Soria 14 21 Sta. Cruz de Tenerife 81 110 Tarragona 61 70 Teruel 26 31 Toledo 60 106 Valencia 220 253 Valladolid 58 57 Vizcaya 98 162 Zamora 37 45 Zaragoza 132 147 Tabla 4.4: Casos de mortalidad seg´un provincias y periodos. Parte 2. 76 An´alisis espacial de riesgos de mortalidad No obstante, hay que calcular la probabilidad de que los riesgos estimados del segundo periodo seg´un las tasas del primer periodo sean mayores que 100 para cada provincia. Para ello ejecutamos la siguiente sintaxis en R: > PRP.p2.resp.p1<-1-resultado.p2.resp.p1$summary.linear.predictor[,"0 cdf"] En la Figura (4.11) se representa la probabilidad de que los riesgos estimados del segundo periodo seg´un las tasas del primer periodo sean mayores que 100 para cada provincia, siguiendo las tonalidades, las segmentaciones y las interpretaciones que en las Figuras (4.5) y (4.8). Se observa que las provincias de las Palmas, Baleares, Castell´on, Huesca, Zaragoza, la Rioja, Valladolid, Asturias, la Coru˜na, Pontevedra y Orense, que en la Figura (4.10) ten´ıan un riesgo suavizado por encima de 100, resultan tener una probabilidad de que el riesgo estimado sea mayor que 100, mayor de 0.9. Es decir, que estas provincias han incrementado su riesgo de que sus habitantes mueran por c´ancer de pr´ostata en el periodo de 1992-2008 frente al periodo 1975-1991. Tambi´en hay que destacar que Gerona, Cantabria, Salamanca y Lugo tienen una probabilidad de que el riesgo estimado sea mayor que 100, mayor que 0.8 y menor que 0.9. Por tanto, estas provincias tambi´en han incrementado su riesgo de morir por c´ancer en el segundo periodo frente al primero. Respecto a las provincias de Navarra, L´erida, Valencia, Albacete y Zamora, que ya hemos visto que ten´ıan un exceso de riesgo en comparaci´on con el riesgo de Espa˜na, vemos en este mapa que tienen una probabilidad baja (menor de 0.8), as´ı que no podemos asegurar que su riesgo sea mayor comparado con el de Espa˜na en global en el periodo anterior. Por lo tanto, las provincias que han incrementado su riesgo de que sus habitantes mueran por c´ancer de pr´ostata en el primer periodo frente al segundo periodo son distintas desde el punto de vista bayesiano que desde el punto de vista frecuentista (calculadas en el cap´ıtulo 3). Concretamente, Santa Cruz de Tenerife, Baleares, Valencia, Albacete, L´erida, Navarra y Zamora tienen una SMR entre 100 y 110, pero no se puede asegurar que sus riesgos suavizados sean mayores que 100 ya que no tienen una probabilidad mayor que 0.8. En resumen, el mapa representado en la Figura (4.11) muestra que el riesgo de mortalidad por c´ancer de pr´ostata se ha incrementado en el segundo periodo respecto del primero en la mitad norte de Espa˜na. Posiblemente sea como resultado de un incremento mayor de la esperanza de vida en la mitad norte que en la mitad sur de Espa˜na en los ´ultimos a˜nos [26] y un aumento de la poblaci´on, ya que como hemos comentado en el primer cap´ıtulo, uno de los principales factores de riesgo de padecer c´ancer de pr´ostata era la edad. En la Figura (4.12) representamos la funci´on de densidad de los riesgos suavizados del periodo 1975-1991 y del periodo 1992-2008, estandarizadas seg´un las tasas espec´ıficas del periodo 1975- 1991 (estandarizaci´on interna). El eje de las y’s est´a en escala logar´ıtmica. En esta gr´afica observamos que, para los riesgos suavizados por encima de 100, la funcion de densidad del primer periodo limita un ´area mayor que la funci´on de densidad del segundo periodo. 77 Paula Navarro Esteban Probabilidad de que los riesgos estimados del periodo 1992−2008 según tasas del periodo 1975−1991 sean superiores a 100 [0.0,0.1) [0.1,0.2) [0.2,0.8) [0.8,0.9) [0.9,1.0] Figura 4.11: Probabilidad de que los riesgos estimados del periodo 1992-2008 seg´un tasas del periodo 1975-1991 sean superiores a 100 78 An´alisis espacial de riesgos de mortalidad 70 80 90 100 110 120 130 0.00 0.01 0.02 0.03 0.04 RMEs Densidad Periodo 1975−1991 1992−2008 Función de densidad de las RMEs Figura 4.12: Funci´on de densidad de las RMEs calculadas durante el periodo 1975-1991 y el periodo 1992-2008 79