Análisis de la siniestralidad vial mediante modelos de datos de conteo

Grado en Matemáticas TRABAJO FIN DE GRADO Análisis de la siniestralidad vial mediante modelos de datos de conteo Estadística e Investigación Operativa Autor: José Manuel López Jiménez Director: José Luis Pino Mejías Sevilla, Junio 2019 Índice general Resumen....................................... iii Abstract....................................... iv Introducción..................................... v ÍndicedeFiguras.................................. vii ÍndicedeCuadros.................................. ix 1. Modelo Lineal Generalizado 1 1.1. Motivación.................................... 1 1.2. Componentes ................................. 2 1.2.1. Componente aleatoria . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2. Componente sistemática . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.3. Funciónlink.............................. 5 1.3. Hipótesis.................................... 5 1.4. EstimaciónenMLG ............................. 6 1.4.1. Método de máxima verosimilitud . . . . . . . . . . . . . . . . . . 6 1.4.2. AplicaciónalMLG.......................... 7 1.5. Inferencia: Ajuste del modelo . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5.1. Bondad de ajuste en MLG . . . . . . . . . . . . . . . . . . . . . . 10 1.5.2. Residuos................................ 13 1.5.3. Interpretación............................. 14 2. Modelos de conteo 15 2.1. Variabledeconteo .............................. 15 2.2. Modelos para variables de conteo . . . . . . . . . . . . . . . . . . . . . . 15 2.2.1. Sobredispersión............................ 15 3. Modelo de Regresión de Poisson 19 3.1. Características ................................ 20 3.2. Propiedades.................................. 20 3.3. Bondaddeajuste................................ 21 3.4. Dispersión.................................... 21 3.5. Modelo de Regresión Bivariante de Poisson . . . . . . . . . . . . . . . . . 22 3.5.1. Distribución de Poisson Bivariada . . . . . . . . . . . . . . . . . . 22 3.5.2. Modelo de Regresión Bivariada de Poisson . . . . . . . . . . . . . 22 4. Modelo de Regresión Binomial Negativa 25 4.1. Sobredispersión constante . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2. Sobredispersión variable . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2.1. Derivación en términos de Poisson-Gamma . . . . . . . . . . . . . 27 i 4.2.2. Derivación en términos de Binomial Negativa . . . . . . . . . . . 28 5. Exceso de ceros en datos de conteo 31 5.1. Modelos de regresión truncados en ceros . . . . . . . . . . . . . . . . . . 32 5.1.1. Poisson cero truncado . . . . . . . . . . . . . . . . . . . . . . . . 32 5.1.2. Binomial negativa cero truncado . . . . . . . . . . . . . . . . . . . 32 5.2. Modelos de regresión Hurdle . . . . . . . . . . . . . . . . . . . . . . . . . 33 5.3. Modelos de regresión cero inflado . . . . . . . . . . . . . . . . . . . . . . 34 5.3.1. Poisson cero inflado . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.3.2. Binomial negativa cero inflado . . . . . . . . . . . . . . . . . . . . 34 6. Criterios de selección de modelos 35 6.1. Criterio de información de Akaike . . . . . . . . . . . . . . . . . . . . . . 35 6.2. Criterio de información Bayesiana . . . . . . . . . . . . . . . . . . . . . . 35 7. Aplicación a los accidentes de tráfico 37 7.1. Estadística en los accidentes . . . . . . . . . . . . . . . . . . . . . . . . . 37 7.2. Aplicación a nuestra base de datos . . . . . . . . . . . . . . . . . . . . . 38 7.2.1. Lectura de los datos . . . . . . . . . . . . . . . . . . . . . . . . . 38 7.2.2. Estudio descriptivo . . . . . . . . . . . . . . . . . . . . . . . . . . 39 7.2.3. Aplicación del modelo de Poisson . . . . . . . . . . . . . . . . . . 45 7.2.3.1. Vías interurbanas . . . . . . . . . . . . . . . . . . . . . . 45 7.2.3.2. Víasurbanas ......................... 51 7.2.4. Aplicación del modelo Binomial Negativa . . . . . . . . . . . . . . 57 7.2.4.1. Vías interurbanas . . . . . . . . . . . . . . . . . . . . . . 57 7.2.4.2. Víasurbanas ........................ 59 7.2.5. Aplicación del modelo de Poisson Bivariante . . . . . . . . . . . . . 61 7.2.6. Comparación y selección de modelos . . . . . . . . . . . . . . . . 63 7.2.6.1. Modelos univariantes . . . . . . . . . . . . . . . . . . . . 63 7.2.6.1.1. Vías interurbanas . . . . . . . . . . . . . . . . . 63 7.2.6.1.2. Vías urbanas . . . . . . . . . . . . . . . . . . . 65 7.2.6.2. Modelos bivariantes . . . . . . . . . . . . . . . . . . . . 66 7.3. Conclusión................................... 67 7.3.1. Víasinterurbanas........................... 67 7.3.2. Víasurbanas ............................. 67 7.3.3. Poisson Bivariante . . . . . . . . . . . . . . . . . . . . . . . . . . 68 A. Apéndice: Código del trabajo en R 69 B. Apéndice: Paquete bivpois 77 C. Apéndice: Depuración de los datos 109 Bibliografía 135 ii Resumen El objetivo de este trabajo es estudiar la siniestralidad vial mediante la aplicación de modelos de datos de conteo, los cuales previamente comentaremos para después poder aplicarlos a nuestra base de datos. Comenzaremos por el estudio del Modelo Lineal Generalizado. Luego nos adentraremos en los modelos de datos de conteo como lo son el de Poisson y el Binomial Negativa. También veremos otros tipos de modelos más complejos que se podrían aplicar (como los truncados o en dos partes). Veremos los posibles problemas que se nos puede plantear debido al uso de dichas variables de conteo, cómo detectarlos y cómo solucionarlos. Para concluir, aplicaremos todo lo estudiado a nuestra base de datos que describiremos detalladamente en el último capítulo del trabajo, siempre con ayuda del software R. iii Abstract The aim of this work is to study the road accident rate by means of the application of models for count of data, which previously will be commented in order to can be applied to our database. We shall start with the so-called Generalized Linear Model. Then we shall enter the models for count of data such as the one due to Poisson and the so-called Negative Binomial. Moreover, we shall deal with other kinds of more complex models that could be applied (as for instance, the truncate model or the two-part model). We will analyze the problems that are likely to arise due to the use of the mentioned variables of count, as well as how to detect them and how to solve them. To finish, and always with the help of the software R, all the items that have been studied along this memory will be applied to our database which will be described in detail in the final chapter of this work. iv Introducción El trabajo titulado “Análisis de la siniestralidad vial mediante datos de conteo” consta de 7 capítulos de los cuales los 6 primeros abordan los diferentes modelos estadísticos y sus características, acabando en el último capítulo con su aplicación a una base de datos. En primer lugar, se comienza con el concepto de modelo y la necesidad de tenerlo para poder realizar estudios estadísticos para comprender el medio. Como los datos de conteo se pueden ver como un caso particular del Modelo Lineal Generalizado (MLG), dicho modelo se estudia en este primer capítulo. Así pues, se estudian sus diferentes componentes (aleatoria, sistemática y función link). Se mencionan las hipótesis realizadas sobre dicho modelo. También se realiza una estimación de los parámetros mediante el Método de Máxima Verosimilitud y se hace inferencia, es decir, se ve como se ajusta nuestro modelo a un conjunto de observaciones y se evalúa la adecuación del modelo mediante diferentes medidas de bondad de ajuste. Como el ajuste no es perfecto, se introduce el concepto de residuos, que nos marca la discrepancia entre el modelo obtenido y los datos. Se finaliza con una interpretación del modelo. En el segundo capítulo, se introduce el concepto de dato y variable de conteo, así como los diferentes problemas que surgen al utilizar dichas variables. Por este motivo se hace una lista con los diferentes modelos que son adecuados para el estudio de dichas variables de conteo. En el tercer capítulo se estudia el modelo de regresión de Poisson, modelo base para el estudio de variables de conteo pero con un problema, la equidispersión, es decir, la igualdad entre media y varianza, cosa que provoca que dicho modelo sea poco aplicable en la práctica ya que, por lo general, los datos poseen una mayor varianza que media (hecho que se conoce como sobredispersión). Es por ello que el trabajo continúa con el estudio de modelos que podamos aplicar. Además, se estudian estadísticos para ver la bondad de ajuste del modelo y se dan criterios para estudiar la dispersión de los datos. Se concluye el capítulo con el estudio del modelo de Regresión Bivariante de Poisson, modelo cuyo vector respuesta es bidimensional. En el cuarto capítulo se estudia el modelo de regresión Binomial Negativa, modelo estrella para el estudio de las variables de conteo y de datos que poseen sobredispersión. En el caso de sobredispersión constante, se parte de un modelo mixto Poisson-Gamma. En el caso de sobredispersión variable, podemos obtener la binomial negativa desde dos caminos. La primera forma es viéndolo como un modelo de Poisson con heterogeneidad Gamma de media uno. La segunda forma es tratándolo como una función de probabilidad propiamente dicha, como un miembro de la familia del MLG. En el quinto capítulo se aborda el problema del exceso de ceros en datos de conteo. Dicha cantidad de ceros no es compatible con los modelos de Poisson y Binomial Negativa. A continuación se introduce el concepto de falsos ceros y ceros auténticos. Se plantean diferentes tipos de modelos: truncados en cero, Hurdle y cero inflado. Se concluye el estudio teórico en este sexto capítulo con dos criterios de selección de modelos, uno basándose en el estadístico AIC y el otro en el estadístico BIC. En el séptimo y último capítulo se encuentra la parte práctica del trabajo, donde se aplican los modelos estudiados a lo largo del trabajo a nuestra base de datos elaborada en v el Apéndice C. Se comparan dichos modelos, seleccionando los más adecuados para así finalizar con las conclusiones. Todo ello se ha realizado con el software estadístico R. La parte del código tanto del trabajo como de la depuración de los datos como la utilizada auxiliarmente se encuentran en los apéndices del trabajo. Para concluir, se muestra la bibliografía utilizada para elaborar dicho trabajo. vi Índice de figuras 7.1. Figura con todos los pares de variables para fallecidos interurbana . . . . 40 7.2. Número de fallecidos en vías interurbanas en Sevilla por año . . . . . . . . 41 7.3. Número de fallecidos en vías interurbanas en Sevilla por población . . . . . 41 7.4. Número de fallecidos en vías interurbanas por año en cada Comunidad Autónoma................................... 42 7.5. Figura con todos los pares de variables en vías urbanas . . . . . . . . . . 43 7.6. Número de fallecidos en vías urbanas en Sevilla por año . . . . . . . . . . 44 7.7. Número de fallecidos en vías urbanas en Sevilla por población . . . . . . 44 7.8. Número de fallecidos en vías urbanas por año en cada Comunidad Autónoma 45 7.9. Residuos y observaciones influyentes para fallecidos en vías interurbanas conPoisson .................................. 50 7.10. Residuos y observaciones influyentes para fallecidos en vías interurbanas conPoisson .................................. 50 7.11. Residuos y observaciones influyentes para fallecidos en vías urbanas con Poisson..................................... 55 7.12. Residuos y observaciones influyentes para fallecidos en vías urbanas con Poisson..................................... 56 vii 1.2. COMPONENTES y por lo tanto µ=E(Y) = b0(θ). Además, sabemos que E(l00(θ, y)) = −E[(l0(θ, y))2], entonces V ar(l0(θ, y)) = E[(l0(θ, y))2] = −E(l00(θ, y)) = b00(θ) a(φ). Por otro lado, V ar(l0(θ, y)) = V ar y−b0(θ) a(φ)!=1 a2(φ)V ar(Y), y en consecuencia V ar(Y) = a(φ)b00(θ). Se comprueba que la varianza de Y se puede escribir como producto de dos funciones, una que depende del parámetro natural y otra que depende del parámetro de dispersión. 4. Elementos de la familia exponencial Cuadro 1.1: Elementos de la familia exponencial Familia exponencial Distribuciones Expresión Rango θ a(φ)b(θ)V(µ) Bernouilli B(p)0,1 ln( p 1−p) 1 ln(1 + eθ)p(1 −p) Binomial Bi(n, p)[0,n] ln( p 1−p) 1 nln(1+eθ)np(1 −p) Normal N(µ, σ2) (−∞,∞)µ σ2θ2 21 Gamma G(µ, v) (0,∞)−1 µ 1 v-ln(-θ)µ2 Poisson P(µ)N∪{0}ln(µ) 1 eθµ Binomial negativa BN(p,r) N∪{0}ln(1-p) 1 −r(ln(1−eθ)) r(1 −p) p2 1.2.2. Componente sistemática La componente sistemática, también llamada predictor lineal, se escribe como combinación lineal de las variables explicativas xi: η=β0+ p X i=1 βixi. 4CAPÍTULO 1. MODELO LINEAL GENERALIZADO 1.3. HIPÓTESIS 1.2.3. Función link Sea µ = E ( Y ). Se denomina función link a la función g que nos relaciona µ con η como: g(µ) = η=β0+ p X i=1 βixi. Es decir, nos relaciona la componente aleatoria con la componente sistemática, o lo que es lo mismo, el valor esperado con las variables explicativas. La función g más simple que podemos tomar es la identidad, la cual da lugar al modelo lineal clásico: µ=η=β0+ p X i=1 βixi. Además, cada elemento de la familia exponencial viene determinado por una función link canónica, que consiste en relacionar el parámetro natural directamente con el predictor lineal: θi=θ(µi) = ηi=β0+ p X j=1 βjxij g(µi) = θ(µi). En el cuadro 1.2 se recogen las funciones link canónicas más importantes que aparecen en el apartado 4.4 de Hilbe (2012). Cuadro 1.2: Funciones link canónicas Funciones Link Distribuciones Expresión Bernouilli g(µi) = θ(µi) =ln( µi 1−µi ) = ηi Binomial g(µi) = θ(µi) =ln( µi 1−µi ) = ηi Gamma g(µi) = θ(µi) = −1 µi =ηi Binomial negativa g(µi) = θ(µi) =ln( αµi 1 + αµi ) = ηi Normal g(µi) = θ(µi) = µi=ηi Poisson g(µi) = θ(µi) =ln(µi) = ηi 1.3. Hipótesis En este apartado se estudian las principales hipótesis que se hacen sobre el MLG. 1. Linealidad de los parámetros. Se supone que se establece una linealidad entre las variables de nuestro modelo. 2. Grados de libertad positivos. CAPÍTULO 1. MODELO LINEAL GENERALIZADO 5 1.4. ESTIMACIÓN EN MLG Esto quiere decir que, como mínimo, necesitamos que haya mayor número de datos que de variables. En caso contrario no existe la inversa de XTX ya que no tendría rango máximo y por tanto nuestro sistema XTXβ =XTY sería indeterminado. 3. Parámetros constantes. Suponemos que los parámetros β0, .., βp∈R. 4. Regresores no estocásticos. Suponemos que nuestra matriz numérica de datos Xes fija. 5. Independecia lineal entre variables explicativas. Esta hipótesis nos dice que cada variable explicativa contiene información adicional de la variable respuesta que no está contenida en otra. 1.4. Estimación en MLG Hay dos formas más comunes de realizar una estimación de los parámetros, mediante mínimos cuadrados o mediante el método de maxima verosimilitud . Utilizamos el segundo método ya que nos aporta mejores propiedades. Consideremos la función de densidad de la familia exponencial como f(y;θ) = exp{a(y)b(θ) + c(θ) + d(y)}. 1.4.1. Método de máxima verosimilitud Sea Y1, ..., Ynvariables aleatorias cuya función de densidad es: f(y1, ..., yn;θ1, ..., θp). Consideremos y = ( y1, ..., yn )y θ = ( θ1, ..., θp ). Sea Ωel espacio paramétrico (conjunto de valores que puede tomar θ ), según el apartado 4.2 de Dobson (2008), se define el estimador de máxima verosimilitud de θ como el valor ˆ θ que maximiza la función de verosimilitud: L(ˆ θ;y)≥L(θ;y)∀θ∈Ω. Equivalentemente, también maximiza la función log-verosimilitud, la cual es el logaritmo de la función de verosimilitud (l(θ;y) =logL(θ;y)), esto es: l(ˆ θ;y)≥l(θ;y)∀θ∈Ω. Este estimador se obtiene haciendo : ∂l(θ;y) ∂θj = 0,∀j∈ {1, ..., p}, 6CAPÍTULO 1. MODELO LINEAL GENERALIZADO 1.4. ESTIMACIÓN EN MLG y luego viendo que la matriz de las segundas derivadas evaluadas en ese valor obtenido de ˆ θes definida negativa: ∂2l(θ;y) ∂θj∂θk . Una propiedad a destacar que pose la estimación por máxima verosimilitud es la propiedad de la invarianza , esto es que si nosotros tenemos una función g ( θ ), entonces su estimador por máxima verosimilitud es g ( ˆ θ ). Otras propiedades a destacar es la consistencia, suficiencia y eficiencia asintótoca. 1.4.2. Aplicación al MLG Se quiere obtener un estimador de β del MLG. La estimación se va a obtener numéricamente mediante un método iterativo. Para ello hemos aplicado el método visto anteriormente y basado en el apartado 3.2 y 4.4 de Dobson (2008). Sea Y1, ..., Yn variables aleatorias independientes, podemos expresar su función loglogaritmo como: l(θ;y) = Xa(yi)b(θi) + Xc(θi) + Xd(yi), donde E ( Yi ) = µi = −c0 ( θi ) /b0 ( θi )y g ( µi ) = xT iβ = ηi , donde g es una función monótona y diferenciable. Una propiedad de la familia exponencial es que tienen la suficiente regularidad como para asegurar que el máximo global de la función log-verosimilitud se alcanza en la solución única de las ecuaciones ∂l/∂β = 0. Primero calculemos la V ar(Y). Para ello, U=∂l ∂θ =a(y)b0(θ) + c0(θ), y U0=∂2l ∂θ2=a(y)b00(θ) + c00(θ). Así que E(U) = b0(θ)E(a(Y)) + c0(θ), pero como E(U)=0, se tiene que E(a(Y)) = −c0(θ)/b0(θ), V ar(U) = (b0(θ))2V ar(a(Y)), E(−U0) = −b00(θ)E(a(Y)) −c00(θ). Luego con todo esto y usando que V ar(U) = E(U2) = E(−U0), tenemos que V ar(Y)=[b00(θ)c0(θ)−c00(θ)b0(θ)]/[b0(θ)]3. Equivalentemente tenemos que V ar(Yi)=[b00(θi)c0(θi)−c00(θi)b0(θi)]/[b0(θi)]3. CAPÍTULO 1. MODELO LINEAL GENERALIZADO 7 1.4. ESTIMACIÓN EN MLG La función score respecto al parámetro βjse define como Uj=∂l(θ;y) ∂βj = n X i=1 ∂li ∂βj , donde li=yib(θi) + c(θi) + d(yi). Para obtener Ujhacemos ∂li ∂βj =∂li ∂θi ∂θi ∂µi ∂µi ∂βj . Derivando liy usando E(Yi)obtenemos ∂li ∂θi =yib0(θi) + c0(θi) = b0(θi)(yi−µi). Derivando µiy usando V ar(Yi) ∂µi ∂θi =−c00(θi) b0(θi)+c0(θi)b00(θi) (b0(θi))2=b0(θi)V ar(Yi). Derivando g(µi) ∂µi ∂βj =∂µi ∂ηi ∂ηi ∂βj =xij ∂µi ηi , donde xij es el elemento j-ésimo de xT i. Por lo que ∂li ∂βj =∂li ∂θi ∂µi ∂βj /∂µi ∂θi =(yi−µi)xij V ar(Yi) ∂µi ∂ηi!. Y así obtenemos Uj= n X i=1 (yi−µi)xij V ar(Yi) ∂µi ∂ηi!. En general, Uj = 0 es no lineal, por lo que tenemos que resolverlo mediante iteración numérica. Para ello se usa el método de Newton-Raphson , para así escribir la m-ésima aproximación como b(m)=b(m−1) − ∂2l ∂βj∂βk!−1 β=b(m−1) U(m−1), donde U(m−1) es el vector Ujevaluado en β=b(m−1). Un procedimiento alternativo a este es el llamado método de score , el cual consiste en sustituir la matriz de las segundas derivadas por E ∂2l ∂βj∂βk!. 8CAPÍTULO 1. MODELO LINEAL GENERALIZADO 1.4. ESTIMACIÓN EN MLG Se define la matriz de información de fisher como F(β) = E[j(β;Y)], donde j ( β ; Y )es la matriz de información observada correspondiente al parámetro β, y viene determinada por j(β;y) = −∂2 ∂β∂βTl(β;y). Ahora, aplicando esto, obtenemos la ecuación b(m)=b(m−1) + (F(m−1))−1U(m−1), donde F(m−1) es F evaluada en b(m−1) . Multiplicando por F(m−1) en ambos lados de la ecuación obtenemos F(m−1)b(m)=F(m−1)b(m−1) +U(m−1).∗ Usando que el elemento (j, k)-ésimo de la matriz F es Fjk =E"∂li ∂βj ∂li ∂βk#=E (yi−µi)2xijxik V ar(Yi)2 ∂µi ∂ηi!2 = n X i=1 xijxik V ar(Yi) ∂µi ∂ηi!2 , escribimos la matriz F como F=XTWX, donde W es una matriz nxn diagonal cuyos elementos son wii =1 V ar(Yi) ∂µi ∂ηi!2 . Usando esto en *, y evaluando en b(m−1), obtenemos F(m−1)b(m−1) +U(m−1) =XTWz, donde los elementos de z son zi=X k xikb(m−1) k+ (yi−µi) ∂ηi ∂µi!, con µiy∂ηi/∂µievaluados en b(m−1). Concluimos entonces que, usando el método de score, podemos aproximar los β mediante la siguiente expresión XTWXb(m−1) =XTWz. CAPÍTULO 1. MODELO LINEAL GENERALIZADO 9 1.5. INFERENCIA: AJUSTE DEL MODELO 1.5. Inferencia: Ajuste del modelo Después de haber estudiado el modelo, nos plateamos ahora ver si se ajusta o no a la realidad. Es decir, ver la relación entre los datos observados y los datos esperados. Aquí vamos a tratar el tema de selección de variables, bondad de ajuste, coste e interpretabilidad. Tenemos que tener en cuenta que el tema de selección de variables es muy importante puesto que tener más variables de la cuenta nos provoca un mayor coste, y si tenemos menos nuestro modelo queda pobre. También afecta a la interpretabilidad del modelo el tener muchas variables ya que no sabremos si el papel de algunas de ellas está enmascarado por otras. Pero si tenemos pocas, quizás su interpretabilidad es fácil, pero puede que nuestro modelo entonces no se ajuste a la realidad. Ahí debemos de saber seleccionar cuales son las variables que mejor ajustan nuestro modelo, mediante procesos de p-valor o de contraste de hipótesis. Tenemos que tener cuidado a la hora de quitar variables puesto que podemos perder informacion sobre nuestro objetivo. Este proceso queda un poco a manos del investigador, que deberá establecer un criterio para quitar o no más variables sabiendo cuáles son las más influyentes en nuestro modelo. En la búsquedad de la estabilidad, siguiendo el apartado 5.7 de Dobson (2008), nos plateamos buscar un modelo de entre los dos siguientes: 1. Modelo saturado: El número de parámetros es igual al número de observaciones, es decir, ˆµi = yi . Este modelo tiene la misma función link que el modelo de interés. 2. Modelo nulo: En este modelo no se representa lo que estamos estudiando, es solamente un modelo de referencia donde consideramos un único parámetro, el valor esperado µ, para todas las observaciones. Lo que buscamos con este proceso es poder expresar nuestro modelo con el menor número de variables posibles (porque así tendremos un menor coste) pero manteniendo la validez del mismo. 1.5.1. Bondad de ajuste en MLG En este apartado se han estudiado diferentes procesos, según el apartado 5.5 de Dobson, que sirven para comparar nuestro modelo con el modelo de interés, para así obtener una visión de la exactitud y validez de nuestro modelo. Es decir, estudiar lo bien que se ajusta nuestro modelo a un conjunto de observaciones. Basándonos en el apartado 1.4 de Atoche (2017), se tiene lo siguiente. - Estadístico de máxima verosimilitud. Denotado por λ , este estadístico se basa en comparar la función de verosimilitud del modelo saturado con la del modelo de interés, obteniendo L ( βmax ; y )y L ( β ; y )respectivamente. Si el modelo de interés describe bien los datos, entonces L ( β ; y )debe ser aproximado a L ( βmax ; y ). Si el modelo es pobre, entonces L(β;y)será mucho menor que L(βmax;y). Así, obtenemos: λ=L(βmax;y) L(β;y). 10 CAPÍTULO 1. MODELO LINEAL GENERALIZADO 1.5. INFERENCIA: AJUSTE DEL MODELO Equivalentemente, podemos fijarnos en la diferencia entre las funciones log-logartimos de los modelos: logλ =l(βmax;y)−l(β;y). Grandes valores de logλ nos indica que el modelo no representa bien los datos. - Estadístico de desviación. Denotada por D ( y, µ ), es una medida de bondad de ajuste que compara el modelo saturado con el modelo de interés. Ésta se define como: D(y, ˆµ)=2φ(l(y, y)−l(ˆµ, y)), donde l ( y ; y )es el logaritmo de verosimilitud del modelo de interés y l ( ˆµ ; y )es el logaritmo de verosimilitud del modelo saturado. Bajo ciertas condiciones en las cuales n es fijo, el parámetro de dispersión es pequeño y conocido y las observaciones individuales convergen a una distribución normal y la desviación satisface: D(y, µ)≈χ2 n−p. Este resultado provee un valor de referencia para determinar la bondad de ajuste del modelo cuando el parámetro de dispersión es conocido. D ( y, ˆµ )es no negativa y un valor grande indica un ajuste pobre del modelo. Cuando el parámetro de dispersión φ es desconocido, entonces el resultado de la desviación no aplica. - Estadístico Chi-cuadrado de Pearson. Denotador por χ2 , este estadístico se define como: χ2=X i (yi−ˆµi)2 V ar(Yi). Cuando el número de observaciones es suficiente grande, ambos estadísticos, D y χ2 , son equivalentes. -Coeficiente de determinación. Denotado por R2, se define como: R2= 1 −D(y;µ) D(y;µ0), verifica que 0 ⩽R2⩽ 1. La interpretación de este coeficiente es la siguiente: cuanto mayor sea el número de variables que incorporemos en nuestro modelo, mayor sera este coeficiente y por tanto, se supone que mejor será el ajuste con el modelo. Hecho que no siempre ocurre y por el cual se define el siguiente coeficiente. - Coeficiente de determinación ajustado. Denotado por R2 adj , se introduce con el fin de corregir el coeficiente de determinación y se define como: R2 adj = 1 −D(y;µ)/(n−k) D(y, µ0)/(n−1). Este coeficiente incorpora una penalización sobre el número de variables en el modelo. Cuantas más variables introduzcamos en el modelo, D ( y ; µ )decrece, pero también lo hace su correspondiente grados de libertad. Por lo tanto, el numerador puede incrementarse si la disminución de la desviación a causa del aumento del número de variables no es compensado por la pérdida de grados de libertad. CAPÍTULO 1. MODELO LINEAL GENERALIZADO 11 1.5. INFERENCIA: AJUSTE DEL MODELO La relación entre R2yR2 adj es la siguiente: R2 adj = 1 −n−1 n−kD(y;µ) D(y;µ0)= 1 −n−1 n−k(1 −R2). Intervalos de confianza y Test de Hipótesis Dos de las herramientas más usadas de la inferencia estadística son los intervalos de confianza y los tests de hipótesis. Por ejemplo, los tests de hipótesis son necesarios para comparar el ajuste de dos modelos ajustados a los datos. Tanto para realizar tests como intervalos de confianza necesitamos las distribuciones muestrales de los estadísticos involucrados. Vamos a realizar un contraste de hipótesis y un intervalo de confianza para nuestros coeficientes βi. 1. Contraste de hipótesis. Una hipótesis se contrasta comparando nuestras predicciones con la realidad. Si coinciden, dentro de un margen de error admisible, mantenemos la hipótesis. En caso contrario, la rechazamos y buscamos otras que sean capaces de explicar los datos. Planteamos una hipótesis nula a priori creíble y sólo la rechazamos cuando existe suficiente evidencia en los datos en contra de la misma. Si rechazamos la nula, implícitamente no rechazamos otra hipótesis llamada alternativa. Vamos a plantearnos el siguiente contraste: H1:η∈L⊂Rq, donde L está parametrizado como η=X1β, y vamos a considerar la hipótesis H0:η∈L0⊂Rs, donde η=X0αys<q, con la hipótesis alternativa H1:η∈L\L0. Para este contraste, usamos el estadístico de razón de verosimilitud , que tiene la forma: −2logl =D(µ(ˆ β);µ(β(ˆα))), donde Des la desviación. Cuando H0 es cierto, la desviación es asintóticamente igual a la distribución χ2 ( k−m ). Planteamos el siguiente contraste: (H0:Aβ =S H1:Aβ 6=S, donde A es una matriz de rango s⩽p+ 1 y S un vector de dimensión s. Para este contraste, se usa el estadístico de Wald: SW= (Aˆ β−S)T(AF−1(ˆ β)AT)(Aˆ β−S). 2. Intervalo de confianza. Para elaborar un IC (intervalo de confianza), necesitamos usar un estadístico. 12 CAPÍTULO 1. MODELO LINEAL GENERALIZADO 1.5. INFERENCIA: AJUSTE DEL MODELO Estadístico razón de verosimilitud para así obtener: IC(βj,1−α) = {βj∈R/l1(βj)⩾l0}, donde l0=l(ˆ β)−0.5χ2 1−α(1), yl1(βj)denota la función log-verosimilitud para βj. Estadístico de Wald para así obtener: IC(β, 1−α) = {β∈Rp+1/(ˆ β−β)(V ar(ˆ β))−1(ˆ β−β)< χ2 p+1,n−1}. 1.5.2. Residuos Los residuos aparecen cuando el modelo que planteamos no coincide exactamente con el problema que queremos modelar. En la realidad esto ocurre siempre, por ello los estudiamos y nos dan una indicación de las diferencias entre los valores ajustados y observados. Estos errores pueden venir a raiz de los datos o por uso de un mal modelo para su aproximación. Un aspecto importante en los MLG es que la ortogonalidad de los residuales y los valores ajustados no se cumple. Es decir, la descomposición datos = ajustados + residuos no se cumple. En general, según el apartado 4.4 de Hilbe (2007), podemos definir tres tipos de residuales: 1. Residual básico. rb i=yi−ˆyi,∀i= 1, ..., n. 2. Residual de Pearson. rp i=yi−ˆµi qV ar(ˆ µi),∀i= 1, ..., n, 3. Residual de Score. rs i=yi−ˆµi V ar(ˆµi) ∂η ∂µ!−1 i 4. Residual de desviación. ed i=qdi∗signo(yi−ˆµi),∀i= 1, ..., n, donde dise define como la componente de desviación: di= 2(l( ˆµi, yi)−l( ˆµi, µi)). Los residuales sirven para detectar datos atípicos y desvíos del modelo con relación a los datos. En MLG se deben usar gráficas similares a las discutidas en modelos lineales. CAPÍTULO 1. MODELO LINEAL GENERALIZADO 13 3.1. CARACTERÍSTICAS A partir de esto, obtenemos lo que se conoce como Modelo de Regresión de Poisson , que se escribe E(Yi|xi) = eβ0+β1xi1+...+βnxin . La función log-verosimilitud es l(µ;y) = n X i=1{yiln(yi)−µi−ln(yi!)}. En el caso particular en el que Y= 0, esta expresión se reduce a li(µ;0) = −µi. Además, se puede escribir la función log-verosimilitud en términos de Xβ l(Xβ;y) = Pn i=1{yiln(exp(Xiβ)) −exp(Xiβ)−ln(yi!)} =Pn i=1{yiXiβ−exp(Xiβ)−ln(yi!)} li(Xβ;yi= 0) = −exp(Xiβ) . 3.1. Características La distribución de Poisson Yposee diversas características: 1. Sean Y1, ..., Yn∈P(µ)independientes, entonces Y= n X i=1 Yi∈P(nµ). 2. Equidispersión. 3. Como consecuencia del teorema central del límite , para valores grandes de µ , una variable Poisson puede aproximarse por una normal de la siguiente forma Y−µ √µ≈N(0,1). 4. La distribución de Poisson es el caso límite de una distribución binomial en el caso en el que n−→ ∞ yp−→ 0de tal forma que µ=pn =cte. 3.2. Propiedades Estudiamos qué propiedades deben cumplir los datos para el modelo de Poisson: 1. La probabilidad de que ocurra un solo suceso sobre un intervalo pequeño es aproximadamente proporcional al tamaño del intervalo. 2. La probabilidad de que ocurran dos sucesos en un mismo intervalo cuya longitud tiende a cero es despreciable. 20 CAPÍTULO 3. MODELO DE REGRESIÓN DE POISSON 3.3. BONDAD DE AJUSTE 3. La probabilidad de que ocurra un suceso en un intervalo de una cierta longitud no varía en cualquier otro intervalo con esa misma longitud. 4. La probabilidad de que ocurra un suceso en un intervalo es independiente de la probabilidad de que ocurra un suceso en otro intervalo que no se superpone al anterior. Cuando las propiedades 3 y 4 no se cumplen, los datos pueden tener mayor varianza que media, que es lo que se define como sobredispersión. 3.3. Bondad de ajuste Una medida de discrepancia entre los valores observados y los predichos es la llamada desviación, la cual se escribe como D= 2Pn i=1{yiln(yi)−yi−yiln(ˆµi) + ˆµi} = 2Pn i=1{yiln yi ˆµi!−(yi−ˆµi)} Di(yi= 0) = 2ˆµi . El estadístico χ2 de Pearson también nos sirve para estudiar esta bondad, el cual se escribe como χ2= n X i=1 (yi−ˆµi)2 ˆµi . Es de destacar el estadístico G2= 2 n X i=1 yiln(yi/ˆµi), el cual nos dice que para G2= 0, el ajuste de nuestro modelo es perfecto. 3.4. Dispersión Cuando se trabaja con una base de datos de conteo puede ocurrir que el modelo Poisson puede parecer sobredisperso y en realidad no lo es, o puede efectivamente presentar sobredispersión. En el primer caso, una simple corrección al modelo puede hacer desaparecer la variabilidad no deseada. Si después de los ajustes realizados sobre el modelo el problema de sobredispersión no desaparece, se presenta el segundo caso, donde se deben buscar modelos alternativos, que puede ser el modelado a partir de la distribución Binomial Negativa (no tiene restricción de igualdad en media y varianza), o el uso de modelos más complejos como son los modelos compuestos y modelos en 2 partes. Podemos usar, según el apartado 13.4 de Hilbe (2012), un test basado en este supuesto siguiendo el siguiente procedimiento: 1. Obtener ˆµi. 2. Calcular z=(yi−ˆµi)2−yi ˆµi√2. CAPÍTULO 3. MODELO DE REGRESIÓN DE POISSON 21 3.5. MODELO DE REGRESIÓN BIVARIANTE DE POISSON 3. Escribir z como un modelo de constante única. El test de hipótesis que vamos a usar es el siguiente (H0:V(y) = E(y) H1:V(y) = E(y) + αg{E(y)}, el cual se lleva a cabo a través del t test y α es el nivel de significación. El estadístico que se suele usar para este contraste de hipótesis es CT = − 2( lPoisson −lBN ). Se tiene que la distribución asintótica del CT es χ2 1. Si existe evidencia suficiente para probar que los datos no siguen una distribución Poisson, entonces será necesario emplear un modelo de conteo alternativo que se ajuste al tipo de supuesto violado en la distribución de los datos, los cuales aparecen nombrados en el apartado de dispersión del capítulo 2. 3.5. Modelo de Regresión Bivariante de Poisson El anterior modelo estudiado hacía referencia al caso en el que solamente teníamos una variable respuesta con una distribución de Poisson. Ahora vamos a tratar el caso en el que tenemos un vector respuesta de dos dimensiones que sigue una distribucion bivariante Poisson. 3.5.1. Distribución de Poisson Bivariada Sean X1, X2, X3 variables aleatorias que siguen una distribución de Poisson cuyos respectivos parámetros son λ1, λ2, λ3> 0. Según el apartado 3.2 de Quinteiro (2014), se dice que X = X1 + X3 e Y = X2 + X3 siguen cojuntamente una distribución Bivariada Poisson (BP(λ1, λ2, λ3)). La función de probabilidad es fBP (x, y) = exp{−(λ1+λ2+λ3)}λx 1 x! λy 2 y! min(x,y) X k=0 x k! y k!k!λ3 λ1λ2k , x, y ∈N∪{0}. Marginalmente cada variable aleatoria sigue una distribución de Poisson con E ( X ) = λ1 + λ3 y E ( Y ) = λ2 + λ3 . Además, cov ( X, Y ) = λ3 . En el caso en el que λ3 = 0 ambas distribuciones son independientes y la distribución Bivariada de Poisson se reduce al producto de dos distribuciones Poisson independientes. 3.5.2. Modelo de Regresión Bivariada de Poisson Sean i = 1 , ..., n el número de observaciones. Consideramos el modelo para la i-ésima observación (Xi, Yi)∼BP (λ1i, λ2i, λ3i), log(λki) = w0 kiβkk= 1,2,3, 22 CAPÍTULO 3. MODELO DE REGRESIÓN DE POISSON 3.5. MODELO DE REGRESIÓN BIVARIANTE DE POISSON donde w0 ki denota el vector de las variables explicativas para la observación i-ésima usada para el modelo λki (donde para k = 1 , 2 , 3nos referimos al modelo para X , Y y el modelo que nos relaciona ambas variables ,respectivamente) y βk vector de coeficientes. Las variables explicativas usadas para modelizar cada parámetro λki no tienen por qué ser las mismas. Se utiliza el algoritmo EM para obtener la estimación por máxima verosimilitud de los coeficientes beta. CAPÍTULO 3. MODELO DE REGRESIÓN DE POISSON 23 Capítulo 4 Modelo de Regresión Binomial Negativa La distribución binomial negativa se utiliza para modelos con datos de conteo. Principalmente para aquellos casos en los que aparece sobredispersión y no podemos abarcarlo con el modelo de Poisson. La distribución binomial negativa es una distribución discreta que toma valores no negativos. Esta distribución modela la probabilidad de observar y fallos antes del r-ésimo éxito en un número finito de sucesos independientes e idénticamente distribuidas (i.i.d.) Bernoulli. Sea Y∼BN ( r, p ), con 0 < p < 1y0 < r . La función de probabilidad del modelo binomial negativo se escribe como P(Y=y) = y+r−1 r−1!pr(1 −p)y. El modelo de regresión binomial negativo comenzó a estudiarse en 1949 por Anscombe. En 1987, Lawless caracterizó la binomial negativa como un modelo mixto, dando fórmulas para su media, varianza, log-verosimilitud y momentos. Más tarde, en 1990, Breslow citó el trabajo de Lawless mientras manipulaba el modelo de Poisson para ajustar los parámetros de la binomial negativa. Desde principios hasta finales de los 80, se contruyó el modelo binomial negativo como un modelo mixto para estudiar aquellos datos con sobredispersión. McCullagh y Nelder (1989) ajustando el valor de τ , podría ser considerado como un MLG. Los autores mencionaron de la existencia de una función canónica link, rompiendo con el concepto de modelo mixto, pero no consiguieron explicarlo. No fue hasta la mitad de los años 90 cuando la binomial negativa fue construida como un elemento de la familia MLG. Hilbe (2011) detalla como los modelos de regresión binomial negativo y geométrico pueden ser derivados de sus respectivas distribuciones de probabilidad. La idea que hay detrás de todo es que el modelo tradicionalmente llamado como binomial negativo no es mas que una log-link binomial negativa. La función canónica link puede tener propiedades de no convergencia y esto ocurre cuando τ aparece en la función link y su inversa. La binomial negativa está basada en el modelo mixto Poisson-Gamma. Hay dos métodos para motivar el modelo de regresión binomial negativo, ambos se estudian a continuación. 25 4.1. SOBREDISPERSIÓN CONSTANTE 4.1. Sobredispersión constante En primer lugar, basado en el apartado 14.1 de Hilbe (2012), consideremos una Poisson- Gamma, donde Y∼P(λ)en la que P[Yi=yi] = e−λiλyi i/yi!, yi∈Z+, λi>0,∀i= 1, ..., n, donde el parámetro λi es una variable aleatoria. Suponemos que su distribución es Ga ( δ, µi ) con función de densidad fλi(λi) = 1 (µi−1)!δµiλµi−1 ie−δλi, donde µi=eXiβ+offseti. Así, podemos expresar la esperanza y varianza de la variable de Poisson como E(λi) = eXiβ+offseti/δ V(λi) = eXiβ+offseti/δ2. Con lo visto hasta ahora, se puede escribir la función de densidad como fYi(yi|xi) = R∞ 0 e−λiλyi i yi! δµi (µi−1)!λµi−1 ie−λiδdλi =δµi yi!(µi−1)! R∞ 0λ(yi+µi)−1 ie−λi(δ+1)dλi =δµi yi!(µi−1)! (yi+µi−1)! (δ+ 1)yi+µiZ∞ 0 (δ+ 1)yi+µi (yi+µi−1)!λ(yi+µi)−1 ie−λi(δ+1)dλi | {z } 1 =(yi+µi−1)! yi!(µi−1)! δ 1 + δ!µi1 1 + δyi . Los momentos vienen dados por E(λi) = eXiβ+offseti/δ, V(λi) = eXiβ+offseti/δ2. Llamamos índice de sobredispersión al cociente entre la varianza y la media como V(Yi) E(Yi)=1 + δ δ, el cual es constante para todas las observaciones. Es decir, la varianza se puede expresar como un múltiplo de la media. Si consideramos ahora α= 1/δ (α > 0), podemos reescribir el modelo como fYi(yi|xi) = (yi+µi−1)! yi!(µi−1)! 1 1 + αµiα 1 + αyi , y V(Yi) E(Yi)= 1 + α. Sin embargo, su distribución no se puede escribir como un elemento de la familia exponencial. 26 CAPÍTULO 4. MODELO DE REGRESIÓN BINOMIAL NEGATIVA 4.2. SOBREDISPERSIÓN VARIABLE 4.2. Sobredispersión variable En segundo lugar, basado en el apartado 14.2 de Hilbe (2012), se puede obtener la binomial negativa desde dos caminos. La primera forma es viéndolo en su sentido tradicional, es decir, como un modelo de Poisson con heterogeneidad Gamma de media uno. La segunda forma es tratándolo como una función de probabilidad propiamente dicha, independientemente de la distribución de Poisson. Lo que hacemos es ver su función de probabilidad como la probabilidad de observar y fallos antes del r-ésimo éxito en un número finito de sucesos independientes e idénticamente distribuidas (i.i.d.) Bernoulli. Ambos métodos convergen a la misma función log-verosimilitud. 4.2.1. Derivación en términos de Poisson-Gamma La primera forma, como dijimos antes, consiste en verlo como un modelo de Poisson con heterogeneidad Gamma de media uno. Primero, tenemos que Yi∼P ( µi )donde µi = λiui de tal forma que lnµi=xiβ+i =lnλi+lnui, donde ui∼G(ν, ν)es el llamado efecto no observado, con ν > 0. La distribución de ycondicionada a xiyuiquedaría f(yi|ui) = e−λiui(λiui)yi yi!. La función de densidad de ycondicionada a xiquedaría f(yi|xi) = Z∞ 0 e−λiui(λiui)yi yi!g(ui)dui, cuyo resultado depende de la función g () que elijamos para modelar ui . Aquí lo hacemos para el caso en el que ui = exp ( i ). Gracias a que la media de Gamma es uno, se tiene el siguiente resultado f(yi|xi) = R∞ 0 e−λiui(λiui)yi yi! νν Γ(ν)uν−1 ie−νuidui =λyi i Γ(yi+ 1) νν Γ(ν)R∞ 0e−(λi+ν)uiu(yi+ν)−1 idui =λyi i Γ(yi+ 1) νν Γ(ν) Γ(yi+ν) (λi+ν)yi+νZ∞ 0 (λi+ν)yi+ν Γ(yi+ν)e−(λi+ν)uiu(yi+ν)−1 idui | {z } 1 =λyi i Γ(yi+ 1) νν Γ(ν) Γ(yi+ν) (λi+ν)yi+ν =λyi i Γ(yi+ 1) νν Γ(ν)Γ(yi+ν)ν λi+νν1 ννλi λi+νyi1 λyi i =Γ(yi+ν) Γ(yi+ 1)Γ(ν)ν λi+ννλi λi+νyi =Γ(yi+ν) Γ(yi+ 1)Γ(ν) 1 1 + λi/ν !ν 1−1 1 + λi/ν !yi . CAPÍTULO 4. MODELO DE REGRESIÓN BINOMIAL NEGATIVA 27 4.2. SOBREDISPERSIÓN VARIABLE Reparametrizando ν= 1/α obtenemos el modelo de Regresión f(yi|xi) = Γ(yi+ 1/α) Γ(yi+ 1)Γ(1/α)1 1 + αµi1/α 1−1 1 + αµiyi , µi=λi=exp(xiβ+offseti). En la parametrización hemos supuesto que α > 0. Para el caso α = 0 obtenemos el modelo de Poisson. A medida que α aumenta, aumenta también la sobredispersión del modelo. Escrito como miembro de la familia exponencial sería f(y;µ, α) = exp yln αµ 1 + αµ+1 αln 1 1 + αµ+lnΓy+1 α−lnΓ(y+ 1) −lnΓ1 α. De aquí, podemos obtener las componentes del MLG binomial negativo θ=ln αµ 1 + αµ, b(θ) = −1 αln 1 1 + αµ, b0(θ) = ∂b ∂µ ∂µ ∂θ =1 1 + αµ{µ(1 + αµ)} =µ, b00(θ) = ∂2b ∂µ2∂µ ∂θ 2+∂b ∂µ ∂2µ ∂θ2 =(−α (1 + αµ)2){µ2(1 + αµ)2}+1 1 + αµ(1 + αµ)(µ+ 2αµ2) =−αµ2+µ+ 2αµ2 =µ+αµ2, V(µ) = µ+αµ2, ∂V (µ) ∂µ = 1 + 2αµ. La función log-verosimilitud y la desviación se obtienen l(µ;y, α) = Pn i=1 yiln αµi 1 + αµi−1 αln(1 + αµi) + lnΓyi+1 α−lnΓ(yi+ 1) −lnΓ1 α, D= 2Pn i=1 yiln yi µi−yi+1 αln 1 + αyi 1 + αµi. 4.2.2. Derivación en términos de Binomial Negativa La segunda forma, como dijimos antes, consiste en ver su función de probabilidad como la probabilidad de observar y fallos antes del r-ésimo éxito en un número finito de sucesos independientes e idénticamente distribuidas (i.i.d.) Bernoulli. Componente aleatoria 28 CAPÍTULO 4. MODELO DE REGRESIÓN BINOMIAL NEGATIVA 4.2. SOBREDISPERSIÓN VARIABLE La función de probabilidad del modelo binomial negativo se escribe como f(y;r, p) = y+r−1 r−1!pr(1 −p)y. Escrito como miembro de la familia exponencial sería f(y;r, p) = exp (yln(1 −p) + rln(p) + ln y+r−1 r−1!). Así tenemos θ=ln(1 −p), eθ= 1 −p, p= 1 −eθ, b(θ) = −rln(p) = −rln(1 −eθ), con a ( φ ) = 1. La primera y segunda derivada se pueden obtener derivando b ( φ )respecto de θ b0(θ) = ∂b ∂p ∂p ∂θ =−r p{−(1 −p)}=r(1 −p) p =µ, b00(θ) = ∂2b ∂p2∂p ∂θ2+∂b ∂p ∂2p ∂θ2 =r p2(1 −p)2+r p(1 −p) =r(1 −p)2+rp(1 −p) p2=r(1 −p) p2. En términos de la media, la varianza se puede escribir como V(µ) = b00(θ) = µ+µ2 r. Reparametrizando α= 1/r, podemos expresar lo escrito anteriormente como p= 1 −eθ=1 1 + αµ,∗ θ=ln(1 −p) = ln αµ 1 + αµ, b(θ) = −1 αln(p) = 1 αln(1 + αµ), b0(θ) = 1−p αp =µ=1 α(e−θ−1), b00(θ) = 1−p αp2=µ+αµ2, g0(θ) = ∂ ∂µln αµ 1 + αµ=1 µ+αµ2, donde la varianza escrita en términos de la media viene dada por V(µ) = b00(θ) = µ+αµ2. CAPÍTULO 4. MODELO DE REGRESIÓN BINOMIAL NEGATIVA 29 Capítulo 7 Aplicación a los accidentes de tráfico En este capítulo se procede a aplicar todos los modelos vistos anteriormente para el estudio de variables de conteo a una base de datos. 7.1. Estadística en los accidentes Cada vez que ocurre algún tipo de accidente, los cuerpos de seguridad del Estado son los encargados de, a parte de solucionar todo lo posible, anotar y dejar constancia de lo ocurrido. Estos datos pasan a digitalizarse y se hacen muchos estudios con ellos. La finalidad de estos puede ser variada, pero siempre con la intención de disminuir el número de fallecidos. Por ejemplo: 1. Dentro de la ciudad, el ayuntamiento puede considerar necesario modificar los límites de velocidad permitidos dentro de ella. 2. Dentro de la ciudad, el ayuntamiento puede considerar necesario incorporar diversos tipos de obstáculos (como resaltos y badenes) o luces (actualmente se colocan en los pasos de peatones y los semáforos) con la intención de disminuir ese número de fallecidos. 3. La DGT (Dirección General de Tráfico) puede considerar necesario modificar los límites de velocidad permitidos, poner zonas de prohibido adelantar, incorporar radares y controles en zonas de riesgo y días de mayor movimiento de personas. 4. Modificar la carretera o redirigir el tráfico en aquellos lugares donde existen puntos negros de las carreteras (puntos donde se concentran gran cantidad de accidentes). 5. También los seguros necesitan de estos datos para subsistir. 6. Incorporación de cada vez más sistemas de seguridad en los vehículos. Como podemos observar, la capacidad de obtener todos estos datos tras los accidentes pueden servir de mucho para poner medidas y evitar o reducir los accidentes y por tanto, las víctimas en las carreteras. Esto es, que gracias a los avances de la tecnología nos permiten hacer estudios para poder sacar conclusiones y extraer medidas más reales y efectivas, las cuales poder aplicar a las carreteras, vehículos y conciencia de las personas para ese número de fallecidos que, desgraciadamente, es positivo. 37 7.2. APLICACIÓN A NUESTRA BASE DE DATOS 7.2. Aplicación a nuestra base de datos Todos los datos necesarios para realizar nuestro estudio han sido cogidos de las paginas de la DGT, del INE y del Ministerio de Fomento. A partir de ellos se ha creado una base de datos en el fichero datos.RData. 7.2.1. Lectura de los datos Nuestra base de datos está formada por 780 observaciones y 14 variables. En nuestro estudio tenemos dos variables objetivo y 9 variables explicativas (3 más que se explicarán más adelante). Las variables objetivo son: número de fallecidos en carreteras interurbanas y urbanas, representadas por fallecidos_interurbana y fallecidos_urbana, respectivamente. Las variables explicativas son: provincia, año, comunidades autónomas, población, número de vehículos, km de carreteras totales y precipitación, temperatura y horas de sol media en un año. Las variables aparecen recogidas por provincia, representadas por provincia, year, CCAA, poblacion, numero_vehiculos, km_carretera, precipitacion_media, temperatura_media y horas_sol_media, respectivamente. Es de notar que cada vez que estudiemos un modelo vamos a hacerlo para cada una de nuestras dos variables objetivo. En primer lugar, cargamos nuestra base de datos ya previamente depurada en el script Base_de_datos.R, y luego vemos un resumen de los valores que toman cada variable. load(file="datos.RData") summary(datos_f) ## provincia fallecidos_interurbana fallecidos_urbana ## Length:780 Min. : 0.0 Min. : 0.00 ## Class :character 1st Qu.: 22.0 1st Qu.: 3.00 ## Mode :character Median : 37.0 Median : 7.00 ## Mean : 48.0 Mean : 12.23 ## 3rd Qu.: 63.5 3rd Qu.: 13.00 ## Max. :253.0 Max. :149.00 ## year numero_vehiculos poblacion precipitacion_media ## Min. :2002 Min. : 39411 Min. : 66526 Min. : 6.658 ## 1st Qu.:2005 1st Qu.: 202288 1st Qu.: 321755 1st Qu.: 26.600 ## Median :2009 Median : 382112 Median : 584128 Median : 32.104 ## Mean :2009 Mean : 569163 Mean : 863809 Mean : 40.236 ## 3rd Qu.:2013 3rd Qu.: 645352 3rd Qu.:1002091 3rd Qu.: 43.019 ## Max. :2016 Max. :4474787 Max. :6466996 Max. :133.567 ## temperatura_media horas_sol_media km_carretera CCAA ## Min. :10.68 Min. :128.4 Min. : 26.83 Length:780 ## 1st Qu.:13.16 1st Qu.:207.4 1st Qu.:2369.99 Class :character ## Median :15.53 Median :236.7 Median :3322.39 Mode :character ## Mean :15.50 Mean :227.0 Mean :3184.52 ## 3rd Qu.:18.25 3rd Qu.:255.7 3rd Qu.:3984.99 ## Max. :21.75 Max. :283.4 Max. :6364.66 38 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS head(datos_f) ## provincia fallecidos_interurbana fallecidos_urbana year ## 1 Araba/Álava 40 9 2002 ## 2 Albacete 61 7 2002 ## 3 Alicante/Alacant 155 8 2002 ## 4 Almería 97 8 2002 ## 5 Ávila 45 3 2002 ## 6 Badajoz 74 19 2002 ## numero_vehiculos poblacion precipitacion_media temperatura_media ## 1 174235 288558 44.375 11.808 ## 2 202325 356661 33.150 14.842 ## 3 998821 1492669 20.933 18.342 ## 4 339226 514237 17.058 18.583 ## 5 93558 163305 28.358 11.258 ## 6 336611 663723 26.525 16.917 ## horas_sol_media km_carretera CCAA ## 1 169.133 1463.33 PVasco ## 2 261.133 3710.45 CLMancha ## 3 265.158 2704.29 CValenciana ## 4 268.508 2368.19 Andalucía ## 5 226.983 2545.57 CLeón ## 6 250.375 4882.83 Extremadura Calculamos la media y la varianza tanto de los fallecidos en carreteras interurbanas como urbanas. mean(datos_f$fallecidos_interurbana) ## [1] 48.00256 var(datos_f$fallecidos_interurbana) ## [1] 1446.339 mean(datos_f$fallecidos_urbana) ## [1] 12.22564 var(datos_f$fallecidos_urbana) ## [1] 347.2712 Se observa que la media y la varianza de ambas variables son bastante diferentes, por lo que ya tenemos que esperar que haya sobredispersión en los datos. 7.2.2. Estudio descriptivo En primer lugar, realizamos el estudio descriptivo con los fallecidos en vías interurbanas. Se estudian las relaciones existentes entre fallecidos_interurbana y las demás variables (menos fallecidos_urbana) con el siguiente gráfico: CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 39 7.2. APLICACIÓN A NUESTRA BASE DE DATOS attach(datos_f) pairs(fallecidos_interurbana~year+numero_vehiculos+poblacion+ precipitacion_media+temperatura_media+horas_sol_media+ km_carretera,panel=panel.smooth) fallecidos_interurbana 2002 2012 0e+00 6e+06 12 18 0 4000 0 200 2002 year numero_vehiculos 0e+00 0e+00 poblacion precipitacion_media 20 12 22 temperatura_media horas_sol_media 150 0 150 0 6000 0e+00 4e+06 20 80 150 250 km_carretera Figura 7.1: Figura con todos los pares de variables para fallecidos interurbana El número de fallecidos en vías interurbanas desciende a lo largo de los años, aumenta con la población y el número de vehículos, aumenta levemente con respecto a los km de carretera y horas media de sol, mientras que la precipitación y temperatura media no destaca. Es de notar que aquellas gráficas donde hay huecos es debido a las dos provincias con más habitantes (Madrid y Barcelona). Además, viendo el cuadro referente a la población y el número de vehículos, podemos intuir que ambas están altamente correladas, cosa que puede afectar a nuestro estudio. Se dibujan, a modo de ejemplo, 2 gráficas referentes a la provincia de Sevilla. library(ggplot2) ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=year, y=fallecidos_interurbana)) + geom_point(size=1)+geom_line() 40 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS 40 80 120 2004 2008 2012 2016 year fallecidos_interurbana Figura 7.2: Número de fallecidos en vías interurbanas en Sevilla por año En la gráfica se ve como el número de fallecidos, que al inicio de nuestro estudio superaba los 110 fallecidos por año, ha disminuido notablemente hasta mantenerse por debajo de los 40 fallecidos por año, menos de la mitad. ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=poblacion, y=fallecidos_interurbana)) + geom_point(size=1)+geom_line() 40 80 120 1700000 1750000 1800000 1850000 1900000 1950000 poblacion fallecidos_interurbana Figura 7.3: Número de fallecidos en vías interurbanas en Sevilla por población CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 41 7.2. APLICACIÓN A NUESTRA BASE DE DATOS Esta gráfica es muy similar a la anterior, cosa que nos inclina a pensar que tanto la población como el año están correladas. Para acabar con el estudio descriptivo de los fallecidos_interurbana, se dibuja: library(lattice) xyplot(fallecidos_interurbana~year|CCAA,col=1) year fallecidos_interurbana 0 50 100 150 200 250 200520102015 Andalucía Aragón 200520102015 Asturias Baleares 200520102015 Canarias Cantabria Cataluña Ceuta CLeón 0 50 100 150 200 250 CLMancha 0 50 100 150 200 250 CValenciana Extremadura Galicia Madrid Melilla Murcia 200520102015 Navarra PVasco 200520102015 0 50 100 150 200 250 Rioja Figura 7.4: Número de fallecidos en vías interurbanas por año en cada Comunidad Autónoma En esta gráfica se muestra por cada comunidad autónoma como el número de fallecidos en vías interurbanas desciende a lo largo de los años, pudiendose notar aquellas comunidades donde hay más y menos provincias. En segundo lugar, realizamos el estudio descriptivo con los fallecidos en vías urbanas. Se estudian las relaciones existentes entre fallecidos_urbana y las demás variables (menos fallecidos_interurbana) con el siguiente gráfico: pairs(fallecidos_urbana~year+numero_vehiculos+poblacion+precipitacion_media+ temperatura_media+horas_sol_media+km_carretera,panel=panel.smooth) 42 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS fallecidos_urbana 2002 2012 0e+00 6e+06 12 18 0 4000 0 150 2002 year numero_vehiculos 0e+00 0e+00 poblacion precipitacion_media 20 12 22 temperatura_media horas_sol_media 150 0 100 0 6000 0e+00 4e+06 20 80 150 250 km_carretera Figura 7.5: Figura con todos los pares de variables en vías urbanas El número de fallecidos en vías urbanas aumenta con la población y el número de vehículos, la interacción con las demás variables no es de destacar. Al igual que ocurría con los fallecidos en vías interurbanas, aquellas gráficas donde hay huecos es debido a las dos provincias con más habitantes (Madrid y Barcelona). Además, viendo el cuadro referente a la población y el número de vehículos, podemos intuir que ambas están altamente correladas, cosa que puede afectar a nuestro estudio. Se dibujan, a modo de ejemplo, 2 gráficas referentes a la provincia de Sevilla. ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=year, y=fallecidos_urbana)) + geom_point(size=1)+geom_line() CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 43 7.2. APLICACIÓN A NUESTRA BASE DE DATOS 10 15 20 25 30 35 2004 2008 2012 2016 year fallecidos_urbana Figura 7.6: Número de fallecidos en vías urbanas en Sevilla por año En la gráfica se ve como el número de fallecidos, excepto por un par de años donde ha habido un pico con mayor número de fallecidos, ha ido descendiendo desde estar entorno a los 33 fallecidos por año hasta los 12. ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=poblacion, y=fallecidos_urbana)) + geom_point(size=1)+geom_line() 10 15 20 25 30 35 1700000 1750000 1800000 1850000 1900000 1950000 poblacion fallecidos_urbana Figura 7.7: Número de fallecidos en vías urbanas en Sevilla por población 44 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS Esta gráfica es muy similar a la anterior, cosa que nos inclina a pensar que tanto la población como el año están correladas. Para acabar con el estudio descriptivo de los fallecidos_urbana, se dibuja: xyplot(fallecidos_urbana~year|CCAA,col=1) year fallecidos_urbana 0 50 100 150 200520102015 Andalucía Aragón 200520102015 Asturias Baleares 200520102015 Canarias Cantabria Cataluña Ceuta CLeón 0 50 100 150 CLMancha 0 50 100 150 CValenciana Extremadura Galicia Madrid Melilla Murcia 200520102015 Navarra PVasco 200520102015 0 50 100 150 Rioja Figura 7.8: Número de fallecidos en vías urbanas por año en cada Comunidad Autónoma En esta gráfica se muestra por cada comunidad autónoma como el número de fallecidos en vías urbanas desciende a lo largo de los años, pudiendose notar aquellas comunidades donde hay más y menos provincias. 7.2.3. Aplicación del modelo de Poisson Comenzamos el estudio con el modelo de Poisson. Se ha ultilizado las siguientes variables explicativas: year, numero_vehiculos, poblacion, precipitacion_media, temperatura_media, horas_sol_media, km_carretera, CCAA. 7.2.3.1. Vías interurbanas modeloTotal_i_p1<-glm(fallecidos_interurbana~year+ numero_vehiculos+poblacion+precipitacion_media+ temperatura_media+horas_sol_media+ km_carretera+CCAA,family="poisson",data=datos_f) summary(modeloTotal_i_p1) ## ## Call: CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 45 7.2. APLICACIÓN A NUESTRA BASE DE DATOS ## km_carretera 1.132e-04 4.006e-05 2.826 0.004710 ** ## CCAAAragón 3.920e-01 1.858e-01 2.110 0.034894 * ## CCAAAsturias 4.126e-01 3.521e-01 1.172 0.241276 ## CCAABaleares 3.975e-01 2.525e-01 1.574 0.115384 ## CCAACanarias -1.995e-01 2.894e-01 -0.689 0.490538 ## CCAACantabria 3.676e-01 3.979e-01 0.924 0.355537 ## CCAACataluña 6.449e-01 1.548e-01 4.166 3.10e-05 *** ## CCAACeuta 3.359e-01 2.774e-01 1.211 0.225907 ## CCAACLeón 9.005e-01 2.521e-01 3.572 0.000355 *** ## CCAACLMancha 1.506e-01 1.760e-01 0.855 0.392291 ## CCAACValenciana 1.885e-01 1.605e-01 1.175 0.239980 ## CCAAExtremadura -3.054e-02 2.010e-01 -0.152 0.879215 ## CCAAGalicia 3.573e-01 2.669e-01 1.338 0.180746 ## CCAAMadrid 2.217e-01 2.568e-01 0.863 0.387953 ## CCAAMelilla 8.206e-01 2.377e-01 3.452 0.000556 *** ## CCAAMurcia -2.309e-02 2.518e-01 -0.092 0.926939 ## CCAANavarra 2.352e-01 3.230e-01 0.728 0.466414 ## CCAAPVasco 5.261e-01 3.530e-01 1.490 0.136177 ## CCAARioja 1.063e+00 2.555e-01 4.161 3.17e-05 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 587.86 on 779 degrees of freedom ## Residual deviance: 397.59 on 755 degrees of freedom ## AIC: 2101.8 ## ## Number of Fisher Scoring iterations: 5 Como podemos observar, un aumento de la variable km_carretera provoca un aumento del número esperado de fallecidos en vías urbanas por cien mil habitantes; mientras que ocurre lo contrario con las variable year, hecho que se puede explicar gracias a que con los años las personas se van concienciando más, hay más avances. Es de destacar un dato bastante relevante y es que Andalucía tiene un número esperado de fallecidos en vías urbana por cien mil habitantes menor que prácticamente todas las comunidades autónomas, excepto en Extremadura, Canarias y Murcia, y estas no son influyentes en nuestro modelo. Comprobamos que no hay gran correlación entre las variables: vif(modeloTotal_u_p) ## GVIF Df GVIF^(1/(2*Df)) ## year 1.554399 1 1.246755 ## numero_vehiculos_milhabr 2.715358 1 1.647834 ## precipitacion_media 4.338545 1 2.082917 ## temperatura_media 15.222225 1 3.901567 ## horas_sol_media 5.578642 1 2.361915 ## km_carretera 3.194961 1 1.787445 52 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS ## CCAA 567.294223 18 1.192600 Mostramos que nuestro modelo se puede reducir más, en el sentido de disminuir el número de variables del estudio, mediante la función: modeloTotal_u_p_reducido=step(modeloTotal_u_p) ## Start: AIC=2101.8 ## fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + precipitacion_media + ## temperatura_media + horas_sol_media + km_carretera + CCAA ## ## Df Deviance AIC ## - precipitacion_media 1 397.59 2099.8 ## - horas_sol_media 1 397.80 2100.0 ## <none> 397.59 2101.8 ## - temperatura_media 1 399.83 2102.1 ## - numero_vehiculos_milhabr 1 401.39 2103.6 ## - km_carretera 1 405.37 2107.6 ## - CCAA 18 457.42 2125.6 ## - year 1 479.78 2182.0 ## ## Step: AIC=2099.81 ## fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + temperatura_media + ## horas_sol_media + km_carretera + CCAA ## ## Df Deviance AIC ## - horas_sol_media 1 397.81 2098.0 ## <none> 397.59 2099.8 ## - temperatura_media 1 399.87 2100.1 ## - numero_vehiculos_milhabr 1 401.39 2101.6 ## - km_carretera 1 405.37 2105.6 ## - CCAA 18 461.23 2127.4 ## - year 1 479.79 2180.0 ## ## Step: AIC=2098.03 ## fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + temperatura_media + ## km_carretera + CCAA ## ## Df Deviance AIC ## <none> 397.81 2098.0 ## - temperatura_media 1 400.82 2099.0 ## - numero_vehiculos_milhabr 1 401.81 2100.0 ## - km_carretera 1 405.48 2103.7 ## - CCAA 18 465.28 2129.5 ## - year 1 480.72 2178.9 Reduzcamos nuestro modelo: modeloTotal_u_p_reducido<-glm(fallecidos_urb_milhabr~year+ numero_vehiculos_milhabr+temperatura_media+ km_carretera+CCAA,family="poisson",data=datos_f) CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 53 7.2. APLICACIÓN A NUESTRA BASE DE DATOS summary(modeloTotal_u_p_reducido) ## ## Call: ## glm(formula = fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + ## temperatura_media + km_carretera + CCAA, family = "poisson", ## data = datos_f) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.26563 -0.41525 -0.01642 0.32266 2.58170 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.525e+02 1.696e+01 8.993 < 2e-16 *** ## year -7.699e-02 8.632e-03 -8.920 < 2e-16 *** ## numero_vehiculos_milhabr 1.214e-05 6.060e-06 2.003 0.045190 * ## temperatura_media 6.375e-02 3.688e-02 1.729 0.083842 . ## km_carretera 1.124e-04 4.007e-05 2.804 0.005040 ** ## CCAAAragón 3.952e-01 1.858e-01 2.128 0.033362 * ## CCAAAsturias 3.405e-01 3.003e-01 1.134 0.256804 ## CCAABaleares 3.667e-01 2.427e-01 1.511 0.130779 ## CCAACanarias -2.531e-01 2.418e-01 -1.047 0.295301 ## CCAACantabria 2.534e-01 2.918e-01 0.868 0.385327 ## CCAACataluña 6.138e-01 1.404e-01 4.371 1.23e-05 *** ## CCAACeuta 3.052e-01 2.692e-01 1.134 0.256909 ## CCAACLeón 9.017e-01 2.517e-01 3.582 0.000341 *** ## CCAACLMancha 1.523e-01 1.763e-01 0.864 0.387826 ## CCAACValenciana 1.819e-01 1.547e-01 1.176 0.239637 ## CCAAExtremadura -2.365e-02 2.002e-01 -0.118 0.905953 ## CCAAGalicia 3.018e-01 2.139e-01 1.411 0.158296 ## CCAAMadrid 2.182e-01 2.562e-01 0.852 0.394442 ## CCAAMelilla 7.959e-01 2.295e-01 3.468 0.000525 *** ## CCAAMurcia -1.856e-02 2.479e-01 -0.075 0.940319 ## CCAANavarra 1.882e-01 2.972e-01 0.633 0.526501 ## CCAAPVasco 4.434e-01 2.754e-01 1.610 0.107376 ## CCAARioja 1.039e+00 2.501e-01 4.155 3.25e-05 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 587.86 on 779 degrees of freedom ## Residual deviance: 397.81 on 757 degrees of freedom ## AIC: 2098 ## ## Number of Fisher Scoring iterations: 5 54 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS Observamos que la interpretación coincide con la del modelo sin reducir, aunque con coeficientes distintos lógicamente. Ahora sí, prácticamente todas las variables son influyentes mientras que antes solo lo eran 2. Aunque luego lo comprobaremos, podemos estudiar si hay sobredispersión viendo si el cociente entre el estadístico de desviación (397.59) y su grados de libertad (755) es mayor que uno, cosa que no se cumple. Dibujamos 2 gráficas relativas a los residuos y observaciones. plot(modeloTotal_u_p_reducido,which=c(1,4)) −0.5 0.0 0.5 1.0 −2 −1 0 1 2 3 Predicted values Residuals glm(fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + temperatura_ ... Residuals vs Fitted 727 269 302 Figura 7.11: Residuos y observaciones influyentes para fallecidos en vías urbanas con Poisson CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 55 7.2. APLICACIÓN A NUESTRA BASE DE DATOS 0 200 400 600 800 0.000 0.010 0.020 0.030 Obs. number Cook's distance glm(fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + temperatura_ ... Cook's distance 727 260 51 Figura 7.12: Residuos y observaciones influyentes para fallecidos en vías urbanas con Poisson En la primera gráfica podemos ver como nuestro modelo aproxima bien a valores entre 0 y 1 fallecidos, pero no ocurre lo mismo una vez que aumenta este número, puesto que como la varianza es mayor, los datos aparecen más dispersos, y de ahí un mayor residuo, lo cual es coherente con el modelo de Poisson. En la segunda podemos ver que hay 3 observaciones influyentes en nuestro modelo, pero sobre todo la 727, observación que corresponde a Ceuta en el año 2015 con 3 fallecidos en vías urbanas. Con la siguiente función vemos si se cumple la hipótesis de la equidispersión, comprobando así también que no hay sobredispersión de los datos mediante la función: dispersiontest(modeloTotal_u_p_reducido,trafo=1) ## ## Overdispersion test ## ## data: modeloTotal_u_p_reducido ## z = -19.116, p-value = 1 ## alternative hypothesis: true alpha is greater than 0 ## sample estimates: ## alpha ## -0.5706904 Nos da un p-valor de uno, por lo que se acepta la hipótesis de la equidispersión, y por tanto aceptamos el estudio mediante el modelo de Poisson. 56 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS 7.2.4. Aplicación del modelo Binomial Negativa Comenzamos el estudio con el modelo Binomial Negativo como miembro de la familia del MLG. Hemos ultilizado las siguientes variables explicativas: year, numero_vehiculos_milhabr, precipitacion_media, temperatura_media, horas_sol_media, km_carretera, CCAA. 7.2.4.1. Vías interurbanas library(MASS) modeloTotal_i_bn<-glm.nb(fallecidos_int_milhabr~year+ numero_vehiculos_milhabr+precipitacion_media+ temperatura_media+horas_sol_media+ km_carretera+CCAA,data=datos_f) summary(modeloTotal_i_bn) ## ## Call: ## glm.nb(formula = fallecidos_int_milhabr ~ year + numero_vehiculos_milhabr + ## precipitacion_media + temperatura_media + horas_sol_media + ## km_carretera + CCAA, data = datos_f, init.theta = 87.47198285, ## link = log) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -3.8288 -0.5985 -0.0789 0.4125 4.1309 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.849e+02 8.437e+00 33.766 < 2e-16 *** ## year -1.412e-01 4.285e-03 -32.947 < 2e-16 *** ## numero_vehiculos_milhabr 2.645e-05 2.665e-06 9.924 < 2e-16 *** ## precipitacion_media -1.826e-03 1.264e-03 -1.445 0.14835 ## temperatura_media -1.249e-01 1.706e-02 -7.321 2.47e-13 *** ## horas_sol_media 4.464e-03 9.276e-04 4.812 1.49e-06 *** ## km_carretera -3.272e-05 1.797e-05 -1.821 0.06865 . ## CCAAAragón 2.028e-01 8.055e-02 2.517 0.01182 * ## CCAAAsturias -3.285e-02 1.656e-01 -0.198 0.84276 ## CCAABaleares -4.211e-01 1.280e-01 -3.290 0.00100 ** ## CCAACanarias -9.556e-02 1.403e-01 -0.681 0.49582 ## CCAACantabria -7.604e-03 1.820e-01 -0.042 0.96667 ## CCAACataluña 8.012e-02 7.139e-02 1.122 0.26178 ## CCAACeuta -3.515e+01 2.475e+06 0.000 0.99999 ## CCAACLeón 5.080e-03 1.065e-01 0.048 0.96197 ## CCAACLMancha 2.010e-01 7.216e-02 2.785 0.00536 ** ## CCAACValenciana -7.815e-02 7.672e-02 -1.019 0.30836 ## CCAAExtremadura 1.089e-01 8.513e-02 1.280 0.20063 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 57 7.2. APLICACIÓN A NUESTRA BASE DE DATOS ## CCAAGalicia 2.131e-01 1.133e-01 1.882 0.05986 . ## CCAAMadrid -1.486e+00 1.900e-01 -7.823 5.14e-15 *** ## CCAAMelilla -3.512e+01 2.517e+06 0.000 0.99999 ## CCAAMurcia 9.972e-03 1.203e-01 0.083 0.93392 ## CCAANavarra -1.116e-01 1.364e-01 -0.819 0.41303 ## CCAAPVasco -3.451e-01 1.499e-01 -2.303 0.02130 * ## CCAARioja 2.217e-01 1.177e-01 1.883 0.05975 . ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(87.472) family taken to be 1) ## ## Null deviance: 3426.70 on 779 degrees of freedom ## Residual deviance: 571.68 on 755 degrees of freedom ## AIC: 3503.6 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 87.5 ## Std. Err.: 37.2 ## ## 2 x log-likelihood: -3451.636 Los coeficientes de este modelo son practicamente iguales que los del modelo de Poisson excepto los referidos a las comunidades autónomas, de los cuáles las más significativas son Baleares, Madrid y País Vasco con menos fallecidos por cien mil habitantes y Aragón y Castilla La Mancha con más (todos respecto a Andalucía). Comprobamos que no hay problemas destacables de las correlaciones entre variables. vif(modeloTotal_i_bn) ## GVIF Df GVIF^(1/(2*Df)) ## year 1.650018 1 1.284530 ## numero_vehiculos_milhabr 2.628555 1 1.621282 ## precipitacion_media 3.846280 1 1.961194 ## temperatura_media 12.101229 1 3.478682 ## horas_sol_media 5.309612 1 2.304259 ## km_carretera 2.129291 1 1.459209 ## CCAA 250.812237 18 1.165866 Comprobemos que nuestro modelo no se puede reducir más, en el sentido de disminuir el número de variables del estudio, mediante la función: modeloTotal_i_bn_reducido=step(modeloTotal_i_bn) ## Start: AIC=3501.64 ## fallecidos_int_milhabr ~ year + numero_vehiculos_milhabr + precipitacion_media + ## temperatura_media + horas_sol_media + km_carretera + CCAA ## ## Df Deviance AIC 58 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS ## <none> 571.68 3501.6 ## - precipitacion_media 1 573.78 3501.7 ## - km_carretera 1 574.99 3502.9 ## - horas_sol_media 1 594.90 3522.8 ## - temperatura_media 1 625.26 3553.2 ## - numero_vehiculos_milhabr 1 668.08 3596.0 ## - CCAA 18 1079.63 3973.6 ## - year 1 1762.78 4690.7 Por último, se compara el modelo de Poisson y el Binomial Negativo. Para ello se usó la siguiente prueba de razón de verosimilitud: tt=2*(logLik(modeloTotal_i_bn)-logLik(modeloTotal_i_p)) pchisq(tt,df=1,lower.tail=FALSE) ## 'log Lik.'0.008846886 (df=26) El hecho que salga un valor del χ2 tan próximo a cero nos indica que es más apropiado utilizar para nuestro estudio el modelo Binomial Negativo que el modelo de Poisson con el mismo número de variables. 7.2.4.2. Vías urbanas library(MASS) modeloTotal_u_bn<-glm.nb(fallecidos_urb_milhabr~year+ numero_vehiculos_milhabr+precipitacion_media+ temperatura_media+horas_sol_media+ km_carretera+CCAA,data=datos_f) summary(modeloTotal_u_bn) ## ## Call: ## glm.nb(formula = fallecidos_urb_milhabr ~ year + numero_vehiculos_milhabr + ## precipitacion_media + temperatura_media + horas_sol_media + ## km_carretera + CCAA, data = datos_f, init.theta = 55483.63471, ## link = log) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.24275 -0.41636 -0.01771 0.32226 2.57908 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.520e+02 1.700e+01 8.939 < 2e-16 *** ## year -7.678e-02 8.643e-03 -8.884 < 2e-16 *** ## numero_vehiculos_milhabr 1.189e-05 6.083e-06 1.954 0.050709 . ## precipitacion_media -1.264e-04 2.839e-03 -0.045 0.964478 ## temperatura_media 5.845e-02 3.908e-02 1.496 0.134712 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 59 7.2. APLICACIÓN A NUESTRA BASE DE DATOS ## horas_sol_media 9.457e-04 2.024e-03 0.467 0.640348 ## km_carretera 1.132e-04 4.006e-05 2.826 0.004711 ** ## CCAAAragón 3.920e-01 1.858e-01 2.110 0.034897 * ## CCAAAsturias 4.126e-01 3.521e-01 1.172 0.241285 ## CCAABaleares 3.975e-01 2.525e-01 1.574 0.115393 ## CCAACanarias -1.995e-01 2.894e-01 -0.689 0.490542 ## CCAACantabria 3.676e-01 3.979e-01 0.924 0.355549 ## CCAACataluña 6.449e-01 1.548e-01 4.166 3.10e-05 *** ## CCAACeuta 3.359e-01 2.774e-01 1.211 0.225917 ## CCAACLeón 9.005e-01 2.521e-01 3.572 0.000355 *** ## CCAACLMancha 1.506e-01 1.760e-01 0.855 0.392295 ## CCAACValenciana 1.885e-01 1.605e-01 1.175 0.239985 ## CCAAExtremadura -3.054e-02 2.010e-01 -0.152 0.879218 ## CCAAGalicia 3.573e-01 2.669e-01 1.338 0.180755 ## CCAAMadrid 2.217e-01 2.568e-01 0.863 0.387960 ## CCAAMelilla 8.206e-01 2.377e-01 3.452 0.000556 *** ## CCAAMurcia -2.309e-02 2.518e-01 -0.092 0.926941 ## CCAANavarra 2.352e-01 3.230e-01 0.728 0.466420 ## CCAAPVasco 5.261e-01 3.531e-01 1.490 0.136186 ## CCAARioja 1.063e+00 2.555e-01 4.161 3.17e-05 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## (Dispersion parameter for Negative Binomial(55483.63) family taken to be 1) ## ## Null deviance: 587.85 on 779 degrees of freedom ## Residual deviance: 397.58 on 755 degrees of freedom ## AIC: 2103.8 ## ## Number of Fisher Scoring iterations: 1 ## ## ## Theta: 55484 ## Std. Err.: 276884 ## Warning while fitting theta: iteration limit reached ## ## 2 x log-likelihood: -2051.817 library(pscl) vuong(modeloTotal_u_p_reducido,modeloTotal_u_bn) ## Vuong Non-Nested Hypothesis Test-Statistic: ## (test-statistic is asymptotically distributed N(0,1) under the ## null that the models are indistinguishible) ## ------------------------------------------------------------- ## Vuong z-statistic H_A p-value ## Raw -0.3406407 model2 > model1 0.36669 ## AIC-corrected 5.9585734 model1 > model2 1.2722e-09 ## BIC-corrected 20.6335184 model1 > model2 < 2.22e-16 60 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 7.2. APLICACIÓN A NUESTRA BASE DE DATOS Concluimos que el modelo de Poisson y el Binomial Negativo son equivalentes (porque el p-valor es muy pequeño tanto en el estadístico AIC como el BIC), por eso hemos descartado este estudio para nuestro modelo. 7.2.5. Aplicación del modelo de Poisson Bivariante Procedemos a aplicar el modelo de Poisson Bivariante a través de la librería bivpois, que actualmente no se encuentra disponible en el software R pero se adjunta el script con las funciones en el anexo. Para ello hacemos uso de la función lm.bp, cuyos argumentos más importantes son: lm.bp(l1, l2, l1l2 = NULL, l3=1, data, common.intercept =FALSE, zeroL3 = FALSE) donde: 1. l1: fórmula de la forma x~X1+... +Xppara los parámetros del logλ1. 2. l2: fórmula de la forma y~X1+... +Xppara los parámetros del logλ2. 3. l1l2: fórmula de la forma ~ X1 + ... + Xp para los parámetros comunes del logλ1 y logλ2 . Si la variable explicativa se encuentra también en l1 y/o l2 entonces el modelo ajusta la interacción entre los parámetros. Se pueden usar aquí términos especiales de la forma c ( X1, X2 ). Estos términos implican parámetros comunes de λ1 y λ2 para distintas variables. 4. l3: fórmula de la forma ~ X1+... +Xppara los parámetros del logλ3. 5. data: Base de datos (data frame) que contiene las variables en el modelo. 6. common.intercept: función lógica que específica si se debe usar un término independiente común entre λ1ylambda2. Por defecto su valor es FALSE. 7. zeroL3: argumento lógico que controla si λ3 debería ser igual a cero (y por lo tanto se ajustaría el modelo de Poisson doble). En nuestro caso consideraremos x=fallecidos_int_milhabr, y=fallecidos_urb_milhabr. Nuestros datos vienen dados del siguiente data frame: cc=cbind(cbind(datos_f$year,datos_f[,7:10]),datos_f[,15:17]) head(cc) ## datos_f$year precipitacion_media temperatura_media horas_sol_media ## 1 2002 44.375 11.808 169.133 ## 2 2002 33.150 14.842 261.133 ## 3 2002 20.933 18.342 265.158 ## 4 2002 17.058 18.583 268.508 ## 5 2002 28.358 11.258 226.983 ## 6 2002 26.525 16.917 250.375 ## km_carretera fallecidos_int_milhabr fallecidos_urb_milhabr ## 1 1463.33 14 3 ## 2 3710.45 17 2 ## 3 2704.29 10 1 ## 4 2368.19 19 2 ## 5 2545.57 28 2 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO 61 7.3. CONCLUSIÓN 7.3.3. Poisson Bivariante A la vista de los resultados, prácticamente todos los coeficientes para el modelo fallecidos_int_milhabr son mayores que los del modelo fallecidos_urb_milhabr, esto quiere decir que un aumento de las variables cuyos coeficientes son positivos provocan un mayor incremento en el número esperado de fallecidos en vías interurbanas por cien mil habitantes que en vías urbanas, mientras que un aumento de las variables cuyos coeficientes son negativos provocan un mayor decrecimiento en el número esperado de fallecidos en vías urbanas por cien mil habitantes que en vías urbanas. El tercer modelo nos indica la relación entre los fallecidos en vías interurbanas y los fallecidos en vías urbanas, el cual es constante. 68 CAPÍTULO 7. APLICACIÓN A LOS ACCIDENTES DE TRÁFICO Apéndice A Apéndice: Código del trabajo en R library(readxl) library(dplyr) library(lmtest) library(MASS) library(sandwich) library(AER) library(car) library(pscl) #Lectura de la base de datos load(file="datos.RData") #Estudio previo summary(datos_f) mean(datos_f$fallecidos_interurbana) var(datos_f$fallecidos_interurbana) mean(datos_f$fallecidos_urbana) var(datos_f$fallecidos_urbana) q=c(0) w=c(0) e=c(0) l=c(0) for(i in (1:52)) { q=c(q,mean(datos_f[datos_f$provincia==provincias[i],]$fallecidos_interurbana)) w=c(w,mean(datos_f[datos_f$provincia==provincias[i],]$fallecidos_urbana)) e=c(e,mean(datos_f[datos_f$provincia==provincias[i],]$numero_vehiculos)) l=c(l,mean(datos_f[datos_f$provincia==provincias[i],]$poblacion)) 69 } q1=c(0) w1=c(0) for(i in (1:52)) { q1=c(q1,sum(datos_f[datos_f$provincia==provincias[i],]$fallecidos_interurbana)) w1=c(w1,sum(datos_f[datos_f$provincia==provincias[i],]$fallecidos_urbana)) } fallecidos_interurbana_total=cbind(provincias,q1[2:53]) fallecidos_interurbana_media=cbind(provincias,q[2:53]) fallecidos_urbana_total=cbind(provincias,w1[2:53]) fallecidos_urbana_media=cbind(provincias,w[2:53]) numero_vehiculos_medio=cbind(provincias,e[2:53]) poblacion_media=cbind(provincias,l[2:53]) heatmap(cor(datos_f[,-c(1,11)]),scale="none") #Estudio descriptivo attach(datos_f) ## interurbana pairs(fallecidos_interurbana~year+numero_vehiculos+poblacion+precipitacion_media+ temperatura_media+horas_sol_media+km_carretera,panel=panel.smooth) library(ggplot2) ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=year, y=fallecidos_interurbana)) + geom_point(size=1)+ geom_line() ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=poblacion, y=fallecidos_interurbana)) + geom_point(size=1)+geom_line() ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=numero_vehiculos, y=fallecidos_interurbana)) + geom_point(size=1)+ geom_line() library(lattice) xyplot(fallecidos_interurbana~year|CCAA,col=1) xyplot(fallecidos_interurbana~poblacion|year,col=1) ## urbana pairs(fallecidos_urbana~year+numero_vehiculos+poblacion+precipitacion_media+ temperatura_media+horas_sol_media+km_carretera,panel=panel.smooth) 70 APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=year, y=fallecidos_urbana)) + geom_point(size=1)+ geom_line() ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=poblacion, y=fallecidos_urbana)) + geom_point(size=1)+ geom_line() ggplot( datos_f[datos_f$provincia=="Sevilla",], aes(x=numero_vehiculos, y=fallecidos_urbana)) + geom_point(size=1)+ geom_line() xyplot(fallecidos_urbana~year|CCAA,col=1) #vemos como en algunas comunidades como catalunya destaca barcelona pero #las demas provincias se mantienen con un numero mas bajo de fallecidos por año. xyplot(fallecidos_urbana~poblacion|year,col=1) #vemos como las cuidades con menor poblacion tienen un menor numero #de fallecidos urbanos y por eso se concentran por debajo de los 35-40, #mientras tanto las que tienen mas poblacion como barcelona y madrid #tienen muchos fallecidos siempre mas de 50 y llegando hasta los 150. #Estudio de los modelos # INTERURBANA #Modelo con todas las variables para fallecidos_interurbana con glm modeloTotal_i_p1<-glm(fallecidos_interurbana~year+numero_vehiculos+ poblacion+precipitacion_media+temperatura_media+horas_sol_media+ km_carretera+CCAA,family="poisson",data=datos_f) summary(modeloTotal_i_p1) # coefconfi<-cbind(Estimate=coef(modeloTotal_i_p1), confint(modeloTotal_i_p1)) # exp(coefconfi) #colinealidad vif(modeloTotal_i_p1) #Como poblacion y numero de vehiculos estan muy correladas, vamos a crear #las siguientes nuevas variables para asi"ocultar" la variable poblacion datos_f$fallecidos_int_milhab<- 100000*datos_f$fallecidos_interurbana/datos_f$poblacion datos_f$fallecidos_urb_milhab<- 100000*datos_f$fallecidos_urbana/datos_f$poblacion datos_f$numero_vehiculos_milhab<- 100000*datos_f$numero_vehiculos/datos_f$poblacion # #Modelo con todas las variables para fallecidos_int_milhab con glm # APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R 71 # modeloTotal_i2<-glm(fallecidos_int_milhab~year+numero_vehiculos_milhab+ # precipitacion_media+temperatura_media+horas_sol_media+km_carretera+ # CCAA,data=datos_f) # # summary(modeloTotal_i2) # # modeloTotal_i2$residuals # # coefconfi2<-cbind(Estimate=coef(modeloTotal_i2), confint(modeloTotal_i2)) # # exp(coefconfi2) # # coeftest(modeloTotal_i2,vcov=sandwich) # # #colinealidad # vif(modeloTotal_i2) #Ya esta solucionado el problema que teniamos antes # de la correlacion entre las variables # # #Reduccion del modelo # modeloTotal_i2_reducido=step(modeloTotal_i2) #Modelo con todas las variables para fallecidos_int_milhabr con poisson datos_f$fallecidos_int_milhabr<-round(datos_f$fallecidos_int_milhab) datos_f$fallecidos_urb_milhabr<-round(datos_f$fallecidos_urb_milhab) datos_f$numero_vehiculos_milhabr<-round(datos_f$numero_vehiculos_milhab) save(datos_f,file="datos2.RData") modeloTotal_i_p<-glm(fallecidos_int_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+CCAA, family="poisson",data=datos_f) modeloTotal_i_p22<-glm(fallecidos_int_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera, family="poisson",data=datos_f) summary(modeloTotal_i_p) # modeloTotal_i_p$residuals # coefconfip<-cbind(Estimate=coef(modeloTotal_i_p), confint(modeloTotal_i_p)) # exp(coefconfip) # coeftest(modeloTotal_i_p,vcov=sandwich) #Graficas plot(modeloTotal_i_p,which=c(1,4)) #colinealidad vif(modeloTotal_i_p) #Reduccion del modelo modeloTotal_i_p_reducido=step(modeloTotal_i_p) 72 APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R #dispersion dispersiontest(modeloTotal_i_p,trafo=1) #Como el pvalor es casi 1, se acepta la hipotesis de que la media y la varianza #sean iguales y por tanto aceptaríamos el estudio con poisson #Modelo con todas las variables para fallecidos_int_milhabr con binomial negativa modeloTotal_i_bn<-glm.nb(fallecidos_int_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+ CCAA,data=datos_f) summary(modeloTotal_i_bn) # modeloTotal_i_bn$residuals # coefconfibn<-cbind(Estimate=coef(modeloTotal_i_bn)) # coefconfibn #colinealidad vif(modeloTotal_i_bn) #Reduccion de modelo modeloTotal_i_bn_reducido=step(modeloTotal_i_bn) #comparacion poisson-bn (se tienen que hacer con dos modelos con las misma #variables) tt=2*(logLik(modeloTotal_i_bn)-logLik(modeloTotal_i_p)) pchisq(tt,df=1,lower.tail=FALSE) #el hecho de que salga un numero proximo #a cero nos dice que es mejor usar el modelo binomial negativo que el de poisson #Como no tenemos exceso de ceros, los dos siguientes modelos no van #a poderse realizar #Modelo con todas las variables para fallecidos_int_milhabr con hurdle modeloTotal_i_h<-hurdle(fallecidos_int_milhabr~numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+year+ CCAA,data=datos_f) #Modelo con todas las variables para fallecidos_int_milhabr con cero inflado modeloTotal_i_z<-zeroinfl(formula=fallecidos_int_milhabr~numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+CCAA|year, data=datos_f) # URBANA #Modelo con todas las variables para fallecidos_urb_milhab con glm modeloTotal_u<-glm(fallecidos_urb_milhab~year+numero_vehiculos_milhab+ APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R 73 precipitacion_media+temperatura_media+horas_sol_media+km_carretera+CCAA, data=datos_f) summary(modeloTotal_u) # modeloTotal_u$residuals # coefconfu<-cbind(Estimate=coef(modeloTotal_u), confint(modeloTotal_u)) # exp(coefconfu) #colinealidad vif(modeloTotal_u) #Reduccion del modelo modeloTotal_u_reducido=step(modeloTotal_u) modeloTotal_u_reducido<-glm(fallecidos_urb_milhab~year+numero_vehiculos_milhab+ temperatura_media+km_carretera+CCAA,data=datos_f) summary(modeloTotal_u_reducido) # modeloTotal_u_reducido$residuals # coefconfured<-cbind(Estimate=coef(modeloTotal_u_reducido), # confint(modeloTotal_u_reducido)) # exp(coefconfured) #Modelo con todas las variables para fallecidos_urb_milhabr con poisson modeloTotal_u_p<-glm(fallecidos_urb_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+CCAA, family="poisson",data=datos_f) modeloTotal_u_p22<-glm(fallecidos_urb_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera, family="poisson",data=datos_f) summary(modeloTotal_u_p) # modeloTotal_u_p$residuals # coefconfup<-cbind(Estimate=coef(modeloTotal_u_p), # confint(modeloTotal_u_p)) # exp(coefconfup) # coeftest(modeloTotal_u_p, vcov = sandwich) #Graficas plot(modeloTotal_u_p,which=c(1,4)) #colinealidad vif(modeloTotal_u_p) #Modelo reducido modeloTotal_u_p_reducido=step(modeloTotal_u_p) 74 APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R modeloTotal_u_p_reducido<-glm(fallecidos_urb_milhabr~year+ numero_vehiculos_milhabr+temperatura_media+km_carretera+CCAA, family="poisson",data=datos_f) summary(modeloTotal_u_p_reducido) # modeloTotal_u_p_reducido$residuals # coefconfup<-cbind(Estimate=coef(modeloTotal_u_p_reducido), # confint(modeloTotal_u_p_reducido)) # exp(coefconfupred) # coeftest(modeloTotal_u_p_reducido, vcov = sandwich) #colinealidad vif(modeloTotal_u_p_reducido) #Graficas plot(modeloTotal_u_p_reducido,which=c(1,4)) #dispersion dispersiontest(modeloTotal_u_p_reducido,trafo=1) #Como el pvalor es 1, se acepta la hipotesis de que la media #y la varianza sean iguales por lo que se acepta el modelo de #Poisson para el estudio de nuestro modelo #Modelo con todas las variables para fallecidos_urb_milhabr con #binomial negativa modeloTotal_u_bn<-glm.nb(fallecidos_urb_milhabr~year+ numero_vehiculos_milhabr+precipitacion_media+temperatura_media+ horas_sol_media+km_carretera+CCAA,data=datos_f) #comparacion poisson-bn vuong(modeloTotal_u_p_reducido,modeloTotal_u_bn) #Como no tenemos exceso de ceros, los dos siguientes modelos #no van a poderse realizar #Modelo con todas las variables para fallecidos_urb_milhabr con hurdle modeloTotal_u_h<-hurdle(fallecidos_urb_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+CCAA, data=datos_f) #Modelo con todas las variables para fallecidos_urb_milhabr con cero inflado modeloTotal_u_z<-zeroinfl(fallecidos_urb_milhabr~year+numero_vehiculos_milhabr+ precipitacion_media+temperatura_media+horas_sol_media+km_carretera+CCAA, APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R 75 data=datos_f,dist = "negbin") #Modelo de Regresion Bivariante X=fallecidos_interurbana Y=fallecidos_urbana cc=cbind(cbind(datos_f$year,datos_f[,7:10]),datos_f[,15:17]) ModeloBivariante1<-lm.bp(fallecidos_int_milhabr~.,fallecidos_urb_milhabr~., data=cc) #Poisson Bivariante con lambda3=cte ModeloBivarianteBic1<-ModeloBivariante1$BIC ModeloBivariante1 ModeloBivarianteAic1<-ModeloBivariante1$AIC ModeloBivariante1$parameters ModeloBivariante2<-lm.bp(fallecidos_int_milhabr~.,fallecidos_urb_milhabr~., l3=~.,data=cc) #Poisson Bivariante ModeloBivarianteBic2<-ModeloBivariante2$BIC ModeloBivariante2 ModeloBivarianteAic2<-ModeloBivariante2$AIC ModeloBivariante2$parameters ModeloBivariante3<-lm.bp(fallecidos_int_milhabr~.,fallecidos_urb_milhabr~., data=cc,zeroL3 = TRUE) #Poisson doble ModeloBivarianteBic3<-ModeloBivariante3$BIC ModeloBivariante3 ModeloBivarianteAic3<-ModeloBivariante3$AIC ModeloBivariante3$parameters #Comparacion de modelos univariantes ## Vías interurbanas (mismo número de variables) tt=2*(logLik(modeloTotal_i_bn)-logLik(modeloTotal_i_p)) pchisq(tt,df=1,lower.tail=FALSE) ## Vías urbanas vuong(modeloTotal_u_p_reducido,modeloTotal_u_bn) 76 APÉNDICE A. APÉNDICE: CÓDIGO DEL TRABAJO EN R Apéndice B Apéndice: Paquete bivpois Aquí se encuentra recogido el código con todas las funciones del paquete bivpois que actualmente no está disponible dentro del software R. "lm.bp" <- function( l1, l2, l1l2=NULL, l3=~1, data, common.intercept=FALSE, zeroL3=FALSE, maxit=300, pres=1e-8, verbose=getOption('verbose') ) # # { options(warn=-1) # # definition of function call templist<-list( l1=l1, l2=l2, l1l2=l1l2, l3=l3, data=substitute(data), common.intercept=common.intercept, zeroL3=zeroL3, maxit=maxit, pres=pres, verbose=verbose) tempcall<-as.call( c(expression(lm.bp), templist)) rm(templist) # # # ------------------------------------------------------------------------------------ # Karlis and Ntzoufras (2003, 2004, 2005) # EM algorithms for Bivariate Poisson Models # ------------------------------------------------------------------------------------ # # l1 : formula for the first linear predictor (of lambda1) # l2 : formula for the second linear predictor (of lambda2) # l1l2 : formula for common variables on both lambda1 and lambda2 # l3 : formula for the third first linear predictor/covariance parameter (lambda3) # common.intercept: logical argument defining whether common intercept should be used for lamdba1,lambda2 # # data : data.frame which contains data {required arguement} # zeroL3 : Logical argument controlling whether lambda3 is zero (DblPoisson) or not 77 # # fit model on lambda3 internal.data2$y3<-s m0<-glm( formula2, family=poisson, data=data2 ) beta3<-m0$coef lambda3<-m0$fitted # # fit model on lambda1 & lambda2 internal.data1$y1y2<-c(x1,x2) m<-glm( formula1, family=poisson, data=data1 ) p3<-length(m$coef) beta<-m$coef # creating names for parameters names(beta)<-newnamesbeta( beta ) # # lambda<-fitted(m) lambda1<-lambda[1:n] lambda2<-lambda[(n+1):(2*n)] ##### end of M step ###### # # detailed or compressed printing during the EM iterations if (verbose) { printvector<-c( i, beta, beta3,loglike[i], difllike ) names(printvector)<-c( 'iter', names(beta), paste('(l3):', names(beta3),sep=''), 'loglike','Rel.Dif.loglike')} else { printvector<-c( i, loglike[i], difllike ) names(printvector)<-c( 'iter','loglike','Rel.Dif.loglike')} # lengthpvec<-length(printvector) print.default( printvector, digits=4 ) } # # calculation of BIC and AIC for bivpoisson model noparams<- m$rank + m0$rank AIC<- -2*loglike[i] + noparams * 2 BIC<- -2*loglike[i] + noparams * log(2*n) # # # Calculation of BIC, AIC of Poisson saturated model x.mean<-x x.mean[x==0]<-1e-12 y.mean<-y y.mean[y==0]<-1e-12 AIC.sat <- sum(log( dpois( x , x.mean ) ) + log( dpois(y,y.mean))) 84 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS BIC.sat <- -2 * AIC.sat + (2*n)* log(2*n) AIC.sat <- -2 * AIC.sat + (2*n)* 2 # # AICtotal<-c(AIC.sat, AIC); BICtotal<-c(BIC.sat, BIC ); names(AICtotal)<-c('Saturated','BivPois') names(BICtotal)<-c('Saturated','BivPois') # # spliting parameter vector betaparameters<-splitbeta( beta ) # # putting all betas in one vector allbeta<-c(betaparameters$beta1,betaparameters$beta2, beta3) names(allbeta)<-c( paste( '(l1):',names(betaparameters$beta1), sep=''),paste('(l2):', names(betaparameters$beta2), sep='' ), paste('(l3):', names(beta3), sep='' ) ) # # Calculation of output result<-list(coefficients=allbeta, fitted.values=data.frame(x=m$fitted[1:n]+ lambda3,y=m$fitted[(n+1):(2*n)]+lambda3), residuals=data.frame(x=x-m$fitted[1:n]-lambda3, y=y-m$fitted[(n+1):(2*n)]-lambda3), beta1=betaparameters$beta1, beta2=betaparameters$beta2, beta3=beta3, lambda1=m$fitted[1:n], lambda2=m$fitted[(n+1):(2*n)], lambda3=lambda3, loglikelihood=loglike[1:i], parameters=noparams, AIC=AICtotal, BIC=BICtotal,iterations=i, call=tempcall ) # # } # end of elseif options(warn=0) # class(result)<-c('lm.bp','lm') result # # } "bivpois.table" <- function(x, y, lambda = c(1, 1, 1)) { # ------------------------------------------------------------------------------ # Karlis and Ntzoufras (2003, 2004) # EM algorithms for Bivariate Poisson Models # ------------------------------------------------------------------------------ # x : 1st count variable APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 85 # y : 2nd count variable # lambda: parameters of the bivariate poisson distribution. # ------------------------------------------------------------------------------ j<-0 n <- length(x) maxy <- c(max(x), max(y)) #Set initial values for parameters lambda1 <- lambda[1] lambda2 <- lambda[2] lambda3 <- lambda[3] if((x == 0) | (y == 0)) { prob <- matrix(NA, nrow = maxy[1] + 1, ncol = maxy[2]+1, byrow = T) prob[maxy[1] + 1, maxy[2] + 1] <- exp( - lambda3) * dpois(x[j], lambda1[j]) * dpois(y[j], lambda2[j]) } else { prob <- matrix(NA, nrow = maxy[1] + 1, ncol = maxy[2]+1, byrow = T) k <- 1 m <- 1 prob[k, m] <- exp( - lambda1 - lambda2 - lambda3) for(i in 2:(maxy[1] + 1)) { prob[i, 1] <- (prob[i - 1, 1] * lambda1)/(i - 1) } for(j in 2:(maxy[2] + 1)) { prob[1, j] <- (prob[1, j - 1] * lambda2)/(j - 1) } for(j in 2:(maxy[2] + 1)) { for(i in 2:(maxy[1] + 1)) { prob[i, j] <- (lambda1 * prob[i - 1, j] + lambda3 * prob[i - 1, j - 1])/(i - 1) } } } result <- prob result } "lm.dibp" <- function ( l1, l2, l1l2=NULL, l3=~1, data, common.intercept=FALSE, zeroL3=FALSE, distribution='discrete', jmax=2,maxit=300, pres=1e-8, verbose=getOption('verbose') ) { options(warn=-1) # # definition of function call templist<-list( l1=l1, l2=l2, l1l2=l1l2, l3=l3, data=substitute(data), common.intercept=common.intercept, zeroL3=zeroL3, 86 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS distribution=distribution,jmax=jmax, maxit=maxit, pres=pres, verbose=verbose) tempcall<-as.call( c(expression(lm.dibp), templist)) rm(templist) # # # # ----------------------------------------------------------------------- # Karlis and Ntzoufras (2003, 2004, 2005) # EM algorithms for Bivariate Poisson Models # ----------------------------------------------------------------------- # # PARAMETERS COMMON WITH lm.bp # l1 : formula for the first linear predictor (of lambda1) # l2 : formula for the second linear predictor (of lambda2) # l1l2 : formula for common variables on both lambda1 and lambda2 # l3 : formula for the third first linear predictor/covariance parameter (lambda3) # common.intercept: logical argument defining whether common intercept should be used for lamdba1,lambda2 # # data : data.frame which contains data {required arguement} # zeroL3 : Logical argument controlling whether lambda3 is zero (DblPoisson) or not # maxit : maximum number of iterations # pres : precision of the relative likelihood difference after which EM stops # verbose : Logical argument controlling whether beta parameters will we # printed while EM runs. Default value is taken options()$verbose value. # # PARAMETERS ADDITIONAL TO lm.bp # distribution : Selection of diagonal inflation distribution. # Three choices are available: # ='discrete': Discrete, jmax is the number of diagonal elements [0,1,...,] # ='poisson': Poisson with mean theta. # ='geometrics': Geometric with success probability theta. # Default is DISCRETE(2). theta[1] and theta[2] stand for theta_1, theta_2 # while theta_0=1-theta[1]-theta[2]. # jmax : Used only for DISCRETE diagonal distribution (distribution='discrete'). # Indicates the number of parameters of the DISCRETE distribution. # ---------------------------------------------------------------------- # # APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 87 # # set common or noncommon intercept if (common.intercept){ formula1.terms<-'1'} else {formula1.terms<-'internal.data1$noncommon'} # # namex<-as.character(l1[2]) namey<-as.character(l2[2]) x<-data[,names(data)==namex] y<-data[,names(data)==namey] # # # Data length n<-length(x) lengthprintvec<-1 # # # # definition of diagonal inflated distribution maxy<-max(c(x,y)) # # changing distribution to codes 1,2,3 dist<-distribution if ( charmatch( dist, 'poisson', nomatch=0) ==1 ) {distribution<-2} else if ( charmatch( dist, 'geometric', nomatch=0) ==1 ) {distribution<-3} else if ( charmatch( dist, 'discrete', nomatch=0) ==1 ) {distribution<-1} if ( distribution==1 ){ dilabel<-paste('Inflation Distribution: Discrete with J=',jmax) if (jmax==0) {theta<-0} else { theta<-1:jmax*0+1/(jmax+1) } di.f<-function (x, theta){ JMAX<-length(theta) if (x>JMAX) { res<-0 } else if (x==0) { res<-1-sum(theta) } else { res<-theta[x] } res } } else if ( distribution==2 ){ dilabel<-'Inflation Distribution: Poisson' theta<-1.0; di.f<-function (x, theta){ if (theta>0) { res<-dpois( x, theta ) } else { if (x==0) { res<-1} else {res<-1e-12} } } 88 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS } else if ( distribution==3 ){ dilabel<-'Inflation Distribution: Geometric' theta<-0.5; di.f<-function (x, theta){ if (theta>0) { if(theta==1) {theta<-0.9999999} res<-dgeom( x, theta ) } else if (theta==1){ if (x==0) { res<-1} else {res<-1e-12} } else {res<-1e-12} } } else { stop(paste(distribution, 'Not known distribution.', sep=':')) } # ------ # setting up data frames, vectors and data # # form dataframes used # data1 includes modelling on lambda1 and lambda2 # data2 includes modelling on lambda3 # internal.data1 and internal.data2 are data frames used for additional internal variables # internal.data1<-data.frame( y1y2=c( x, y) ) internal.data2<-data.frame( y3 = rep(0, n ) ) # p<-length(as.data.frame(data)) data1<-rbind(data, data) names(data1)<-names(data) # # removing x and y data1<-data1[ , names(data1)!=namex] data1<-data1[ , names(data1)!=namey] # # # # define full model if (as.character(l1[3])=='.') {l1<-formula( paste(as.character(l1[2]), paste( names(data1),'',collapse='+',sep='' ), sep='~') ) } if (as.character(l2[3])=='.') {l2<-formula( paste(as.character(l2[2]), paste( names(data1),'',collapse='+',sep='' ), sep='~') ) } if (as.character(l3[2])=='.') {l3<-formula( paste( '', paste( names(data1),'',collapse='+',sep='' ) , sep='~') ) } # APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 89 # define the formula used for covariance term formula2<-formula(paste('internal.data2$y3~',as.character(l3[2]),sep='')) # internal.data1$noncommon<- as.factor(c(1:n*0,1:n*0+1)) contrasts(internal.data1$noncommon)<-contr.treatment(2, base=1) internal.data1$indct1<-c(1:n*0+1,1:n*0 ) internal.data1$indct2<-c(1:n*0 ,1:n*0+1) # # if (!zeroL3){ data2<-data1[1:n,] names(data2)<-names(data1) } ##### # # add the common terms # if ( !is.null(l1l2) ) { formula1.terms<-paste( formula1.terms, as.character(l1l2[2]),sep='+') } # # add the special common terms (if any) # # # # in this section we identify non-common parameters # if a variable X is common in all formulas the we use term x*noncommon to include x+x:noncommon terms # otherwise use I(internal.data1$indct1*x) to add sepererate parameter on lambda1 # templ1<- labels(terms(l1)) # # run this only if there are terms in l1 formula if (length( templ1 )>0){ for ( k1 in 1:length( templ1 ) ){ if (!is.null(l1l2)){checkvar1<-sum(labels(terms(l1l2))==templ1[k1])==1} else{ checkvar1<-FALSE } checkvar2<-sum(labels(terms(l2))==templ1[k1] )==1 if (checkvar1&checkvar2) {formula1.terms<-paste(formula1.terms, paste('internal.data1$noncommon*',templ1[k1],sep=''), sep='+')} else{ formula1.terms<-paste(formula1.terms, paste('+I(internal.data1$indct1*',templ1[k1],sep=''), sep='') formula1.terms<-paste(formula1.terms, ')',sep='') } } } 90 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS # # if a variable X is not common st # otherwise use I(internal.data1$indct1*x) to add sepererate parameter on lambda1 # templ2<- labels(terms(l2)) # # run this only if there are terms in l1 formula if (length( templ2 )>0){ for ( k1 in 1:length( templ2 ) ){ if ( !is.null(l1l2) ) {checkvar1<-(sum(labels(terms(l1l2))==templ2[k1] )+ sum(labels(terms(l1))==templ2[k1] ))!=2 } else{ checkvar1<-TRUE } if ( checkvar1 ) { formula1.terms<-paste(formula1.terms, paste('+I(internal.data1$indct2*',templ2[k1],sep=''), sep='') formula1.terms<-paste(formula1.terms, ')',sep='') } } } # rm(templ1) rm(templ2) rm(Checkvar1) rm(Checkvar2) # # # # # # # This bit creates labels for special terms of type c(x1,x2) used in l1l2 # # formula1<-formula(paste('internal.data1$y1y2~',formula1.terms,sep='')) tmpform1<-as.character(formula1[3]) newformula<-formula1 while( regexpr('c\\(',tmpform1) != -1) { temppos1<-regexpr('c\\(',tmpform1)[1] tempfor <-substring( tmpform1, first = temppos1+ 2 ) temppos2<-regexpr('\\)', tempfor)[1] tempvar <-substring( tempfor , first = 1, last = temppos2-1 ) temppos3<-regexpr(',', tempvar)[1] tempname1<-substring(tempfor , first = 1, last = temppos3-1 ) tempname2<-substring(tempfor , first = temppos3+2, last=temppos2-1) tempname2<-sub( '\\)','', tempname2 ) APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 91 tempvar1<-data[, names(data)==tempname1] tempvar2<-data[, names(data)==tempname2] data1$newvar1<-c(tempvar1, tempvar2) # if( is.factor(tempvar1)& is.factor(tempvar2) ){ data1$newvar1<-as.factor(data1$newvar1) if (all(levels(tempvar1)==levels(tempvar2))){ attributes(data1$newvar1)<-attributes(tempvar1)} } tempvar<-sub( ',','..', tempvar ) names(data1)[names(data1)=='newvar1']<-tempvar newformula<-sub( 'c\\(','', tmpform1 ) newformula<-sub( '\\)','', newformula ) newformula<-sub( ',','..', newformula ) tmpform1<-newformula formula1<-formula(paste('internal.data1$y1y2~',newformula,sep='')) } ##### rm(temppos1) rm(temppos2) rm(temppos3) rm(tmpform1) rm(tempfor) rm(tempvar) rm(tempvar1) rm(tempvar2) rm(tempname1) rm(tempname2) # ------ # initial values for parameters prob<-0.20 s<-rep(0,n) vi<-1:n*0 v1<-1-c(vi,vi) like<-1:n*0 zero<- ( x==0 )|( y==0 ) if (zeroL3) { lambda3<-rep(0,n) } else { lambda3<-rep( max(0.1, cov(x,y,use='complete.obs')), n) } # # # # # Initial values for lambda internal.data1$v1<-1-c(vi,vi); lambda<-glm( formula1, family=poisson, data=data1, weights=internal.data1$v1, maxit=100)$fitted 92 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS # lambda1<-lambda[1:n] lambda2<-lambda[(n+1):(2*n)] # difllike<-100.0 loglike0<-1000.0 i<-0 ii<-0 if (zeroL3) { # # fit double poisson diagonal inflated model loglike<-rep(0, maxit) lambda3<-1:n*0 while ( (difllike>pres) && (i <= maxit) ) { i<-i+1 ##### E step ###### for (j in 1:n) { if (zero[j]) { s[j]<-0; # calculation of log-likelihood if (x[j]==y[j]) { density.di<-di.f( 0.0, theta ) like[j]<-log( (1-prob)*exp(-lambda1[j]- lambda2[j])+prob*density.di ); vi[j]<-prob*density.di*exp(-like[j]) } else{ like[j]<-log(1-prob)+log(dpois(x[j],lambda1[j]))+ log(dpois(y[j],lambda2[j])); vi[j]<-0.0 ;} } else { if (x[j]==y[j]) { density.di<-di.f( x[j],theta ); like[j]<-log( (1-prob)*dpois( x[j],lambda1[j] )* dpois( y[j],lambda2[j] ) + prob*density.di ); vi[j] <- prob*density.di*exp( -like[j] ) } else { vi[j]<-0.0; like[j]<-log(1-prob)+log( dpois(x[j],lambda1[j])* dpois(y[j],lambda2[j]) )} } } #### end of E-step ######### x1<-x; x2<-y; loglike[i]<-sum( like ) ; difllike<-abs( (loglike0-loglike[i])/loglike0 ) APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 93 else if (distribution==2) { meandiag<-theta } else if (distribution==3) { meandiag<- (1-theta)/theta } # fittedval1[x==y]<-prob*meandiag + fittedval1[x==y] fittedval2[x==y]<-prob*meandiag + fittedval2[x==y] # # # saving output result<-list(coefficients=allparameters, fitted.values= data.frame(x=fittedval1,y=fittedval2), residuals= data.frame(x=x-fittedval1,y=y-fittedval2), beta1=betaparameters$beta1, beta2=betaparameters$beta2, beta3=beta3, p=prob, theta=theta, diagonal.distribution=dilabel, lambda1=m$fitted[1:n], lambda2=m$fitted[(n+1):(2*n)], lambda3=lambda3, loglikelihood=loglike[1:i], parameters=noparams, AIC=AICtotal, BIC=BICtotal, iterations=i , call=tempcall) # } # end of elseif # options(warn=0) class(result)<-c('lm.dibp','lm') # result # # } "newnamesbeta" <- function( bvec ) { # # ------------------------------------------------------------------------ # Karlis and Ntzoufras (2004) # EM algorithms for Bivariate Poisson Models # ------------------------------------------------------------------------ # Internal function for renaming parameters according to their interpretation # # (c) June 2004 I. Ntzoufras, Chios, Greece # (c) Revised on May 2005 by I. Ntzoufras, Athens, Greece # ------------------------------------------------------------------------ # # (l1): is parameter rererring to the log(lamdba1) # (l2): is parameter rererring to the log(lamdba2) # (l2-l1): this parameter is added to the current parameter of l1 # no label: the parameter is common for both log(lambda1) and log(lambda2) # ------------------------------------------------------------------------ 100 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS names(bvec)<-sub('\\)','',names(bvec)) #remove right parenthesis names(bvec)<-sub('\\(Intercept','(Intercept)',names(bvec)) # replace "(Intercept" with "(Intercept)" names(bvec)[pmatch('internal.data1$noncommon2', names(bvec))]<-'(l2-l1):(Intercept)' # replace 'internal.data1$noncommon2'with 'l2-l1'for intercept names(bvec)<-sub('internal.data1\\$noncommon2:','(l2-l1):',names(bvec)) # the same for the rest of parameters names(bvec)<-sub('internal.data1\\$noncommon0:','(l1):',names(bvec)) # replace 'internal.data1\\$noncommon0:'by '(l1)' names(bvec)<-sub('internal.data1\\$noncommon1:','(l2):',names(bvec)) # replace 'internal.data1\\$noncommon1:'by '(l2)' names(bvec)<-sub(':internal.data1\\$noncommon2','(l2-l1):',names(bvec)) # same as above with ":" in front of expressions names(bvec)<-sub(':internal.data1\\$noncommon0','(l1):',names(bvec)) names(bvec)<-sub(':internal.data1\\$noncommon1','(l2):',names(bvec)) names(bvec)<-sub('I\\(internal.data1\\$indct1 \\* ','(l1):',names(bvec)) # replace 'I(internal.data1$indct1 * 'with '(l1):' names(bvec)<-sub('I\\(internal.data1\\$indct2 \\* ','(l2):',names(bvec)) # replace 'I(internal.data1$indct2 * 'with '(l2):' names(bvec) } "pbivpois" <- function(x, y=NULL, lambda = c(1, 1, 1), log=FALSE) { # -------------------------------------------------------------------------- # Karlis and Ntzoufras (2003, 2004) # EM algorithms for Bivariate Poisson Models # -------------------------------------------------------------------------- # x : matrix or vector of length n # y : vector of length n. If x is matrix then it is not used # lambda : parameters of the bivariate poisson distribution # log : argument controlling the calculation of the log-probability # or the probability function. # -------------------------------------------------------------------------- # if ( is.matrix(x) ) { var1<-x[,1] var2<-x[,2] } else if (is.vector(x)&is.vector(y)){ if (length(x)==length(y)){ var1<-x APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 101 var2<-y } else{ stop('lengths of x and y are not equal') } } else{ stop('x is not a matrix or x and y are not vectors') } n <- length(var1) logbp<-vector(length=n) # for (k in 1:n){ x0<-var1[k] y0<-var2[k] xymin<-min( x0,y0 ) lambdaratio<-lambda[3]/(lambda[1]*lambda[2]) # i<-0:xymin sums<- -lgamma(var1[k]-i+1)-lgamma(i+1)-lgamma(var2[k]-i+1)+i* log(lambdaratio) maxsums <- max(sums) sums<- sums - maxsums logsummation<- log( sum(exp(sums)) ) + maxsums logbp[k]<- -sum(lambda) + var1[k] * log( lambda[1] ) + var2[k] * log( lambda[2] ) + logsummation } if (log) { result<- logbp } else { result<-exp(logbp) } result # end of function bivpois } "simple.bp" <- function(x, y, ini3=1.0, maxit=300, pres=1e-8) { # # -------------------------------------------------------------------------- # Karlis and Ntzoufras (2003, 2004) # (last revision 25/8/2005) # Athens University of Economics and Business # # EM algorithms for Bivariate Poisson Models # -------------------------------------------------------------------------- # # x : matrix or vector of length n # y : vector of length n. If x is matrix then it is not used # ini3 : initial value for lambda3 102 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS # maxit : maximum number of iterations # pres : precision of the relative likelihood difference after which EM stops # -------------------------------------------------------------------------- # Data length # # if ( is.matrix(x) ) { var1<-x[,1] var2<-x[,2] } else if (is.vector(x)&is.vector(y)){ if (length(x)==length(y)){ var1<-x var2<-y } else{ stop('lengths of x and y are not equal') } } else{ stop('x is not a matrix or x and y are not vectors') } # # # n<-length(var1) # # initial values s<-rep(0,n) like<-1:n*0 zero<- ( var1==0 )|( var2==0 ) # # # # Initial values for lambda lambda3<- ini3 lambda1<- max( 0.1, mean(var1)-lambda3 ) lambda2<- max( 0.1, mean(var2)-lambda3 ) # # difllike<-1000.0 loglike0<-1000.0 i<-0 loglike<-rep(0,maxit) while ( (difllike>pres) && (i <= maxit) ) { APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 103 i<-i+1 ##### E step ###### for (j in 1:n) { if (zero[j]) { s[j]<-0; like[j]<- log(dpois(var1[j], lambda1)) + log(dpois(var2[j], lambda2))-lambda3; } else { lbp1<- pbivpois( var1[j]-1, var2[j]-1, lambda=c(lambda1,lambda2,lambda3), log=TRUE ); lbp2<- pbivpois( var1[j] , var2[j] , lambda=c(lambda1,lambda2,lambda3) , log=TRUE ); s[j]<-exp( log(lambda3) + lbp1 - lbp2 ); like[j]<-lbp2; } } ##### end of E step ###### x1<-var1-s x2<-var2-s loglike[i]<-sum(like) difllike<-abs( (loglike0-loglike[i])/loglike0 ) loglike0<-loglike[i] # # ##### M step ###### # # fit model on lambda3 lambda1<-mean(x1) lambda2<-mean(x2) lambda3<-mean(s) ##### end of M step ###### printpars<-c(i,lambda1, lambda2, lambda3, loglike[i] ) names(printpars)<-c('Iter.','lambda1','lambda2','lambda3','loglike') print( round(printpars ,3 ) ) cat( 'Relative Difference in Loglike:', difllike, '\n') } # # calculation of BIC and AIC of Bivariate Poisson model noparams<- 3 AIC<- -2*loglike[i] + noparams * 2 BIC<- -2*loglike[i] + noparams * log(2*n) # # # Calculation of BIC, AIC of Poisson saturated model x.mean<-var1 x.mean[var1==0]<-1e-12 104 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS y.mean<-var2 y.mean[var2==0]<-1e-12 AIC.sat <- sum(log(dpois(var1,x.mean)) + log( dpois(var2,y.mean))) BIC.sat <- -2 * AIC.sat + (2*n)* log(2*n) AIC.sat <- -2 * AIC.sat + (2*n)* 2 # # # Calculation of BIC, AIC of simple Poisson model x.mean<-mean(var1) y.mean<-mean(var2) AIC.pois <- sum(log( dpois(var1,x.mean)) + log(dpois(var2,y.mean))) BIC.pois <- -2 * AIC.pois + 2* log(2*n) AIC.pois <- -2 * AIC.pois + 2* 2 AICtotal<-c(AIC.sat, AIC.pois, AIC) BICtotal<-c(BIC.sat, BIC.pois, BIC ) names(AICtotal)<- c( 'Saturated','DblPois','BivPois') names(BICtotal)<- c( 'Saturated','DblPois','BivPois') # # Calculation of fitted values result<-list(lambda=c(lambda1, lambda2, lambda3),loglikelihood=loglike[1:i], parameters=noparams, AIC=AICtotal, BIC=BICtotal ,iterations=i ) # result # # } "splitbeta" <- function( bvec ){ # # ------------------------------------------------------------------------ # Karlis and Ntzoufras (2004) # EM algorithms for Bivariate Poisson Models # ------------------------------------------------------------------------ # Internal function for spliting beta parameters according to their interpretation # # (c) June 2004 I. Ntzoufras, Chios, Greece # ------------------------------------------------------------------------ # p3<-length(bvec) indx1<-grep( '\\(l1\\):', names(bvec) ) # identify parameters for lambda1 indx2<-grep( '\\(l2\\):', names(bvec) ) # identify parameters for lambda2 indx3<-grep( '\\(l2-l1\\):', names(bvec) ) # identify difference parameters for lambda2 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 105 # # create temporary labels to identify common parameters tempnames<-sub( '\\(l2-l1)\\:','k', names(bvec) ) tempnames<-sub( '\\(l2)\\:','k', tempnames ) tempnames<-sub( '\\(l1)\\:','k', tempnames ) indx4<-tempnames%in%names(bvec) # common parameters are identified as TRUE # beta1<-c(bvec[indx4],bvec[indx1]) beta2<-c(bvec[indx4],bvec[indx3],bvec[indx2]) indexbeta2<-c( rep(0,sum(indx4)), rep(1,length(indx3)), rep(2,length(indx2)) ) names(beta1)<-sub('\\(l1\\):','',names(beta1)) names(beta2)<-sub('\\(l2\\):','',names(beta2)) names(beta2)<-sub('\\(l2-l1\\):','',names(beta2)) beta1<-beta1[order(names(beta1))] indexbeta2<-indexbeta2[order(names(beta2))] beta2<-beta2[order(names(beta2))] ii<-1:length(beta2) ii<-ii[indexbeta2==0] for ( i in ii ) { # beta2[i]<-sum( beta2[ grep( names(beta2)[i], names(beta2) ) ] ) beta2[i]<-sum( beta2[ names(beta2)[i]==names(beta2) ] ) } beta2<-beta2[indexbeta2%in%c(0,2)] btemp<-list(beta1=beta1,beta2=beta2) btemp } # -------------------------------------------- # COMMANDS FOR FITTING CONSTANT MODEL FOR EX1 # SECTION 4.1.1 # -------------------------------------------- # # # library(bivpois) # load bivpois library # data(ex1.sim) # load data of example 1 # # ---------------------------------------------------------------------------- # # xtemp<-readline(prompt = "Press Enter to Continue") # # # Simple Bivariate Poisson Model # ex1.simple<-simple.bp( ex1.sim$x, ex1.sim$y ) # fit simple model of # section 4.1.1 # # names(ex1.simple) # monitor output variables 106 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS # ex1.simple$lambda # view lambda1 # ex1.simple$BIC # view BIC # ex1.simple # view all results of the model # xtemp<-readline(prompt = "Press Enter to Continue") # # # # plot of loglikelihood vs. iterations # win.graph() # plot( 1:ex1.simple$iterations, ex1.simple$loglikelihood, xlab='Iterations', # ylab='Log-likelihood', type='l') # # # ------------------------------------- # # COMMANDS FOR FITTING MODELS FOR EX1 # # SECTION 4.1.2 # # ------------------------------------- # ex1.m2<-lm.bp(x~1 , y~1 , data=ex1.sim, zeroL3=TRUE) # Model 2: DblPoisson(l1, l2) # ex1.m3<-lm.bp(x~1 , y~1 , data=ex1.sim) # Model 3: BivPoisson(l1, l2, l3) # ex1.m4<-lm.bp(x~. , y~. , data=ex1.sim, zeroL3=TRUE) # Model 4: DblPoisson (l1=Full, l2=Full) # ex1.m5<-lm.bp(x~. , y~. , data=ex1.sim) # Model 5: BivPoisson(l1=full, l2=full, l3=constant) # ex1.m6<-lm.bp(x~z1 , y~z1+z5 , l1l2=~z3, data=ex1.sim, zeroL3=TRUE) # Model 6: DblPois(l1,l2) # ex1.m7<-lm.bp(x~z1 , y~z1+z5 , l1l2=~z3, data=ex1.sim) # Model 7: BivPois(l1,l2,l3=constant) # ex1.m8<-lm.bp(x~. , y~. , l3=~., data=ex1.sim) # Model 8: BivPoisson(l1=full, l2=full, l3=full) # ex1.m9<-lm.bp(x~. , y~. , l3=~.-z5, data=ex1.sim) # Model 9: BivPoisson(l1=full, l2=full, l3=z1+z2+z3+z4) # ex1.m10<-lm.bp(x~z1 , y~z1+z5 , l1l2=~z3, l3=~., data=ex1.sim) # Model 10: BivPoisson(l1, l2, l3=full) # ex1.m11<-lm.bp(x~z1 , y~z1+z5 , l1l2=~z3, l3=~.-z5, data=ex1.sim) # Model 11: BivPoisson(l1, l2, l3= z1+z2+z3+z4) # ex1.m11$coef # monitor all beta parameters of model 11 # # # ex1.m11$beta1 # monitor all beta parameters of lambda1 of model 11 # ex1.m11$beta2 # monitor all beta parameters of lambda2 of model 11 # ex1.m11$beta3 # monitor all beta parameters of lambda3 of model 11 APÉNDICE B. APÉNDICE: PAQUETE BIVPOIS 107 Apéndice C Apéndice: Depuración de los datos Los datos que provienen de la base de datos de la DGT (referentes al número de fallecidos en vías urbanas e interurbanas por ciudad en cada año) fueron modificados con el mismo nombre de la provincia para poder hacer el estudio. La variable km_carretera ha sido creada a partir de los datos tomados del Ministerio de Fomento, además ha sido considerada constante a lo largo de los años para cada provincia. Las variables creadas a partir de los datos tomados del INE (precipitación, temperatura y horas de sol media) han sido consideradas como media de los datos tomados por cada mes en cada provincia. Además, se considera como constante la variación de los km de carretera en cada provincia por años. Por último, una vez añadidas las variables CCAA (comunidades autónomas) y poblacion (población de la provincia por año), se ha creado una tabla con todos los datos, pasando las variables a variables de tipo numérica, y se concluye generando la tabla (datos_f) con la que trabajamos que se encuentra en el archivo datos.RData. library(readxl) library(dplyr) library(lmtest) library(MASS) library(sandwich) library(AER) library(car) library(pscl) t2016<-read_xlsx("TG16.xlsx") t2016a<-t2016[-c(1,3),c(1,4,9)] t2016b<-t2016a[-c(1),] t2016b$year<-2016 t2016pv<-read_xlsx("PV16.xlsx") t2016PV<-t2016pv[-c(1,3),] t2016bb=cbind(t2016b,t2016PV[,c(9)]) t2016b=cbind(t2016bb,t2016PV[,c(10)]) t2016b=t2016b[-c(53),] #datost<-rbind(t2002b,t2016b) 109 if(t2005b$provincia[i]=="LAS PALMAS" ) t2005b$provincia[i]<-"Palmas (Las)" if(t2005b$provincia[i]=="PONTEVEDRA" ) t2005b$provincia[i]<-"Pontevedra" if(t2005b$provincia[i]=="SALAMANCA" ) t2005b$provincia[i]<-"Salamanca" if(t2005b$provincia[i]=="STA. C. TENERIFE" ) t2005b$provincia[i]<-"Santa Cruz de Tenerife" if(t2005b$provincia[i]=="CANTABRIA" ) t2005b$provincia[i]<-"Cantabria" if(t2005b$provincia[i]=="SEGOVIA" ) t2005b$provincia[i]<-"Segovia" if(t2005b$provincia[i]=="SEVILLA" ) t2005b$provincia[i]<-"Sevilla" if(t2005b$provincia[i]=="SORIA" ) t2005b$provincia[i]<-"Soria" if(t2005b$provincia[i]=="TARRAGONA" ) t2005b$provincia[i]<-"Tarragona" if(t2005b$provincia[i]=="TERUEL" ) t2005b$provincia[i]<-"Teruel" if(t2005b$provincia[i]=="TOLEDO" ) t2005b$provincia[i]<-"Toledo" if(t2005b$provincia[i]=="VALENCIA" ) t2005b$provincia[i]<-"Valencia/València" if(t2005b$provincia[i]=="VALLADOLID" ) t2005b$provincia[i]<-"Valladolid" if(t2005b$provincia[i]=="VIZCAYA" ) t2005b$provincia[i]<-"Bizkaia" if(t2005b$provincia[i]=="ZAMORA" ) t2005b$provincia[i]<-"Zamora" if(t2005b$provincia[i]=="ZARAGOZA" ) t2005b$provincia[i]<-"Zaragoza" if(t2005b$provincia[i]=="CEUTA" ) t2005b$provincia[i]<-"Ceuta" if(t2005b$provincia[i]=="MELILLA" ) t2005b$provincia[i]<-"Melilla" } t2006<-read_xls("TG06.xls") t2006a<-t2006[-c(1,3),c(1,10,15)] t2006b<-t2006a[-c(1),] t2006b$year<-2006 t2006pv<-read_xlsx("PV06.xlsx") t2006PV<-t2006pv[-c(1,3),] t2006bb=cbind(t2006b,t2006PV[,c(8)]) t2006b=cbind(t2006bb,t2006PV[,c(13)]) t2006b=t2006b[-c(53),] nombresantiguos<-unique(t2006b$provincia) colnames(t2006b)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "año","numero_vehiculos","poblacion") t2006b$numero_vehiculos[1]<-"181038" t2006b$poblacion[1]<-299921 for(i in (1:length(t2006b$provincia))) {cat(t2006b$provincia[i],"\n") if(t2006b$provincia[i]=="ÁLAVA") t2006b$provincia[i]<-"Araba/Álava" if(t2006b$provincia[i]=="ALBACETE") t2006b$provincia[i]<-"Albacete" if(t2006b$provincia[i]=="ALICANTE" ) t2006b$provincia[i]<-"Alicante/Alacant" if(t2006b$provincia[i]=="ALMERÍA" ) t2006b$provincia[i]<-"Almería" if(t2006b$provincia[i]=="ÁVILA" ) t2006b$provincia[i]<-"Ávila" if(t2006b$provincia[i]=="BADAJOZ" ) t2006b$provincia[i]<-"Badajoz" if(t2006b$provincia[i]=="ILLES BALEARS" ) t2006b$provincia[i]<-"Balears (Illes)" 116 APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS if(t2006b$provincia[i]=="BARCELONA" ) t2006b$provincia[i]<-"Barcelona" if(t2006b$provincia[i]=="BURGOS" ) t2006b$provincia[i]<-"Burgos" if(t2006b$provincia[i]=="CÁCERES" ) t2006b$provincia[i]<-"Cáceres" if(t2006b$provincia[i]=="CÁDIZ" ) t2006b$provincia[i]<-"Cádiz" if(t2006b$provincia[i]=="CASTELLÓN" ) t2006b$provincia[i]<-"Castellón/Castelló" if(t2006b$provincia[i]=="CIUDAD REAL" ) t2006b$provincia[i]<-"Ciudad Real" if(t2006b$provincia[i]=="CÓRDOBA" ) t2006b$provincia[i]<-"Córdoba" if(t2006b$provincia[i]=="A CORUÑA" ) t2006b$provincia[i]<-"Coruña (A)" if(t2006b$provincia[i]=="CUENCA" ) t2006b$provincia[i]<-"Cuenca" if(t2006b$provincia[i]=="GIRONA" ) t2006b$provincia[i]<-"Girona" if(t2006b$provincia[i]=="GRANADA" ) t2006b$provincia[i]<-"Granada" if(t2006b$provincia[i]=="GUADALAJARA" ) t2006b$provincia[i]<-"Guadalajara" if(t2006b$provincia[i]=="GUIPÚZCOA" ) t2006b$provincia[i]<-"Gipuzkoa" if(t2006b$provincia[i]=="HUELVA" ) t2006b$provincia[i]<-"Huelva" if(t2006b$provincia[i]=="HUESCA" ) t2006b$provincia[i]<-"Huesca" if(t2006b$provincia[i]=="JAÉN" ) t2006b$provincia[i]<-"Jaén" if(t2006b$provincia[i]=="LEÓN" ) t2006b$provincia[i]<-"León" if(t2006b$provincia[i]=="LLEIDA" ) t2006b$provincia[i]<-"Lleida" if(t2006b$provincia[i]=="LA RIOJA" ) t2006b$provincia[i]<-"Rioja (La)" if(t2006b$provincia[i]=="LUGO" ) t2006b$provincia[i]<-"Lugo" if(t2006b$provincia[i]=="MADRID" ) t2006b$provincia[i]<-"Madrid" if(t2006b$provincia[i]=="MÁLAGA" ) t2006b$provincia[i]<-"Málaga" if(t2006b$provincia[i]=="MURCIA" ) t2006b$provincia[i]<-"Murcia" if(t2006b$provincia[i]=="NAVARRA" ) t2006b$provincia[i]<-"Navarra" if(t2006b$provincia[i]=="OURENSE" ) t2006b$provincia[i]<-"Ourense" if(t2006b$provincia[i]=="ASTURIAS" ) t2006b$provincia[i]<-"Asturias" if(t2006b$provincia[i]=="PALENCIA" ) t2006b$provincia[i]<-"Palencia" if(t2006b$provincia[i]=="LAS PALMAS" ) t2006b$provincia[i]<-"Palmas (Las)" if(t2006b$provincia[i]=="PONTEVEDRA" ) t2006b$provincia[i]<-"Pontevedra" if(t2006b$provincia[i]=="SALAMANCA" ) t2006b$provincia[i]<-"Salamanca" if(t2006b$provincia[i]=="STA. C. TENERIFE" ) t2006b$provincia[i]<-"Santa Cruz de Tenerife" if(t2006b$provincia[i]=="CANTABRIA" ) t2006b$provincia[i]<-"Cantabria" if(t2006b$provincia[i]=="SEGOVIA" ) t2006b$provincia[i]<-"Segovia" if(t2006b$provincia[i]=="SEVILLA" ) t2006b$provincia[i]<-"Sevilla" if(t2006b$provincia[i]=="SORIA" ) t2006b$provincia[i]<-"Soria" if(t2006b$provincia[i]=="TARRAGONA" ) t2006b$provincia[i]<-"Tarragona" if(t2006b$provincia[i]=="TERUEL" ) t2006b$provincia[i]<-"Teruel" if(t2006b$provincia[i]=="TOLEDO" ) t2006b$provincia[i]<-"Toledo" if(t2006b$provincia[i]=="VALENCIA" ) t2006b$provincia[i]<-"Valencia/València" if(t2006b$provincia[i]=="VALLADOLID" ) t2006b$provincia[i]<-"Valladolid" if(t2006b$provincia[i]=="VIZCAYA" ) t2006b$provincia[i]<-"Bizkaia" if(t2006b$provincia[i]=="ZAMORA" ) t2006b$provincia[i]<-"Zamora" if(t2006b$provincia[i]=="ZARAGOZA" ) t2006b$provincia[i]<-"Zaragoza" if(t2006b$provincia[i]=="CEUTA" ) t2006b$provincia[i]<-"Ceuta" if(t2006b$provincia[i]=="MELILLA" ) t2006b$provincia[i]<-"Melilla" } APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS 117 t2007<-read_xls("TG07.xls") t2007a<-t2007[-c(1,3),c(1,10,15)] t2007b<-t2007a[-c(1),] t2007b$year<-2007 t2007pv<-read_xls("PV07.xls") t2007PV<-t2007pv[-c(1,3),] t2007bb=cbind(t2007b,t2007PV[,c(8)]) t2007b=cbind(t2007bb,t2007PV[,c(13)]) t2007b=t2007b[-c(53),] nombresantiguos<-unique(t2007b$provincia) colnames(t2007b)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "año","numero_vehiculos","poblacion") if(t2007b$numero_vehiculos[1]=="TOTAL") t2007b$numero_vehiculos[1]<-"191501" t2007b$poblacion[1]<-301548 for(i in (1:length(t2007b$provincia))) {cat(t2007b$provincia[i],"\n") if(t2007b$provincia[i]=="ÁLAVA") t2007b$provincia[i]<-"Araba/Álava" if(t2007b$provincia[i]=="ALBACETE") t2007b$provincia[i]<-"Albacete" if(t2007b$provincia[i]=="ALICANTE" ) t2007b$provincia[i]<-"Alicante/Alacant" if(t2007b$provincia[i]=="ALMERÍA" ) t2007b$provincia[i]<-"Almería" if(t2007b$provincia[i]=="ÁVILA" ) t2007b$provincia[i]<-"Ávila" if(t2007b$provincia[i]=="BADAJOZ" ) t2007b$provincia[i]<-"Badajoz" if(t2007b$provincia[i]=="ILLES BALEARS" ) t2007b$provincia[i]<-"Balears (Illes)" if(t2007b$provincia[i]=="BARCELONA" ) t2007b$provincia[i]<-"Barcelona" if(t2007b$provincia[i]=="BURGOS" ) t2007b$provincia[i]<-"Burgos" if(t2007b$provincia[i]=="CÁCERES" ) t2007b$provincia[i]<-"Cáceres" if(t2007b$provincia[i]=="CÁDIZ" ) t2007b$provincia[i]<-"Cádiz" if(t2007b$provincia[i]=="CASTELLÓN" ) t2007b$provincia[i]<-"Castellón/Castelló" if(t2007b$provincia[i]=="CIUDAD REAL" ) t2007b$provincia[i]<-"Ciudad Real" if(t2007b$provincia[i]=="CÓRDOBA" ) t2007b$provincia[i]<-"Córdoba" if(t2007b$provincia[i]=="A CORUÑA" ) t2007b$provincia[i]<-"Coruña (A)" if(t2007b$provincia[i]=="CUENCA" ) t2007b$provincia[i]<-"Cuenca" if(t2007b$provincia[i]=="GIRONA" ) t2007b$provincia[i]<-"Girona" if(t2007b$provincia[i]=="GRANADA" ) t2007b$provincia[i]<-"Granada" if(t2007b$provincia[i]=="GUADALAJARA" ) t2007b$provincia[i]<-"Guadalajara" if(t2007b$provincia[i]=="GUIPÚZCOA" ) t2007b$provincia[i]<-"Gipuzkoa" if(t2007b$provincia[i]=="HUELVA" ) t2007b$provincia[i]<-"Huelva" if(t2007b$provincia[i]=="HUESCA" ) t2007b$provincia[i]<-"Huesca" if(t2007b$provincia[i]=="JAÉN" ) t2007b$provincia[i]<-"Jaén" if(t2007b$provincia[i]=="LEÓN" ) t2007b$provincia[i]<-"León" if(t2007b$provincia[i]=="LLEIDA" ) t2007b$provincia[i]<-"Lleida" if(t2007b$provincia[i]=="LA RIOJA" ) t2007b$provincia[i]<-"Rioja (La)" if(t2007b$provincia[i]=="LUGO" ) t2007b$provincia[i]<-"Lugo" if(t2007b$provincia[i]=="MADRID" ) t2007b$provincia[i]<-"Madrid" 118 APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS if(t2007b$provincia[i]=="MÁLAGA" ) t2007b$provincia[i]<-"Málaga" if(t2007b$provincia[i]=="MURCIA" ) t2007b$provincia[i]<-"Murcia" if(t2007b$provincia[i]=="NAVARRA" ) t2007b$provincia[i]<-"Navarra" if(t2007b$provincia[i]=="OURENSE" ) t2007b$provincia[i]<-"Ourense" if(t2007b$provincia[i]=="ASTURIAS" ) t2007b$provincia[i]<-"Asturias" if(t2007b$provincia[i]=="PALENCIA" ) t2007b$provincia[i]<-"Palencia" if(t2007b$provincia[i]=="LAS PALMAS" ) t2007b$provincia[i]<-"Palmas (Las)" if(t2007b$provincia[i]=="PONTEVEDRA" ) t2007b$provincia[i]<-"Pontevedra" if(t2007b$provincia[i]=="SALAMANCA" ) t2007b$provincia[i]<-"Salamanca" if(t2007b$provincia[i]=="STA. C. TENERIFE" ) t2007b$provincia[i]<-"Santa Cruz de Tenerife" if(t2007b$provincia[i]=="CANTABRIA" ) t2007b$provincia[i]<-"Cantabria" if(t2007b$provincia[i]=="SEGOVIA" ) t2007b$provincia[i]<-"Segovia" if(t2007b$provincia[i]=="SEVILLA" ) t2007b$provincia[i]<-"Sevilla" if(t2007b$provincia[i]=="SORIA" ) t2007b$provincia[i]<-"Soria" if(t2007b$provincia[i]=="TARRAGONA" ) t2007b$provincia[i]<-"Tarragona" if(t2007b$provincia[i]=="TERUEL" ) t2007b$provincia[i]<-"Teruel" if(t2007b$provincia[i]=="TOLEDO" ) t2007b$provincia[i]<-"Toledo" if(t2007b$provincia[i]=="VALENCIA" ) t2007b$provincia[i]<-"Valencia/València" if(t2007b$provincia[i]=="VALLADOLID" ) t2007b$provincia[i]<-"Valladolid" if(t2007b$provincia[i]=="VIZCAYA" ) t2007b$provincia[i]<-"Bizkaia" if(t2007b$provincia[i]=="ZAMORA" ) t2007b$provincia[i]<-"Zamora" if(t2007b$provincia[i]=="ZARAGOZA" ) t2007b$provincia[i]<-"Zaragoza" if(t2007b$provincia[i]=="CEUTA" ) t2007b$provincia[i]<-"Ceuta" if(t2007b$provincia[i]=="MELILLA" ) t2007b$provincia[i]<-"Melilla" } t2008<-read_xls("TG08.xls") t2008a<-t2008[-c(1,3),c(1,10,15)] t2008b<-t2008a[-c(1),] t2008b$year<-2008 t2008pv<-read_xls("PV08.xls") t2008PV<-t2008pv[-c(1,3),] t2008bb=cbind(t2008b,t2008PV[-c(54,55),c(8)]) t2008b=cbind(t2008bb,t2008PV[-c(54,55),c(13)]) t2008b=t2008b[-c(53),] nombresantiguos<-unique(t2008b$provincia) colnames(t2008b)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "año","numero_vehiculos","poblacion") if(t2008b$numero_vehiculos[1]=="TOTAL") t2008b$numero_vehiculos[1]<-"196087" t2008b$poblacion[1]<-306527 for(i in (1:length(t2008b$provincia))) {cat(t2008b$provincia[i],"\n") if(t2008b$provincia[i]=="ÁLAVA") t2008b$provincia[i]<-"Araba/Álava" APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS 119 if(t2008b$provincia[i]=="ALBACETE") t2008b$provincia[i]<-"Albacete" if(t2008b$provincia[i]=="ALICANTE" ) t2008b$provincia[i]<-"Alicante/Alacant" if(t2008b$provincia[i]=="ALMERÍA" ) t2008b$provincia[i]<-"Almería" if(t2008b$provincia[i]=="ÁVILA" ) t2008b$provincia[i]<-"Ávila" if(t2008b$provincia[i]=="BADAJOZ" ) t2008b$provincia[i]<-"Badajoz" if(t2008b$provincia[i]=="ILLES BALEARS" ) t2008b$provincia[i]<-"Balears (Illes)" if(t2008b$provincia[i]=="BARCELONA" ) t2008b$provincia[i]<-"Barcelona" if(t2008b$provincia[i]=="BURGOS" ) t2008b$provincia[i]<-"Burgos" if(t2008b$provincia[i]=="CÁCERES" ) t2008b$provincia[i]<-"Cáceres" if(t2008b$provincia[i]=="CÁDIZ" ) t2008b$provincia[i]<-"Cádiz" if(t2008b$provincia[i]=="CASTELLÓN" ) t2008b$provincia[i]<-"Castellón/Castelló" if(t2008b$provincia[i]=="CIUDAD REAL" ) t2008b$provincia[i]<-"Ciudad Real" if(t2008b$provincia[i]=="CÓRDOBA" ) t2008b$provincia[i]<-"Córdoba" if(t2008b$provincia[i]=="A CORUÑA" ) t2008b$provincia[i]<-"Coruña (A)" if(t2008b$provincia[i]=="CUENCA" ) t2008b$provincia[i]<-"Cuenca" if(t2008b$provincia[i]=="GIRONA" ) t2008b$provincia[i]<-"Girona" if(t2008b$provincia[i]=="GRANADA" ) t2008b$provincia[i]<-"Granada" if(t2008b$provincia[i]=="GUADALAJARA" ) t2008b$provincia[i]<-"Guadalajara" if(t2008b$provincia[i]=="GUIPÚZCOA" ) t2008b$provincia[i]<-"Gipuzkoa" if(t2008b$provincia[i]=="HUELVA" ) t2008b$provincia[i]<-"Huelva" if(t2008b$provincia[i]=="HUESCA" ) t2008b$provincia[i]<-"Huesca" if(t2008b$provincia[i]=="JAÉN" ) t2008b$provincia[i]<-"Jaén" if(t2008b$provincia[i]=="LEÓN" ) t2008b$provincia[i]<-"León" if(t2008b$provincia[i]=="LLEIDA" ) t2008b$provincia[i]<-"Lleida" if(t2008b$provincia[i]=="LA RIOJA" ) t2008b$provincia[i]<-"Rioja (La)" if(t2008b$provincia[i]=="LUGO" ) t2008b$provincia[i]<-"Lugo" if(t2008b$provincia[i]=="MADRID" ) t2008b$provincia[i]<-"Madrid" if(t2008b$provincia[i]=="MÁLAGA" ) t2008b$provincia[i]<-"Málaga" if(t2008b$provincia[i]=="MURCIA" ) t2008b$provincia[i]<-"Murcia" if(t2008b$provincia[i]=="NAVARRA" ) t2008b$provincia[i]<-"Navarra" if(t2008b$provincia[i]=="OURENSE" ) t2008b$provincia[i]<-"Ourense" if(t2008b$provincia[i]=="ASTURIAS" ) t2008b$provincia[i]<-"Asturias" if(t2008b$provincia[i]=="PALENCIA" ) t2008b$provincia[i]<-"Palencia" if(t2008b$provincia[i]=="LAS PALMAS" ) t2008b$provincia[i]<-"Palmas (Las)" if(t2008b$provincia[i]=="PONTEVEDRA" ) t2008b$provincia[i]<-"Pontevedra" if(t2008b$provincia[i]=="SALAMANCA" ) t2008b$provincia[i]<-"Salamanca" if(t2008b$provincia[i]=="STA. C. TENERIFE" ) t2008b$provincia[i]<-"Santa Cruz de Tenerife" if(t2008b$provincia[i]=="CANTABRIA" ) t2008b$provincia[i]<-"Cantabria" if(t2008b$provincia[i]=="SEGOVIA" ) t2008b$provincia[i]<-"Segovia" if(t2008b$provincia[i]=="SEVILLA" ) t2008b$provincia[i]<-"Sevilla" if(t2008b$provincia[i]=="SORIA" ) t2008b$provincia[i]<-"Soria" if(t2008b$provincia[i]=="TARRAGONA" ) t2008b$provincia[i]<-"Tarragona" if(t2008b$provincia[i]=="TERUEL" ) t2008b$provincia[i]<-"Teruel" if(t2008b$provincia[i]=="TOLEDO" ) t2008b$provincia[i]<-"Toledo" if(t2008b$provincia[i]=="VALENCIA" ) t2008b$provincia[i]<-"Valencia/València" if(t2008b$provincia[i]=="VALLADOLID" ) t2008b$provincia[i]<-"Valladolid" if(t2008b$provincia[i]=="VIZCAYA" ) t2008b$provincia[i]<-"Bizkaia" 120 APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS if(t2008b$provincia[i]=="ZAMORA" ) t2008b$provincia[i]<-"Zamora" if(t2008b$provincia[i]=="ZARAGOZA" ) t2008b$provincia[i]<-"Zaragoza" if(t2008b$provincia[i]=="CEUTA" ) t2008b$provincia[i]<-"Ceuta" if(t2008b$provincia[i]=="MELILLA" ) t2008b$provincia[i]<-"Melilla" } t2009<-read_xls("TG09.xls") t2009a<-t2009[-c(1,3),c(1,10,15)] t2009b<-t2009a[-c(1),] t2009b$year<-2009 t2009pv<-read_xlsx("PV09.xlsx") t2009PV<-t2009pv[-c(1,3),] t2009bb=cbind(t2009b,t2009PV[-c(54,55),c(9)]) t2009b=cbind(t2009bb,t2009PV[-c(54,55),c(14)]) t2009b=t2009b[-c(53),] nombresantiguos<-unique(t2009b$provincia) colnames(t2009b)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "año","numero_vehiculos","poblacion") t2009b$numero_vehiculos[1]<-"199308" t2009b$poblacion[1]<-307656 for(i in (1:length(t2009b$provincia))) {cat(t2009b$provincia[i],"\n") if(t2009b$provincia[i]=="ÁLAVA") t2009b$provincia[i]<-"Araba/Álava" if(t2009b$provincia[i]=="ALBACETE") t2009b$provincia[i]<-"Albacete" if(t2009b$provincia[i]=="ALICANTE" ) t2009b$provincia[i]<-"Alicante/Alacant" if(t2009b$provincia[i]=="ALMERÍA" ) t2009b$provincia[i]<-"Almería" if(t2009b$provincia[i]=="ÁVILA" ) t2009b$provincia[i]<-"Ávila" if(t2009b$provincia[i]=="BADAJOZ" ) t2009b$provincia[i]<-"Badajoz" if(t2009b$provincia[i]=="ILLES BALEARS" ) t2009b$provincia[i]<-"Balears (Illes)" if(t2009b$provincia[i]=="BARCELONA" ) t2009b$provincia[i]<-"Barcelona" if(t2009b$provincia[i]=="BURGOS" ) t2009b$provincia[i]<-"Burgos" if(t2009b$provincia[i]=="CÁCERES" ) t2009b$provincia[i]<-"Cáceres" if(t2009b$provincia[i]=="CÁDIZ" ) t2009b$provincia[i]<-"Cádiz" if(t2009b$provincia[i]=="CASTELLÓN" ) t2009b$provincia[i]<-"Castellón/Castelló" if(t2009b$provincia[i]=="CIUDAD REAL" ) t2009b$provincia[i]<-"Ciudad Real" if(t2009b$provincia[i]=="CÓRDOBA" ) t2009b$provincia[i]<-"Córdoba" if(t2009b$provincia[i]=="A CORUÑA" ) t2009b$provincia[i]<-"Coruña (A)" if(t2009b$provincia[i]=="CUENCA" ) t2009b$provincia[i]<-"Cuenca" if(t2009b$provincia[i]=="GIRONA" ) t2009b$provincia[i]<-"Girona" if(t2009b$provincia[i]=="GRANADA" ) t2009b$provincia[i]<-"Granada" if(t2009b$provincia[i]=="GUADALAJARA" ) t2009b$provincia[i]<-"Guadalajara" if(t2009b$provincia[i]=="GUIPÚZCOA" ) t2009b$provincia[i]<-"Gipuzkoa" if(t2009b$provincia[i]=="HUELVA" ) t2009b$provincia[i]<-"Huelva" if(t2009b$provincia[i]=="HUESCA" ) t2009b$provincia[i]<-"Huesca" APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS 121 if(t2009b$provincia[i]=="JAÉN" ) t2009b$provincia[i]<-"Jaén" if(t2009b$provincia[i]=="LEÓN" ) t2009b$provincia[i]<-"León" if(t2009b$provincia[i]=="LLEIDA" ) t2009b$provincia[i]<-"Lleida" if(t2009b$provincia[i]=="LA RIOJA" ) t2009b$provincia[i]<-"Rioja (La)" if(t2009b$provincia[i]=="LUGO" ) t2009b$provincia[i]<-"Lugo" if(t2009b$provincia[i]=="MADRID" ) t2009b$provincia[i]<-"Madrid" if(t2009b$provincia[i]=="MÁLAGA" ) t2009b$provincia[i]<-"Málaga" if(t2009b$provincia[i]=="MURCIA" ) t2009b$provincia[i]<-"Murcia" if(t2009b$provincia[i]=="NAVARRA" ) t2009b$provincia[i]<-"Navarra" if(t2009b$provincia[i]=="OURENSE" ) t2009b$provincia[i]<-"Ourense" if(t2009b$provincia[i]=="ASTURIAS" ) t2009b$provincia[i]<-"Asturias" if(t2009b$provincia[i]=="PALENCIA" ) t2009b$provincia[i]<-"Palencia" if(t2009b$provincia[i]=="LAS PALMAS" ) t2009b$provincia[i]<-"Palmas (Las)" if(t2009b$provincia[i]=="PONTEVEDRA" ) t2009b$provincia[i]<-"Pontevedra" if(t2009b$provincia[i]=="SALAMANCA" ) t2009b$provincia[i]<-"Salamanca" if(t2009b$provincia[i]=="STA. C. TENERIFE" ) t2009b$provincia[i]<-"Santa Cruz de Tenerife" if(t2009b$provincia[i]=="CANTABRIA" ) t2009b$provincia[i]<-"Cantabria" if(t2009b$provincia[i]=="SEGOVIA" ) t2009b$provincia[i]<-"Segovia" if(t2009b$provincia[i]=="SEVILLA" ) t2009b$provincia[i]<-"Sevilla" if(t2009b$provincia[i]=="SORIA" ) t2009b$provincia[i]<-"Soria" if(t2009b$provincia[i]=="TARRAGONA" ) t2009b$provincia[i]<-"Tarragona" if(t2009b$provincia[i]=="TERUEL" ) t2009b$provincia[i]<-"Teruel" if(t2009b$provincia[i]=="TOLEDO" ) t2009b$provincia[i]<-"Toledo" if(t2009b$provincia[i]=="VALENCIA" ) t2009b$provincia[i]<-"Valencia/València" if(t2009b$provincia[i]=="VALLADOLID" ) t2009b$provincia[i]<-"Valladolid" if(t2009b$provincia[i]=="VIZCAYA" ) t2009b$provincia[i]<-"Bizkaia" if(t2009b$provincia[i]=="ZAMORA" ) t2009b$provincia[i]<-"Zamora" if(t2009b$provincia[i]=="ZARAGOZA" ) t2009b$provincia[i]<-"Zaragoza" if(t2009b$provincia[i]=="CEUTA" ) t2009b$provincia[i]<-"Ceuta" if(t2009b$provincia[i]=="MELILLA" ) t2009b$provincia[i]<-"Melilla" } t2010<-read_xls("TG10.xls") t2010a<-t2010[-c(1,3),c(1,10,15)] t2010b<-t2010a[-c(1),] t2010b$year<-2010 t2010pv<-read_xlsx("PV10.xlsx") t2010PV<-t2010pv[-c(1,3),] t2010bb=cbind(t2010b,t2010PV[-c(54,55),c(9)]) t2010b=cbind(t2010bb,t2010PV[-c(54,55),c(14)]) t2010b=t2010b[-c(53),] nombresantiguos<-unique(t2010b$provincia) colnames(t2010b)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "año","numero_vehiculos","poblacion") 122 APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS t2010b$numero_vehiculos[1]<-"202542" t2010b$poblacion[1]<-310562 for(i in (1:length(t2010b$provincia))) {cat(t2010b$provincia[i],"\n") if(t2010b$provincia[i]=="ARABA/ÁLAVA") t2010b$provincia[i]<-"Araba/Álava" if(t2010b$provincia[i]=="ALBACETE") t2010b$provincia[i]<-"Albacete" if(t2010b$provincia[i]=="ALICANTE" ) t2010b$provincia[i]<-"Alicante/Alacant" if(t2010b$provincia[i]=="ALMERÍA" ) t2010b$provincia[i]<-"Almería" if(t2010b$provincia[i]=="ÁVILA" ) t2010b$provincia[i]<-"Ávila" if(t2010b$provincia[i]=="BADAJOZ" ) t2010b$provincia[i]<-"Badajoz" if(t2010b$provincia[i]=="ILLES BALEARS" ) t2010b$provincia[i]<-"Balears (Illes)" if(t2010b$provincia[i]=="BARCELONA" ) t2010b$provincia[i]<-"Barcelona" if(t2010b$provincia[i]=="BURGOS" ) t2010b$provincia[i]<-"Burgos" if(t2010b$provincia[i]=="CÁCERES" ) t2010b$provincia[i]<-"Cáceres" if(t2010b$provincia[i]=="CÁDIZ" ) t2010b$provincia[i]<-"Cádiz" if(t2010b$provincia[i]=="CASTELLÓN" ) t2010b$provincia[i]<-"Castellón/Castelló" if(t2010b$provincia[i]=="CIUDAD REAL" ) t2010b$provincia[i]<-"Ciudad Real" if(t2010b$provincia[i]=="CÓRDOBA" ) t2010b$provincia[i]<-"Córdoba" if(t2010b$provincia[i]=="A CORUÑA" ) t2010b$provincia[i]<-"Coruña (A)" if(t2010b$provincia[i]=="CUENCA" ) t2010b$provincia[i]<-"Cuenca" if(t2010b$provincia[i]=="GIRONA" ) t2010b$provincia[i]<-"Girona" if(t2010b$provincia[i]=="GRANADA" ) t2010b$provincia[i]<-"Granada" if(t2010b$provincia[i]=="GUADALAJARA" ) t2010b$provincia[i]<-"Guadalajara" if(t2010b$provincia[i]=="GIPUZKOA" ) t2010b$provincia[i]<-"Gipuzkoa" if(t2010b$provincia[i]=="HUELVA" ) t2010b$provincia[i]<-"Huelva" if(t2010b$provincia[i]=="HUESCA" ) t2010b$provincia[i]<-"Huesca" if(t2010b$provincia[i]=="JAÉN" ) t2010b$provincia[i]<-"Jaén" if(t2010b$provincia[i]=="LEÓN" ) t2010b$provincia[i]<-"León" if(t2010b$provincia[i]=="LLEIDA" ) t2010b$provincia[i]<-"Lleida" if(t2010b$provincia[i]=="LA RIOJA" ) t2010b$provincia[i]<-"Rioja (La)" if(t2010b$provincia[i]=="LUGO" ) t2010b$provincia[i]<-"Lugo" if(t2010b$provincia[i]=="MADRID" ) t2010b$provincia[i]<-"Madrid" if(t2010b$provincia[i]=="MÁLAGA" ) t2010b$provincia[i]<-"Málaga" if(t2010b$provincia[i]=="MURCIA" ) t2010b$provincia[i]<-"Murcia" if(t2010b$provincia[i]=="NAVARRA" ) t2010b$provincia[i]<-"Navarra" if(t2010b$provincia[i]=="OURENSE" ) t2010b$provincia[i]<-"Ourense" if(t2010b$provincia[i]=="ASTURIAS" ) t2010b$provincia[i]<-"Asturias" if(t2010b$provincia[i]=="PALENCIA" ) t2010b$provincia[i]<-"Palencia" if(t2010b$provincia[i]=="LAS PALMAS" ) t2010b$provincia[i]<-"Palmas (Las)" if(t2010b$provincia[i]=="PONTEVEDRA" ) t2010b$provincia[i]<-"Pontevedra" if(t2010b$provincia[i]=="SALAMANCA" ) t2010b$provincia[i]<-"Salamanca" if(t2010b$provincia[i]=="STA. C. TENERIFE" ) t2010b$provincia[i]<-"Santa Cruz de Tenerife" if(t2010b$provincia[i]=="CANTABRIA" ) t2010b$provincia[i]<-"Cantabria" if(t2010b$provincia[i]=="SEGOVIA" ) t2010b$provincia[i]<-"Segovia" if(t2010b$provincia[i]=="SEVILLA" ) t2010b$provincia[i]<-"Sevilla" if(t2010b$provincia[i]=="SORIA" ) t2010b$provincia[i]<-"Soria" APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS 123 if(t2010b$provincia[i]=="TARRAGONA" ) t2010b$provincia[i]<-"Tarragona" if(t2010b$provincia[i]=="TERUEL" ) t2010b$provincia[i]<-"Teruel" if(t2010b$provincia[i]=="TOLEDO" ) t2010b$provincia[i]<-"Toledo" if(t2010b$provincia[i]=="VALENCIA" ) t2010b$provincia[i]<-"Valencia/València" if(t2010b$provincia[i]=="VALLADOLID" ) t2010b$provincia[i]<-"Valladolid" if(t2010b$provincia[i]=="BIZKAIA" ) t2010b$provincia[i]<-"Bizkaia" if(t2010b$provincia[i]=="ZAMORA" ) t2010b$provincia[i]<-"Zamora" if(t2010b$provincia[i]=="ZARAGOZA" ) t2010b$provincia[i]<-"Zaragoza" if(t2010b$provincia[i]=="CEUTA" ) t2010b$provincia[i]<-"Ceuta" if(t2010b$provincia[i]=="MELILLA" ) t2010b$provincia[i]<-"Melilla" } t2011<-read_xls("TG11.xls") t2011a<-t2011[-c(1,3),c(1,10,15)] t2011b<-t2011a[-c(1),] t2011b$year<-2011 t2011pv<-read_xlsx("PV11.xlsx") t2011PV<-t2011pv[-c(1,3),] t2011bb=cbind(t2011b,t2011PV[-c(54,55),c(9)]) t2011b=cbind(t2011bb,t2011PV[-c(54,55),c(14)]) t2011b=t2011b[-c(53),] nombresantiguos<-unique(t2011b$provincia) colnames(t2011b)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "año","numero_vehiculos","poblacion") t2011b$numero_vehiculos[1]<-"203328" t2011b$poblacion[1]<-312763 for(i in (1:length(t2011b$provincia))) {cat(t2011b$provincia[i],"\n") if(t2011b$provincia[i]=="ARABA/ÁLAVA") t2011b$provincia[i]<-"Araba/Álava" if(t2011b$provincia[i]=="ALBACETE") t2011b$provincia[i]<-"Albacete" if(t2011b$provincia[i]=="ALICANTE" ) t2011b$provincia[i]<-"Alicante/Alacant" if(t2011b$provincia[i]=="ALMERÍA" ) t2011b$provincia[i]<-"Almería" if(t2011b$provincia[i]=="ÁVILA" ) t2011b$provincia[i]<-"Ávila" if(t2011b$provincia[i]=="BADAJOZ" ) t2011b$provincia[i]<-"Badajoz" if(t2011b$provincia[i]=="ILLES BALEARS" ) t2011b$provincia[i]<-"Balears (Illes)" if(t2011b$provincia[i]=="BARCELONA" ) t2011b$provincia[i]<-"Barcelona" if(t2011b$provincia[i]=="BURGOS" ) t2011b$provincia[i]<-"Burgos" if(t2011b$provincia[i]=="CÁCERES" ) t2011b$provincia[i]<-"Cáceres" if(t2011b$provincia[i]=="CÁDIZ" ) t2011b$provincia[i]<-"Cádiz" if(t2011b$provincia[i]=="CASTELLÓN" ) t2011b$provincia[i]<-"Castellón/Castelló" if(t2011b$provincia[i]=="CIUDAD REAL" ) t2011b$provincia[i]<-"Ciudad Real" if(t2011b$provincia[i]=="CÓRDOBA" ) t2011b$provincia[i]<-"Córdoba" if(t2011b$provincia[i]=="A CORUÑA" ) t2011b$provincia[i]<-"Coruña (A)" if(t2011b$provincia[i]=="CUENCA" ) t2011b$provincia[i]<-"Cuenca" 124 APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS if(t2011b$provincia[i]=="GIRONA" ) t2011b$provincia[i]<-"Girona" if(t2011b$provincia[i]=="GRANADA" ) t2011b$provincia[i]<-"Granada" if(t2011b$provincia[i]=="GUADALAJARA" ) t2011b$provincia[i]<-"Guadalajara" if(t2011b$provincia[i]=="GIPUZKOA" ) t2011b$provincia[i]<-"Gipuzkoa" if(t2011b$provincia[i]=="HUELVA" ) t2011b$provincia[i]<-"Huelva" if(t2011b$provincia[i]=="HUESCA" ) t2011b$provincia[i]<-"Huesca" if(t2011b$provincia[i]=="JAÉN" ) t2011b$provincia[i]<-"Jaén" if(t2011b$provincia[i]=="LEÓN" ) t2011b$provincia[i]<-"León" if(t2011b$provincia[i]=="LLEIDA" ) t2011b$provincia[i]<-"Lleida" if(t2011b$provincia[i]=="LA RIOJA" ) t2011b$provincia[i]<-"Rioja (La)" if(t2011b$provincia[i]=="LUGO" ) t2011b$provincia[i]<-"Lugo" if(t2011b$provincia[i]=="MADRID" ) t2011b$provincia[i]<-"Madrid" if(t2011b$provincia[i]=="MÁLAGA" ) t2011b$provincia[i]<-"Málaga" if(t2011b$provincia[i]=="MURCIA" ) t2011b$provincia[i]<-"Murcia" if(t2011b$provincia[i]=="NAVARRA" ) t2011b$provincia[i]<-"Navarra" if(t2011b$provincia[i]=="OURENSE" ) t2011b$provincia[i]<-"Ourense" if(t2011b$provincia[i]=="ASTURIAS" ) t2011b$provincia[i]<-"Asturias" if(t2011b$provincia[i]=="PALENCIA" ) t2011b$provincia[i]<-"Palencia" if(t2011b$provincia[i]=="LAS PALMAS" ) t2011b$provincia[i]<-"Palmas (Las)" if(t2011b$provincia[i]=="PONTEVEDRA" ) t2011b$provincia[i]<-"Pontevedra" if(t2011b$provincia[i]=="SALAMANCA" ) t2011b$provincia[i]<-"Salamanca" if(t2011b$provincia[i]=="STA. C. TENERIFE" ) t2011b$provincia[i]<-"Santa Cruz de Tenerife" if(t2011b$provincia[i]=="CANTABRIA" ) t2011b$provincia[i]<-"Cantabria" if(t2011b$provincia[i]=="SEGOVIA" ) t2011b$provincia[i]<-"Segovia" if(t2011b$provincia[i]=="SEVILLA" ) t2011b$provincia[i]<-"Sevilla" if(t2011b$provincia[i]=="SORIA" ) t2011b$provincia[i]<-"Soria" if(t2011b$provincia[i]=="TARRAGONA" ) t2011b$provincia[i]<-"Tarragona" if(t2011b$provincia[i]=="TERUEL" ) t2011b$provincia[i]<-"Teruel" if(t2011b$provincia[i]=="TOLEDO" ) t2011b$provincia[i]<-"Toledo" if(t2011b$provincia[i]=="VALENCIA" ) t2011b$provincia[i]<-"Valencia/València" if(t2011b$provincia[i]=="VALLADOLID" ) t2011b$provincia[i]<-"Valladolid" if(t2011b$provincia[i]=="BIZKAIA" ) t2011b$provincia[i]<-"Bizkaia" if(t2011b$provincia[i]=="ZAMORA" ) t2011b$provincia[i]<-"Zamora" if(t2011b$provincia[i]=="ZARAGOZA" ) t2011b$provincia[i]<-"Zaragoza" if(t2011b$provincia[i]=="CEUTA" ) t2011b$provincia[i]<-"Ceuta" if(t2011b$provincia[i]=="MELILLA" ) t2011b$provincia[i]<-"Melilla" } t2012<-read_xls("TG12.xls") t2012a<-t2012[-c(1,3),c(1,10,15)] t2012b<-t2012a[-c(1),] t2012b$year<-2012 t2012pv<-read_xlsx("PV12.xlsx") t2012PV<-t2012pv[-c(1,3),] t2012bb=cbind(t2012b,t2012PV[-c(54,55),c(9)]) t2012b=cbind(t2012bb,t2012PV[-c(54,55),c(14)]) APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS 125 if(t2015b$provincia[i]=="CEUTA" ) t2015b$provincia[i]<-"Ceuta" if(t2015b$provincia[i]=="MELILLA" ) t2015b$provincia[i]<-"Melilla" } datos_f<-rbind(t2002b,t2003b,t2004b,t2005b,t2006b,t2007b,t2008b,t2009b, t2010b,t2011b,t2012b,t2013b,t2014b,t2015b,t2016b) x<-read_xlsx("P_media.xlsx") x<-rbind(c("Araba/Álava",44.375),x) x<-x[,2] x<-rbind(x,x,x,x,x,x,x,x,x,x,x,x,x,x,x) datos_f=cbind(datos_f,x) y<-read_xlsx("T_media.xlsx") y<-rbind(c("Araba/Álava",11.808),y) y<-y[,2] y<-rbind(y,y,y,y,y,y,y,y,y,y,y,y,y,y,y) datos_f=cbind(datos_f,y) z<-read_xlsx("H_sol.xlsx") z<-rbind(c("Araba/Álava",169.133),z) z<-z[,2] z<-rbind(z,z,z,z,z,z,z,z,z,z,z,z,z,z,z) datos_f=cbind(datos_f,z) h<-read_xlsx("km_carretera.xlsx") h<-rbind(c("Araba/Álava",1463.33),h) h<-h[,2] h<-rbind(h,h,h,h,h,h,h,h,h,h,h,h,h,h,h) datos_f=cbind(datos_f,h) r<-read_xlsx("CCAA.xlsx") r<-rbind(c("Araba/Álava","PVasco"),r) r<-r[,2] r<-rbind(r,r,r,r,r,r,r,r,r,r,r,r,r,r,r) datos_f=cbind(datos_f,r) colnames(datos_f)<-c("provincia","fallecidos_interurbana","fallecidos_urbana", "year","numero_vehiculos","poblacion","precipitacion_media", "temperatura_media","horas_sol_media","km_carretera","CCAA") #Lectura de los datos str(datos_f) datos_f$fallecidos_interurbana<-as.numeric(datos_f$fallecidos_interurbana) datos_f$fallecidos_urbana<-as.numeric(datos_f$fallecidos_urbana) datos_f$numero_vehiculos<-as.numeric(datos_f$numero_vehiculos) 132 APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS datos_f$poblacion<-round(as.numeric(datos_f$poblacion),digits = 0) datos_f$precipitacion_media<-as.numeric(datos_f$precipitacion_media) datos_f$temperatura_media<-as.numeric(datos_f$temperatura_media) datos_f$horas_sol_media<-as.numeric(datos_f$horas_sol_media) datos_f$km_carretera<-as.numeric(datos_f$km_carretera) unique(datos_f$provincia)->provincias table(datos_f$provincia) # Tabla save(datos_f,file = "datos.RData") APÉNDICE C. APÉNDICE: DEPURACIÓN DE LOS DATOS 133 Bibliografía [1] Anscombe, F.J. 1948. The transformation of poisson, binomial and negative-binomial data. Biometrika. 35, 3/4 (1948), 246–254. [2] Atoche Calzada, P. 2017. Modelos de regresión con datos de conteo: Aplicación a competiciones deportivas. (2017). [3] Calcaterra, E.M. 2017. Una revisión de los modelos de conteo con excesos de ceros. Disponible en http://www.iesta.edu.uy/wp-content/uploads/2018/01/pasantia_ martinez_voucher.pdf. [4] Dobson, A.J. and Barnett, A. 2008. An introduction to generalized linear models. Chapman; Hall/CRC. [5] Hachuel, L.S., Boggio, G.S. and Harvey, G. 2010. Modelos alternativos para el análisis de datos de conteo con exceso de ceros. (2010). [6] Hilbe, J.M. 2011. Negative binomial regression. [7] Lawless, J.F. 1987. Negative binomial and mixed poisson regression. Canadian Journal of Statistics. 15, 3 (1987), 209–225. [8] López-González, E. and Ruiz-Soler, M. 2011. Análisis de datos con el modelo lineal generalizado. una aplicación con r. Revista española de pedagogía. 248 (2011), 59–80. [9] Luque-Calvo, P.L. 2017. Escribir un trabajo fin de estudios con r markdown. Disponible en http://destio.us.es/calvo. [10] Madsen, H. and Thyregod, P. 2010. Introduction to general and generalized linear models. CRC Press. [11] McCullagh, N., P. 1989. Generalized linear models. Disponible en http://www. utstat.toronto.edu/~brunner/oldclass/2201s11/readings/glmbook.pdf. [12] Pino-Mejías, J.L. 2017. Apuntes del máster data science & big data: Modelo de regresión de poisson y modelo de regresión con exceso de ceros. [13] Quinteiro, E.M.G. 2014. Aplicación de modelos de regresión de poisson bivariados a los resultados de los partidos de la liga española de fútbol. [14] R Core Team 2013. R: A language and environment for statistical computing. R Foundation for Statistical Computing. [15] William, J., Hilbe, J.M. and Hilbe 2012. Generalized linear models and extensions. [16] Zeileis, A., Kleiber, C. and Jackman, S. 2008. Regression models for count data in r. Journal of statistical software. 27, 8 (2008), 1–25. 135