Modelo de regresión binominal negativa

FACULTAD DE MATEMÁTICAS DEPARTAMENTO DE ESTADÍSTICA E INVESTIGACIÓN OPERATIVA Trabajo Fin de Grado: MODELO DE REGRESIÓN BINOMIAL NEGATIVA Mario Alcaide Delgado - Grado en Matemáticas Junio 2015 Dirigido por: Joaquín García de las Heras Índice general Introducción 1 Summary 3 1. Modelo Lineal Generalizado. 5 1.1. Componentes del Modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.1. Componente aleatoria. . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.2. Componente sistemática. . . . . . . . . . . . . . . . . . . . . . . . . 8 1.1.3. Funciónenlace. ............................. 8 1.2. Estimación de los parámetros. . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3. Adecuación de los modelos. . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.1. Bondaddeajuste............................. 14 1.3.2. Residuos................................. 17 1.4. Interpretación.................................. 18 2. Modelos para Datos de Recuento 19 2.1. Datosderecuento................................ 19 2.2. Modelos específicos para variables de recuento . . . . . . . . . . . . . . . . 20 2.2.1. Modelo de Regresión de Poisson . . . . . . . . . . . . . . . . . . . . 21 2.2.2. Modelo de Regresión Binomial Negativa . . . . . . . . . . . . . . . 23 2.2.3. Excesodeceros............................. 24 2.3. Problema de sobredispersión . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3. Modelos de Regresión Binomial Negativa 29 3.1. DerivacióndelModelo ............................. 30 3.1.1. Derivación del modelo a partir de una distribución Poisson . . . . . 30 3.1.2. Derivación del modelo como modelo lineal generalizado. . . . . . . . 37 3.2. Estimación.................................... 39 3.3. Adecuacióndelmodelo............................. 43 3.4. Interpretación del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4. Aplicaciones 47 4.1. Simulación.................................... 47 4.2. Aplicación a una base de datos particular. . . . . . . . . . . . . . . . . . . 50 Bibliografía 57 iii Introducción En este trabajo nos centraremos, principalmente, en el estudio del modelo de regresión binomial negativa. Aunque este modelo se puede derivar de un modelo de regresión Poisson en el que, tradicionalmente, se ha introducido un nuevo término o componente aleatoria, también puede ser pensado como un miembro de la familia de modelos lineales generalizados. El planteamiento teórico como miembro de esta familia, permitirá aplicar al modelo binomial negativo los distintos test de bondad de ajuste, análisis de residuos y demás resultados desarrollados para los modelos lineales generalizados. Si bien, el modelo de regresión de Poisson es, en general, el método más útil para modelizar datos de conteo, el requisito de igualdad de media y varianza de esta distribución dificulta su aplicabilidad, ya que en numerosas ocasiones nos encontramos con experimentos cuyos datos presentan mayor varianza que media debido a múltiples causas entre las que destaca la frecuencia alta de ceros. Estos datos son denominados Poisson sobredispersos, pero son más comúnmente designados como, simplemente, sobredispersos. Esta problemática motiva el uso de modelos alternativos que presenten una mayor flexibilidad y que se adapten mejor a este tipo de datos. El presente trabajo se ha estructurado en cuatro capítulos. En el primer capítulo recogeremos diferentes resultados para el análisis estadístico de los modelos lineales generalizados, ya que el modelo de regresión binomial negativo puede considerarse como un caso particular del mismo. Se definen las distintas componentes del modelo lineal generalizado: componente aleatoria, componente sistemática y la función enlace. Trataremos la estimación de los parámetros del modelo mediante el método de máxima verosimilitud. Asimismo, recogemos las distintas técnicas para la evaluación del modelo, contraste de bondad de ajuste, diferentes tipos de residuos,... . Con el fin de analizar la aplicabilidad del modelo de regresión binomial negativa, en el segundo capítulo introduciremos las variables de conteo o recuento, que cuentan con una amplia presencia en diversos ámbitos de la investigación aplicada. Nos centraremos en una breve descripción de diversas alternativas estadísticas para el estudio de este tipo de datos, algunas de ellas son casos particulares de la modelización lineal (modelo de regresión Poisson, modelo de regresión binomial negativa) que nos permite considerar y analizar el comportamiento de las variables de conteo frente a los valores del conjunto de variables explicativas. Finalizamos tratando el problema de sobredispersión que, en general, pueden presentar este tipo de datos. En el tercer capítulo nos centraremos en el análisis y estudio del modelo de regresión 1 2ÍNDICE GENERAL binomial negativa, tratando las diferentes derivaciones del mismo. Una vez formulado, realizaremos un estudio teórico de dicho modelo (estimación de los parámetros, medida de bondad de ajuste, interpretación, etc.) con el fin de su evaluación. El siguiente capítulo está dedicado a dos casos prácticos, donde en ambos se ha utilizado el software R. En primer lugar se ha realizado una breve simulación de un conjunto de datos para el estudio de este modelo y en el segundo caso, se ha aplicado la técnica a un conjunto de datos dados. En ambos casos, se ha aplicado la regresión Poisson y la regresión binomial negativa. En cada una de las aplicaciones vemos cómo la regresión binomial negativa puede ser usada para modelar una variable respuesta de conteo, caracterizada por la presencia de sobredispersión en los datos. Finalizamos con una recopilación de aquellos textos (libros, artículos, etc.) que nos han permitido realizar este trabajo. Introduction In this text we are going to make a study about negative binomial regression model. Although this model can be traditionally derived from a Poisson regression model in which a new term or random component is introduced, this model can be also thought of as a member of the family of generalized linear models. The theoretical approach as a member of this family, will let us apply the different goodness-of-fit tests, analysis of residuals and other results developed for generalized linear models to the negative binomial model. Although the Poisson regression model is the standard method used to model count response data, the equality of mean and variance in this distribution makes difficult its application, because we found too many situations in which we come up with experiments whom data presents much bigger variance tan mean due to many different causes in between we find the high frequency of zero data. This data are termed Poisson overdispersed, but are more commonly simply designated as overdispersed. This problematic motivates the use of alternative models which introduces a greater flexibility and adapt well to this type of data. The present text is structured in four chapters. In the first chapter we collect different results for statistical analysis of generalized linear models, as the negative binomial regression model can be considered as a particular case of them. We defined the different components of generalized linear model: Random component, systematic component and the link function. We will try to estimate the parameters of the model using maximum likelihood method. We also collect the different techniques for the evaluation of the model, goodness-of-fit test, different types of residual, ... . In order to analyze the applicability of the negative binomial regression model, in the second chapter we will introduce the count variables, that have a strong presence in various fields of applied investigation. We will focus on a brief description of various statistical alternatives for the study of this type of data, some of which are special cases of the linear modeling (Poisson regression model, negative binomial regression model) that allows us to consider and analyze the behavior of count variables against the values of the set of covariables. We finished addressing the problem of overdispersion which generally can show this type of data. In the third chapter we focus on the analysis and study of negative binomial regression model, studying its different derivations. Also, once formulated, we will make a theoretical study of the model (estimation of parameters, goodness-of-fit measure, interpretation, etc.) with the purpose of its evaluation. 3 4ÍNDICE GENERAL The next chapter is dedicated to two practical cases where the software R is used. Firstly we make a brief simulation of a set of data for the study of this model and in the second case, the technique has been applied to a data set. In both cases, Poisson regression and negative binomial regression have been applied. In each the applications we see how the negative binomial regression can be used to model a response count variable, characterized by the presence of overdispersion in the data. We conclude with a compilation of those texts (books, articles, etc.) that we have used to make this text. Capítulo 1 Modelo Lineal Generalizado. En la mayoría de los estudios aplicados, el investigador (economista, sociólogo, psicólogo,...) a partir de los resultados del fenómeno real bajo estudio, trata de obtener un modelo simulado de la realidad que permita reconstruir el mecanismo subyacente en el fenómeno real. Un modelo es el resultado de un proceso de abstracción desde la realidad al sistema mátematico. Lo podemos reflejar mediante el siguiente diagrama: A través de la literatura estadística muchos autores, entre los que destacan R.A.Fisher, J.Neyman y D.R.Cox, han tratado el problema de la modelización, es decir, el problema de representar la realidad y su variabilidad e incertidumbre a través de un modelo matemático que permita el estudio, el análisis y la compresión de la misma, con el objetivo de tratar de transformarla, de predecir su futuro o, simplemente, de conocerla. Entre los diversos criterios de clasificación de modelos existentes, resulta de especial interés la clasificación en el ámbito del modelado estadístico, la distinción entre modelos deterministas y no deterministas [Lindsey, 1997]. Los modelos deterministas son aquellos con capacidad para predecir o explicar un fenómeno determinado, de forma que contienen información suficiente para representar sin error la realidad que modelan. Por el contrario, los modelos probabilísticos o estocásticos tratan con la variedad y la incertidumbre. Como ocurre en los sistemas científicos, socioeconómicos o sanitarios , en los que el fenómeno real bajo estudio presenta diversas fuentes de incertidumbre (aleatoriedad). Estos modelos 5 12 CAPÍTULO 1. MODELO LINEAL GENERALIZADO. si(β) = ∂ ∂β li(β, φ, yi) = yixiDi(β)σ−2 i(β)−µi(β)xiDi(β)σ−2 i(β) = xiDi(β)σ−2 i(β)[yi−µi(β)] Otros conceptos importantes que aparecen en la estimación máximo-verosímil del vector de parámetros son: Matriz de información de Fisher esperada: F(β) = Cov s(β) = PiFi(β) Fi(β) = xixt iwi(β):wi(β) = D2 i(β)σ−2 i(β) Matriz información de Fisher observada: Fobs(β) = −∂2l(β) ∂β∂βt Se puede comprobar que F(β) = E(Fobs(β)). Para las funciones enlace canónicas θ(µi) = xt iβ, se simplifica la forma de las matrices de información y la función marcador: s(β) = 1 a(φ)Pixi[yi−µi(β)] F(β) = 1 a(φ)PiV(µi(β))xixt i;F(β) = Fobs(β) La obtención de la estimación máxima-verosímil no se plantea, generalmente, como el cálculo de un máximo global, sino como las soluciones de las ecuaciones de verosimilitud: s(ˆ β)=0 lo que corresponde a un máximo local, es decir, con la matriz de segundas derivadas Fobs(ˆ β)definida positiva. Estas ecuaciones no son, generalmente, lineales por lo que ha de ser resueltas a través de métodos iterativos. En muchos casos se utiliza el método iterativo Fisher scoring o mínimos cuadrados ponderados iterados, cuyas iteraciones vienen definidas, a partir de un estimador inicial ˆ β0, por: ˆ β(k+1) =ˆ β(k)+F−1(ˆ β(k))s(ˆ β(k))k= 0,1,2, ... Sobre este método debemos realizar las siguientes consideraciones: 1. El parámetro escala φse cancela en el término F−1(ˆ β(k))s(ˆ β(k))por lo que se explica el comentario inicial hecho para dicho parámetro. 2. Como punto inicial del proceso iterativo ˆ β(0) se puede utilizar el estimador mínimocuadrático de los puntos (g(yi), xi). 3. El proceso iterativo puede terminar con el criterio 1.3. ADECUACIÓN DE LOS MODELOS. 13 kˆ β(k+1) −ˆ β(k)k kˆ β(k)k< ε para un ε > 0fijado. Otros procedimientos alternativos son el método de Newton-Raphson y métodos Quasi- Newton. Obtenemos a través de este método las estimaciones de los parámetros del modelo ˆ β. Estas estimaciones máximo-verosímiles presentan propiedades de consistencia, eficiencia asintótica y distribución normal asintótica. En caso de que el parámetro dispersión sea desconocido, puede considerarse el siguiente estimador consistente: ˆ φ=1 n−pX i [yi−µi(ˆ β)]2 v(µi(ˆ β)) 1.3. Adecuación de los modelos. Una vez estimados los parámetros, se debe valorar la magnitud de discrepancia entre los datos observados y los esperados por el modelo. Una discrepancia pequeña entre los datos observados y el conjunto de valores estimados ˆµpuede ser tolerable, en cuanto que una discrepancia grande no. De esta manera, si se admite una combinación satisfactoria de la distribución de la componente aleatoria y de la función enlace, el objetivo es determinar cuantos términos son necesarios en la estructura lineal para una descripción razonable de los datos. Un número grande de variables explicativas puede llevar a que un modelo explique bien los datos pero con un aumento de complejidad en su interpretación. Por otro lado, un número pequeño de variables explicativas puede llevar a un modelo de fácil interpretación pero que se ajuste pobremente a los datos. Lo que se desea en realidad es un modelo intermedio. En el proceso del ajuste del modelo se evalúan generalmente un conjunto de modelos que constituyen aproximaciones de los datos observados. Tratamos de construir un modelo intermedio entre el modelo saturado y el modelo nulo, los cuales se caracterizan por : Modelo saturado: En este modelo el número de parámetros estimados es igual al número de observaciones. En estos datos individuales, este modelo no constituye parámetros a estimar, ninguna simplificación de los datos puesto que sólo reproduce lo que está ocurriendo. Modelo nulo: Este es un modelo muy simple, el cual se utiliza como modelo de referencia. Contiene como único parametro al valor esperado µpara todas las observaciones. Habitualmente es incapaz de representar adecuadamente la estructura de los datos, asume un efecto nulo de las variables explicativas El concepto de ”mejor modelo”depende de la finalidad que se persiga. Cuando la finalidad del modelado es del tipo predictivo se seleccionan las variables que expliquen el mayor porcentaje de variabilidad de la respuesta, y para ello se emplean fundamentalmente los criterios estadísticos que veremos a continuación. Por otro lado, cuando la finalidad es explicativa, son los argumentos teóricos los que deben tomar un mayor protagonismo, por 14 CAPÍTULO 1. MODELO LINEAL GENERALIZADO. lo que el proceso de selección debe ser guiado por el investigador y se basa en la especificación de un modelo máximo inicial y de un conjunto sucesivo de modelos restringidos que se comparan mediante un ajuste condicional. A partir de ahí, el proceso de adecuación seguiría, a grandes rasgos, los siguientes pasos: Evaluar los términos de interacción a partir de su significación estadística Analizar la necesidad de mantener variables de control en el modelo. Haciendo uso de criterios de relevancia práctica más que en base a criterios estadísticos, se debe evaluar si: 1. La eliminación de variables confusas sesgará la estimación de los parámetros de interés 2. O bien, en el caso de variables de ajuste, si su supresión del modelo implicaría una pérdida de precisión de las estimaciones. Valorar si las variables explicativas de interés deben permanecer en el modelo. Para ello se deben emplear tanto criterios estadísticos como de relevancia substantiva. Tal como se ha expuesto anteriormente, el objetivo del proceso de modelado es la obtención de un modelo que sea capaz de representar los datos y, al mismo tiempo reducir la complejidad, es decir, atender a los criterios de bondad de ajuste. 1.3.1. Bondad de ajuste. En un modelo lineal generalizado, la bondad del ajuste se puede evaluar de distintas formas, entre las que destacan: La función ó estadístico desviación: D(y; ˆµ) = 2{l(y;y)−l(ˆµ;y)} Que es la distancia entre el logaritmo de la función verosimilitud del modelo saturado (con n parámetros, uno por observación) y el modelo que se está investigando. Un valor pequeño de la desviación indica que para un número menor de parámetros, se obtiene un ajuste tan bueno como cuando se ajusta el modelo saturado. Para probar la adecuación de un Modelo Lineal Generalizado, el valor de la desviación debe ser comparado con el percentil de alguna distribución de probabilidad referente. Si el modelo es correcto el estadístico desviación se distribuye asitóticamente según una X2 n−pcon n−pgrados de libertad [McCullagh y Nelder, 1989]. D(y, ˆµ)∼ X2 n−p El coeficiente de determinación R2: La medida R2es definida como la reducción proporcional en la incertidumbre, debido a la inclusión de los regresores. Bajo ciertas condiciones también puede ser explicada como la varianza explicada por el modelo ajustado. 1.3. ADECUACIÓN DE LOS MODELOS. 15 Este coeficiente viene dado por: R2= 1 −D(y; ˆµ) D(y; ˆµ0) donde D(y; ˆµ)yD(y; ˆµ0)son las funciones de desviación del modelo ajustado y nulo (modelo simple que se usa como referencia) , respectivamente y se verifica que 0≤R2≤1 El estadístico Chi-cuadrado de Pearson: X2= n X i=1 (yi−ˆµi)2 V( ˆµi) donde V(ˆµ)es la función varianza estimada para la distribución de la variable respuesta. En cuanto a la selección del modelo, hacemos uso de los criterios de información para seleccionar entre diferentes modelos. Estos se basan en la comparación de log-verosimilitudes pero penalizando a aquellos modelos con más variables explicativas. Estos criterios son: El Criterio de información Akaike (AIC): En el caso general el AIC es: AIC =k−2ln(ˆ L) donde kes el número de parámetros en el modelo estadístico y ˆ Les el máximo valor de la función verosimilitud para el modelo estimado. Dado un conjunto de modelos candidatos para los datos, el modelo preferido es el que tiene el valor mínimo en el AIC. Este criterio no sólo recompensa la bondad de ajuste, sino también incluye una penalidad, que es una función creciente del número de parámetros estimados. El Criterio de información bayesiano (BIC): BIC =−2lnˆ L+kln(n) donde kes el número de parámetros libres a ser estimados, nes el tamaño de la muestra y ˆ Les el máximo valor de la función de verosimilitud. Está estrechamente relacionado con el AIC y al igual que él resuelve el problema de selección de modelos mediante la introducción de un término penalización para el número de parámetros en el modelo, el término de penalización es mayor en el BIC que el en AIC. Se selecciona aquel modelo con menor valor en el criterio que se utilize. 16 CAPÍTULO 1. MODELO LINEAL GENERALIZADO. Análisis de las variables explicativas A continuación vamos a realizar inferencias sobre el vector de parámetros desconocidos β de dimensión p, dependiendo de los objetivos o de la metodología a usar para el análisis de datos. La mayoría de las cuestiones que se plantean en la realidad pueden ser formuladas a través de una hipótesis lineal de la forma Cβ, siendo Cuna matriz de rango total s≤p siendo pla dimensión del vector de parámetros βyξun vector de constantes conocido de dimensión s H0:Cβ=ξ H1:Cβ6=ξ En general, se recogen los siguientes procedimientos para este contraste: Estadístico de razón de Verosimilitud. Definido por: ΛRV =−2{l(˜ β;y)−l(ˆ β;y)} que compara el máximo del logaritmo de la verosimilitud con el máximo obtenido bajo la restricción definida por H0 Estadístico de Wald, que se basa en la distibución normal asintótica del vector ˆ β. Se define por : ξW= [Cˆ β−ξ]T[CF−1(ˆ β)C0][Cˆ β−ξ] Determina la distancia ponderada entre el estimador Cˆ βyCβy su valor determinado por la hipótesis nula, donde F−1(ˆ β)denota la estimación de la matriz de información de Fisher de ˆ β. Estadístico score, obtenida a partir de la función score. Se basa en el hecho de que la función score se anula en el estimador de máxima verosimilitud, por lo que la evaluación de ésta en el estimador obtenido bajo restricción lineal, el resultado será significativamente diferente de cero si la hipótesis nula no es cierta. Así se utiliza la distancia poderada de s(˜ β)a cero, es decir: ξSR = [s(˜ β)]TF−1(˜ β)s(˜ β) donde F(ˆ β)es la matriz de covarianza asintótica estimada bajo H0:Cβ=ξ Asintóticamente y bajo hipótesis nula, los tres estadísticos definidos anteriormente se distribuyen según una ley Chi-cuadrado con sgrados de libertad X2 s. Sea el caso especial de testar la significación de un subconjunto de covariantes: H0:βr= 0 H1:βr6= 0 1.3. ADECUACIÓN DE LOS MODELOS. 17 donde βres un subvector de β. Los estadísticos anteriores se reducen a: ξW= (ˆ β)TF−1 r(ˆ β)ˆ βr ξSR =sT r(˜ β)F−1 r(˜ β)sr(˜ β) Donde Fres la submatriz de Fcorrespondiente a los elementos de βrysrel subvector de scorrespondiente a dichos elementos. Para las hipótesis relativas a un único coeficiente βi, el estadístico de Wald es el más utilizado. Éste coincide con el cuadrado del estadístico t2: tj=ˆ β ˆajj donde ˆajj es el elemento j-ésimo diagonal de la matriz de covarianzas asintóticas F(ˆ β)de ˆ β. Para hipótesis relativas a varios coeficientes, la razón de Máxima Verosimilitud es preferida por ser un test uniformemente más potente [Codeiro, 2000]. Regiones de confianza para β Los intervalos de confianza y/o regiones de confianza pueden contruirse usando cualquiera de estos estadísticos. Usando,por ejemplo el estadístico de Wald, una región de confianza para βcon un nivel de confianza 100(1 −α)% viene dada por: {β∈Rp|(ˆ β−β)T[ˆ V ar(ˆ β)]−1(ˆ β−β)<X2 p,1−α} 1.3.2. Residuos En la práctica, puede ocurrir que aún escogiendo cuidadosamente un modelo y después al ajustarlo a un conjunto de datos, el resultado sea insatisfactorio. Las desviaciones sistemáticas se originan por haber escogido inadecuadamente la función de enlace o las variables explicativas incluídas en el modelo. Las discrepancias aisladas pueden ocurrir debido a puntos extremos, o porque estos realmente son erróneos como resultado de una mala transcripción de los datos o por factores no controlados en el momento de la toma de datos. La verificación de la adecuación del modelo es un requisito fundamental que se realiza sobre el conjunto de datos para analizar posibles desviaciones de las suposiciones hechas para el modelo, así como la existencia de observaciones extremas con alguna interferencia desproporcionada en los resultados del ajuste. Como en la regresión lineal, los residuos son los utilizados para verificar dicha adecuación del modelo. Expresan la discrepancia entre una observación y su valor ajustado, y también pueden indicar la presencia de valores anómalos o discordantes que puedan requerir de una investigación más detallada. Entre otros residuos los más destacados son: El residuo básico: Definido como la diferencia entre el valor observado, yi, de la variable respuesta y el valor ajustado, ˆyi, por el modelo. rb i=yi−ˆyicon i= 1, ..., n 18 CAPÍTULO 1. MODELO LINEAL GENERALIZADO. EL residuo de Pearson: Es la contribución individual al estadístico X2de Pearson, se define como: rp i=yi−ˆµi qˆ φV ar(ˆµi)i= 1,2, ..., n siendo ˆ φun estimador consistente del parámetro escala φ. Y su versión estudentadizada viene dada por: rp i=yi−ˆµi qˆ φV ar(ˆµi)(1 −hi)i= 1,2, ..., n siendo hiel elemento diagonal de la matriz de proyección, H , donde : H=W1/2X(XTWX)−1XTW1/2 con W una matriz diagonal, cuyos elementos de la diagonal principal vienen dados por: wi=1 V ar(µi)(∂µi ∂η )2 La ventaja de usar este residuo estudentadizado frente al anterior es que capta mejor la variabilidad de los datos, debido a que usa el valor de hi, el cual es útil para medir la influencia de la i-ésima observación. EL residuo desviación, que se define como: rD i=sign(yi−ˆµi)qdii= 1,2, ..., n donde dies llamado la componente desviación, di= 2(l(yi, yi)−l(ˆµi, yi)) Y su versión estudentadizada: rD0 i=rD i qˆ φ(1 −hi) donde hies el i-ésimo elemento de la diagonal de la matriz H y ˆ φes la estimación del parámetro escala φ. Es importante destacar que los errores de especificación mencionados pueden provocar una violación del supuesto distribucional de relación media-varianza de diversos modelos lineales generalizados, como es el caso del modelo de regresión de Poisson. 1.4. Interpretación Una vez se ha obtenido el modelo adecuado haciendo uso de los criterios antes mencionados de bondad de ajuste, estudio de residuos,... el proceso de modelado se cierra con la interpretación del modelo. Recordando que la transformación provocada por la aplicación de una función enlace (excepto en la función de enlace identidad) da lugar en la mayoría de los casos a una ecuación del modelo expresada en términos multiplicativos, en la que la interpretación de los parámetros se realiza en términos del factor de cambio en el valor esperado para un incremento unitario de las variables explicativas. Capítulo 2 Modelos para Datos de Recuento Las variables de recuento cuentan con una amplia presencia en diversos ámbitos de la investigación aplicada tanto en las Ciencias Sociales como en las Ciencias de la Salud. En este sentido, encontramos ejemplos de investigaciones aplicadas con variables de recuento en disciplinas como Demografía [Merklesson y Roth, 2000], Farmacología [Lindsey, 2001], Relaciones Laborales [Sturman, 1999]ó Criminología [Osgood, 2000], por citar algunos. Mucho más numerosas resultan las investigaciones que contienen variables de recuento en Medicina [Böhning, 1994], Ciencias Políticas [King, G., 1989] y Ciencias Económicas [Meliciani, 2000]. De hecho, estas tres disciplinas merecen ser consideradas aparte, puesto que no sólo cuentan con una extensa aplicación de investigaciones con variables de recuento, sino que han hecho valiosas aportaciones en el tratamiento estadístico de este tipo de variables [Cameron y Trivedi, 1986], [Winkelmann, 2000]. En este capítulo nos centraremos en los modelos de datos de recuento o conteo, un caso particular de modelización lineal que nos permite considerar y analizar el comportamien- to de las variables de conteo frente a los valores del conjunto de variables explicativas. Introduciremos los distintos tipos de modelos específicos para este tipo de variables y las distintas características que estos presentan. 2.1. Datos de recuento Se denominan variables de recuento (count data) a aquellas variables que determinan el número de sucesos o eventos que ocurren en una misma unidad de observación en un intervalo espacial o temporal definido. A partir de esta definición propuesta por [Lindsey, 1995], se derivan dos características principales de una variable recuento, que la diferencia de una variable cuantitativa continua, estas son, su naturaleza discreta y no negativa. La variable a tratar Y, que toma los valores 0,1,2, ..., se caracteriza por tomar infinitos números de valores que podemos ordenar en orden creciente, y cuya probabilidad va en descenso a medida que sea mayor el valor de la variable. Ejemplos de tales eventos o conteos pueden ser: 19 20 CAPÍTULO 2. MODELOS PARA DATOS DE RECUENTO Conteos en el tiempo Número de accidentes de tráfico en un tramo de cierta carretera en un mes. Número del registro de partículas de una desintegración radioactiva por segundo. Número de patentes solicitadas por una empresa en 1 año. Conteos en el espacio Número de células sanguíneas en una muestra de sangre. Número de árboles infectados por hectárea en un bosque. Número de artículos publicados por una revista. Los modelos de datos de conteo se caracterizan porque no tienen un límite superior natural, toman el valor cero ( en un porcentaje no despreciable) para algunos miembros de la población y suelen tomar pocos valores distintos. El objetivo,como hemos señalado, es analizar Ycomo función de las variables explicativas X1, ..., Xp, a través de: E(Y/X=x) podiendo ser cada una de las variables explicativas de cualquier tipo. 2.2. Modelos específicos para variables de recuento No todos los modelos de predicción son aplicables a este tipo de variables, pues pueden originarse problemas como pérdida de información o de inconsistencia. La problemática que surgiría al construir un modelo cualquiera para este tipo de variables sería: Modelo de regresión lineal 1. Las predicciones de Ypueden salirse del rango de valores en el que está definido. 2. Las estimaciones pueden ser inconsistentes. 3. Puede tener validez para hacer una exploración previa de las relaciones. 2.2. MODELOS ESPECÍFICOS PARA VARIABLES DE RECUENTO 21 Modelo de elección binaria Si la variable Ytoma muchos valores, plantear un modelo de elección binaria nos conduce a una pérdida de eficiencia (porque perdemos información) ya que agregamos todos los valores mayores que uno en un solo valor. Modelos ordenados 1. Si la variable Ytoma muchos valores o si tiene pocas observaciones en algunas de las variables, es necesario agrupar si queremos estimar un modelo ordenado. Esto en determinados contextos, puede suponer pérdida de información. 2. Un aspecto positivo de los modelos ordenados es que podemos utilizarlos cuando queremos analizar variables que toman valores enteros negativos A continuación presentamos un breve resumen de algunos de los modelos propuestos que ofrecen una mejor aproximación para este tipo de datos, son los Modelos para variables de recuento y estos son : →Modelo de regresión de Poisson →Modelo de regresión Binomial Negativa →Modelo en dos partes →Modelo con exceso de ceros. Donde los modelos de regresión de Poisson y de regresión binomial negativa, a diferencia del resto de los mencionados, pertenecen a la familia de modelos lineales generalizados. Los modelos para variables de conteo se encuentran en la confluencia entre el modelo lineal generalizado y el estudio de las variables de recuento. De entre estos modelos destaca, por su papel como modelo de referencia en el estudio de las variables de recuento, el modelo de regresión de Poisson. 2.2.1. Modelo de Regresión de Poisson La distribución habitual para procesos de conteo es la Distribución de Poisson Sea Y∼P(µ) con función de probabilidad: P(Y=y) = e−µµy y!y= 0,1,2, .. de parámetro µ > 0 la media y la varianza viene dada por: E(Y) = µ V ar(Y) = µ 28 CAPÍTULO 2. MODELOS PARA DATOS DE RECUENTO Capítulo 3 Modelos de Regresión Binomial Negativa Como hemos visto en el capítulo anterior, a pesar de que el modelo de referencia en estudios de variables de recuento es el modelo de regresión de Poisson (MRP), éste presenta varios problemas a la hora de tratar con datos en los que la media y la varianza condicionadas no coinciden, en concreto, con datos que presentan sobredispersión. Una forma de relajar esta restricción de igualdad media-varianza del Modelo de Regresión de Poisson es especificar una distribución que permita un modelado más flexible. En este sentido, el modelo paramétrico estandar para datos de recuento con presencia de sobredispersión es el Modelo de Regresión Binomial Negativa (MRBN). A pesar de las ventajas potenciales de este modelo para datos de conteo, son pocos los estudios publicados que han usado esta distribución. Esos estudios se han centrado mayormente en los recuentos de aves raras, accidentes de tráfico, la utilización de servicios de salud, lesiones neurológicas y leucocitos [Welsh, Cunningham y Chambers, 2000], [Abdel y Radwan, 2000], [Sormani, Bruzzi y Miller, 1999] y [Finch y Chen, 1999]. En el artículo sobre un ensayo clínico realizado por [Byers, Allore, Gill y Peduzzi 2003] podemos ver un estudio practico más detallado del modelo de regresión binomial negativa. El modelo de regresión binomial negativa es un modelo estadístico atípico, en el sentido de que los investigadores se refieren a éste como un único modelo, como, por ejemplo, el modelo de regresión Poisson, la regresión logística,... cuando en realidad existen distintos modelos binomiales negativos, que dependerán del tipo de problema de fondo que se esté abordando. [Boswell y Patil, 1970] identificaron 13 tipos distintos que derivan de la distribución binomial negativa mientras que otros autores argumentan que existen más. Podemos describir la distribucíon binomial negativa como aquella variable que estudia la probabilidad de observar un número determinado de fracasos, y, antes del r-ésimo éxito en una serie de experimentos Bernoulli independientes. Bajo tal descripción de r, este sería un entero positivo, sin embargo, no hay razón matemática de peso para limitar este parámetro para números enteros; solo limitar r como positivo. Diremos que la variable aleatoria de conteo Yicon i= 1, ..., n sigue una distribución Binomial Negativa de parámetros ryp,Yi∼BN(r, p), con función de probabilidad dada 29 30 CAPÍTULO 3. MODELOS DE REGRESIÓN BINOMIAL NEGATIVA por: P(Yi=yi) = yi+r−1 r−1!pr(1 −p)yi(3.1) donde 0<p<1,r > 0,yi= 0,1,2.... El valor esperado y la varianza vienen dados por: E(Yi) = r(1 −p) p(3.2) V(Yi) = r(1 −p) p2(3.3) estableciéndose entre ambos la siguiente relación: V(Yi) = 1 pE(Yi) Como 0< p < 1, se verifica que V(Yi)> E(Yi)lo que justifica la aptitud natural de esta distribución para modelar datos que se caracterizan por la existencia de sobredispersión. En este capítulo nos centraremos principalmente en la derivación del modelo de regresión binomial negativo, trataremos los distintos métodos para la estimación de los parámetros del mismo y estudiaremos la bondad de ajuste, los residuos asociados,... del modelo con el fin de evaluarlo. 3.1. Derivación del Modelo Podemos recoger las distintas formas de derivar el modelo de regresión binomial negativo (MRBN) en dos orígenes diferenciados. Por una parte el MRBN se puede derivar como una distribución Poisson compuesta con una Gamma, en el cual la distribución Gamma es usada para ajustar los datos Poisson que presentan sobredispersión. De esta forma se deriva el modelo tradicional binomial negativo, el cual es comunmente simbolizado como modelo de regresión BN2, [Cameron y Trivedi, 1998]. Por otra parte el MRBN puede ser pensado como miembro de la familia exponencial de distribuciones y, por tanto, ser considerado como un modelo lineal generalizado. Este enfoque sólo se tiene si el parámetro de dispersión o heterogeneidad es introducido en la distribución como una constante [Hilbe 1993]. Tal interpretación permite a los investigadores aplicar al MRBN los test de bondad de ajuste, análisis de residuos y cualquier otro estudio que hemos desarrollado para los Modelos Lineales Generalizados. 3.1.1. Derivación del modelo a partir de una distribución Poisson Partiendo de un modelo de regresión Poisson (MRP), podemos derivar el modelo de regresión binomial negativo (MRBN) siguiendo dos enfoques distintos. 3.1. DERIVACIÓN DEL MODELO 31 Un primer enfoque, el más común, introducido por [GTM, 1984], en el que la variable respuesta sigue una distribución Poisson, cuya media está especificada de forma incompleta debido a una situación de heterogeneidad no observada, para suplir dicha situación se introduce un nuevo término error. Un segundo enfoque, [Hausman, Hall y Griliches, 1984], que supone que la variable respuesta sigue una distribución Poisson en la que su media no se considera como un parámetro fijo, sino que se interpreta como un parámetro estocástico que varía aleatoriamente como una distribución Gamma. En cualquiera de los dos casos, el modelo Poisson es insuficiente para la modelización de los datos, debido a la heterogeneidad que estos presentan. Ambas motivaciones conducen a una distribución binomial negativa. A continuación recogemos la función de densidad de la variable Gamma: T∼G(τ, ω)⇒f(t) = 1 wτΓ(τ)tτ−1et ω con t, ω, τ > 0 y siendo la función Gamma: Γ(v) = Z∞ 0tv−1e−tdt 1er Caso Mientras que en el MRP la media condicional de Yies: µi=exp(xiβ) en el modelo que se va a realizar esta media es reemplazada por la variable aleatoria ˜µi, verificándose la siguiente ecuación estocástica: ˜µi=exp(xiβ+εi) bajo la hipótesis de que el término error considerado, εi, está incorrelacionado con las variables explicativas. El término de error εipuede ser el resultado del efecto conjunto de variables no incluidas en el modelo o bien una fuente de aleatoriedad intrínseca. Sea cual sea su origen, εirepresenta la heterogeneidad no observada de los datos. En el MRP, la variación en µies introducida a través de la heterogeneidad observada, de forma que diferentes valores de xidan diferentes valores de µi. Así, todos los individuos con los valores xitienen la misma µi. En este modelo propuesto, la variación en ˜µes debida tanto a la variación en xentre los individuos, como a la heterogeneidad no observada introducida a través de ε. Para una combinación de valores en las variables independientes, existe una distribución de diversas ˜µien lugar de una µiúnica. La relación entre ambas medias viene dada por: 32 CAPÍTULO 3. MODELOS DE REGRESIÓN BINOMIAL NEGATIVA ˜µi=exp(xiβ)exp(εi) = µi·exp(εi) = µiδi siendo δiequivalente a exp(εi). La aproximación del MRBN dependerá de la asignación ó hipótesis acerca de la media del término de error. La asignación más conveniente y usual es que: E(δi) = 1 V(δi) = ωi que implicará que el recuento esperado después de añadir la nueva fuente de variación sea la misma que para el MRP: E(˜µi) = E(µiδi) = µi·E(δi) = µi Bajo estas suposiciones se verifica que: P(Yi=yi|xi, δ) = exp(−˜µi)˜µyi i yi!=exp(−µiδi)(µiδi)yi yi! por lo que se vuelve a obtener una distribución Poisson. Sin embargo, puesto que el error δies desconocido no podemos calcular P(yi|xi, δi). Para calcular P(yi|xi)sin tener en cuenta δi, promediamos P(yi|xiδi)por la probabilidad de cada valor de δi. Si ges la función de densidad de probabilidad de δi, entonces la densidad marginal de Yide i-ésimo individuo puede ser obtenida integrando con respecto a δi: P(Yi=yi|xi) = Z∞ 0[P(Yi=yi|xiδi)]g(δi)dδi=Z∞ 0 e−exp(xiβi)δiexp(xiβi)yiδyi yi!g(δi)dδi (3.4) Esta expresión define la distribución de Poisson compuesta [Cameron y Trivedi, 1986]. Tal como indican estos mismos autores, las distribuciones de Poisson compuestas proporcionan una generalización natural de los modelos de Poisson básicos y su aplicación obedece generalmente a una necesidad de mayor flexibilidad, especialmente en situaciones de sobredispersión. La ecuación de la distribución de Poisson compuesta (3.4) calcula la probabilidad de la variable respuesta Yicomo una mixtura de dos distribuciones de probabilidad. Asimismo, la forma de (3.4) depende de la selección de g(δi), es decir, de la función de densidad de probabilidad que se asuma para δi. En este sentido, suponiendo que la variable error δisigue una distribución G(τ, ω), con función de densidad: g(δi) = 1 ωτΓ(τ)δτ−1exp(−δi ω)para τ, ω > 0 3.1. DERIVACIÓN DEL MODELO 33 El resultado de la ecuación de la regresión de Poisson compuesta (3.4) conduce a una distribución binomial negativa, es decir, tal como indica [Poortema, 1999], la distribución binomial negativa es la distribución compuesta resultante si la conjugada de la distribución Poisson, es decir, la distribución Gamma, es utilizada para la composición. En resumen, bajo estas consideraciones y partiendo de que δi∼G(τ, ω)podemos probar que Yi|xi∼BN τ, 1 1 + ωµi! con función de probabilidad: P(Yi=yi/xi) = Γ(yi+τ) Γ(yi+ 1)Γ(τ) 1/ω 1/ω +µi!τ µi 1/ω +µi!yi para yi= 0,1,2... (3.5) y cuyo valor esperado es: E(Yi|xi) = τωµi y la varianza condicional: V ar(Yi|xi) = τωµi(1 + ωµi)(3.6) 2oCaso En este enfoque se supone que el parámetro µino es un parámetro fijo, sino que este se distribuye de forma aleatoria según una ley Gamma. Se supone que: Yi|xi, µi∼P(µi) µi|xi∼G(τi, ω), τi, ω > 0 y sabiendo que P(Yi=yi/xi) = Z∞ 0P(Yi=yi|xi, µi)f(µi|xi)dµi demostramos que Yi|xisigue una distribución binomial negativa: P(Yi=yi|xi) = Z∞ 0 e−µiµyi i yi! 1 wτiΓ(τi)µτi−1 ieµi ωdµi =1 ωτiΓ(τi)yi!Z∞ 0e−(1+ 1 ω)µiµyi+τi−1 idµi Esta última integral corresponde a Γ(yi+τi) 1 + 1 ωyi+τiyyi! = Γ(yi+ 1), por lo que: P(Yi=yi|xi) = Γ(yi+τi) Γ(yi+ 1)Γ(τi)    1 ω1 + 1 ω    τi    1 1 + 1 ω     yi =Γ(yi+τi) Γ(yi+ 1)Γ(τi)1 1 + ωτiω 1 + ωyi 34 CAPÍTULO 3. MODELOS DE REGRESIÓN BINOMIAL NEGATIVA Comparándola con (3.1), se verifica que Yi|xi∼BN τi,1 1 + ω, y haciendolo con (3.2) y (3.3) concluimos que: E(Yi) = ωτi(3.7) V(Yi) = ω(1 + ω)τi(3.8) El modelo de regresión se establece como: E(Yi) = exp(xiβ) La relación entre la varianza y el valor esperado condicional es lineal y viene dada por la expresión: V(Yi|xi) = (1 + ω)E(Yi|xi) Parametrizaciones más comunes A continuación, presentaremos las parametrizaciones más comunes del modelo de regresión binomial negativa, dependiendo del enfoque que sigamos. Siguiendo el primer enfoque señalado, [Long, 1997] afirma que la asignación más común es que δisigue una distribución G(1 ω, ω), con función de densidad: g(δi) = 1 ω1/ωΓ(1/ω)δ 1 ω−1 iexp(−δi ω)para ω > 0 [Johnson,Kotz y Balakrishnan 1994] demuestran que si δisigue una distribución Gamma, entonces se verifica también que E(δi)=1yV ar(δi) = ω, como habiamos considerado previamente. En resumen, bajo estas consideraciones y partiendo de que δi∼G(1 ω, ω)podemos probar que Yi|xi∼BN 1 ω,1 1 + ωµi! con función de probabilidad: P(Yi=yi/xi) = Γ(yi+ 1/ω) Γ(yi+ 1)Γ(1/ω) 1/ω 1/ω +µi!1/ω µi 1/ω +µi!yi para yi= 0,1,2... (3.9) y cuyo valor esperado es el mismo que para la distribución de Poisson: E(Yi|xi) = exp(xiβ) = µi sin embargo, la varianza condicional si difiere en relación a la de la distribución de Poisson: V ar(Yi|xi) = µi(1 + ωµi) = exp(xiβ) + ωexp(xiβ)2(3.10) Puesto que µi>0yω > 0, la varianza condicional de Yidel individuo i-ésimo en el MRBN será mayor que la media condicional exp(xiβ). Se observa que a medida que ωtiende a 3.1. DERIVACIÓN DEL MODELO 35 cero la varianza tiende a µiV ar(Yi|xi)→µi, obteniéndose la hipótesis de equidispersión. Por otro lado, el incremento de la frecuencia relativa de valores de recuento altos y bajos ocasiona una varianza condicional elevada en Yi. De esta forma, en una situación de sobredispersión, la distribución binomial negativa corrige, especialmente, la probabilidad asociada a valores bajos de recuento que habitualmente presentan un ajuste deficiente a través del MRP. Siguiendo el segundo enfoque señalado, [Cameron y Trivedi, 1986] construyó un modelo más general que contempla varias relaciones posibles entre la media y varianza condicionada. Suponiendo que: Yi|xi, µi∼P(µi) µi|xi∼G(vi,µi vi ), vi, µi>0 Y usando un razonamiento similar al anterior se obtiene que: Yi|xi∼BN vi,vi µi+vi! con media y varianza condicional dada por: E(Yi|xi) = µi(3.11) V(Yi|xi) = µi+1 vi µ2 i(3.12) Los modelos de regresión se obtienen definiendo: µi=exp(xiβ)(3.13) vi=1 α[exp(xiβ)]a, α > 0, a = 0,1,2, ... (3.14) Sustituyendo en el valor de vien (1.10) se obtiene: V ar(Yi|xi) = exp(xiβ) + α[exp(xiβ)]−a[exp(xiβ)]2 =exp(xiβ) + α[exp(xiβ)]2−a =µi+αµ2−a i Esta nueva expresión origina una amplia variedad de formas para relacionar la varianza y la media. Bajo esta expresión ya se pueden empezar a considerar dos modelos principales del tipo Binomial Negativa que recogemos a continuación. 36 CAPÍTULO 3. MODELOS DE REGRESIÓN BINOMIAL NEGATIVA El modelo Binomial Negativa I (BN1): (relación lineal entre la varianza y la media) Es el modelo resultante si consideramos a= 1, entonces vi=µi α. Luego: Yi|xi∼BN µi α,1 1 + α con una función de probabilidad: P(yi/xi) = Γ(yi+α−1µi) Γ(yi+ 1)Γ(α−1µi)1 1 + α µi αα 1 + αyi para yi= 0,1,2... y su varianza condicional es: V ar(Yi|xi) = µi+αµi= (1 + α)µi El modelo Binomial Negativa II (BN2): (relación cuadrática entre varianza y media) Es el modelo resultante si consideramos a= 0, entonces vi=1 α. Luego: Yi|xi∼BN 1 α,1 1 + αµi! con función de probabilidad: P(yi/xi) = Γ(yi+α−1) Γ(yi+ 1)Γ(α−1) α−1 α−1+µi!α−1 µi α−1+µi!yi para yi= 0,1,2... (3.15) Al hacer efectiva la igualdad vi=α−1, se hace evidente que al incrementar α, que es conocido como el parámetro de dispersión, aumenta la varianza condicional de Y: V ar(Yi|xi) = µi1 + µi α−1=µi(1 + αµi) = µi+αµ2 i(3.16) Se observa que si el parámetro de dispersión α= 0 habría equidispersión, pues V ar(Yi|xi) = µi+αµ2 i=µi La función de probabilidad (3.15) y la varianza (3.16) caracterizan al modelo de regresión binomial negativo de tipo 2 (NB2) al que se conoce como el modelo tradicional binomial negativo. De la forma en que hemos deducido el modelo binomial negativo, se es bastante exigente en términos de supuestos, en particular al exigir que el error de especificación siga una distribución Gamma. Es por ello que, por lo general, se comienza estimando el modelo más sencillo el de Poisson, sujeto a una serie de pruebas y, sólo si se considera conveniente, 3.1. DERIVACIÓN DEL MODELO 37 analizar el resto de modelos, que también deberán ser evaluados. Además de los modelos BN1 y BN2 considerados anteriormente, algunos autores como [Cameron y Trivedi, 1986], o [Winkelmann y Zimmermann, 1991] proponen un modelo binomial negativo más general, el denominado hipermodelo Binomial Negativo K(BN k), en el cual V ar(Yi|xi) = µi+αµ2−k i. Consideramos que esta extensión queda fuera de nuestro estudio. 3.1.2. Derivación del modelo como modelo lineal generalizado. Él método que abordamos a continuación, es el método más usual para derivar cualquier miembro de la familia de modelos lineales generalizados. La distribución binomial negativa es un miembro de la familia exponencial siempre que el parámetro de dispersión o heterogeneidad sea introducido en la distribución como una constante, por tanto la podemos usar como distribución de la componente aleatoria del modelo lineal generalizado. Como se ha comentado anteriormente, otra decisión importante es la selección del enlace y por lo tanto, según la elección de éste saldrán a su vez diferentes modelos. En la literatura se suelen elegir dos enlaces diferentes el enlace canónico y el enlace logarítmico. La ventaja de utilizar éste último se debe a que nos permite hacer una comparación con el modelo de regresión Poisson. Componente aleatoria La forma de la función de probabilidad binomial negativa con la que comenzamos se expresa como: Y∼BN(r, p) P(Y=y) = y+r−1 r−1!pr(1 −p)y(3.17) con y= 0,1, ... que se puede expresar como miembro de la familia exponencial negativa con la siguiente estructura: P(Y=y) = exp (yln(1 −p) + r(ln(p)) + ln y+r−1 r−1!!) (3.18) Aplicando los resultados obtenidos en el capítulo 1, se obtiene que: θ=ln(1 −p)⇒p= 1 −exp(θ)(3.19) b(θ) = −rln(p)⇒ −r(1 −exp(θ)) (3.20) a(φ),es el escalar que tomamos como 1 (3.21) La primera y la segunda derivada de b(θ), con respecto de θ, respectivamente dan lugar a la función media y varianza: 44 CAPÍTULO 3. MODELOS DE REGRESIÓN BINOMIAL NEGATIVA El residuo de Pearson: rp i=(yi−ˆµi) qˆµi+ ˆαˆµ2 i siendo su versión estudentadizado: rpt i=(yi−ˆµi) ˆ φqˆµi+ ˆαˆµ2 i siendo ˆ φun estimador consistente del parámetro escala φ. Su versión estandarizada viene dada por: rpst i=(yi−ˆµi) q(1 −hi)(ˆµi+ ˆαˆµ2 i) siendo hiel i-ésimo elemento de la matriz de proyección. El residuo desviación: rd i=sgn(yi−ˆµi)√D siendo Del estadístico desviación para el MRBN anteriormente señalado. Su versión estudentadizado viene dada por: rdt i=sgn(yi−ˆµi)√D ˆ φ Y su versión estandarizada: rdst i=sgn(yi−ˆµi)√D q(1 −hi) Además, consideraremos el residuo Anscombe, desarrollado para la binomial negativa NB2 [Hilbe, 1994], [Hardin y Hilbe, 2001]. El residuo Anscombe [Ascombe, 1972] toma valores cercanos a los de la desviación estandarizada. Aunque hay veces en la que este no es el caso, y el residuo Ascombe realiza una mejor interpretación que el residuo desviación. El residuo Ascombe trata de normalizar la diferencia entre valores observados y ajustados de manera que la heterogeneidad en los datos llega a ser fácilmente identificable. Los residuos Anscombe usan la función de varianza del modelo Poisson :V(µ) = µ B2 : V(µ) = µ(1 + αµ) Este residuo se define de la siguiente forma: rA i=A(yi)−A(ˆµi) A0(ˆµi)qV(ˆµi) 3.4. INTERPRETACIÓN DEL MODELO 45 donde A(.) = Zdµi/V 1/3(µi) El residual Anscombe calculado para un modelo de regresión binomial negativa viene dado por: rA={3/ˆα{(1 + ˆαyi)2/3−(1 + ˆαˆµi)2/3}+ 3(y2/3 i−ˆµ2/3 i)} 2(ˆαˆµ2 i+ ˆµi)1/6 3.4. Interpretación del modelo Una vez se ha obtenido el modelo adecuado de regresión binomial negativa haciendo uso de los criterios antes mencionados de bondad de ajuste, estudio de residuos,... el proceso de modelado se cierra con la interpretación del modelo. Es importante recordar que la transformación provocada por la aplicación de una función enlace da lugar, en la mayoría de los casos, a una ecuación del modelo expresada en términos multiplicativos, en la que, como señalamos en el primer capítulo, la interpretación de los parámetros se realiza en términos del factor de cambio en el valor esperado para un incremento unitario de las variables explicativas. Particularmente el modelo tradicional de regresión binomial negativa viene dado por: ln(E[Y]) = ˆ β0+X1ˆ β1+X2ˆ β2+.... +Xpˆ βp cuya función enlace es la logarítmica. Dependiendo del valor del βestimado, cada variable explicativa puede influir de forma positiva o negativa a la variable de interés. Al igual que en la regresión logística si una variable explicativa es cualitativa, se deben introducir las variables ”dummys”. En este caso, si la variable explicativa Xkes una variable dummy, entonces se tiene la siguiente relación: E(Yi/Xik = 1) E(Yi/Xik = 0) =exp(ˆ βk) la media condicional es exp(ˆ βk)veces mayor si Xktoma valor 1 en vez de 0. En el caso de que sea una variable explicativa Xkcon k= 1...p una variable continua cada unidad de aumento de esta variable produce un aumento de ˆ βjdel log valor esperado de la variable respuesta. Esta interpretación proporciona una restricción a valores positivos de la variable respuesta, adecuándose al tipo de datos de recuento. 46 CAPÍTULO 3. MODELOS DE REGRESIÓN BINOMIAL NEGATIVA Capítulo 4 Aplicaciones En este capítulo trataremos dos ejemplos prácticos trabajados en R, teniendo como fin el correcto uso del modelo de regresión binomial negativa (MRBN). En cada uno de ellos veremos como la regresión binomial negativa puede ser usada para modelar una variable respuesta de conteo, caracterizada por la presencia de sobredispersión en los datos. En el primer caso se ha realizado una simulación de datos que se ajustan a una distribución binomial negativa en la que estudiaremos las distintas mejoras que presenta el modelo MRBN frente al modelo de regresión Poisson (MRP). Y en el segundo caso aplicamos la regresión binomial negativa a una base de datos dada con el posterior análisis de los resultados obtenidos. En ambos ejemplos haremos uso de los distintos comandos implementados en R para poder realizar tales estudios, interpretando posteriormente los resultados obtenidos. 4.1. Simulación En este apartado vamos a realizar una breve simulación a través de R, comenzaremos construyendo un conjunto simulado de datos que satisfaga ciertas especificaciones para demostrar el efecto del ajuste del modelo de regresión binomial negativa. Implementando en R las librerías apropiadas generamos, a través de la función ”nb2_syn”, un conjunto de datos aleatorios con nobservaciones de forma que la variables respuesta y explicativas quedan definida de acuerdo al modelo tradicional de regresión binomial negativa con parámetro de dispersión αy tal que los coeficientes de las variables explicativas del modelo sean (β0, β1, ..., βp). En nuestro estudio particular generaremos n= 500 observaciones, 3 variables explicativas, siendo los coeficientes (β0, β1, β2, β3) =(2,0.75,-0.3,2), y una variable respuesta, con parámetro de dispersión α= 0,5. Las instrucciones vienen dadas por: 47 48 CAPÍTULO 4. APLICACIONES Una vez generado el conjunto de datos, a continuación abordamos su estudio como un MRBN y como un MRP. Notamos que el lenguaje R permite un estudio del MRBN desde los dos enfoques planteados en el tercer capítulo. Comenzamos analizando los datos como un modelo de regresión binomial negativa enfocado como una Poisson compuesta, usando el método de máxima verosimilitud. A continuación establecemos las instrucciones y las salidas obtenidas: En esta salida se obtiene para cada βisu estimación puntual, el intervalo de confianza correspondiente a un nivel de confianza de un 95%, además del error estandar de estimación. Dado a que generamos los datos para ser modelados usando la regresión binomial negativa, no nos sorprende el buen ajuste del modelo. Notamos que mediante la función ”ml.nb2”, R realiza una estimación de los parámetros del MRBN, tanto de los parámetros βcomo del parámetro de dispersión α. Hay que hacer constar que en esta salida el valor 1,929 no es la estimación de αsino de su inversa, siendo su valor 0,5181. A continuación vamos a analizar el mismo conjunto de datos pero enfocando el modelo de regresión binomial negativa (MRBN) como un miembro de la familia de modelos lineales generalizados. Recogemos las siguientes instrucciones con sus salidas correspondientes: 4.1. SIMULACIÓN 49 mediante la función ”glm.nb”se obtiene una salida similar a la anterior salvo por el parámetro de dispersión (definido por theta, θ, en este algoritmo), el cual R lo estima previamente mediante el método de máxima verosimilitud para su posterior inclusión en el algoritmo de modelos lineales generalizados (GLM). Estos resultados aparecen compañados de unos estadísticos adicionales como el estadístico desviación, el criterio AIC,... . El estadístico desviación, bajo las hipótesis del modelo correcto, sigue una distribución chi-cuadrado con 496 grados de libertad. Este estadístico tiene un valor de 528.3, que evaluando la siguiente relación para detectar la sobredispersión en los datos, como señalamos en el segundo capítulo D gl ⇒528,3 496 = 1,065 observamos como el modelo acomoda la sobredispersión de los datos, es decir, el modelo se adapta a la heterogeneidad adicional que podría existir en los datos. A continuación analizamos el conjunto de datos binomiales negativos a través de la regresión Poisson que, como caso particular del GLM, vienen dados a través de la función ”glm[..., family = ”poisson”, ..]”. Obtenemos: 50 CAPÍTULO 4. APLICACIONES Entre los resultados obtenidos aparecen las estimaciones de los coeficientes para el MRP. También obtenemos el estadístico desviación que, a diferencia del caso anterior, éste indica que los datos presentan una fuerte sobredispersión, pues en este caso: D gl ⇒38379 496 = 77,377 Los modelos de regresión Poisson sobredispersos algunas veces conducen a la confusión de que algunas variables explicativas contribuyen significativamente al modelo, cuando no esto no ocurre. Debemos tener en cuenta que cuando nos referimos de forma general a un conjunto de datos sobredispersos Poisson no tiene por qué ser sinónimo de datos binomiales negativos, pues existen distintos tipos de datos de conteo sobredispersos los cuales no pueden ser explicados mediante un modelo de regresión binomial negativo. 4.2. Aplicación a una base de datos particular. En el siguiente ejemplo abordaremos de una forma más detallada, a través de R, un modelo de ajuste que se adecue al siguiente caso: Estudio, realizado por los administradores escolares, del comportamiento de asistencia de los estudiantes de tercero de secundaria. Las variables explicativas, que usaremos para el estudio de la cantidad de días de ausencia, vienen dadas por el tipo de programa en el que está matriculado el estudiante y la nota alcanzada en un examen de matemáticas. Descripción de los datos Disponemos de 314 observaciones de estudiantes de secundaria en el archivo ”nb_data.dta”. La variable respuesta de interés es la cantidad de días de ausencia, daysabs, que viene dada en función de las variables math, que da la puntuación estandarizada en un examen de matemáticas realizado por cada estudiante, y la variable prog, que es una variable nominal de tres niveles que indica el tipo de programa de instrucción en el que está matriculado el estudiante, estos son el general, vocacional y académico. Instrucciones y resumen descriptivo de la base de datos: 4.2. APLICACIÓN A UNA BASE DE DATOS PARTICULAR. 51 Observamos como la varianza muestral (49.52) es bastante mayor a la media (5.955), lo que nos hace empezar a cuestionarnos sobre el problema de sobredispersión en los datos. Realizamos el siguiente gráfico que representa, dependiendo del tipo de programa en el que está matriculado el estudiante, el número de estudiantes que se ausentan una determinada cantidad de días. 52 CAPÍTULO 4. APLICACIONES La siguiente tabla muestra las medias y varianzas condicionadas de los dias de ausencia por cada tipo de programa. Se aprecia que el valor medio del número de días de ausencia, nuestra variable respuesta, parece variar según la variable prog. Y vemos que las varianzas dentro de cada nivel de prog son mayores a la media. Esta diferencia sugiere que el exceso de dispersión está presente y que un modelo binomial negativo sería apropiado para su análisis. Como hemos señalado en el capítulo anterior, podemos hacer uso del modelo de regresión binomial negativo para el estudio de datos de recuento que presentan sobredispersión, cosa que sospechamos que ocurra en estos datos debido a las claras diferencias media-varianza obtenidas. El análisis de regresión binomial negativo A continuación se utiliza la función glm.nb de la librería MASS para estimar los parámetros del modelo de regresión binomial negativa tradicional, es decir, los parámetros del modelo que obtendremos serán los del modelo binomial negativo NB2 que hemos estudiado. Mediante esta función R realiza una estimación previa del parámetro dispersión θmediante el método de máxima verosimilitud, para posteriormente incluirlo en el algoritmo GLM. 4.2. APLICACIÓN A UNA BASE DE DATOS PARTICULAR. 53 En esta salida se recoge la estimación de los coeficientes de regresión para cada una de las variables, junto a los errores estándar y p-valores. A través de los p-valores vemos que todos los βison no nulos, es decir, que todas las variables explicativas (math, progAcademic yprogVocational) influyen en la variable objetivo. A nivel interpretativo, tenemos que un aumento de la variable math produce una disminución del valor esperado del número de días de ausencia, de forma análoga ocurre para las variables dummy definidas (progAcademic yprogVocational), que tienen como referencia el tipo de instrucción categórica General. Nos muestra también la desviación nula calculada a partir del modelo con solo el coeficiente β0, la desviación del modelo completo, el resíduo desviación, y nos muestra el criterio AIC y la 2∗log verosimilitud. Hay que hacer constar que en esta salida el valor 1,0327 no es la estimación de θsino de su inversa, siendo su valor 0,968. Comprobación del modelo seleccionado Como hemos mencionado anteriormente, los modelos binomiales negativos consideran que la medias condicionadas no son iguales a las varianzas condicionadas. Esta desigualdad media-varianza es capturada por la estimación del parámetro de dispersión. Al contrario de la binomial negativa,el modelo de Poisson el cual si establece la igualdad, conocida como equidispersión. Luego, el modelo Poisson puede relacionarse con el modelo binomial negativo. Podemos usar una prueba de razón de verosimilitud para comparar estos dos modelos y probar el modelo seleccionado. El valor del chi-cuadrado asociado es de 2,157 ·10−203, con 5 grados de libertad. Esto sugiere fuertemente que el modelo binomial negativo, con la estimación del parámetro de dispersión,es más apropiado que el modelo de Poisson. Podemos obtener los intervalos de confianza para los coeficientes a través de la función de verosimilitud.