Modelos de regresión de Poisson

Traballo Fin de Grao Modelos de regresión de Poisson María García García 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Modelos de regresión de Poisson María García García Julio 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Estadística e Investigación Operativa Título: Modelos de regresión de Poisson Breve descrición do contido Un modelo de regresión nos permite establecer la relación entre una variable respuesta y una o varias variables explicativas. Habitualmente se considera que la variable respuesta es una variable continua y se utiliza el método de mínimos cuadrados para obtener estimaciones del modelo de regresión. Sin embargo, en muchas situaciones prácticas nos encontramos que la variable respuesta es el resultado de un recuento con lo cual no podíamos utilizar los métodos clásicos. Surgen así los modelos de regresión de Poisson. El trabajo se organizará en las siguientes secciones: Presentación del modelo de regresión de Poisson. Cálculo de los estimadores de un modelo de Poisson. Inferencia sobre los parámetros de un modelo de Poisson. Diagnosis y validación del modelo de regresión de Poisson. Sobre-dispersión en el modelo de Poisson. Además, presentaremos diferentes modelos de regresión de Poisson aplicados a conjuntos de datos. Para ello utilizaremos el software estadístico libre (https://www.r-project.org/). iii iv Recomendacións Outras observacións Índice general Resumen v 1. Introducción 1 2. Modelo de Poisson 9 2.1. LavariabledePoisson.............................. 9 2.2. El modelo de Poisson simple . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3. El modelo de Poisson múltiple . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3. Inferencia sobre los parámetros 27 3.1. Intervalos de conanza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2. Contrastes de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4. Diagnosis y validación del modelo 39 4.1. Validación de un modelo de Poisson . . . . . . . . . . . . . . . . . . . . . . 40 v vi ÍNDICE GENERAL 4.2. Bondad del modelo ajustado . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5. Sobre-dispersión 47 5.1. Métodos para determinar la sobre-dispersión . . . . . . . . . . . . . . . . . . 47 5.1.1. Método basado en la deviance ...................... 48 5.1.2. Método basado en los residuos de Pearson . . . . . . . . . . . . . . . 48 5.2. Corrección de la sobre-dispersión . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2.1. Modelo de regresión Binomial Negativo . . . . . . . . . . . . . . . . 49 5.2.2. Modelo de regresión Quasi-Poisson . . . . . . . . . . . . . . . . . . . 56 Anexo A: Código 59 Bibliografía 67 Resumo Os modelos de regresión de Poisson permiten representar a dependencia dunha variable resposta resultado dun reconto, con respecto a unha ou varias variables explicativas, aproximando a variable resposta discreta a partir das variables explicativas cun certo erro. O obxectivo deste traballo é estudar estes modelos en profundidade: realizar a estimación dos parámetros mediante o método de máxima verosimilitude e levar a cabo a Inferencia sobre eles expoñendo diversas metodoloxías. Co propósito de ilustrar todos os conceptos desenvolvidos empregaremos unha aplicación a datos reais. Para os datos analizados, deberemos comprobar que as hipótesis básicas do modelo se verican, pois senón as conclusións extraídas poderían non ser certas; e en caso de que non se cumpran, estudaremos posibles melloras do modelo. Ademais, mediremos a bondade de axuste do modelo, é dicir, a discrepancia entre os valores observados e os valores esperados. Unha das hipóteses máis importante e restritiva destes modelos, é a igualdade entre a media e a varianza da variable resposta, por seguir esta unha distribución de Poisson. Porén, na práctica, podemos atopar diferentes situacións nas cales a varianza sexa maior que a media, este fenómeno é o que se coñece como sobre-dispersión. Consideraremos varios métodos que permiten a identicación de datos con sobre-dispersión e estudaremos dous procedementos para corrixir este problema. Palabras chave: Modelo de regresión de Poisson; método de máxima verosimilitude; sobre-dispersión; modelo de regresión Binomial Negativo; modelo de regresión QuasiPoisson. vii 4 INTRODUCCIÓN donde x=1 n n P i=1 xi es la media muestral de la variable explicativa. Y=1 n n P i=1 Yi se corresponde con la media muestral de la variable respuesta. SxY =1 n n P i=1 (xi−x)(Yi−Y) es la covarianza muestral. S2 x=1 n n P i=1 (xi−x)2 se trata de la varianza muestral de la variable explicativa. Observamos que la recta de regresión aproximada por el método de mínimos cuadrados pasa por el vector de medias (x, Y ) y tiene pendiente b β1=SxY S2 x . Por otro lado, la varianza del error se estima mediante la varianza de los residuos. Es decir, consideraremos el estimador: bσ2=1 n−2 n X i=1 bεi2=1 n−2 n X i=1 (Yi−b β0−b β1xi)2. Veamos ahora las propiedades de los estimadores que acabamos de calcular. Con respecto a b β1 , es un estimador insesgado, es decir, E(b β1) = β1 , y V ar(b β1) = σ2 nS2 x . Además, como b β1 es combinación lineal de las variables Y1, ..., Yn que son normales e independientes (por serlo ε1, ..., εn y estar trabajando bajo diseño jo), entonces b β1 también tiene distribución Normal, esto es: b β1∈Nβ1,σ2 nS2 x. Por otro lado, b β0 es también un estimador insesgado, E(b β0) = β0 , y su varianza es V ar(b β0) = σ21 n+x2 nS2 x . Como pasaba antes, al ser b β0 combinación lineal de las variables Y1, ..., Yn que tienen distribución normal y son independientes, entonces b β0 también es Normal, esto es: b β0∈Nβ0, σ21 n+x2 nS2 x. INTRODUCCIÓN 5 Y nalmente, como consecuencia del Teorema de Fisher 3 , bσ2 es un estimador insesgado de σ2 y tiene una distribución ji-cuadrado 4 (n−2)bσ2 σ2∈χ2 n−2. Este modelo de regresión lineal simple se puede generalizar para casos en los que haya varias variables explicativas X1, X2, ..., Xp−1 mediante un modelo de regresión lineal múltiple de la siguiente forma: Y=β0+β1X1+... +βp−1Xp−1+ε. donde Y es la variable respuesta; X1, ..., Xp−1 las variables explicativas; β0, β1, ..., βp−1 los coecientes; y ε el error del modelo, que debe vericar que E(ε|X1, ..., Xp−1) = 0 . Si consideramos una muestra {xi,1, xi,2, ..., xi,p−1, Yi} con i= 1, ...n se tiene que Yi=β0+β1xi,1+... +βp−1xi,p−1+εi. Normalmente se usa la notación vectorial, xi= (1, xi,1, ..., xi,p−1) denota el vector la asociado a los valores que toman las variables explicativas para el i -ésimo individuo y β= (β0, β1, ..., βp−1)0 es el vector de coecientes. 3 Este teorema es el siguiente: Teorema 1.3. Sean X1, ..., Xn∈N(µ, σ2) independientes. Entonces se verican: I. X∈Nµ, σ2 n . II. nS2 σ2=(n−1)S2 c σ2∈χ2 n−1 . III. X y S2 (o S2 c ) son independientes. donde X=1 n n P i=1 Xi , S2=1 n n P i=1 (Xi−X)2 y S2 c=1 n−1 n P i=1 (Xi−X)2 . 4 Esto signica que: Denición 1.4. Sean Z1, ..., Zm variables aleatorias Normales estándar independientes. Diremos que la variable aleatoria W=Z2 1+... +Z2 m∈χ2 m sigue una distribución ji-cuadrado con m grados de libertad. Sus principales características son las siguientes: Sus valores posibles son todos no negativos, esto es, W≥0 . E(W) = m . V ar(W) = 2m . 6 INTRODUCCIÓN Habitualmente expresaremos, el modelo de regresión múltiple en forma matricial. Es decir, podemos escribir el modelo como     Y1 . . . Yn    =    1x11 ··· x1,p−1 . . .. . ..... . . 1xn1··· xn,p−1    ·       β0 β1 . . . βp−1        +    ε1 . . . εn    , esto es: Y=Xβ +ε siendo Y el vector de respuestas, X se denomina matriz de diseño y es una matriz n×p donde cada la representa a uno de los n individuos y cada columna una de las p−1 variables, β el vector de los parámetros y ε el vector de los errores que cumple que ε∈Nn(0, σ2In) siendo In la matriz identidad de dimensión n . Del mismo modo que en el caso del modelo lineal simple, la estimación del vector de coecientes β se hace mediante el método de mínimos cuadrados, es decir, el estimador b β será el que cumpla n X i=1 (Yi−xib β)2= m´ın β n X i=1 (Yi−xiβ)2 donde xi es la i -ésima la de la matriz de diseño X . Esto expresado matricialmente es m´ın β(Y−Xβ)0(Y−Xβ) = m´ın βφ(β). Luego, derivando la función φ respecto a β , igualando a cero y despejando obtenemos el estimador b β= (X0X)−1X0Y . La estimación de la varianza del error σ2 vendría dada por bσ2=1 n−p n X i=1 bεi2=1 n−p n X i=1 (Yi−xib β)2. Veamos las propiedades de estos estimadores. Con respecto a b β , suponiendo que E(ε) = 0 , b β es un estimador insesgado E(b β) = E((X0X)−1X0Y) = (X0X)−1X0E(Y)=(X0X)−1X0Xβ =β. INTRODUCCIÓN 7 Calculamos la matriz de covarianzas del vector aleatorio b β aplicando la hipótesis de homocedasticidad Cov(b β, b β) = Cov((X0X)−1X0Y, (X0X)−1X0Y)=(X0X)−1X0Cov(Y, Y ) ((X0X)−1X0)0 = (X0X)−1X0σ2In(X0)0((X0X)−1)0= (X0X)−1X0σ2InX((X0X)0)−1 = (X0X)−1X0σ2InX(X0X)−1=σ2(X0X)−1. El estimador b β tiene una distribución normal con el vector de medias y la matriz de covarianzas que hemos calculado anteriormente, es decir, b β∈Npβ, σ2(X0X)−1. Mientras que, de nuevo como consecuencia del Teorema de Fisher, se tiene que (n−p)bσ2 σ2∈χ2 n−p. Como hemos visto, usualmente se considera como variable respuesta una variable continua y se emplea el método de mínimos cuadrados para obtener las estimaciones del modelo. Pero, a menudo la variable respuesta es una variable discreta por lo que no podemos utilizar los métodos clásicos, y es así como aparecen los modelos de regresión de Poisson . Cuando la variable respuesta es un recuento 5 ilimitado {0,1,2,3, ...} podemos usar un modelo de regresión para datos de conteo. Por ejemplo, si queremos estudiar la utilización de la sanidad por los adultos de entre 18-65 años, la variable respuesta sería el número de visitas al médico. Como podemos observar se trata de un recuento del número total de visitas al médico que un adulto de entre 18-65 años ha hecho durante un año (por ejemplo), y sólo son posibles números no negativos y discretos (incluimos el cero porque una de las posibilidades es que el número de visitas en un año fuese cero). En algunos casos el conteo es sucientemente grande y se puede usar un modelo de regresión lineal como los que hemos descrito anteriormente, pero, no es lo habitual. Por lo tanto, dedicaremos este trabajo al estudio de los modelos de regresión para datos de conteo de Poisson. 5 Se dene un recuento como: Denición 1.5. Una variable de conteo o recuento es el número de sucesos o eventos que ocurren en una misma unidad de observación en un intervalo de espacio o tiempo denido. Más detalles de esta denición pueden encontrarse en [7]. 8 INTRODUCCIÓN El modelo de regresión de Poisson es un caso particular de modelos más generales que son los modelos lineales generalizados o también conocidos como modelos GLM (de las siglas en inglés Generalized Lineal Models ). Capítulo 2 Modelo de Poisson A lo largo de este Capítulo vamos a presentar el modelo de regresión de Poisson. En primer lugar, vamos a recordar la denición y principales características de la distribución de Poisson. En segundo lugar, introduciremos el modelo de Poisson simple y veremos cómo estimar los coecientes asociados al mismo. Por último, presentaremos el modelo de Poisson múltiple y la estimación de sus parámetros. Además, para ilustrar los conceptos que vamos a ir desarrollando a lo largo del Capítulo utilizaremos una aplicación a datos reales que usaremos a lo largo de todo el trabajo. 2.1. La variable de Poisson Recordemos la denición de la distribución de Poisson. Si Y es una variable que sigue una distribución de Poisson de parámetro λ > 0 , luego, por denición: P(Y=y) = e−λλy y!, para todo y= 0,1,2, ... (2.1) siendo y el número de veces que ocurre un evento y λ un parámetro positivo que representa el número de veces que se espera que ocurra el evento en un período determinado. A la vista de la denición anterior, una variable de Poisson verica que su esperanza es E(Y) = λ y su varianza es V ar(Y) = λ . El hecho de que E(Y) = V ar(Y) jugará un papel 9 10 CAPÍTULO 2. MODELO DE POISSON λ = 1 y Masa de probabilidad P(Y=y) 0.0 0.1 0.2 0.3 0.4 λ = 4 y Masa de probabilidad P(Y=y) 0.0 0.1 0.2 0.3 0.4 λ = 20 y Masa de probabilidad P(Y=y) 0.0 0.1 0.2 0.3 0.4 λ = 40 y Masa de probabilidad P(Y=y) 0.0 0.1 0.2 0.3 0.4 Figura 2.1: Representación gráca de las funciones de masa de probabilidad de una distribución de Poisson para diversos valores del parámetro. importante a la hora de analizar un modelo de regresión de Poisson. En la Figura 2.1 podemos observar la forma que presenta la función de masa de probabilidad de esta distribución para distintos valores del parámetro λ . Fijémonos que a medida que los valores de λ aumentan, una variable aleatoria que sigue una distribución de Poisson se aproxima a una Normal. Esta distribución de Poisson es la probabilidad de que un determinado número de eventos ocurran durante un intervalo de tiempo especíco, suponiendo que la probabilidad de que suceda este evento en un intervalo de tiempo dado es proporcional a la longitud de dicho intervalo e independiente de que sucedan otros eventos. Por ejemplo, podemos usarla para modelar el número de llamadas telefónicas entrantes a un servicio técnico o el número de terremotos en un periodo de tiempo jo. Aunque en la práctica las hipótesis no siempre se cumplen, por ejemplo, la tasa de llamadas telefó- 2.2. EL MODELO DE POISSON SIMPLE 11 nicas entrantes varía dependiendo de la hora del día y el ritmo de los terremotos no es completamente independiente. Una propiedad importante de esta distribución es que la suma de variables aleatorias de Poisson es también una variable de Poisson. Además, supongamos que Yi∈Pois(λi) para i= 1,2, ... y son independientes, entonces, PiYi∈Pois(Piλi) . Esta distribución presenta muchas otras propiedades que no citaremos aquí, pero, que pueden encontrarse, por ejemplo, en [4]. 2.2. El modelo de Poisson simple Una vez denida la distribución de Poisson, puede resultar interesante estudiar la relación que existe entre una variable que sigue esta distribución y una o varias variables explicativas. Empezaremos por el caso más sencillo, el modelo de regresión de Poisson simple, en el cual consideraremos una única variable explicativa. Puede profundizarse más en este tipo de modelos en [16]. Sea Y∈Pois(λ) una variable respuesta resultado de un conteo que toma valores en el conjunto {0,1,2,3, ...} y X una variable explicativa, queremos ver qué relación hay entre ellas mediante un modelo de regresión. Sabemos que el concepto de regresión se formaliza como la función: λ(x) = E(Y|X=x) para todo x . Es decir, la función de regresión es la media condicionada de la variable respuesta en función de cada valor de la variable explicativa. Entonces, si Y∈Pois(λ) (toma valores en {0,1,2,3, ...} ), nunca toma valores negativos, y por tanto, no es coherente aplicar un modelo lineal directo. Pues si λ(x) se expresa como una función lineal, la recta puede cruzar el eje X y dar predicciones negativas, lo cual se contradice con el hecho de que Y es un recuento. Luego, necesitamos una función de enlace o función link previa a la aplicación de cualquier modelo lineal. La función de enlace es una función g de la esperanza de Y , E(Y) 12 CAPÍTULO 2. MODELO DE POISSON (en este caso, E(Y) = λ(x, β) ), que relaciona λ con el predictor lineal: g(λ(x, β)) = β0+β1x. Como función link parece razonable elegir el logaritmo porque la función de regresión está en el intervalo (0,+∞) , solo toma valores no negativos, entonces log(λ(x, β)) = β0+β1x, de modo que, aplicando exponenciales, la función de regresión del modelo queda expresada de la siguiente forma λ(x, β) = eβ0+β1x=eβ0eβ1x=eβ0eβ1x. (2.2) Por ello, a este modelo se le llama a menudo modelo log-lineal . Por otra parte, a la vista de la expresión (2.2), podemos interpretar los parámetros del modelo de regresión de Poisson de la siguiente manera: eβ0 : Valor inicial de la variable respuesta, es decir, el valor de la función de regresión cuando x= 0 . eβ1 : Tasa de incremento de la respuesta esperada al incrementar una unidad la variable explicativa, esto es, pasando de x a x+ 1 tenemos que λ(x+ 1, β) = eβ0eβ1x+1 =eβ0eβ1xeβ1=λ(x, β)eβ1. Para realizar la estimación de los parámetros, consideramos una muestra aleatoria simple de tamaño n (x1, Y1), ..., (xn, Yn) donde x1, ..., xn son las realizaciones de la variable explicativa e Y1, ..., Yn son las realizaciones de la variable respuesta y además Yi∈Pois(λ(xi, β)) . Entonces, en este caso: log(λ(xi, β)) = β0+β1xi y λ(xi, β) = eβ0+β1xi. Por lo tanto, para el valor de la muestra xi , la predicción del Yi sería b Yi=λ(xi,b β) = eb β0+b β1xi (que es un valor de la curva de regresión). Luego, los errores que estamos cometiendo en la predicción son bεi=Yi−b Yi=Yi−eb β0+b β1xi , denominados residuos de la regresión. 2.2. EL MODELO DE POISSON SIMPLE 13 Para estimar los coecientes de este modelo β0 y β1 vamos a usar el método de máxima verosimilitud . El estimador de máxima verosimilitud (EMV) es aquel valor o valores que maximizan la masa de probabilidad (o densidad) de la muestra, es decir, será el valor que maximiza la función de verosimilitud. En nuestro caso, teniendo en cuenta la función de masa de probabilidad dada en (2.1), la función de verosimilitud es: L(β0, β1) = n Y i=1 "e−λ(xi,β)λ(xi, β)Yi Yi!#= e − n P i=1 λ(xi,β)n Q i=1 λ(xi, β)Yi n Q i=1 Yi! . Ahora, como la función logaritmo es monótona creciente y podemos suponer que la función L(β0, β1) es positiva, podemos aplicar la función logaritmo a L(β0, β1) porque sus máximos van a coincidir. Por lo tanto, aplicamos la función logaritmo a la expresión anterior pues nos va a facilitar los cálculos de derivación, usamos las propiedades de los logaritmos y que log(λ(xi, β)) = β0+β1xi⇒elog(λ(xi,β)) =eβ0+β1xi⇒λ(xi, β) = eβ0+β1xi y podemos concluir que: l(β0, β1) = log(L(β0, β1)) = log       e − n P i=1 λ(xi,β)n Q i=1 λ(xi, β)Yi n Q i=1 Yi!       =log e − n P i=1 λ(xi,β)!+log n Y i=1 λ(xi, β)Yi!−log n Y i=1 Yi!! =− n X i=1 λ(xi, β)log(e) + n X i=1 log(λ(xi, β)Yi)− n X i=1 log(Yi!) =− n X i=1 λ(xi, β) + n X i=1 Yilog(λ(xi, β)) − n X i=1 log(Yi!) = n X i=1 [−λ(xi, β) + Yilog(λ(xi, β)) −log(Yi!)] = n X i=1 [−eβ0+β1xi+Yi(β0+β1xi)−log(Yi!)] = n X i=1 [Yi(β0+β1xi)−eβ0+β1xi−log(Yi!)]. 20 CAPÍTULO 2. MODELO DE POISSON es decir, tenemos que logλ =Xβ ⇒λ=eXβ. De manera análoga a lo que vimos en la Sección 2.2, vamos a usar el método de máxima verosimilitud para realizar la estimación de los coecientes de este modelo (el vector de parámetros). La función de verosimilitud es: L(β) = n Y i=1 "e−λ(xi,β)λ(xi, β)Yi Yi!#= e − n P i=1 λ(xi,β)n Q i=1 λ(xi, β)Yi n Q i=1 Yi! . Ahora, aplicando la función logaritmo a la expresión anterior y usando que log(λ(xi, β)) = xiβ⇒elog(λ(xi,β)) =exiβ⇒λ(xi, β) = exiβ tenemos: l(β) = log       e − n P i=1 λ(xi,β)n Q i=1 λ(xi, β)Yi n Q i=1 Yi!       =log e − n P i=1 λ(xi,β)!+log n Y i=1 λ(xi, β)Yi!−log n Y i=1 Yi!! =− n X i=1 λ(xi, β)log(e) + n X i=1 log(λ(xi, β)Yi)− n X i=1 log(Yi!) =− n X i=1 λ(xi, β) + n X i=1 Yilog(λ(xi, β)) − n X i=1 log(Yi!) = n X i=1 [−λ(xi, β) + Yilog(λ(xi, β)) −log(Yi!)] = n X i=1 [−exiβ+Yixiβ−log(Yi!)] = n X i=1 [Yixiβ−exiβ−log(Yi!)] (2.3) En consecuencia, la función l(β) es: 2.3. EL MODELO DE POISSON MÚLTIPLE 21 l(β) = n X i=1 [Yixiβ−exiβ−log(Yi!)] = n X i=1 [Yi(β0+xi,1β1+... +xi,p−1βp−1)−e(β0+xi,1β1+...+xi,p−1βp−1)−log(Yi!)]. Derivamos la expresión anterior con respecto a βm y teniendo en cuenta que xi= (1, xi,1, ..., xi,p−1) , luego xiβ=β0+xi,1β1+... +xi,p−1βp−1, y entonces: ∂ ∂βm l(β) = n X i=1 [Yixi,m −exiβxi,m] = n X i=1 [Yi−exiβ]xi,m = n X i=1 [Yi−λ(xi, β)]xi,m para todo m. Por consiguiente, deducimos también que ∂ ∂β l(β) = n X i=1 [Yi−λ(xi, β)]xi porque la coordenada m -ésima de este vector gradiente acabamos de ver que era ∂ ∂βm l(β) = n X i=1 [Yi−λ(xi, β)]xi,m. Obtenemos las ecuaciones de verosimilitud igualando a cero estas expresiones: n X i=1 [Yi−λ(xi, β)]xi= 0. Nota 2.2 . Sabemos que la coordenada i -ésima del producto de una matriz S∈ Mm×n por un vector u∈Rn es [Su]i= n P j=1 Si,juj y que Si,j = (Sj,i)0 . Entonces, utilizando la Nota 2.2 anterior, tenemos que: 22 CAPÍTULO 2. MODELO DE POISSON n X i=1 [Yi−λ(xi, β)]xi,m = 0 ⇒ n X i=1 Yixi,m = n X i=1 λ(xi, β)xi,m ⇒ n X i=1 (xm,i)0Yi= n X i=1 (xm,i)0λ(xi, β) ⇒[X0Y]m= [X0λ]m⇒X0Y=X0λ es decir que, escrito matricialmente, podemos expresar las ecuaciones de verosimilitud de la siguiente forma: X0Y=X0λ⇒X0Y−X0λ= 0 ⇒X0(Y−λ) = 0. Resolviendo estas ecuaciones tendríamos que obtener el estimador de máxima verosimilitud b β (EMV). Sin embargo, igual que ocurría en el caso del modelo de Poisson simple, no hay una fórmula explícita para el estimador de máxima verosimilitud b β de la regresión de Poisson, porque se obtiene un sistema de ecuaciones implícitas que no tiene solución explícita. Y entonces, debemos recurrir a métodos numéricos para encontrar una solución. Otra vez, un ejemplo de método numérico empleado es el método de Newton-Raphson, pero, en la práctica este proceso no se acostumbra hacer a mano, se suele usar un software estadístico como . Para poder aplicar el método de Newton-Raphson vamos a necesitar la forma de la matriz hessiana, (véase [14]): ∂2 ∂β ∂β0l(β) = ∂ ∂β0∂ ∂β l(β)=∂ ∂β0 n X i=1 [Yi−λ(xi, β)]xi!. Para hacer esto, vamos a hacerlo componente a componente como antes: ∂2 ∂βm∂βl l(β) = ∂ ∂βl∂ ∂βm l(β)=∂ ∂βl n X i=1 [Yi−λ(xi, β)]xi,m! =∂ ∂βl n X i=1 [Yi−exiβ]xi,m!=∂ ∂βl n X i=1 [Yi−eβ0+β1xi,1+...+βp−1xp−1]xi,m! = n X i=1 [−exiβxi,l xi,m] = − n X i=1 [λ(xi, β)xi,l xi,m] = − n X i=1 [λ(xi, β) (xl,i)0xi,m]. 2.3. EL MODELO DE POISSON MÚLTIPLE 23 Y ahora sí, teniendo en cuenta que xi es la la i -ésima de X y entonces se tiene que x0 i es la i -ésima columna de X0 : ∂2 ∂β ∂β0l(β) = − n X i=1 x0 ixiλ(xi, β). Y podemos expresar esta matriz como ∂2 ∂β ∂β0l(β) = − n X i=1 x0 ixiλ(xi, β) = −X0V X (2.4) donde la matriz V viene dada por V=    λ(x1, β)··· 0 . . ..... . . 0··· λ(xn, β)    . Es decir: −X0V X =−    . . .. . .. . . 1x1··· xp−1 . . .. . .. . .    ·    λ1··· 0 . . ..... . . 0··· λn    ·       ··· 1··· ··· x1··· . . . ··· xp−1···        . Ahora ya podemos aplicar el método de Newton-Raphson para varias variables (explicado exhaustivamente en [15]). Su expresión general es: βk+1 =βk− H (βk)−1∇(βk). Y luego, se concluye que la k -ésima iteración, siendo β0 un iterante inicial (no confundir con el intercepto), es: βk+1 =βk+ (X0VkX)−1X0(Y−λ). (2.5) Para k sucientemente grande b β≈βk+1 , y entonces βk+1 es un estimador razonable del vector de coecientes β . Vamos a ver un ejemplo de la estimación de los parámetros del modelo de Poisson múltiple utilizando la misma base de datos que hemos empleado para el modelo de Poisson simple. 24 CAPÍTULO 2. MODELO DE POISSON Ejemplo 2.3 ( Modelo de Poisson múltiple ) . Emplearemos otra vez como ejemplo, los datos recogidos en la base de datos esdcomp . En este caso, el objetivo es estudiar si el número de visitas de pacientes, la residencia médica, el género, los ingresos y las horas trabajadas afectan al número de quejas recibidas. Vamos a presentar un modelo de regresión que nos permita explicar el número de quejas recibidas por un determinado doctor/a en función de visits , residency , gender , revenue y hours . Es decir, tenemos 5 variables explicativas que son: visits , residency , gender , revenue y hours , mientras que la variable respuesta es complaints . Para tener una idea sobre la relación entre la variable respuesta y las variables explicativas podemos realizar los diagramas de dispersión que se muestran en la Figura 2.4 obtenidos con el siguiente comando de : > pairs(datos) visits 048 ●● ● ● ● ● ● ● ● ●●● ● ● ● ● ●● ● ● ● ● ● ●●● ● ● ●● ● ● ●● ● ● ● ●● ●● ●● ● ● ● ● ● ●● ● ● ● ● ●● ●●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● 1.0 1.4 1.8 ● ●● ● ●●● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ● ●● ● ●● ● ●● ●● ● ● ●●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● 1000 2000 3000 600 1200 ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ●● ●● ● 0 2 4 6 8 10 ● ● ● ● ● ● ● ●● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ● ●● ● ● ● ● ● ● complaints ● ● ● ● ●● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ●● ● ●● ●● ● ● ● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ●● ●●● ● ●● ● ● ●● ●●● ● ●● ●● ● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ● ●● ● ●● ●● ● ● ● ● ● residency ● ● ●●●● ●● ● ● ● ●● ● ● ●●● ●● ● ●● ● ●● ●●● ● ●●●●● ● ●●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● 1.0 1.4 1.8 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● 1.0 1.4 1.8 ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ●● ● ● ●● ● ● ● ●● ● ● ●● ● ● ● ● ●●● ● ●● ● ● ●●● ● ●● ● ● ●● ●● ● ● ● ●● ● ● ● ●●●● ● ● ●● ●● ● ● ● gender ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ●●●● ● ●● ● ● ●● ●● ● ● ● ● ● ● ●● ●●● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ●● ● ● ● ●● ● ● ●●●● ● ● ●● ● ● ● ● ●● ● ●●●● ● ● ●●● ● ●● ● ● revenue 220 260 300 340 ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● 600 1000 1400 1800 1000 2500 ● ● ● ● ● ● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ●● ●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ●●● ● ● ● ● ●● ● ● ● ● ● ●●● ● ● ●● ● ● ●● ● ● ● ●● ●●●● ● 1.0 1.4 1.8 ● ● ● ● ●● ● ● ● ● ●● ● ●● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ● ●● ● ●● ● ●● ●● ● ● ●●●● ● ● ● ● 220 280 340 ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ●● ● ●● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● hours Figura 2.4: Diagramas de dispersión. Ofrece una matriz de diagramas de dispersión simples de todos los pares y todas las 2.3. EL MODELO DE POISSON MÚLTIPLE 25 ordenaciones posibles que se pueden formar con las variables, cada variable frente a las otras. Se puede observar que el diagrama de dispersión de la segunda la y quinta columna se corresponde con el de la Figura 2.2 (quejas frente a ingresos). Mostramos a continuación el ajuste del modelo de Poisson múltiple obtenido gracias al software estadístico : > mod_poisson_multiple = glm(complaints~visits+residency+gender+revenue+ hours, family = poisson(link = log), data = datos) > > mod_poisson_multiple$coefficients (Intercept) visits residencyY genderM revenue -0.0803447860 0.0009499373 -0.2319740072 0.1122391151 -0.0033827203 hours -0.0001569430 Podemos ver que la estimación del intercepto es −0.0803 . Mientras que, el resto de coecientes estimados son: b β1 = 0.0009 b β2 = −0.2320 b β3 = 0.1122 b β4 = −0.0034 b β5 = −0.0002 . A modo de ejemplo, sabemos que eb β4 es lo que disminuye la variable respuesta cuando la variable revenue se incrementa una unidad y las demás variables se mantienen constantes. Interpretaciones análogas podrían detallarse para las restantes variables explicativas. 26 CAPÍTULO 2. MODELO DE POISSON Capítulo 3 Inferencia sobre los parámetros asociados a un modelo de regresión de Poisson En este Capítulo, vamos a construir intervalos de conanza y efectuar contrastes de hipótesis para los parámetros del modelo de regresión de Poisson. En ambos casos, tendremos varias formas de realizar estas tareas de Inferencia sobre los parámetros. Y una vez más, vamos a emplear el mismo ejemplo del Capítulo 2 que nos permitirá ilustrar estos conceptos. 3.1. Intervalos de conanza En primer lugar vamos a recordar la denición de intervalo de conanza. Empleamos los intervalos de conanza para saber qué seguridad tenemos de que la estimación puntual obtenida se aproxime al verdadero valor del parámetro, nos van a permitir precisar la incertidumbre existente en la estimación. Denición 3.1. Un intervalo de conanza es un intervalo construido en base a la muestra, y por tanto aleatorio, que contiene al parámetro con una cierta probabilidad, denominada nivel de conanza (a menudo expresado en porcentaje). Diremos que (a, b) 27 28 CAPÍTULO 3. INFERENCIA SOBRE LOS PARÁMETROS es un intervalo de conanza para un parámetro γ con un nivel de conanza 1−α (donde α∈[0,1] ), si P(a<γ<b)≥1−α . Con el propósito de aplicar este tipo de técnicas de Inferencia, nos encontramos con que tenemos varias opciones para proceder. La primera de ellas es que, igual que hacíamos en regresión lineal, podemos recurrir a la distribución del estimador, pero en este caso, acudiremos a la distribución asintótica . Podemos aproximar la distribución del estimador de máxima verosimilitud b β dado en (2.5) mediante una distribución Normal usando el Teorema Central del Límite. Esto es: √n(b β−β)∼Np0, I1(β)−1(a) ⇒√n(b β−β)∼Np0, I1(b β)−1 (b) ⇒b β−β∼Np 0,I1(b β)−1 n! (c) ⇒b β−β∼Np0, In(b β)−1 (d) ⇒b β−β∼Np0,(X0b V X)−1 ⇒b β∼Npβ, (X0b V X)−1. donde En ( a ) hemos aproximado la matriz de información de Fisher , I1(β) , sustituyendo el verdadero valor del parámetro (que es desconocido) por su estimación b β . En ( b ) usamos que V ar b β−β=V ar 1 √n√n(b β−β) =1 √n2 V ar √n(b β−β) =1 nV ar √n(b β−β) =1 nI1(b β)−1=I1(b β)−1 n. En ( c ) hemos empleado que si I1(β) es la información que aporta una muestra de tamaño 1 tal que In(β) = n I1(β)⇒(In(β))−1= (n I1(β))−1⇒(In(β))−1=1 nI1(β)−1 ⇒(In(β))−1=I1(β)−1 n. 3.1. INTERVALOS DE CONFIANZA 29 En ( d ) utilizamos que la matriz de información de Fisher, denotada por In(β) , es una matriz denida como: In(β) = E−∂2 ∂β ∂β0logL(β) y como por (2.4) tenemos que ∂2 ∂β ∂β0logL(β) = ∂2 ∂β ∂β0l(β) = −X0V X y E(X0V X) = X0V X (dado que estamos trabajando bajo diseño jo), luego: In(β) = X0V X. Y por consiguiente, hemos llegado a que el estimador de los parámetros es asintóticamente insesgado y su matriz de covarianzas asintótica es (X0V X)−1 . Así pues, ya podemos construir los intervalos de conanza para cada parámetro βk basándonos en el pivote: b βk−βk q(X0b V X)−1 k,k ∼N(0,1). (3.1) De este modo, los intervalos de conanza asintóticos para un parámetro βk de nivel 1−α serían de la forma: b βk−zα/2q(X0b V X)−1 k,k ,b βk+zα/2q(X0b V X)−1 k,k  donde zα/2 representa el cuantil 1 de orden 1−α/2 de una distribución Normal estándar. Otra opción para construir intervalos de conanza para los parámetros asociados a un modelo de Poisson, es utilizar el perl de verosimilitud o prole likelihood , para investigar un poco más sobre este método se propone consultar [9]. El perl de verosimilitud de un parámetro βk se dene como el máximo de la función de verosimilitud con respecto a los demás parámetros, es decir: PL(βk) = m´ax β1,...,βk−1,βk+1,...,βpL(β) = m´ax β1,...,βk−1,βk+1,...,βpL(β1, ..., βk−1, βk, βk+1, ..., βp). 1 Entendemos por cuantil: Denición 3.2. El cuantil de orden 1−p (con 0< p < 1 ) de una distribución Z es el valor de la variable zp de modo que la probabilidad de coger un valor mayor que él es p . Es decir, P(Z > zp) = p (o equivalentemente, P(Z≤zp) = 1 −p ). Si se desea profundizar en esta denición se recomienda consultar [2] 36 CAPÍTULO 3. INFERENCIA SOBRE LOS PARÁMETROS > summary(mod_poisson_multiple) Call: glm(formula = complaints ~ visits + residency + gender + revenue + hours, family = poisson(link = log), data = datos) Deviance Residuals: Min 1Q Median 3Q Max -1.8989 -0.9193 -0.3835 0.4981 1.8221 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.0803448 1.1542122 -0.070 0.94450 visits 0.0009499 0.0003386 2.806 0.00502 ** residencyY -0.2319740 0.2029388 -1.143 0.25301 genderM 0.1122391 0.2235043 0.502 0.61554 revenue -0.0033827 0.0041553 -0.814 0.41560 hours -0.0001569 0.0006634 -0.237 0.81298 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 89.447 on 43 degrees of freedom Residual deviance: 49.995 on 38 degrees of freedom AIC: 184.77 Number of Fisher Scoring iterations: 5 Comentemos la salida de este comando. En la columna Estimate aparecen los coe- cientes estimados. En la segunda columna, Std. Error , tenemos los errores típicos estimados en base a la distribución asintótica, es decir, las desviaciones típicas estimadas de los coecientes estimados 3 . En la columna z value el cociente entre las estimaciones de los coecientes y los errores típicos, es decir, el estadístico de contraste para la hipótesis 3 Empleando Estimate y Std. Error se podría calcular el intervalo de conanza para los parámetros sin más que calcular los cuantiles correspondientes de la Normal estándar. 3.2. CONTRASTES DE HIPÓTESIS 37 nula de que el coeciente correspondiente sea igual a cero . Y en la última columna, los Pr(>|z|) son los niveles críticos para el contraste de que el coeciente vale cero. 4 Si lo hacemos con la primera opción que hemos señalado antes (es decir, utilizando la distribución asintótica), tenemos que observar estos valores Pr(>|z|) de los que acabamos de hablar. Notemos que todas las variables explicativas, excepto visits , son poco signicativas porque sus niveles críticos son superiores a los niveles de signicación usuales ( 10 % , 5 % , 1 % ). Y, por lo tanto, rechazamos la hipótesis nula, hay evidencias de que el coeciente asociado a la variable visits es distinto de cero (hay evidencias para rechazar la hipótesis nula). Así pues, si tenemos que descartar una variable, descartamos hours porque es la variable con mayor nivel crítico (el menos signicativo). La segunda forma de hacerlo (utilizando el perl de verosimilitud), que tenemos, es usando la función drop1 . Esta función plantea el contraste de hipótesis mediante la diferencia de las deviances que hemos comentado antes, estas deviances se corresponden al modelo más general y el modelo simplicado: > #--CONTRASTES DE HIPÓTESIS: > > #--Segunda opción: > drop1(mod_poisson_multiple, test = "Chi") Single term deletions Model: complaints ~ visits + residency + gender + revenue + hours Df Deviance AIC LRT Pr(>Chi) <none> 49.995 184.78 visits 1 57.568 190.35 7.5730 0.005925 ** residency 1 51.319 184.10 1.3237 0.249933 gender 1 50.251 183.03 0.2558 0.613035 revenue 1 50.665 183.44 0.6703 0.412964 hours 1 50.051 182.83 0.0559 0.813085 --- 4 En la última la de la salida del comando summary , podemos ver el número de iteraciones que se llevaron a cabo para resolver el sistema y obtener las estimaciones de los coecientes, que en este caso es igual a 5 (el método iterativo que aplica es un algoritmo llamado IRLS (iteratively reweighted least squares)). 38 CAPÍTULO 3. INFERENCIA SOBRE LOS PARÁMETROS Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 El modelo con todas las variables explicativas tiene una deviance de 49.995 . Si quitamos, por ejemplo, la variable hours , la deviance es 50.051 : una diferencia de 0.0559 (como vemos en la penúltima columna de la la correspondiente a esa variable). El estadístico de contraste χ2 = 0.0559 sigue (aproximadamente) una distribución ji-cuadrado con 6−5 = 1 grado de libertad, el cual da un p -valor de 0.8131 (esto puede ser comprobado con el comando pchisq de ). Por tanto, observemos que como decíamos antes, la variable hours es la candidata a ser descartada del modelo. Esto se debe a que, la menor diferencia de las deviances del modelo con todas las variables y el modelo sin una de ellas, se da con el modelo al que se le ha quitado la variable hours . Ahora, lo hacemos de la última forma que mencionamos, para realizar el contraste usamos la función anova y la opción test = Chi . Consideramos dos modelos, por ejemplo: el modelo de regresión de Poisson múltiple del Ejemplo 2.3 y ese mismo modelo sin la variable hours . > #--Tercera opción: > mod_poisson_multiple_sin_hours <- glm(complaints~visits+residency+gender+ revenue, family = poisson(link = log), data = datos) > anova(mod_poisson_multiple_sin_hours, mod_poisson_multiple, test="Chi") Analysis of Deviance Table Model 1: complaints ~ visits + residency + gender + revenue Model 2: complaints ~ visits + residency + gender + revenue + hours Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 39 50.051 2 38 49.995 1 0.055908 0.8131 La diferencia de las deviances es 0.0559 y sigue aproximadamente una distribución jicuadrado con 1 grado de libertad. El nivel crítico de 0.8131 no es inferior a los niveles de signicación habituales, y entonces no tenemos evidencias para rechazar la hipótesis nula, por tanto, no hay pruebas que demuestran que es mejor el modelo con todas las variables explicativas que el modelo sin la variable explicativa hours . Capítulo 4 Diagnosis y validación del modelo En los Capítulos 2 y 3 hemos presentado y analizado un modelo de regresión de Poisson suponiendo que las hipótesis del modelo son ciertas. En la práctica, debemos asegurarnos de que realmente estas hipótesis se verican para nuestros datos, puesto que en caso contrario, las conclusiones extraídas del modelo podrían no ser ciertas. Y en caso de que dichas hipótesis no se cumplan, estudiaremos posibles mejoras del modelo de regresión de Poisson. La validación y diagnosis de un modelo consiste precisamente en estudiar si las hipótesis básicas del modelo se verican en un conjunto de observaciones y en medir el ajuste del modelo de regresión de Poisson. A estas tareas dedicaremos este Capítulo 4. Recordemos las hipótesis que los modelos de regresión de Poisson deben cumplir para ajustarse a los datos adecuadamente, podemos encontrarlas en [6], y son las siguientes: Respuesta Poisson: La variable respuesta sigue una distribución de Poisson, es decir, es un recuento. Independencia: Las observaciones del error deben ser independientes entre sí. Media igual a varianza: La media de una variable aleatoria de Poisson debe ser igual a su varianza. Log-linealidad: El log(λ(x, β)) debe ser una función lineal de x . 39 40 CAPÍTULO 4. DIAGNOSIS Y VALIDACIÓN DEL MODELO 4.1. Validación de un modelo de Poisson Cuando seleccionamos un modelo para nuestros datos, en un principio, no podemos estar seguros de que ese modelo sea adecuado (es decir, si cumple las hipótesis). Para intentar solucionar este problema, lo que se hace es un proceso de validación del modelo de regresión de Poisson, y este procedimiento se lleva a cabo mediante un análisis de los residuos y sus representaciones grácas, que van a ser herramientas importantes en este proceso. Vamos a considerar residuos de distinta índole, cuyas deniciones son las siguientes: Residuos brutos. Diferencia entre el valor de la variable respuesta observado y la predicción para ese valor del modelo (es la distancia vertical entre la observación y la curva del ajuste): bεi=Yi−b Yi=Yi−λ(xi,b β) = Yi−exib β. Extienden la idea de residuos de la regresión clásica, pero, no nos van a resultar demasiado útiles en el caso del modelo de Poisson, porque no tienen por qué ser homocedásticos, ni presentar distribución simétrica entorno a cero. Residuos de Pearson. Estandarización de los residuos brutos, es decir, dividir los residuos brutos entre la raíz cuadrada de la varianza de Yi (tenemos en cuenta que, en la distribución de Poisson, la media coincide con la varianza): Yi−b Yi pV ar(Yi)=Yi−b Yi qb Yi . Este nombre es debido a que la suma de los cuadrados de estos residuos se corresponde con el estadístico de Pearson. Residuos de la deviance . Raíz cuadrada de los sumandos de la expresión de la deviance multiplicada por el signo del residuo bruto: signo(Yi−b Yi)s2·Yilog Yi b Yi−(Yi−b Yi). Es decir, la deviance es la suma de los cuadrados de los residuos de la deviance . Vamos a calcular estos tres tipos de residuos para el Ejemplo 2.1 , con este n, usaremos la función residuals y las opciones type = response , type = pearson y 4.1. VALIDACIÓN DE UN MODELO DE POISSON 41 type = deviance para calcular los residuos brutos, de Pearson y de la deviance, respectivamente. Necesitamos también la función predict y la opción type = response , que nos ofrecen las predicciones de la variable respuesta. Y así, obtenemos los diagramas de dispersión de estos tres tipos de residuos frente a sus correspondientes predicciones de la variable respuesta b Yi tal y como puede verse en la Figura 4.1. Además, vamos a añadir un gráco más a esta Figura 4.1, que va a ser el diagrama de dispersión de los residuos de la deviance frente a los predictores lineales xib β , para hallarlos, usamos otra vez la función predict , con la opción type = link . Hemos usado los predictores lineales, en lugar de las predicciones, para solucionar los problemas de escala en el eje horizontal. ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● 3.0 3.5 4.0 4.5 −2 0 2 4 6 Res. brutos vs. Predicciones Predicciones Res. brutos ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● 3.0 3.5 4.0 4.5 −2 −1 0 1 2 3 4 Res. Pearson vs. Predicciones Predicciones Res. Pearson ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● 3.0 3.5 4.0 4.5 −2 −1 0 1 2 3 Res. deviance vs. Predicciones Predicciones Res. deviance ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● 1.0 1.1 1.2 1.3 1.4 1.5 −2 −1 0 1 2 3 Res. deviance vs. Pred. lineales Predictores lineales Res. deviance Figura 4.1: Diagramas de dispersión de los residuos asociados al modelo de regresión de Poisson analizado en el Ejemplo 2.1 . Los residuos brutos son heterocedásticos, pues en el diagrama de dispersión de los valores 42 CAPÍTULO 4. DIAGNOSIS Y VALIDACIÓN DEL MODELO observados frente a los valores de la variable explicativa (Figura 2.3) ya veíamos que la desviación respecto a la curva del ajuste era mucho mayor en los valores grandes de la variable explicativa (las desviaciones son más acusadas en unas zonas que en otras, la varianza del error no es constante). Esto entra dentro de las hipótesis del modelo de regresión de Poisson, y es por eso que se consideran los otros dos tipos de residuos, que como podemos observar en estos grácos, son más homocedásticos. La homocedasticidad mejora al considerar estos residuos (jémonos en la escala del eje vertical), podríamos decir que los residuos son homocedásticos. En denitiva, como los residuos no presentan ningún patrón podemos decir que el modelo es adecuado para estos datos. Si por el contrario, hubiese patrones en estos grácos, esto es un indicador de que puede haber sobre-dispersión (concepto que veremos en el Capítulo 5). Y ahora, podemos obtener con la función plot del modelo de regresión de Poisson los grácos de diagnosis que nos da (véase la Figura 4.2): El primero de ellos es el mismo que el último de la Figura 4.1, un diagrama de dispersión de los residuos de la deviance frente a los predictores lineales. El situado en la primera la a la derecha, se corresponde con un gráco QQ de los residuos, que representa los cuantiles muestrales de los residuos de Pearson estandarizados frente a los cuantiles teóricos de una normal estándar, observemos que los puntos correspondientes a cada par cuantil-cuantil no caen sobre la diagonal de la gráca, por lo tanto no presentan normalidad. El de la segunda la a la izquierda, es un diagrama de dispersión de las raíces cuadradas de los residuos de Pearson estandarizados frente a los predictores lineales. Y en el último de los grácos, aparecen los conceptos de estadísticos de apalancamiento y de distancia de Cook. Esta última distancia se usa para saber si un dato es una observación inuyente , es decir, un punto que tiene impacto en las estimaciones del modelo. La distancia de Cook es una medida de cómo inuye la observación i -ésima sobre la estimación de β al ser retirada del conjunto de datos, una distancia de Cook grande signica que una observación tiene un peso grande en la estimación de β . Podemos calcularlas para nuestro ejemplo: > cooks.distance(mod_poisson_simple) 4.1. VALIDACIÓN DE UN MODELO DE POISSON 43 1.0 1.1 1.2 1.3 1.4 1.5 −2 −1 0 1 2 3 4 Predicted values Residuals ●● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● Residuals vs Fitted 5 944 ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ●● ●● ● ● ● −2 −1 0 1 2 −2 −1 0 1 2 3 4 Theoretical Quantiles Std. Pearson resid. Normal Q−Q 5 9 44 1.0 1.1 1.2 1.3 1.4 1.5 0.0 0.5 1.0 1.5 2.0 Predicted values Std. Pearson resid. ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● Scale−Location 5 944 0.00 0.05 0.10 0.15 −2 −1 0 1 2 3 4 Leverage Std. Pearson resid. ●● ● ● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ●● ● ● ● ● ● ●● ●● ● Cook's distance 0.5 0.5 1 Residuals vs Leverage 44 5 9 Figura 4.2: Grácos de validación y diagnosis del modelo desarrollado en el Ejemplo 2.1 . 123456 6.515855e-03 6.046649e-02 4.579574e-02 3.145628e-02 3.048989e-01 2.584071e-02 7 8 9 10 11 12 2.388998e-02 8.976202e-02 1.475720e-01 6.809779e-04 7.027896e-03 7.403240e-03 13 14 15 16 17 18 7.416322e-02 7.396199e-03 5.178853e-02 1.125854e-01 9.586121e-03 1.184714e-03 19 20 21 22 23 24 6.727843e-03 8.450078e-02 2.263467e-02 4.043870e-02 1.590412e-03 6.423751e-03 25 26 27 28 29 30 3.662856e-02 2.032518e-02 3.236417e-02 7.812881e-03 2.105833e-02 2.744964e-02 31 32 33 34 35 36 1.934804e-02 3.977915e-02 6.277584e-03 7.320781e-03 1.398771e-02 4.033456e-05 37 38 39 40 41 42 1.622868e-02 1.185866e-01 8.012412e-02 4.368192e-02 3.973153e-02 1.983878e-02 44 CAPÍTULO 4. DIAGNOSIS Y VALIDACIÓN DEL MODELO 43 44 2.661837e-02 4.050698e-01 Observamos que el cambio más grande ocurriría omitiendo la 44 -ésima observación. Y una observación se dice que es una observación atípica (outlier) si es numéricamente distante del resto de los datos. En la Figura 4.2, vemos que las observaciones 5,9 y 44 -ésimas tienen valores residuales grandes. Es importante recalcar que las observaciones atípicas no se deben sacar inmediatamente del modelo, antes se deben estudiar para ver si hay algo raro con ellas, y si es así, se sacan de la base y se ajusta nuevamente el modelo. 4.2. Bondad del modelo ajustado En el modelo de regresión de Poisson no se dispone de un coeciente de determinación R2 , a través del que podríamos medir el ajuste del modelo. Lo máximo que nos podemos aproximar a este coeciente es la variabilidad explicada (también llamada pseudo R2 ) que es la parte de la variabilidad que podemos explicar en base al modelo, que queda justicada por la inuencia de las variables explicativas y no al error del modelo. Esta, se calcula como sigue: Pseudo R2= 100 ×null deviance −deviance null deviance donde la null deviance es la deviance de un modelo que solo contenga el intercepto, y deviance es la deviance del modelo que estamos considerando. Si este valor es grande (próximo a 100 ), signica que el modelo se ajusta bien los datos y es muy útil para realizar predicciones, dado que la variable explicativa explica gran parte de la variabilidad de la variable respuesta y la variabilidad del error es baja. En el Ejemplo 2.1 , podemos observar en la salida del comando summary , los valores de la null deviance y la deviance, y así obtener la variabilidad explicada: Pseudo R2= 100 ×null deviance −deviance null deviance = 100 ×89.447 −86.929 89.447 = 2.82 % De aquí deducimos que la variable explicativa revenue explica el 2.82% de la variabilidad de las quejas recibidas. Entonces, como este valor es muy pequeño, el modelo no ajusta bien 4.2. BONDAD DEL MODELO AJUSTADO 45 los datos (las observaciones/datos están lejos de la curva de ajuste) y no es útil para realizar predicciones, pues la variable explicativa revenue explica poca parte de la variabilidad de la variable respuesta complaints . 52 CAPÍTULO 5. SOBRE-DISPERSIÓN Luego, necesitamos una función de enlace previa a cualquier modelo lineal. Como función link parece coherente escoger el logaritmo porque la función de regresión está en el intervalo (0,+∞) , solo toma valores no negativos, entonces log(µ(xi, β)) = xiβ⇒µ(xi, β) = exiβ siendo xi es la i -ésima la de la matriz X , β un vector columna de los coecientes del modelo y µ(xi, β) es una de las componentes de vector columna µ . Advirtamos que el parámetro de sobre-dispersión no depende de las variables explicativas (es constante). Vamos a usar el método de máxima verosimilitud para realizar la estimación de los coecientes de este modelo (el vector de parámetros β y θ ). La función de máxima verosimilitud es: L(β, θ) = n Y i=1 f(yi;µ(xi, β); θ) = n Y i=1 Γ(yi+θ) yi! Γ(θ) θθµ(xi, β)yi (θ+µ(xi, β))θ+yi = n Q i=1 Γ(yi+θ) n Q i=1 yi! n Q i=1 Γ(θ) θ n P i=1 θn Q i=1 µ(xi, β)yi n Q i=1 (θ+µ(xi, β))θ+yi . Aplicando la función logaritmo a la expresión anterior se tiene que: l(β, θ) = log(L(β, θ)) = log       n Q i=1 Γ(yi+θ) n Q i=1 yi! n Q i=1 Γ(θ) θ n P i=1 θn Q i=1 µ(xi, β))yi n Q i=1 (θ+µ(xi, β))θ+yi       =log n Y i=1 Γ(yi+θ)!−log n Y i=1 yi!!−log n Y i=1 Γ(θ)!+log θ n P i=1 θ! +log n Y i=1 µ(xi, β)yi!−log n Y i=1 (θ+µ(xi, β))θ+yi! = n X i=1 log (Γ(yi+θ)) − n X i=1 log (yi!) − n X i=1 log (Γ(θ)) + n X i=1 θ·log (θ) + n X i=1 yi·log (µ(xi, β)) − n X i=1 ((θ+yi)·log(θ+µ(xi, β))) = n X i=1 log(Γ(yi+θ)) −log(Γ(yi+ 1)) −log(Γ(θ)) + θ·log(θ) +yi·log(µ(xi, β)) −θ·log(θ+µ(xi, β)) −yi·log(θ+µ(xi, β)) 5.2. CORRECCIÓN DE LA SOBRE-DISPERSIÓN 53 = n X i=1 log(Γ(yi+θ)) −log(Γ(yi+ 1)) −log(Γ(θ)) +θ·[log(θ)−log(θ+µ(xi, β))] + yi·[log(µ(xi, β)) −log(θ+µ(xi, β))] = n X i=1 log(Γ(yi+θ)) −log(Γ(yi+ 1)) −log(Γ(θ)) + θ·log θ θ+µ(xi, β) +yi·log µ(xi, β) θ+µ(xi, β). Las estimaciones de los coecientes se obtendrían de una forma muy similar a la que ya vimos para el caso del modelo de regresión de Poisson. Ejemplo 5.2 ( Datos esdcomp analizados en Faraway (2016) ) . Para ilustrar el modelo de regresión Binomial Negativo, vamos a emplear el mismo ejemplo que hemos utilizado anteriormente. Aunque, esta vez, utilizaremos dos variables explicativas: revenue y hours . Vamos a presentar un modelo de regresión que nos permita explicar el número de quejas recibidas por un determinado doctor/a en función de sus ingresos y del número total de horas trabajadas. La variable respuesta es complaints y las variables explicativas son: revenue (medida en dólares por hora) y hours . Para poder aplicar este modelo en necesitamos la función glm.nb del paquete MASS (para más información ver [12]): > library(MASS) > BN_GLM <- glm.nb(complaints~hours+revenue, link = log, data = datos) > summary(BN_GLM) Call: glm.nb(formula = complaints ~ hours + revenue, data = datos, link = log, init.theta = 7.016018421) Deviance Residuals: Min 1Q Median 3Q Max -1.6337 -0.9254 -0.2507 0.7011 1.6847 Coefficients: 54 CAPÍTULO 5. SOBRE-DISPERSIÓN Estimate Std. Error z value Pr(>|z|) (Intercept) -2.1784963 1.0291337 -2.117 0.0343 * hours 0.0014919 0.0003743 3.986 6.73e-05 *** revenue 0.0044594 0.0030849 1.446 0.1483 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for Negative Binomial(7.016) family taken to be 1) Null deviance: 59.936 on 43 degrees of freedom Residual deviance: 40.453 on 41 degrees of freedom AIC: 187.32 Number of Fisher Scoring iterations: 1 Theta: 7.02 Std. Err.: 4.46 2 x log-likelihood: -179.315 Y como podemos observar, la salida del comando summary es similar a la del modelo de regresión de Poisson, la única diferencia es que además nos ofrece: la 2×verosimilitud , la estimación del parámetro de sobre-dispersión b θ= 7.02 y de su error típico (b θ) = 4.46 . Además, como alguno de los parámetros es no signicativo al nivel del 5 % , deberíamos volver a hacer una selección del modelo (de igual forma que se hacía para el modelo de Poisson). Así mismo, podemos obtener, en este caso también, los grácos de validación y diagnosis para el modelo con la función plot de que se muestra en la Figura 5.1. Adicionalmente, podemos realizar un contraste de la sobre-dispersión, es decir, un contraste entre un modelo de regresión de Poisson y un modelo de regresión Binomial Negativa. Para llevar a cabo este contraste, utilizamos la función lrtest del paquete lmtest (para 5.2. CORRECCIÓN DE LA SOBRE-DISPERSIÓN 55 0.0 0.5 1.0 1.5 −1 0 1 2 Predicted values Residuals ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● Residuals vs Fitted 5 15 9 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ● ● −2 −1 0 1 2 −1 0 1 2 Theoretical Quantiles Std. Pearson resid. Normal Q−Q 5 15 9 0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 Predicted values Std. Pearson resid. ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● Scale−Location 5 15 9 0.00 0.05 0.10 0.15 −1 0 1 2 Leverage Std. Pearson resid. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● Cook's distance 0.5 Residuals vs Leverage 5 44 25 Figura 5.1: Grácos de validación y diagnosis del modelo de regresión Binomial Negativo del Ejemplo 5.2 . más información ver[1]). > #--Contraste: > install.packages("lmtest") > library(lmtest) > mod_poisson_2 <- glm(complaints~hours+revenue, family=poisson(link=log), + data=datos) > > lrtest(mod_poisson_2,BN_GLM) Likelihood ratio test Model 1: complaints ~ hours + revenue Model 2: complaints ~ hours + revenue #Df LogLik Df Chisq Pr(>Chisq) 1 3 -91.932 2 4 -89.658 1 4.5484 0.03295 * 56 CAPÍTULO 5. SOBRE-DISPERSIÓN --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Chisq se corresponde con dos veces la diferencia entre las log-verosimilitudes de los dos modelos. Como ya hemos visto, el nivel crítico se obtiene de una distribución ji-cuadrado con Df grados de libertad siendo Df la diferencia entre número de parámetros de ambos modelos. En este caso, el nivel crítico es aproximadamente 0.033 . Fijado un nivel de signicación α= 5 % , entonces podemos armar la hipótesis alternativa porque tenemos pruebas para rechazar la hipótesis nula. Es decir, se rechaza el modelo de Poisson en favor del modelo Binomial Negativo (lo cual signica también que, el parámetro de sobre-dispersión no está próximo a 1 , este es el contraste del que hablábamos en la Subsección 5.1.1 ). 5.2.2. Modelo de regresión Quasi-Poisson Otra alternativa para corregir la sobre-dispersión (válida también para la infra-dispersión) es, usar un modelo de regresión Quasi-Poisson . En este caso, se intenta corregir la sobre-dispersión introduciendo un nuevo parámetro en el modelo de Poisson clásico. En esto se puede hacer mediante la opción family = quasipoisson de la función glm , y así se puede estimar el parámetro de sobre-dispersión. El modelo de regresión Quasi-Poisson es una generalización del modelo de regresión de Poisson (si el parámetro de dispersión es φ= 1 , tenemos el modelo de Poisson) y se usa para modelar una variable de conteo con sobre-dispersión. Es una ligera recticación del modelo de Poisson. Por ejemplo, en el modelo de regresión de Poisson, podemos aproximar el parámetro de sobre-dispersión mediante la expresión (5.1). La diferencia entre ambos es, que el modelo de Poisson supone la hipótesis de que la varianza y la media son iguales, en cambio, el modelo Quasi-Poisson supone que la varianza es una función lineal de la media. Es decir, en el modelo Quasi-Poisson (modelo de Poisson con sobre-dispersión) tenemos que E(Yi) = λ(xi, β) y V ar(Yi) = φ·λ(xi, β) . Por otra parte, en ambos modelos utilizamos la misma función de enlace. El precio que tenemos que pagar por introducir un parámetro de sobre-dispersión en 5.2. CORRECCIÓN DE LA SOBRE-DISPERSIÓN 57 el modelo de Poisson, es que los errores típicos de las estimaciones de los parámetros se multiplican por la raíz cuadrada de b φ , y así los parámetros se vuelven menos signicativos (es decir, el nivel crítico asociado se hace más grande). Sin embargo, debemos destacar que las estimaciones de los coecientes β del modelo no cambian (a diferencia de lo que pasaba en el modelo de regresión Binomial Negativa, en el que sí se modicaban las estimaciones de los coecientes con respecto a las del modelo de Poisson). Es muy importante recalcar que no existe ninguna distribución llamada Quasi-Poisson, no podemos hablar del modelo Quasi-Poisson como lo hacemos con el modelo de Poisson o el modelo Binomial Negativo. Todo lo que hacemos aquí es especicar la relación entre la media y la varianza (función lineal) y la función de enlace (logaritmo). Veamos la aplicación con de este modelo Quasi-Poisson al Ejemplo 5.2 : > mod_quasipoisson_2 <- glm(complaints~hours+revenue, family = quasipoisson(link = log)) > summary(mod_quasipoisson_2) Call: glm(formula = complaints ~ hours + revenue, family=quasipoisson(link = log)) Deviance Residuals: Min 1Q Median 3Q Max -1.9325 -1.1943 -0.2919 0.8651 2.3242 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.2492667 1.0256759 -2.193 0.034040 * hours 0.0015014 0.0003866 3.884 0.000367 *** revenue 0.0046756 0.0029506 1.585 0.120738 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for quasipoisson family taken to be 1.458444) Null deviance: 89.447 on 43 degrees of freedom 58 CAPÍTULO 5. SOBRE-DISPERSIÓN Residual deviance: 61.084 on 41 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 Observamos que la estimación del parámetro de sobre-dispersión es b φ= 1.4584 . Entonces, los errores típicos asociados a los coecientes del modelo han sido multiplicados por √1.4584 = 1.2077 , y así, la mayoría de los parámetros ya no son signicativos. Este es un problema muy habitual que nos encontramos al trabajar con modelos Quasi-Poisson. Nótese que la interpretación de las estimaciones de los parámetros del modelo es análoga a la que habíamos visto para el modelo de Poisson. Por otra parte, este tipo de modelos tienen la ventaja de que podrían tratar la infradispersión, lo que no es posible si utilizamos el modelo de regresión Binomial Negativo. Anexo A: Código Figura 2.1 (Representación gráca de las funciones de masa de probabilidad de una distribución de Poisson para diversos valores del parámetro). > par(mfrow=c(2,2)) > > barplot(dpois(0:57,lambda = 1), # Dibujamos una Poisson(1) + space=0, # No dejamos espacio entre las + # barras + type="h", # Dibujamos puntos + pch=16, + col="blue", # Pintamos en color azul + xlab="y", # Ponemos nombre al eje X + ylab="Masa de probabilidad P(Y=y)", # Ponemos nombre al eje Y + main=expression(paste(lambda, " = 1")), # Ponemos título a la + # gráfica + las=1, # Ponemos los números de los + # ejes en horizontal + bty="l", # Elegimos que solo me pinte la + # línea horizontal abajo y la + # vertical de la izquerda + cex.axis=1.5, # Cambiamos el tamaño de la + # letra de los ejes + cex.lab=1.2, # Cambiamos el tamaño de la + # letra de las descripciones de 59 60 ANEXO A: CÓDIGO + # los ejes + cex.main=2, # Cambiamos el tamaño de la + # letra del título + font.main=4, # Ponemos el título en negrita + # y cursiva + ylim = c(0,0.4)) > > #--Añadimos la Poisson(4), la Poisson(20) y la Poisson(40): > barplot(dpois(0:57,lambda = 4), space=0, type="h",pch=16, col="red", + xlab="y", ylab="Masa de probabilidad P(Y=y)", + main=expression(paste(lambda, " = 4")), las=1, bty="l", cex.axis=1.5, + cex.lab=1.2, cex.main=2, font.main=4, ylim = c(0,0.4)) > > barplot(dpois(0:57,lambda = 20), space=0, type="h", pch=16, col="green", + xlab="y", ylab="Masa de probabilidad P(Y=y)", + main=expression(paste(lambda, " = 20")), las=1, bty="l", cex.axis=1.5, + cex.lab=1.2, cex.main=2, font.main=4, ylim = c(0,0.4)) > > > barplot(dpois(0:57,lambda = 40), space=0, type="h", pch=16, col="yellow", + xlab="y", ylab="Masa de probabilidad P(Y=y)", + main=expression(paste(lambda, " = 40")), las=1, bty="l", cex.axis=1.5, + cex.lab=1.2, cex.main=2, font.main=4, ylim = c(0,0.4)) Ejemplo 2.1 (Regresión de Poisson simple) Lectura de datos > #--Leemos los datos: > install.packages("faraway") > datos <- faraway::esdcomp > View(datos) > #--Hacemos visibles los 6 primeros registros: > head(datos) 61 > attach(datos) Figura 2.3 (Diagrama de dispersión para las quejas recibidas frente a los ingresos) > #--Diagrama de dispersión: > plot(datos$revenue, datos$complaints, type="p", pch=16, cex=1, + col="blue", xlab="Ingresos", ylab="Quejas", + main="Diagrama de dispersión", las=1, cex.axis=1.5, + cex.lab=1.5, cex.main=2) Ajuste del modelo de Poisson simple > mod_poisson_simple = glm(complaints~revenue, family = poisson(link = log), data = datos) > > #--Obtenemos los coeficientes del modelo: > mod_poisson_simple$coefficients > > #--Otra forma: > coef(mod_poisson_simple) > > #--Calculamos las exponenciales de estos coeficientes: > exp(coef(mod_poisson_simple)) Figura 1.3 (Representación del modelo (con variable explicativa revenue y variable respuesta complaints ) ajustado sobre el diagrama de dispersión.) > #--Añadimos el modelo ajustado: > curve(exp(mod_poisson_simple$coefficients[1] + + mod_poisson_simple$coefficients[2] * x), + add = TRUE, lwd=3) 68 BIBLIOGRAFÍA https://dialnet.unirioja.es/servlet/articulo?codigo=4770351, Dialnet. [Consulta: 5 noviembre 2020]. [11] Sheather, S.J. (2009). A modern approach to regression with R , 1st ed., Springer, New York. [12] Venables, W.N. y Ripley, B.D. (2002). Modern Applied statistics with S , 4th ed., Springer, New York, http://www.stats.ox.ac.uk/pub/MASS4. [13] Venables, W. N. y Ripley, B. D. (2010). Modern applied statistics with S , 4th ed., Springer. [14] Vives Brosa, J. El diagnóstico de la sobredispersión en modelos de análisis de datos de recuento (2002). Disponible en: http://hdl.handle.net/10803/5422, Dialnet. [Consulta: 25 enero 2021]. [15] Willis, B. H.; Baragilly, M. y Coomar, D. Maximum likelihood estimation based on NewtonRaphson iteration for the bivariate random eects model in test accuracy meta-analysis , Statistical Methods in Medical Research, 29 (2020), 11971211. Disponible en: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7221455/, Pubmed. [Consulta: 26 noviembre 2020]. [16] Zuur, A.F.; Ieno, E.N.; Walker, N.J. ; Saveliev, A.A. y Smith, G.M. (2009). Mixed Eects Models and Extensions in Ecology with R , 1st ed., Springer-Verlag New York Inc.