La evolución en el estudio de los modelos lineales. El Modelo Lineal Generalizado y los Modelos Mixtos
Abstract
Este proyecto es un trabajo teórico que presenta la evolución natural de ciertos problemas en el campo de la estadística. Supone una continuación del curso de Modelos Lineales y Diseño de Experimentos. No solo estudiaremos resultados y pruebas que no se cubrieron en la sección de Regresión Lineal, sino que también profundizaremos en el Modelo Lineal Generalizado como una extensión natural de la Regresión Lineal. Además, exploraremos a fondo otra forma de extender esta, los Modelos Mixtos. Ambos modelos pueden fusionarse para dar lugar a los Modelos Lineales Mixtos Generalizados. En cuanto a este último, solo discutiremos su funcionalidad, sin profundizar, ya que utiliza las mismas herramientas que los Modelos Mixtos y los Modelos Lineales Generalizados pero con expresiones más extensas y tediosas. Por tanto, no será de gran interés para el proyecto, ya que implicaría repetir cálculos y argumentos ya vistos en otros modelos. Finalmente, veremos cómo podemos aplicar esta potencia teórica en problemas cotidianos.
Full text
Universidad de Sevilla Facultad de Matemáticas Trabajo Fin de Estudios Grado en Matemáticas La evolución en el estudio de los modelos lineales. El Modelo Lineal Generalizado y los Modelos Mixtos. Raúl Vicente Muñoz Tutor Nombre: María José Zafra Garrido Dpto. Estadística e Investigación Operativa
Resumen Este proyecto es un trabajo teórico que presenta la evolución natural de ciertos problemas en el campo de la estadística. Supone una continuación del curso de Modelos Lineales y Diseño de Experimentos. No solo estudiaremos resultados y pruebas que no se cubrieron en la sección de Regresión Lineal, sino que también profundizaremos en el Modelo Lineal Generalizado como una extensión natural de la Regresión Lineal. Además, exploraremos a fondo otra forma de extender esta, los Modelos Mixtos. Ambos modelos pueden fusionarse para dar lugar a los Modelos Lineales Mixtos Generalizados. En cuanto a este último, solo discutiremos su funcionalidad, sin profundizar, ya que utiliza las mismas herramientas que los Modelos Mixtos y los Modelos Lineales Generalizados pero con expresiones más extensas y tediosas. Por tanto, no será de gran interés para el proyecto, ya que implicaría repetir cálculos y argumentos ya vistos en otros modelos. Finalmente, veremos cómo podemos aplicar esta potencia teórica en problemas cotidianos. iii
iv
Abstract In this project, we will observe the natural evolution of modeling certain problems in the eld of statistics. It can be viewed as a theoretical endeavor. It is assumed to be a continuation of the Linear Models and Experimental Design course. Not only will we examine results and proofs that were not covered in the Linear Regression section, but we will also delve into the Generalized Linear Model as a natural extension of linear regression. Furthermore, we will thoroughly explore another way to extend the Linear Regression, Mixed Models. Both models can be merged to give rise to Generalized Mixed Linear Models. Regarding the latter, we will only discuss its functionality without delving deeper, as it utilizes the same tools as Mixed Models and Generalized Linear Models but with longer expressions. Therefore, it would not be of great interest to the project, as it would involve repeating calculations and arguments already seen in other models. Finally, we will see how we can apply all this theoretical power to everyday problems. v
vi
Índice general Resumen III Abstract V Index VII Prefacio XI 1. Preliminares 1 1.1. Distribuciones y propiedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1. Deniciones y propiedades básicas de los vectores aleatorios . . . . . . . 1 1.1.2. Distribuciones fundamentales . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2. La inversa generalizada y formas cuadráticas . . . . . . . . . . . . . . . . . . . . 4 1.2.1. La inversa generalizada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2. Formascuadráticas.............................. 5 1.3. Otros conceptos de interés para los modelos . . . . . . . . . . . . . . . . . . . . 8 1.3.1. Propiedades del estimador de máxima verosimilitud . . . . . . . . . . . . 8 1.3.2. Criterios para la comparación de modelos . . . . . . . . . . . . . . . . . . 9 2. Regresión Lineal 11 2.1. Objetivo........................................ 11 2.2. Variables........................................ 12 2.2.1. Numéricas................................... 12 2.2.2. Categóricas .................................. 12 2.3. Hipótesisdelmodelo ................................. 12 2.4. Estimación de los coecientes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4.1. Ecuaciones normales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4.2. Propiedades de los EMC y una estimación de σ2 .............. 14 2.4.3. Tabla ANOVA y coeciente de determinación . . . . . . . . . . . . . . . 15 2.5. Regresión Lineal múltiple bajo normalidad . . . . . . . . . . . . . . . . . . . . . 16 2.5.1. Estimación de máxima verosimilitud . . . . . . . . . . . . . . . . . . . . 17 2.5.2. Contraste fundamental . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5.3. Contrastes individuales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 vii
viii Sección Índice general 2.5.4. Intervalo de conanza para la respuesta media . . . . . . . . . . . . . . . 20 2.5.5. Predicción para una nueva observación . . . . . . . . . . . . . . . . . . . 21 3. Modelo Lineal Generalizado 23 3.1. Objetivo........................................ 23 3.2. Breve introducción a la familia exponencial simple . . . . . . . . . . . . . . . . . 23 3.3. Hipótesisdelmodelo ................................. 24 3.4. Condiciones de regularidad y conceptos de interés para el modelo . . . . . . . . 26 3.5. Ecuaciones de máxima verosimilitud y la matriz de información de Fisher del modelo......................................... 28 3.5.1. Ecuaciones de máxima verosimilitud . . . . . . . . . . . . . . . . . . . . 28 3.5.2. La matriz de información de Fisher del modelo . . . . . . . . . . . . . . . 30 3.5.3. La simplicación de la función link canónica . . . . . . . . . . . . . . . . 30 3.6. Resolución de las ecuaciones de máxima verosimilitud . . . . . . . . . . . . . . . 31 3.6.1. El método de Newton Raphson . . . . . . . . . . . . . . . . . . . . . . . 31 3.6.2. El método de Fisher Scoring . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.6.3. El método de los Mínimos Cuadrados Ponderados Iterados (IRLS) . . . . 35 3.6.4. Valor inicial para los IRLS . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.6.5. Propiedades de ˆ β , el EMV de β ....................... 36 3.7. Inferencia en el modelo lineal generalizado . . . . . . . . . . . . . . . . . . . . . 38 3.7.1. El estadístico de Wald, un estadístico para la resolución de contrastes de hipótesislineales ............................... 39 3.7.2. El estadístico Score . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.7.3. Un estadístico para medir la bondad del ajuste: el estadístico Desviación 43 3.7.4. Contrastes anidados. La equivalencia asintótica de los estadísticos . . . . 47 3.7.5. Intervalo de conanza para el predictor lineal . . . . . . . . . . . . . . . 49 3.7.6. Residuos. Modicación de los estadísticos cuando ϕ es desconocido . . . . 50 4. Modelos Mixtos 55 4.1. Objetivo........................................ 55 4.2. Hipótesisdelmodelo ................................. 55 4.2.1. Algunas cuestiones sobre las hipótesis del modelo . . . . . . . . . . . . . 57 4.3. Estimaciones de máxima verosimilitud . . . . . . . . . . . . . . . . . . . . . . . 57 4.3.1. Ecuaciones del Modelo Mixto de Henderson . . . . . . . . . . . . . . . . 58 4.3.2. La estimación de máxima verosimilitud de σ2 ............... 59 4.3.3. Máxima verosimilitud restringida . . . . . . . . . . . . . . . . . . . . . . 62 4.3.4. Una propiedad interesante cuando V es desconocido . . . . . . . . . . . . 64 4.4. Inferencia en el vector paramétrico . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.4.1. Test de Wald para el vector paramétrico . . . . . . . . . . . . . . . . . . 65 4.4.2. Test de razón de verosimilitudes del vector paramétrico . . . . . . . . . . 67 4.5. Inferencia en los elementos de la varianza . . . . . . . . . . . . . . . . . . . . . . 67 4.5.1. Test de Wald para componentes de la varianza . . . . . . . . . . . . . . . 68
4.5.2. Test de razón de verosimilitudes para componentes de la varianza . . . . 68 4.6. Inferencia en los efectos aleatorios . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.6.1. Mejor predicción lineal insesgada (BLUP) . . . . . . . . . . . . . . . . . 70 4.6.2. Shrinkage ................................... 70 4.7. Residuos........................................ 71 4.8. Modelos Lineales Mixtos Generalizados . . . . . . . . . . . . . . . . . . . . . . . 72 Conclusiones 73 Bibliografía ......................................... 75 ix
2 Sección 1.1. Distribuciones y propiedades 1.1.2. Distribuciones fundamentales En este apartado se estudiarán las distribuciones más relevantes en los modelos que se abordarán, así como algunas propiedades. Denición 1.3. Una variable aleatoria Y se dice que sigue una normal de media µ y varianza σ2 si su función de densidad viene dada por fY(y) = 1 √2πσ2exp −1 2σ2(y−µ)2 con y∈R y se denota Y∼N(µ, σ2) . Su función generatriz de momentos asociada es MY(t) = exp tµ +1 2σ2t2 Denición 1.4. Dada Y∼N(µ, σ2) se llama tipicación de Y a Z=Y−µ σ vericando que Z∼N(0,1) Denición 1.5. Se dice que un vector aleatorio n -dimensional Y sigue una normal multivariante de parámetros µ y Σ , una matriz simétrica y semidenida positiva, si su función generatriz de momentos es MY(t) = exp t′µ+1 2t′Σt y se denota Y∼Nn(µ, Σ) Daremos algunas propiedades que se deducen directamente de que la función generatriz de momentos caracteriza a la distribución. Proposición 1.3. Toda transformación lineal de una normal es normal. Es decir si Y∼Nn(µ, Σ) y B es una matriz m×n=⇒ B Y∼Nm( B µ, B Σ B ′) Corolario 1.1. Dada Y∼Nn(µ, A ) , se verica que ∀i= 1, . . . , n , Yi∼N(µi, aii) donde aii denota el elemento (i, i) de A . Lema 1.1. Una matriz Σ es simétrica y semidenida positiva ⇐⇒ ∃A vericando que Σ=AA′ . Más aún Σ es invertible ⇐⇒ A lo es. Demostración: Por ser Σ simétrica, ∃ Q ortogonal vericando Σ= QDQ ′ con D matriz diagonal de los autovalores de Σ , que al ser semidenida postiva, son no negativos y por tanto D se puede descomponer como √ D √ D = D . Finalmente, tomando como A = Q √ D se tiene el resultado. Proposición 1.4. Se verica que Y∼Nn(µ, Σ)⇐⇒ Y=µ+ A Z siendo Z=hZ1. . . Zni′ con cada Zj∼N(0,1) independientes y Σ=AA′ . Corolario 1.2. Si Y∼Nn(µ, Σ) =⇒E[Y] = µ y V ar(Y) = Σ
Capítulo 1. Preliminares 3 Proposición 1.5. Para dos distribuciones normales, se dirá que son incorreladas ⇐⇒ son independientes (en general independencia implica incorrelación simplemente). Ahora veremos la función de densidad para distribuciones normales multivariantes. La prueba de esto se basa en hacer la demostración cuando Σ= I n que por independencia es obvio y para el caso general usar el teorema de cambio de variable. Proposición 1.6. La función de densidad de un vector aleatorio n -dimensional con distribución normal multivariante Y de media µ y matriz de covarianzas no singular Σ es fY(y) = 1 p(2π)ndet(Σ)exp((y−µ)′Σ−1(y−µ)) Denición 1.6. Se dene la distribución chi-cuadrado de n grados de libertad como la suma de n variables aleatorias independientes N(0,1) . Es decir si Y∼χ2 n , entonces Y=Pn i=1 Z2 i donde cada Zi∼N(0,1) Proposición 1.7. Su función generatriz de momentos es MY(t) = (1 −2t)−n/2 y la esperanza coincide con el número de grados de libertad. Proposición 1.8. La distribución chi-cuadrado es reproductiva, es decir, si X∼χ2 p e Y∼χ2 q independientes =⇒X±Y∼χ2 p±q Ahora daremos una propiedad que relaciona la normal multivariante con la distribución χ2 . Proposición 1.9. Sea Y un vector aleatorio n -dimensional con densidad, cuyas componentes no son necesariamente independientes. Entonces Y∼Nn(µ, Σ)⇐⇒ (Y−µ)′Σ−1(Y−µ)∼χ2 n Demostración: Se sabe que Y∼Nn(µ, Σ)⇐⇒ Y=µ+ A Z donde Z∼Nn(0, I n) y Σ= AA ′ . Despejando obtenemos que Z= A −1(Y−µ) . Por la denición de la distribución χ2 , Z′Z∼χ2 n , pero gracias a que ( A −1)′= ( A ′)−1 , Z′Z= (Y−µ)′( A −1)′ A −1(Y−µ) = (Y−µ)′( AA ′)−1(Y−µ) = (Y−µ)′Σ−1(Y−µ) y por tanto se tiene el resultado. Denición 1.7. Se dice que una variable T sigue una distribución t -student de n grados de libertad si dadas Z∼N(0,1) e Y∼χ2 n independientes entonces T=Z pY/n y se denota T∼tn
4 Sección 1.2. La inversa generalizada y formas cuadráticas Denición 1.8. Se dice que una variable F sigue una distribución F de Snedecor de n y m grados de libertad si dados Y1∼χ2 n y Y2∼χ2 m independientes, entonces F=Y1/n Y2/m y se denota F∼Fn,m 1.2. La inversa generalizada y formas cuadráticas 1.2.1. La inversa generalizada Para esta parte nos basamos en Morales et al. (2009). Primero se pondrá un poco en contexto sobre el motivo de la nueva herramienta. En un sistema A x=b , si A es una matriz invertible, entonces x= A −1b . La inversa generalizada se introduce en caso de que dicha inversa no exista. Denición 1.9. Dada una matriz A cuadrada de orden m si rg( A ) = r≤m , se dice que A − es una g -inversa de A si AA − A = A Se verá que A −b es solución del sistema A x=b , pero sin embargo, AA − no es la identidad, ni dicha g -inversa es única. Proposición 1.10. Si el sistema A x=b no es incompatible, entonces A −b es una solución del mismo donde A − es una g -inversa de A . Demostración: Si el sistema no es incompatible, ∃x0 con A x0=b , luego dada una g -inversa de A , A ( A −b) = AA − A x0= A x0=b de donde se deduce el resultado. También podemos expresar todas las soluciones del sistema en función de una inversa generalizada. Proposición 1.11. Si el sistema A x=b no es incompatible, todas las soluciones del mismo vienen dadas por ˜x= A −b+ ( A − A − I m)z donde z es un vector de dimensión adecuada y A − una inversa generalizada de A . Demostración: En primer lugar se verica que A ˜x= AA −b+ AA − A z− A z=b+ A z− A z=b
Capítulo 1. Preliminares 5 donde hemos usado el resultado anterior y la denición de A − . Por otro lado dado x0 solución del sistema, tomando z=−x0 , obtenemos que ˜x=x0 probando lo que se quería. Ahora damos un algoritmo para construir una inversa generalizada de A . Algoritmo 1: Construcción de la inversa generalizada. 1: Encontrar una submatriz de A cuadrada de rango igual al de A . Dicha matriz será W . 2: Invertir y trasponer W . 3: Regresar ( W −1)′ a A en las posiciones correspondientes. En el resto de elementos poner ceros. 4: Trasponer la matriz del paso anterior. Denición 1.10. Dada una matriz A se dene la inversa de Moore-Penrose como la única matriz A + que verica AA + A = A , A + AA += A + y tanto AA + como A + A son simétricas. Pruebas de la existencia, la unicidad y la propiedad que nos dice que si A es simétrica entonces su inversa de Moore-Penrose también se pueden ver en López Michelini (2019). Además de la propia denición se puede deducir el siguiente corolario: Corolario 1.3. Dada una matriz A y sea A + su inversa de Moore-Penrose, se verica que ( A +)+= A Esta sección nos ha servido para resolver el problema de que una matriz no tenga inversa. En ese caso se usa una inversa generalizada de dicha matriz, o también podemos usar la inversa de Moore-Penrose, que es única. En los modelos que se estudiarán, resuelve el problema de que no exista alguna de las inversas de ciertas matrices que aparezcan en dichos modelos. 1.2.2. Formas cuadráticas Denición 1.11. Sean A ∈Rn×n simétrica e Y∈Rn se dene la forma cuadrática Q=Y′ A Y=tr(Y′ A Y) Proposición 1.12. En las condiciones anteriores E[Q] = tr( A V ar(Y)) + E[Y]′ A E[Y] Demostración: Tomando µ=E[Y] , obtenemos que manipulando la expresión, E[Q] = E[Y′ A Y] = E[(Y−µ)′ A (Y−µ)] + E[µ′ A µ] + E[(Y−µ)′ A µ] + E[µ′ A (Y−µ)] =
6 Sección 1.2. La inversa generalizada y formas cuadráticas E[(Y−µ)′ A (Y−µ)] + µ′ A µ Ahora vamos a estudiar el primer término pues el segundo nos sirve para lo que queremos probar. Usando la denición de las formas cuadráticas, la linealidad de la esperanza y que tr( TS ) = tr( ST ) , E[(Y−µ)′ A (Y−µ)] = E[tr((Y−µ)′ A (Y−µ))] = tr( A E[(Y−µ)(Y−µ)′]) = tr( A V ar(Y)) que es justamente lo que quería probarse. Ahora demostraremos un lema que nos será útil para un resultado muy importante de formas cuadráticas. Lema 1.2. Sea A matriz simétrica e Y∼Nn(0, I n) entonces la forma cuadrática Q=Y′ A Y∼χ2 k⇐⇒ A es idempotente de rango k . Demostración: Si A es simétrica con A 2= A y rg( A ) = k . Gracias a que es idempotente, sus autovalores son o 0 o 1 , y por tanto, admite una descomposición a través de una cierta matriz ortogonal P de modo que P ′ AP =" I k 0 0 0 # donde las matrices nulas tienen las dimensiones necesarias para que tengan sentido. Ahora bien tomando Z= P Y , como Z es una transformación lineal de una normal multivariante, Z también es normal, además, E[Z] = E[ P Y] = P E[Y] = 0 y V ar(Z) = V ar( P Y) = P V ar(Y) P ′= PP ′= I n , pues P es ortogonal. Por tanto Z∼Nn(0, I n) , es decir, las marginales son independientes y Zi∼N(0,1) para cada i . Por último Y= P Z por lo que Y′ A Y=Z′ P ′ AP Z=Z′" I k 0 0 0 #Z= k X i=1 Z2 i concluyendo que Y′ A Y∼χ2 k . Recíprocamente, si Q=Y′ A Y∼χ2 k entonces MQ(t) = (1 −2t)−k/2 . Como A es simétrica ∃ P matriz ortogonal vericando que P ′ AP = D donde D es una matriz diagonal, es decir D =diag(d1, . . . , dn) . Para probar el resultado basta ver que hay k unos y el resto son ceros. Tomamos W= P ′Y , luego Y= P W por lo que Y′ A Y=W′ P ′ AP W=W′ D W= n X i=1 W2 idi (1.1) Además las variables Yi son independientes entre sí, pues en distribuciones normales, incorrelación equivale a independencia y en este caso V ar(Y) = I n . Luego las variables Wi
Capítulo 1. Preliminares 7 son normales, pues W es una transformación lineal de una normal multivariante y por un argumento completamente análogo al realizado en la otra implicación W∼Nn(0, I n) por lo que las Wi∼N(0,1) , independientes entre sí. Gracias a la independencia y que W2 i∼χ2 1 , MQ(t) = E[exp(Y′ A Y)] =(1,1) E"exp t n X i=1 W2 idi!#= n Y i=1 E[exp(tdiW2 i)] = n Y i=1 MW2 i(tdi) esto último para cada i es (1 −2tdi)−1/2 y nalmente uniendo esta última línea y la hipótesis inicial, Qn i=1(1 −2tdi)−1/2= (1 −2t)−k/2 por lo que hay k de los di que son unos y el resto deben ser ceros. En conclusión A es una matriz idempotente de rango k como se quería probar. Proposición 1.13. Sea A matriz simétrica e Y∼Nn(0,Σ) con Σ invertible, entonces la forma cuadrática Q=Y′ A Y∼χ2 k⇐⇒ A Σ es idempotente de rango k . Demostración: Gracias a la Proposición 1.4, Y= B Z con Z∼Nn(0, I n) y Σ= BB ′ , donde B es invertible. Usando el lema 1.2, gracias a que B ′ AB es simétrica por serlo A , Y′ A Y∼χ2 k⇐⇒ Z′( B ′ AB )Z∼χ2 k⇐⇒ B ′ AB es idempotente de rango k⇐⇒ B ′ AB = B ′ ABB ′ AB = B ′ A Σ AB y como tanto B como Σ son invertibles, lo anterior es equivalente a A = A Σ A ⇐⇒ A Σ= A Σ A Σ . Es decir, siguiendo la cadena de equivalencias, Y′ A Y∼χ2 k⇐⇒ ( A Σ)2= A Σ , es decir A Σ es idempotente. Finalmente veamos su rango. Gracias a que los autovalores de las matrices idempotentes son 0 o 1 y como el rango de una matriz diagonalizable es igual al rango de la matriz diagonal asociada, el rango de una matriz idempotente será igual al número de unos que tenga como autovalor, es decir su rango coincide con la traza. Por tanto, gracias a la conmutatividad de la traza, rg( A Σ) = tr( A Σ) = tr( ABB ′) = tr( B ′ AB ) = rg( B ′ AB ) de donde se deduce el resultado. Observación 1.1. Si Y∼Nn(0,Σ) , gracias al resultado anterior se puede probar parte de la proposición 1.9 pues entonces Y−µ∼Nn(0,Σ) y por tanto (Y−µ)′Σ−1(Y−µ)∼χ2 n⇐⇒ Σ−1Σ= I n es idempotente de rango n lo cual es trivial. Observación 1.2. Este resultado es una generalización del caso unidimensional, pues si Y∼N(µ, σ2) entonces (y−µ)2 σ2∼χ2 1 Vamos a ver otra consecuencia directa de estos resultados de formas cuadráticas aplicada a la inversa generalizada.
8 Sección 1.3. Otros conceptos de interés para los modelos Corolario 1.4. Si Y∼Nn(µ, Σ) con rg(Σ) = q < n y sea Σ+ su inversa de Moore-Penrose =⇒(Y−µ)′Σ+(Y−µ)∼χ2 q Demostración: En primer lugar, sabemos que como Σ es simétrica, Σ+ lo es. Además, Y∼Nn(µ, Σ)⇐⇒ Y=µ+ A Z donde Z∼Nn(0, I n) y Σ= AA ′ . Gracias al lema 1.2, como Z∼Nn(0, I n) y además A Z=Y−µ∼Nn(0,Σ) , entonces como Σ+ es simétrica, (Y−µ)′Σ+(Y−µ)∼χ2 q⇐⇒ Z′( A ′Σ+ A )Z∼χ2 q⇐⇒ A ′Σ+ A es idempotente de rango q , pero por ser Σ+ la inversa de Moore-Penrose, A ′Σ+ AA ′Σ+ A = A ′Σ+ΣΣ+ A = A ′Σ+ A y Σ+Σ también es idempotente con rg( A ′Σ+ A ) = tr( A ′Σ+ A ) = tr(Σ+ AA ′) = tr(Σ+Σ) = rg(Σ+Σ) = rg(Σ) = q La penúltima igualdad se debe a que por un lado, se tiene que rg(Σ+Σ)≤rg(Σ) y por otro, rg(Σ) = rg(ΣΣ+Σ)≤rg(Σ+Σ) . Por tanto, rg( A ′Σ+ A ) = q y el resultado queda probado. 1.3. Otros conceptos de interés para los modelos 1.3.1. Propiedades del estimador de máxima verosimilitud En esta parte se dará un resultado que será de gran utilidad para la teoría que se desarrollará. Apoyándonos en Cox and Hinkley (1979), tenemos el siguiente teorema. Teorema 1.1. Bajo ciertas condiciones de regularidad generales, sea ˆ θ estimador de máxima verosimilitud p -dimensional de θ vericando que es único y Σ , su matriz de covarianzas, entonces, ˆ θ es consistente de θ (al menos débilmente) y si Σ es no singular, ˆ θ∼aNp(θ, Σ)o(ˆ θ−θ)′Σ−1(ˆ θ−θ)∼aχ2 p Observación 1.3. En el caso en que Σ sea singular ( rg(Σ) = q < p ), podemos obtener Σ− una inversa generalizada de Σ y entonces (ˆ θ−θ)′Σ−(ˆ θ−θ)∼aχ2 q , ver Dobson and Barnett (2018). Otra opción es aplicar el corolario 1.4.
Capítulo 1. Preliminares 9 1.3.2. Criterios para la comparación de modelos Basándonos en Posada and Rosero Noguera (2007), se verán dos criterios usados en la comparación de modelos estadísticos. Denición 1.12. Se dene el coeciente de información de Akaike (AIC) como AIC = 2K−2 log Lik donde K es el número de parámetros inuyentes en el modelo y log Lik hace referencia al logaritmo de la función de máxima verosimilitud. Este criterio tiene en cuenta los cambios en la bondad del ajuste y las diferencias en el número de parámetros entre dos modelos. Los mejores modelos son aquellos que presenten el menor valor del AIC . Puede ser utilizado en todos los modelos que se van a presentar. Sin embargo, cuando se utiliza teoría Bayesiana (en nuestro caso, para los efectos aleatorios, introducidos en modelos mixtos), se suele usar el criterio de información bayesiano. Denición 1.13. Se dene el coeciente de información bayesiano (BIC) como BIC = log(n)K−2 log Lik donde n es el número total de datos, K es el número de parámetros inuyentes en el modelo y log Lik hace referencia al logaritmo de la función de máxima verosimilitud. Nuevamente, el modelo con el más bajo valor de BIC es considerado el mejor en explicar los datos con el mínimo número de parámetros. Se pueden implementar fácilmente algoritmos (de paso atrás o paso adelante entre otros) que clasiquen los posibles modelos en función de cualquiera de estos coecientes con el objetivo de encontrar el modelo con el valor más bajo de AIC o BIC .
10 Sección 1.3. Otros conceptos de interés para los modelos
Capítulo 2 Regresión Lineal Para este capítulo nos basaremos en los apuntes de la asignatura de Modelos Lineales y Diseño de Experimentos del curso 2022 - 2023 . 2.1. Objetivo La Regresión Lineal se plantea para estudiar una variable continua Y en función de ciertas p variables explicativas X1, . . . , Xp de forma que el modelo establece la siguiente relación Y=β0+β1X1+. . . +βpXp+ε donde en este caso los βj son escalares ∀j= 0, . . . , p y ε es el error producido en el modelo del que se supone E[ε] = 0 . En este capítulo se estudiará cómo podemos elegir dichos coecientes de forma que se minimice el error producido. Una vez hecho esto, podremos estudiar cómo es la bondad de nuestro ajuste, es decir, si nuestros datos son buenos para utilizar dicho modelo como aproximación. Otro de nuestros objetivos será inferir sobre aproximaciones de nuevas observaciones en función de las variables explicativas que dispongamos así como resolver algunos contrastes sobre parámetros del modelo y dar intervalos de conanza para los mismos. 11
18 Sección 2.5. Regresión Lineal múltiple bajo normalidad pues X ′ H = X ′ . Esto prueba la independencia de ˆ β y e concluyendo el resultado. Gracias a este teorema podemos obtener fácilmente un intervalo de conanza a nivel 1−α para σ2 . Este es el siguiente IC(1 −α, σ2) = SCE χ2 n−p−1,1−α/2 ,SCE χ2 n−p−1,α/2! 2.5.2. Contraste fundamental El contraste fundamental es el siguiente H0:β1=. . . =βp= 0 H1:βj= 0 para algún j Resolveremos dicho contraste con el test de razón de verosimilitudes (TRV). Para ello es necesario en primer lugar obtener los EMV bajo H0 , en cuyo caso yi∼N(β0, σ2) para cada i= 1, . . . , n , por tanto, ahora L(y;β0, σ2) = 1 (2π)n/2(σ2)n/2exp −1 2σ2Pn i=1(yi−β0)2 Tomando logaritmo de la verosimilitud de nuevo, se obtiene que ln L(y;β0, σ2) = −n 2ln 2π−n 2ln σ2−1 2σ2(yi−β0)2 Para obtener ˜ ˜ β0 , el EMV bajo H0 , basta con derivar respecto a β0 e igualar a 0 la expresión anterior, es decir ∂ln L(y;β0, σ2) ∂β0 = 0 ⇐⇒ Pn i=1(yi−β0) = 0 ⇐⇒ ny =nβ0 de donde se deduce que ˜ ˜ β0=y De forma análoga se hace la derivada del logaritmo de la verosimilitud respecto a σ2 , ∂ln L(y;˜ ˜ β0, σ2) ∂σ2= 0 ⇐⇒ − n 2σ2+1 2σ4Pn i=1(yi−y)2= 0 y despejando igual que en el subapartado anterior se obtiene que ˜ ˜σ2=SCT n Ahora podemos calcular el coeciente de razón de verosimilitudes y emplear un estadístico para resolver el contraste λ= Ly1, . . . , yn;y, SCT n Ly1, . . . , yn;ˆ β, SCE n= 1 (2π)n/2SCT nn/2exp −1 2SCT n SCT 1 (2π)n/2SCE nn/2exp −1 2SCE n SCE =SCE SCT n/2 = SCE SCR +SCE n/2 = 1 1 + SCR SCE n/2 Emplearemos primero resultados teóricos sobre poblaciones normales.
Capítulo 2. Regresión Lineal 19 Lema 2.1. Sea y1, . . . , yn m.a.s de una población con distribución N(µ, σ2) =⇒SCT σ2∼χ2 n−1 Proposición 2.2. Bajo H0 se verica que SCR σ2∼χ2 p independiente de SCE . Demostración: Podemos aplicar el lema anterior, pues bajo H0 , cada yi está idénticamente distribuida como una N(β0, σ2) y son independientes entre sí, por tanto SCT σ2∼χ2 n−1 , pero como SCE σ2∼χ2 n−p−1 , gracias al teorema 2.1, tras dividir por σ2 , SCT σ2=SCE σ2+SCR σ2 . En vista a lo anterior, necesariamente SCR es independiente de SCE , por ser una descomposición como suma de formas cuadráticas y SCR σ2∼χ2 p , probando el resultado. Usando lo anterior y el teorema 2.2 podemos construir el siguiente estadístico T=SCR/p SCE/n −p−1∼H0Fp,n−p−1 Así que mediante transformaciones de λ podemos llegar fácilmente al estadístico T y viceversa, por tanto la siguiente región crítica para el estadístico T nos resuelve el contraste fundamental, Rc :T > Fp,n−p−1,1−α que es una región crítica para un test de nivel α . Observación 2.5. Se puede resolver este contraste a través de otro estadístico. Como ˆ β∼Np+1(β, σ2 C ) y esto es equivalente a que (ˆ β−β)′ C −1(ˆ β−β) σ2∼χ2 p+1 Tomando ahora β= 0 como en la hipótesis nula y gracias al teorema 2.2, ˆ β′ C −1ˆ β σ2(p+ 1) SCE σ2(n−p−1) ∼H0Fp+1,n−p−1⇐⇒ ˜ T= ˆ β′ C −1ˆ β ˆσ2(p+ 1) ∼H0Fp+1,n−p−1 Resolvemos el contraste fundamental a través de una región crítica, Rc :˜ T > Fp+1,n−p−1,1−α 2.5.3. Contrastes individuales Este contraste se realiza sobre uno de los coecientes del modelo, es decir, dado j∈ {0,1, . . . , p} se plantea el siguiente contraste H0:βj= 0 H1:βj= 0 Para ello, recordemos que ˆ βj∼N(βj, σ2cjj) . Por propiedades de la normal, ˆ βj−βj pσ2cjj ∼N(0,1)
20 Sección 2.5. Regresión Lineal múltiple bajo normalidad y de nuevo gracias al teorema 2.2 se obtiene que ˆ βj−βj pσ2cjj rSCE σ2(n−p−1) ∼tn−p−1⇐⇒ ˆ βj−βj s.e.(ˆ βj) ∼tn−p−1 que es una cantidad pivotal para hallar un intervalo de conanza para βj a nivel 1−α , esto es, IC(1 −α, βj) = ˆ βj∓s.e.(ˆ βj)tn−p−1,1−α/2 lo cual nos resuelve el contraste, pues Tj=ˆ βj s.e.(ˆ βj) ∼H0tn−p−1 que nos permite obtener la siguiente región crítica a partir de la cual podemos hallar un p -valor, Rc :|Tj|> tn−p−1,1−α/2 El contraste fundamental sirve para saber si las covariables explican en cierta medida Y . Si se acepta H0 no quiere decir que βj= 0 ∀j , sino que aunque realmente dichos coecientes no sean nulos, no nos ayudan a explicar Y . Para el contraste individual, es un razonamiento análogo, si no se rechaza H0 no quiere decir que βj= 0 , sino que la variable Xj no es realmente signicativa en la regresión y quizá podamos eliminar Xj del modelo mediante métodos de selección de variables como por ejemplo uno basado en obtener el modelo de Regresión Lineal con menor AIC . 2.5.4. Intervalo de conanza para la respuesta media Dado un individuo cuyos valores son x′ 0=hx01 . . . x0pi , sabemos que E[Y|X=x0] = β0+β1x01 +. . . +βpx0p=z′ 0β=m(x0) . Consideramos ahora el estadístico ˆm(x0) = z′ 0ˆ β . Como ˆm(x0) es una transformación lineal de ˆ β y este último sigue una normal, entonces, ˆm(x0) también sigue una normal. Calculemos ahora su esperanza y varianza. E[ ˆm(x0)] = E[z′ 0ˆ β] = z′ 0β=m(x0) y V ar( ˆm(x0)) = V ar(z′ 0ˆ β) = z′ 0V ar(ˆ β)z0=σ2z′ 0 C z0 Por tanto, ˆm(x0)∼N(m(x0), σ2z′ 0 C z0) , por lo que apoyándonos nuevamente en el teorema 2.2 y haciendo cálculos análogos a los del apartado anterior, obtenemos que ˆm(x0)−m(x0) s.e.( ˆm(x0)) ∼tn−p−1 Finalmente, IC(1 −α, m(x0)) = ˆm(x0)∓s.e.( ˆm(x0))tn−p−1,1−α/2 es un intervalo de conanza para la respuesta media a nivel 1−α . Cuando se quiere dar una estimación puntual de la misma se tomará como es natural el punto medio de dicho intervalo, es decir ˆm(x0) , que tiene sentido pues estaríamos dando como aproximación cómo ajusta nuestro modelo los coecientes de la regresión.
Capítulo 2. Regresión Lineal 21 2.5.5. Predicción para una nueva observación Dado x′ 0 como en el apartado anterior y por la hipótesis distribucional, sabemos que y0=Y|X=x0 ∼N(m(x0), σ2) y procediendo igual que anteriormente obtenemos que ˆm(x0)∼N(m(x0), σ2z′ 0 C z0) . Gracias a propiedades de la diferencia de normales, tenemos que y0−ˆm(x0)∼N(0, σ2(1 + z′ 0 C z0)) y de igual forma, aplicando exactamente la misma idea que en el apartado anterior se obtiene una cantidad pivotal para y0 , esta es y0−ˆm(x0) s.e.(y0−ˆm(x0)) ∼tn−p−1 Todo esto, nos permite hallar un intervalo de predicción para una nueva observación a nivel 1−α , IP(1 −α, y0) = ˆm(x0)∓s.e.(y0−ˆm(x0))tn−p−1,1−α/2 . Este intervalo sirve para predecir cómo será la nueva observación en función de los valores que toman las covariables. Como anteriormente, para dar una estimación puntual se suele tomar el punto medio del intervalo que es ˆm(x0) . Observación 2.6. Si observamos los dos últimos intervalos obtenidos podemos ver que son bastante parecidos, de hecho, el último de ellos tiene amplitud mayor que el anterior. Esto tiene sentido, pues es más fácil predecir la media que tendrá y0 , que el propio valor que tomará.
22 Sección 2.5. Regresión Lineal múltiple bajo normalidad
Capítulo 3 Modelo Lineal Generalizado 3.1. Objetivo El objetivo del Mdelo Lineal Generalizado es: Dada una variable Y que sigue una cierta distribución que veremos a continuación, estudiaremos ciertas propiedades de dicha variable a partir de los valores de p variables explicativas X1, . . . , Xp . En muchas distribuciones se buscará predecir el valor esperado de la variable (como en el caso de la Bernoulli). En este capítulo se verá que en realidad, el Modelo Lineal Generalizado es una generalización de la Regresión Lineal. 3.2. Breve introducción a la familia exponencial simple Denición 3.1. Dada una variable aleatoria Y , se dice que esta pertenece a la familia exponencial si su función de densidad o de probabilidad es de la forma f(y;θ, ϕ) = exp yθ −b(θ) a(ϕ)+c(y;ϕ) donde θ es el denominado parámetro natural, ϕ su parámetro de escala y las funciones a , b y c son funciones especícas de la distribución. En particular, θ es una función de la esperanza y la varianza. Nota 3.1. Usualmente se toma a(ϕ) = ϕ ω donde w es un cierto peso asociado y ϕ > 0 . -Si los datos son agrupados, ω= 1 -Si la variable respuesta expresa promedio, ω=n -Si la variable es la suma de n respuestas individuales, ω=1 n 23
24 Sección 3.3. Hipótesis del modelo Proposición 3.1. Se verica que MY(t) = exp b(θ+a(ϕ)t)−b(θ) a(ϕ) y además E[Y] = b′(θ) V ar(Y) = a(ϕ)b′′(θ) Demostración: MY(t) = E[exp(tY )] = Zexp yθ −b(θ) a(ϕ)+c(y;ϕ) + tydy = exp b(θ+a(ϕ)t)−b(θ) a(ϕ)Zexp y(θ+a(ϕ)t)−b(θ+a(ϕ)t) a(ϕ)+c(y;ϕ)dy Ahora bien, la integral última es 1 en base a que estamos integrando la función de densidad. Esto concluye la prueba de la primera parte. Para la segunda usaremos las propiedades de la función generatriz, ya calculada, M′ Y(t) = MY(t)b′(θ+a(ϕ)t)⇒E[Y] = M′ Y(0) = MY(0)b′(θ) = b′(θ) y M′′ Y(t) = MY(t)(b′(θ+a(ϕ)t))2+MY(t)b′′(θ+a(ϕ)t)⇒E[Y2]=(b′(θ))2+a(ϕ)b′′(θ) Finalmente se tiene que V ar(Y) = E[Y2]−(E[Y])2=a(ϕ)b′′(θ) concluyendo la prueba. Una vez demostradas estas propiedades, estamos en condiciones de plantear las hipótesis del modelo. 3.3. Hipótesis del modelo Dadas n observaciones independientes, (yi;xi1, . . . , xip) con i= 1, . . . , n , se plantean las siguientes hipótesis: -Hipótesis distribucional: las variables yi=Y|X=xi siguen una distribución de la familia exponencial simple de media µi y parámetro de escala ϕ , que es común. Nota 3.2. Para pertenecer a dicha familia, es además necesario que esta cumpla una serie de condiciones de regularidad, que serán citadas al terminar las hipótesis del modelo. -Hipótesis estructural: E[yi] = µi=h(ηi) ó ηi=g(µi) con i= 1, . . . , n , donde ηi=x′ iβ y recordemos que x′ iβ=β1xi1+. . . +βpxip i= 1, . . . , p h es una función biyectiva conocida como función respuesta. g es función link, vínculo o enlace e inversa de h . β es vector paramétrico desconocido de dimensión p . xi es el vector de diseño de dimensión p , función de las covariables Xi .
Capítulo 3. Modelo Lineal Generalizado 25 Nota 3.3. Se ha introducido una pequeña modicación respecto al modelo de Regresión Lineal y es que en este caso no se considerará el término independiente β0 pues cuando la función h no es la identidad (que es lo que ocurre en el modelo de Regresión Lineal y por eso es una generalización), es más natural considerar los parámetros de dicha forma. Si quisiera considerarse dicho término, los resultados que se obtendrán a lo largo del capítulo serían análogos salvo un pequeño ajuste. En ocasiones es complicado encontrar la función h o g . Sin embargo, para cada elemento de la familia exponencial, existe una función vínculo canónica que viene dada por θi=θ(µi) = ηi . Recordemos que θ es el parámetro natural asociado y era función de la esperanza µi . Con esto tenemos que g(µi) = θ(µi) (3.1) es la conocida función link canónica. Como θi=ηi y sabemos que µi=b′(θi) = b′(ηi) , h(ηi) = b′(ηi) (3.2) es la función respuesta canónica. Ejemplo 3.1. Vamos a ver un ejemplo en el que comprobemos que una distribución sea de la familia exponencial y obtengamos la función link canónica. Tomaremos Y∼Po(λ) . Sabemos que la función de probabilidad de una Poisson es la siguiente P[Y=y] = exp(−λ)λy y! y aplicando logaritmos, ln(P[Y=y]) = −λ+yln(λ)−ln(y!) Para la familia exponencial, f(y;θ, ϕ) = exp yθ −b(θ) a(ϕ)+c(y;ϕ) Así que aplicando logaritmos nuevamente, se obtiene que ln(f(y;θ, ϕ)) = yθ −b(θ) a(ϕ)+c(y;ϕ) Ahora basta tomar θ= ln(λ) , b(θ) = exp(θ) , ϕ= 1 ( a=ϕ ) y c(y;ϕ) = −ln(y!) . Esto prueba que la Poisson pertenece a la familia exponencial simple. Obtengamos la función link canónica. Si Y|X=xi sigue una Poisson de media µi=λi . Por los cálculos anteriores se tiene que x′ iβ=θi= ln(λi) = ln(µi) . Luego g(µi) = ln(µi) = ηi es su función link canónica, cuya inversa es la función respuesta asociada de la cual obtenemos la hipótesis estructural del modelo, que en este caso es E[Y|X=xi] = µi=λi=h(ηi) = exp(ηi) Antes de pasar a la siguiente sección, se realizará una tabla en la que se verán ejemplos de
26 Sección 3.4. Condiciones de regularidad y conceptos de interés para el modelo distribuciones de la familia exponencial junto a sus funciones canónicas asociadas. Distribución b(θ)ϕ c(y;ϕ)θ(µ)µ(θ) N(µ, σ2)θ2 2σ2−1 2y2 ϕ+ ln(2πϕ)µ θ Be(π) ln(1 + exp(θ)) 1 0 ln π 1−πexp(θ) 1 + exp(θ) Po(λ) exp(θ) 1 −ln(y!) ln(λ) exp(θ) Γ(µ, ν)−ln(−θ)ν−1νln(νy)−ln(y)−ln(Γ(ν)) −1 µ−1 θ IG(µ, σ2)−(−2θ)1/2σ2−1 2ln(2πϕy3) + 1 ϕy−1 2µ2(−2θ)−1/2 Cuadro 3.1: Parámetros de distribuciones que siguen una familia exponencial Nota 3.4. En realidad las dos últimas columnas de la tabla anterior son la función link canónica y la función respuesta canónica respectivamente. Nótese que para el caso de la normal, las funciones canónicas son la identidad, lo cual tiene sentido, pues es la regresión lineal, estudiada anteriormente, se predecía cada yi como una combinación lineal de parámetros desconocidos en los que afectaba un cierto error de esperanza nula. 3.4. Condiciones de regularidad y conceptos de interés para el modelo Proposición 3.2. (Condiciones de regularidad de Cramer-Rao). Sea Y v.a. con función de densidad f(y;θ) (para el caso en que Y no sea continua las condiciones serían similares) con θ∈Θ e y1, . . . , yn una muestra independiente de Y . Las condiciones de regularidad son: -La variable Y tiene un campo de variación independiente de parámetros. - Θ ha de ser un intervalo abierto de R . - ∃∂f(y;θ) ∂θ salvo quizá en un conjunto de medida nula. -Se cumple que ∂ ∂θ R. . . Rf(y1, . . . , yn, θ)dy1. . . dyn=R. . . R∂ ∂θf(y1, . . . , yn, θ)dy1. . . dyn - ∃E"∂ ∂θ ln(f(y;θ))2#>0∀θ
Capítulo 3. Modelo Lineal Generalizado 27 Nota 3.5. La proposición anterior ha sido enunciada para el caso en que el espacio paramétrico es 1 -dimensional, sin embargo, esta se puede extender fácilmente a mayor dimensión. Notación: por simplicidad denotaremos por l(θ) = ln(f(y;θ, ϕ)) Denición 3.2. En estas familias se dene la función score como s(θ) = ∂l(θ) ∂θ Ahora se enunciará una propiedad que aunque se cumpla para las familias de distribuciones que veriquen las condiciones de regularidad, solo se demostrará para el caso particular de la familia exponencial simple, que es el que nos interesa. Proposición 3.3. Para una distribución que pertenezca a la familia exponencial simple se verica que E[s(θ)] = 0 y E[(s(θ))2] = −E[s′(θ)] Demostración: Recordemos que f(y;θ, ϕ). =f(y;θ) = exp yθ −b(θ) a(ϕ)+c(y;ϕ) y por tanto l(θ) = yθ −b(θ) a(ϕ)+c(y;ϕ) . Ahora derivamos respecto a θ obteniendo que s(θ) = y−b′(θ) a(ϕ)=1 a(ϕ)(y−E[Y]) y aplicando esperanza se obtiene la primera parte del resultado. Ahora bien, derivando de nuevo respecto a θ , s′(θ) = ∂ ∂θ y−b′(θ) a(ϕ)=−b′′(θ) a(ϕ) que al aplicarle la esperanza se queda igual y nalmente E[(s(θ))2] = 1 a(ϕ)2E[(y−E[Y])2] = 1 a(ϕ)2V ar(Y) = b′′(θ) a(ϕ) , ya que se sabe que V ar(Y) = a(ϕ)b′′(θ) . Esto prueba lo que se quería demostrar. Observación 3.1. El resultado anterior se ha enunciado en dimensión 1 , pues antes denimos la familia exponencial simple que es una variable aleatoria. Sin embargo, un cálculo completamente análogo nos dice que E[s(θ)s(θ)′] = −E∂2l(θ) ∂θ∂θ′ , donde en este caso s(θ) es un vector, ya que suponemos θ como un vector con varias componentes. La esperanza sigue siendo el vector 0 . La demostración se puede realizar de forma similar a la anterior viendo qué ocurre componente a componente. Denición 3.3. Se dene la matriz de información de Fisher esperada como F(θ) = V ar(s(θ)) Denición 3.4. Se dene la matriz de información de Fisher observada como Fobs(θ) = −∂2l(θ) ∂θ∂θ′
34 Sección 3.6. Resolución de las ecuaciones de máxima verosimilitud probando lo que se quería y concluyendo el teorema. Una vez hemos probado este resultado general, vamos a emplearlo en nuestro caso. Corolario 3.1. Una vez elegida una buena aproximación inicial β0 y bajo las condiciones de que H (βk) sea invertible y denida negativa, el método iterativo construido en (3.3) es convergente a ˆ β que es un máximo local de la función de máxima verosimilitud. Demostración: La convergencia está asegurada gracias al teorema anterior tomando como g(·) = s(·) , es decir ˆ β= l´ımk→∞ βk . Además como nos movemos en puntos próximos, gracias a la hipótesis de que la matriz de segundas derivadas es denida negativa, H (ˆ β) es denida negativa y como por construcción s(ˆ β) = 0 , ˆ β es un máximo local de la función de máxima verosimilitud, concluyendo el resultado. Observación 3.4. Las iteraciones que construye el método de Newton no dependen del término ai(ϕ) , pues observando el teorema 3.3 y la expresión (3.3), al efectuar F−1(βk)s(βk) , dicho término se cancela. Nota 3.7. Si la aproximación inicial es lo sucientemente buena o la función a maximizar cóncava, podemos conseguir que el punto al que converjamos sea, de hecho, máximo global de la función de máxima verosimilitud. Por tanto, hemos resuelto el problema que se quería y hemos obtenido a ˆ β como EMV de β . Posteriormente se verán propiedades de dicho estadístico y se realizarán intervalos de conanza y test de hipótesis gracias a dichas propiedades. Veamos un método que se obtiene a partir de este que acabamos de ver. 3.6.2. El método de Fisher Scoring Este método es un derivado del de Newton pues sustituye el Hessiano de l(·) por su valor esperado. Como bien sabemos por los cálculos realizados hasta ahora E[ H (·)] = −F(·) . La ventaja que tiene esta matriz es que por denición es la matriz de covarianzas del vector s(·) , por lo que es semidenida positiva y por tanto las iteraciones del método al realizar la sustitución nos irá maximizando l(·) sin problemas. En este caso la única hipótesis del método es que la matriz de información de Fisher sea invertible. Finalmente este resulta así, dado β0 , βk+1 =βk+F−1(βk)s(βk)∀k≥0 (3.5)
Capítulo 3. Modelo Lineal Generalizado 35 Ahora estudiaremos un método que en realidad no es más que una reescritura de este. 3.6.3. El método de los Mínimos Cuadrados Ponderados Iterados (IRLS) El método que se plantea puede interpretarse como una linealización del modelo. En la Regresión Lineal, se toma y= X β+ε . Veamos como se adapta esto al caso general. Por los teoremas 3.1 y 3.2 s(β) = n X i=1 yi−µi(β) V ar(yi) ∂µi ∂ηi xi= X ′ WD −1(y−µ(β)) F(β) = X ′ WX donde W y D son diagonales vericando que wii =1 V ar(yi)∂µi ∂ηi2 y dii =∂µi ∂ηi Ahora vamos a usar la expresión de (3.5) pero multiplicando a la izquierda por F(βk) , F(βk)βk+1 =F(βk)βk+s(βk) = X ′ W k X βk+ X ′ W k ( Dk)−1(y−µ(βk)) = X ′ W k[ X βk+ ( Dk)−1(y−µ(βk))] = X ′ W kzk (3.6) donde zk tiene los elementos zk i= n X j=1 xijβk j+ (yi−µi(βk))∂ηk i µk i =ηk i+ (yi−µi(βk))∂ηk i µk i (3.7) Ahora bien, dado β0 y despejando en (3.6) obtenemos la expresión de este método βk+1 = ( X ′ W k X )−1 X ′ W kzk∀k≥0 (3.8) En realidad el vector zk es en esta formulación una linealización (es la generalización de la que se habló en el comienzo del método) estimada de la función link g(·) evaluada en y , ya que por una aproximación de Taylor de orden 1 y gracias a (3.7), g(yi)≈g(µi(βk)) + (yi−µi(βk))g′(µi(βk)) = ηk i+ (yi−µi(βk))∂ηk i ∂µk i =zk i (3.9)
36 Sección 3.6. Resolución de las ecuaciones de máxima verosimilitud En la penúltima igualdad se ha usado que por la hipótesis estructural ηi=g(µi) . Obsérvese que en este proceso la variable z ha sido ajustada y aproximada por zk . Con el esquema del método, podemos obtener una nueva iteración βk+1 y a su vez una nueva aproximación zk+1 para el siguiente ciclo. Cuando el método converge lo hace a la mejor aproximación (o una muy buena), ya que se sabe que ˆ β es el EMV de β por lo que se ha visto en los apartados anteriores. Gracias a (3.8), ˆ β= ( X ′ˆ WX )−1 X ′ˆ W ˆz , con la respuesta ajustada estimada dada por (3.6), ˆz= X ˆ β+ ( ˆ D )−1(y−µ(ˆ β)) que es una aproximación buena de z cuando el modelo es bueno y/o está bien planteado y tiene sentido. Observación 3.5. En la ultima expresión se puede ver la linealización de la que hablamos. En el caso de la normal, la expresión queda tal cual como al principio. En otros casos, está ponderada. 3.6.4. Valor inicial para los IRLS Por como está construido, para obtener una aproximación inicial del método, basta dar una aproximación para los µi que es mucho más fácil de implementar en vista de la expresión (3.9) y con este último método podremos hacer iteraciones y ya tendríamos todo ganado para hallar ˆ β . Es decir, lo razonable sería empezar con la inversa de la función link de la media de la variable a predecir. Debemos asegurarnos de que los valores iniciales de partida están dentro del rango implícito de la función de la varianza especíca del Modelo Lineal Generalizado en el que estemos. En conclusión, tomamos µi(β0) = yi , pues en tal caso gracias a lo dicho anteriormente, z0 i=g(yi) y a raíz de aquí, por (3.8), podremos hallar β1 y seguir iterando. 3.6.5. Propiedades de ˆ β , el EMV de β Estudiaremos las propiedades del estimador cuando la función enlace empleada es la función link canónica. Después se comentará el caso general. Corolario 3.2. Usando como función link la canónica, los métodos iterativos estudiados anteriormente son, en realidad, el mismo. Demostración: Por construcción, el método de los IRLS es equivalente al método de Fisher. Ahora bien, gracias al teorema 3.3, cuando usamos la función link canónica, F≡Fobs , luego H =−Fobs =−F de donde se deduce lo que se quería probar.
Capítulo 3. Modelo Lineal Generalizado 37 Por esto, la única hipótesis necesaria para los métodos iterativos es que la matriz de información de Fisher sea invertible (pues H (·) es semidenida negativa). Ahora se probará un resultado sobre propiedades asintóticas para ˆ β construido con el método de Newton y que es EMV de β . Antes de enunciar el teorema, debemos asumir que se verican algunas propiedades de regularidad, ver Fahrmeir and Kaufmann (1985). Entre ellas, Pn i=1 xix′ i tiene rango máximo (ya supusimos rg( X ) = p ) y ∂µi ∂ηi= 0 ∀i= 1, . . . , n . Estas condiciones no son estrictas y, gracias a ellas, la matriz de información de Fisher es invertible para cada β del espacio paramétrico, propiedad que nos será de utilidad para probar la unicidad del estimador de máxima verosimilitud. Para hallar la esperanza y la varianza del estimador nos basaremos en Dobson and Barnett (2018). Teorema 3.5. Bajo las condiciones de regularidad del teorema 1.1, las que acabamos de enunciar y las hipótesis del teorema 3.4, usando como función link la canónica y dado ˆ β construido con el método de Newton, se verica que dicho estimador de β es consistente y además ˆ β∼aNp(β, F−1(ˆ β)) Demostración: Para probar la consistencia y que es asintóticamente normal, vamos a usar el teorema 1.1. Nos basta con ver la unicidad del EMV. Ahora bien, por construcción, ˆ β es EMV de β . Sin embargo, la función de máxima verosimilitud es estrictamente cóncava, ya que el logaritmo de la verosimilitud lo es (y dicha transformación es monótona creciente). Esto último se debe a que la matriz de segundas derivadas, H ≡ −F y F(·) es matriz de covarianzas y por tanto semidenida positiva. Además, gracias a las condiciones de regularidad enunciadas, es invertible. Por todo esto, la función l(·) es estrictamente cóncava y dichas funciones verican que si se alcanza un máximo, este es único. Entonces ˆ β es EMV de β que de hecho es único, por tanto ˆ β es consistente probando la primera parte. Para la segunda parte, nos basta con calcular la esperanza y la varianza del estadístico, pues gracias a la tesis del teorema 1.1, ˆ β es asintóticamente normal. Para un n lo sucientemente grande, F(ˆ β) = F(β) y será denotada por J . Esto se debe a que como ˆ β es consistente (al menos débilmente), ˆ β−→ β . Además gracias al uso de la función link canónica, H (ˆ β). =−F(ˆ β) . Por tanto, como H (·) es continua por ser matriz de segundas derivadas (es parte de las condiciones de regularidad), gracias a la propiedad de consistencia, como la convergencia en probabilidad o casi segura se conserva por transformaciones continuas, queda probada la expresión de J . Intuitivamente, esto tiene sentido, pues dado n lo sucientemente grande, podemos tomar ˆ β=β (produciendo un error poco signicativo) dentro de F , pues esta última es continua e invertible, por lo que si el error entre ˆ β y β es pequeño, también lo será al aplicarles F . Por
38 Sección 3.7. Inferencia en el modelo lineal generalizado otro lado, con un desarrollo de Taylor de segundo orden, como ˆ β es cercano a β l(β) = l(ˆ β)+(β−ˆ β)s(ˆ β) + 1 2(β−ˆ β)′ H (ˆ β)(β−ˆ β) = l(ˆ β)+(β−ˆ β)s(ˆ β)−1 2(β−ˆ β)′ J (β−ˆ β) Ahora podemos derivar la expresión anterior obteniendo que s(β) = s(ˆ β)− J (β−ˆ β) . Del hecho de que s(ˆ β) = 0 y J es invertible, llegamos a que ˆ β−β= J −1s(β) (3.10) Tomando esperanza como E[ J −1s(β)] = J −1E[s(β)] = 0 , se tiene que asintóticamente E[ˆ β] = β . Finalmente, gracias a la expresión de J y que esta es simétrica, V ar(ˆ β) = E[(ˆ β−β)(ˆ β−β)′] = E[ J −1s(β)s(β)′ J −1] = J −1E[s(β)s(β)′] J −1= J −1 Esto prueba que ˆ β∼aNp(β, F−1(ˆ β)) , concluyendo el teorema. Observación 3.6. Usando la función link canónica, como l(·) es cóncava, de hecho todo punto ˆ β vericando que s(ˆ β)=0 es máximo global. Las condiciones impuestas anteriormente nos aseguran que además l(·) es estrictamente cóncava y por ello ˆ β es único. Nota 3.8. Los resultados anteriores siguen siendo ciertos cuando la función link no es la canónica. Las pruebas de esto y unas diferentes de las anteriores se pueden consultar en Fahrmeir and Kaufmann (1985). 3.7. Inferencia en el modelo lineal generalizado Una vez resuelto el problema que presentaban las ecuaciones de máxima verosimilitud y obtenidas las propiedades del estadístico ˆ β , en esta sección estudiaremos algunos estadísticos de interés para el modelo, así como intervalos de conanza sobre los parámetros y la respuesta media en función de unos datos de las covariables. Además veremos un estadístico de bondad de ajuste y su uso para comparar modelos, los modelos anidados. Empezaremos viendo tres estadísticos que serán los protagonistas de la sección. Nota 3.9. En primer lugar vamos a realizar todos los cálculos cuando el parámetro de escala es conocido (por ejemplo en la Poisson o si lo conociéramos directamente). Después de ver los residuos del modelo, habrá un apartado en el que estudiaremos cómo se alteran estos estadísticos si ϕ es desconocido.
Capítulo 3. Modelo Lineal Generalizado 39 3.7.1. El estadístico de Wald, un estadístico para la resolución de contrastes de hipótesis lineales Como vimos en el teorema 3.5, bajo ciertas condiciones de regularidad, ˆ β∼aNp(β, F−1(ˆ β)) o lo que es lo mismo w = (ˆ β−β)′F(ˆ β)(ˆ β−β)∼aχ2 p . Este último es conocido como el estadístico de Wald. Su aplicación más útil es que permite resolver contrastes de hipótesis sobre los coecientes. Además, puede dar intervalos de conanza asintóticos para cada βj . Vamos a dar un primer uso al estadístico de Wald resolviendo contrastes individuales, esto es para j∈ {1, . . . , p} se plantea H0:βj= 0 H1:βj= 0 Aplicando propiedades de la normal, usando la expresión obtenida anteriormente y que V ar(ˆ βj) = vjj donde vjj denotará el elemento (j, j) de F−1(ˆ β) se tiene que ˆ βj∼aN(βj, vjj)⇐⇒ ˆ βj−βj √vjj ∼aN(0,1) que es una cantidad pivotal asintótica para βj que nos permite hallar el siguiente intervalo de conanza asintótico a nivel 1−α , IC(1 −α, βj) = ˆ βj∓√vjjZ1−α/2 Para resolver el contraste individual basta con ver si 0∈IC(1 −α, βj) , pues rechazamos si no está. En cambio si el 0 está, no tendríamos evidencia estadística para rechazar βj= 0 . El signicado de este contraste (muy similar al modelo de regresión lineal) es que si aceptamos H0 querría decir que la variable Xj no es muy signicativa en nuestro modelo, es decir no es signicativa para predecir µ y quizás podría construirse otro modelo en el que se prescinda de dicha variable, reduciendo la dimensión del vector paramétrico y quizá con un menor AIC . Observación 3.7. Podríamos haber tomado otra cantidad pivotal asíntótica para el intervalo anterior, ya que ˆ βj−βj √vjj ∼aN(0,1) ⇐⇒ (ˆ βj−βj)2 vjj ∼aχ2 1 . Sin embargo, dicha cantidad pivotal nos proporcionaría una región de conanza distinta, ya que estamos elevando al cuadrado y no tenemos por qué tener un intervalo. Se ha dado dicha cantidad pivotal asintótica para que se corresponda con la forma de resolver el contraste que viene a continuación. El estadístico de Wald no solo sirve para resolver contrastes individuales, sino que puede ser usado para múltiples contrastes, gracias a que la normalidad respeta las transformaciones lineales, de hecho el contraste individual se puede encuadrar en estos. La mayoría de las cuestiones que se plantean, pueden ser formuladas a través de una hipótesis lineal de la siguiente
40 Sección 3.7. Inferencia en el modelo lineal generalizado forma. H0: C β=d H1: C β=d (3.11) donde, sin pérdida de generalidad, C es una matriz q×p de rango máximo q≤p y d∈Rq . Si tuviera rango menor podríamos eliminar una la sin perder información y acabaríamos reduciéndonos al caso que se menciona. Para resolver este contraste daremos una región crítica usando un estadístico basado en el estadístico de Wald. Una vez resuelto, veremos un caso muy importante del contraste lineal. En primer lugar ˆ β es asintóticamente normal, entonces como C es una matriz de constantes, necesariamente C ˆ β debe ser asintóticamente normal, vericando que para un n sucientemente grande, E[ C ˆ β] = C E[ˆ β] = C β y V ar( C ˆ β) = C V ar(ˆ β) C ′= C F−1(ˆ β) C ′ , luego bajo H0 , la esperanza es d y de forma asintótica C ˆ β∼Nq(d, C F−1(ˆ β) C ′) o lo que es lo mismo w C ,d = ( C ˆ β−d)′( C F−1(ˆ β) C ′)−1( C ˆ β−d)∼H0χ2 q (3.12) pues rg( C F−1(ˆ β) C ′) = q y por tanto, ∃( C F−1(ˆ β) C ′)−1 . Esto último se debe a que por ser F−1(ˆ β) denida positiva y rg( C ) = q , podemos aplicar el siguiente resultado del álgebra lineal. Lema 3.1. Dada C , matriz q×p con rango máximo q≤p y A , simétrica y denida positiva de tamaño p×p=⇒ CAC ′ es simétrica y denida positiva, y por tanto, rg( CAC ′) = q Demostración: Dado w∈Rq no nulo, w′ CAC ′w= ( C ′w)′ A ( C ′w) . Ahora, como A es denida positiva, basta ver que C ′w es no nulo para que la expresión anterior sea siempre positiva. Por ello veamos que C ′ , vista como aplicación lineal es inyectiva, pero esto es cierto, pues el rango de C ′ coincide con el número columnas de C ′ y por tanto ker( C ′) = {0} . Esto prueba el resultado. La expresión dada en (3.12) es una cantidad pivotal asintótica que nos permite resolver el contraste lineal dando la siguiente región crítica, Rc : w C ,d > χ2 q,1−α . En este caso el estadístico de Wald puede ser visto como la distancia ponderada entre el estimador C ˆ β de C β y su valor determinado por la hipótesis nula, donde los pesos están dados por la inversa de la matriz de covarianzas del estimador. Nota 3.10. La expresión de (3.12) solo es cierta bajo H0 , en general se tiene que ( C ˆ β− C β)′( C F−1(ˆ β) C ′)−1( C ˆ β− C β)∼aχ2 q Observación 3.8. El contraste individual no es más que un caso particular de lo hecho
Capítulo 3. Modelo Lineal Generalizado 41 anteriormente, tomando C como un vector la de todo 0 salvo con un 1 en la componente j - ésima y siendo d= 0 un escalar. El estadístico asíntotico que se obtiene es el de la observación 3.7 que tiene sentido pues en este caso q= 1 . Debemos señalar un caso muy relevante de los contrastes lineales, el conocido contraste fundamental, este se plantea como H0:β1=. . . =βp= 0 H1:βj= 0 para algún j y es quivalente a un contraste lineal en el que q=p , C = I p y d es el vector de ceros de longitud p . Gracias a la expresión de (3.12) aplicada a este caso, se verica que w I p,0=ˆ β′F(ˆ β)ˆ β∼H0χ2 p de forma asintótica. Esto último nos resuelve el contraste fundamental con la siguiente región crítica, Rc : w I p,0> χ2 p,1−α . Este contraste sirve para ver si las covariables son signicativas para explicar la esperanza de la variable a predecir. Si aceptamos H0 no quiere decir que todas las covariables valgan 0 , sino que su impacto es pequeño y se pueden despreciar. Lo ideal sería que rechazáramos H0 , pues eso nos querría decir que nuestras covariables explican en cierta forma la esperanza de la variable a predecir, que es indicio de un buen modelo (como se puede comprobar es una interpretación análoga a la Regresión Lineal). 3.7.2. El estadístico Score Este estadístico se basa en el hecho de que la función score se anula en el estadístico de máxima verosimilitud, luego se usa en contrastes de hipótesis, pues la evaluación de la misma en el estimador obtenido bajo la hipótesis nula, será signicativamente diferente de 0 si dicha hipótesis no es cierta. Así se utiliza la distancia ponderada de s(˜ β) (estimador de máxima verosimilitud obtenido bajo H0 ) a 0 u =s(˜ β)′F−1(˜ β)s(˜ β) (3.13) Como E[s(˜ β)] = 0 y V ar(s(˜ β)) = F(˜ β) , nuevamente, Fahrmeir and Kaufmann (1985) prueban que bajo ciertas condiciones de regularidad u ∼aχ2 p . Este estadístico también nos sirve para resolver contrastes individuales a través de una región crítica y se demostrará que no solo para contrastes individuales, sino de forma general, cuando queremos contrastar que un subvector
42 Sección 3.7. Inferencia en el modelo lineal generalizado de β es nulo. Es decir, se plantea lo siguiente H0:βq= 0 H1:βq= 0 (3.14) donde como hemos dicho βq es un subvector de β . Para resolver el contraste mediante el estadístico Score, usaremos lo anterior. Tomamos A q , la submatriz de F−1(β) asociada al subvector de β . Esta siempre es invertible, ya que por ser F−1(β) denida positiva, gracias a que el contraste anterior es de la forma de 3.11, debido a que βq es una transformación lineal de β , podemos aplicar nuevamente el lema 3.1, pues C tendrá q vectores linealmente independientes de la matriz identidad q -dimensional y el resto el vector nulo. Por las propiedades de la normal multivariante, de la expresión de (3.13) deducimos que asintóticamente ˜ u =sq(˜ βq)′˜ A q(˜ β)sq(˜ βq)∼H0χ2 q (3.15) De modo que para resolver el contraste planteado nos basta con dar una región crítica. Esta será la siguiente, Rc :˜ u > χ2 q,1−α . Como dijimos anteriormente, este contraste se puede reescribir como un contraste lineal. Esto nos da pie al siguiente resultado. Proposición 3.5. Para la resolución del contraste planteado en (3.14), el estadístico de Wald y el estadístico Score son asíntóticamente equivalentes bajo H0 a una χ2 q . Demostración: En primer lugar podemos suponer sin pérdida de generalidad que el vector βq es el subvector de las primeras q componentes de β (si no lo fuera reordenando y cambiando índices se tiene) por tanto en realidad el contraste de hipótesis planteado es un contraste de la forma de (3.11) donde C = [ I q 0 ] , 0 es la matriz nula de dimensión q×(p−q) y d es el vector nulo. Gracias a (3.12), el contraste planteado en (3.14) se resuelve con un estadístico de Wald de la siguiente forma, w C , 0 =ˆ β′ q˜ A −1 q(ˆ β)ˆ βq∼H0χ2 q que es justo la misma distribución de la expresión de (3.15), una χ2 q . Esto prueba el resultado. Observación 3.9. La matriz ˜ A q no es más que C F−1(·) C ′ ( C la matriz de la demostración anterior). Por lo ya dicho, se observa que sq(˜ βq) = C s(˜ β) . Su esperanza es nula y V ar(sq(˜ βq)) = C E[s(˜ β)s(˜ β)′] C ′= C F(˜ β) C ′ . Debido a que en la expresión del estadistico que sigue una χ2 , la varianza aparece invertida y gracias a la forma de C , se tiene lo que se quería demostrar. Esta es otra forma de obtener la expresión (3.15).
Capítulo 3. Modelo Lineal Generalizado 43 3.7.3. Un estadístico para medir la bondad del ajuste: el estadístico Desviación El estadístico que se va a presentar a continuación se diferencia de los otros dos pues se puede emplear para medir la bondad del ajuste, es decir, nos permite ver si nuestro modelo es bueno y se ajusta a lo que queremos. Aún así, veremos que planteando un cierto contraste (los contrastes anidados), asintóticamente se puede relacionar con los otros dos estadísticos. Para sus propiedades y los contrastes anidados nos apoyaremos en Lindsey (2000). Antes de dar su expresión debemos ver qué es el modelo saturado. Denición 3.5. El modelo saturado tiene tantos parámetros como observaciones, es decir n parámetros linealmente independientes. Por lo tanto reproduce exactamente los datos pero sin simplicaciones, por lo que tiene poco uso para la interpretación, es un modelo sin patrones ni relaciones. El estadístico desviación se basa es una comparación entre el modelo y el modelo saturado. Su expresión es la siguiente D=−2 ln "L(ˆ β;y) L(ˆ βsat;y)# (3.16) Observación 3.10. Recordemos que L(·) es la función de máxima verosimilitud de las observaciones. Se cumple que ln(L(·)) = l(·) Antes de ver el contraste que nos sirve para ver si el ajuste es bueno o no, usaremos un lema que nos será de utilidad. Lema 3.2. Bajo las condiciones del teorema 1.1, dado ˆ β estimador de máxima verosimilitud de β (por tanto consistente) con longitud ja p , se verica que T=−2(ln(L(β;y)) −ln(L(ˆ β;y))) ∼aχ2 p (3.17) Demostración: Dado n grande y ˆ β único, ˆ β es cercano a β , por tanto, gracias a un desarrollo de Taylor de segundo orden donde aproximaremos H por su valor esperado, ver Lindsey (2000), y usaremos que s(ˆ β) = 0 , obtenemos que ln(L(β;y)) = ln(L(ˆ β;y))+(β−ˆ β)′s(ˆ β)+ 1 2(β−ˆ β)′ H (ˆ β)(β−ˆ β) = ln(L(ˆ β;y))−1 2(β−ˆ β)′ J (β−ˆ β)
50 Sección 3.7. Inferencia en el modelo lineal generalizado 3.7.6. Residuos. Modicación de los estadísticos cuando ϕ es desconocido Esta parte está basada en Mayoral and Socuéllamos (2001) y Fahrmeir et al. (2013). Los residuos se pueden utilizar para ver si el modelo se adecua respecto a la elección de la varianza, función link y los términos del predictor lineal. También se pueden utilizar para detectar observaciones inuyentes o valores anómalos que requieran mayor observación, pero de forma general, el residuo de cada dato, mide la discrepancia entre el valor observado y el pronosticado por el modelo. Alguna de las razones por las que un modelo puede no ser bueno son las siguientes: -Elección incorrecta de la distribución de probabilidad. -Especicación incorrecta de la forma en que la media cambia con las variables explicativas, ya sea porque h(·) o g(·) no sean apropiadas. -Variables explicativas no incorporadas al modelo. -Dependencia entre las observaciones, como por ejemplo, a lo largo del tiempo. A continuación se denirán dos tipos de residuos, los residuos de Pearson y los residuos de la devianza. Para ello, tomaremos ai(ϕ) = ϕ ωi como se vio en la Nota 3.1. Denición 3.6. Se dene el i -ésimo error de Pearson como rP i=√ωi yi−ˆµi qb′′(ˆ θi) Observación 3.14. Se puede considerar una generalización inmediata de los residuos para datos normales (función identidad como link canónico). Ahora daremos otra expresión de los residuos más estandarizada. Basándonos en la expresion de ˆ β de los IRLS, podemos construir una matriz de predicción generalizada (la matriz de predicción de la Regresión Lineal). Esta debe tener la forma ˆ H =ˆ W 1/2 X ( X ′ˆ WX )−1 X ′ˆ W 1/2 Puesto que la variabilidad de los residuos depende de las variables explicativas, podemos corregir su efecto dividiendo por p1−ˆ hii donde ˆ hii denota la componente (i, i) de la matriz de predicción. Ahora podemos denir lo siguiente: Denición 3.7. Se dene el i -ésimo residuo de Pearson estandarizado como rPS i=rP i p1−ˆ hii Ahora deniremos los residuos de la devianza, pero antes debemos saber que es la devianza no escalada. En primer lugar, debido a la expresión que tiene el estadístico desviación, está escalado, es decir depende del parámetro de escala. Sin embargo tiene para cada i al factor
Capítulo 3. Modelo Lineal Generalizado 51 ai(ϕ) dividiendo. Como ai(ϕ) = ϕ ωi , podemos sacar de dicha expresión el parámetro de escala dividiendo, obteniendo como resultado D=Du ϕ donde Du es la devianza no escalada. Antes de pasar a la denición, nótese que Du=Pn i=1 di (simplemente tomamos di como cada uno de los sumandos de la verosimilitud). Denición 3.8. Se dene el i -ésimo residuo de la devianza como di,r =signo(yi−ˆµi)√di vericando que Du=Pn i=1 d2 i,r Ahora veremos dos formas de estimar ϕ , pues necesitamos ver cómo calculamos los intervalos de conanza y cómo resolvemos los contrastes anteriores cuando ϕ es desconocido. Denición 3.9. Se dene el estadístico de Pearson como χ2= n X i=1 ωi(yi−ˆµi)2 b′′(ˆ θi) Se cumple que χ2=Pn i=1(rP i)2 Además, dividido por ϕ , se sabe que dicho estadístico se distribuye como una χ2 n−p para un n grande. La primera estimación que daremos de ϕ , será usando dicho estadístico, del que se puede probar consistencia. ˆ ϕχ2=χ2 n−p Observación 3.15. Intuitivamente, para un n muy grande ( n−p≈n ), como h(·) es continua, por ser ˆ β consistente, ˆµi≈µi y por tanto, bajo las condiciones de regularidad, ˆ θi≈θi , por lo que b′′(ˆ θi)≈b′′(θi) . Gracias a esto, E"ωi(yi−ˆµi)2 b′′(ˆ θi)#=ωi b′′(θi)E[(yi−µi)2] = ωi b′′(θi)V ar(yi) = ϕ de donde se tiene que E[ˆ ϕ] = ϕ , es decir, ˆ ϕ es asintóticamente insesgado que es una condición necesaria de que sea consistente. Ahora daremos otra forma de estimar ϕ , que aunque no tenga por qué ser consistente es muy útil y ya lo conocemos, este estadístico se basa en la devianza no escalada. ˆ ϕDu=Du n−p
52 Sección 3.7. Inferencia en el modelo lineal generalizado Observación 3.16. Aunque no podamos garantizar la consistencia, es asintóticamente insesgado cuando nuestro modelo es bueno, pues dado un n grande E[ˆ ϕDu] = E[Du] n−p=ϕE[D] n−p=ϕ basándonos en que D∼χ2 n−p cuando el modelo es bueno. Una vez visto todo esto, para terminar el capítulo veremos como quedan los estadísticos que empleamos en la sección anterior cuando ϕ es desconocido. Se supondrá la independencia entre los estadísticos que se relacionan para obtener distribuciones como la t -student o la F de Snedecor. Nota 3.14. Se llegará hasta las cantidades pivotales asintóticas, pues obtener intervalos de conanza o regiones críticas a partir de las mismas es análogo a como se hizo, pero cambiando la distribución (pues estas cambiarán). En primer lugar una vez sabemos que nuestro modelo se ajusta bien a los datos, dado un n grande, podemos asumir D∼χ2 n−p y podemos modicar el resto de estadísticos para que no dependan de ϕ . Nótese que igual que se hizo con D , tanto s(·) como F(·) tienen a ϕ en el denominador luego tomaremos ˜s(·) y ˜ F(·) que no dependen de ϕ . Para dar un intervalo de conanza asintótico para cada βj con el estadístico de Wald, como vjj era el elemento (j, j) de F−1(·) , tomamos vjj =ϕ˜vjj de modo que como ˆ βj−βj √vjj ∼aN(0,1) , gracias a la independencia, ˆ βj−βj pϕ˜vjj pDu/ϕ(n−p)⇐⇒ ˆ βj−βj qˆ ϕDu˜vjj ∼tn−p Obteniendo una cantidad pivotal asintótica que no depende de ϕ . Ahora pasaremos al contraste (3.11). Como en la expresión aparece la matriz de información de Fisher inversa y después se vuelve a invertir la expresión, podemos descomponer w C ,d =˜ w C ,d ϕ∼H0χ2 q y usando de nuevo la hipótesis de independencia y la distribución de D , dado un n grande, ˜ w C ,d ϕq Du ϕ(n−p) ∼H0Fq,n−p⇐⇒ ˜ w C ,d qˆ ϕDu ∼H0Fq,n−p cantidad pivotal asintótica que nos resuelve el contraste con ϕ desconocido. Haciendo un cálculo análogo podemos resolver (3.14) teniendo en cuenta que en ˜ u aparece ϕ
Capítulo 3. Modelo Lineal Generalizado 53 en el denominador y podemos tomar ˜ u =˜ ˜ u ϕ donde ˜ ˜ u no depende del parámetro de dispersión. Gracias a esto y los cálculos ya hechos podemos resolver el contraste con la siguiente cantidad pivotal asintótica ˜ ˜ u qˆ ϕDu ∼H0Fq,n−p Para los contrastes anidados de la forma de (3.20), como ocurre con D , podemos descomponer D(p, t) = Du(p, t) ϕ∼H0χ2 p−q de modo que volviendo a hacer lo mismo que hasta ahora, dado un n lo sucientemente grande, Du(p, t) ϕ(p−t) Du ϕ(n−p) ∼H0Fp−t,n−p⇐⇒ Du(p, t) (p−t)ˆ ϕDu ∼H0Fp−t,n−p Proporcionando una cantidad pivotal asintótica que nos permite resolver el contraste. Por último haremos lo mismo para el intervalo de conanza del predictor lineal. Como vimos en dicho apartado, tomamos p =ˆη0−η0 qx′ 0F−1(ˆ β)x0 ∼aN(0,1) y podemos obtener ϕ del denominador de la expresión anterior (teniendo en cuenta que está afectado por la raíz) de modo que p =˜ p √ϕ y actuando de forma análoga a anteriormente ˜ p √ϕ pDu/ϕ(n−p) ∼atn−p⇐⇒ ˜ p qˆ ϕDu ∼atn−p cantidad pivotal asintótica con la que volvemos a poder dar un intervalo de conanza asintótico para el predictor lineal, esta vez sin que dependa del parámetro de escala. Esto concluye el apartado. Observación 3.17. Antes de terminar el capítulo, nótese que los resultados asintóticos son exactos cuando la variable a predecir es una normal (visto en el capítulo anterior) y se usa como función enlace la canónica.
54 Sección 3.7. Inferencia en el modelo lineal generalizado
Capítulo 4 Modelos Mixtos 4.1. Objetivo El objetivo de los Modelo Mixtos es el estudio de una variable numérica en función de otras a través de unas observaciones que se van tomando a lo largo del tiempo, espacio, etc. Al igual que la Regresión Lineal, este modelo tiene efectos jos. Sin embargo, en este caso también debemos incluir efectos aleatorios que tendrán que ver con la correlación de los datos, pues al tomar varias observaciones de un mismo individuo, estas no tienen por qué ser independientes. Los efectos aleatorios incluidos modelarán la variabilidad de la variable a predecir. Nuestro objetivo será estimar los efectos jos, predecir los efectos aleatorios e inferir en algunas características del modelo. Los Modelo Mixtos se usan por ejemplo cuando disponemos de datos longitudinales que son aquellos en los que la respuesta se observa varias veces en diferentes tiempos para cada sujeto o unidad involucrada en el estudio. No se observan necesariamente el mismo número de veces (datos desbalanceados). 4.2. Hipótesis del modelo Verbeke et al. (1997) modela la situación que se ha descrito, se tomarán datos de m individuos y se hará un determinado número ni de veces con i= 1, . . . , m Las hipótesis del modelo son las siguientes: -Hipótesis estructural: E[yij] = x′ ijβ+z′ ijbi o yij =x′ ijβ+z′ ijbi+εij con i= 1, . . . , m y j= 1, . . . , ni donde β es un vector paramétrico p -dimensional de efectos jos y bi es un vector 55
56 Sección 4.2. Hipótesis del modelo de efectos aleatorios q -dimensional. Los vectores xij y zij son vectores de diseño. Hemos notado yij =Y|X=xij . Además podemos escribir la expresión anterior como yi= X iβ+ Z ibi+εii= 1, . . . , m (4.1) donde X i y Z i son matrices ni×p y ni×q respectivamente que presentan en la la j a x′ ij o z′ ij con j= 1, . . . , ni . Los vectores yi y εi son ni -dimensionales con yij y εij en la componente j . -Hipótesis distribucional: bi∼Nq(0, D ) y εi∼Nni(0,Σi) con i= 1, . . . , m -Hipótesis de incorrelación o independencia: bi y εi son independientes ∀i (recordemos que bajo normalidad independencia es equivalente a incorrelación). Además tiene sentido que si jamos un j , yij es independiente de yi′j siempre que i=i′ , pues medimos individuos distintos. En este modelo hay algunas simplicaciones para la matriz D , será desestructurada si Cov(bij, bi,j′) = σjj′ o de estructura diagonal si la expresión anterior es 0 . Lo mismo ocurre para Σi , tendrá una estructura de independencia si Cov(εij, εij′)=0 , una estructura de simetría compuesta si Cov(εij, εij′) = σε o desestructurada si Cov(εij, εij′) = σεjεj′ Por otro lado podemos reescribir el modelo una vez más, esta vez agrupando todo con la siguiente expresión matricial Y = X β+ Z b+ε (4.2) donde X es un vector por bloques que presenta en la i -ésima componente a X i , Z es una matriz diagonal por bloques con Z i en la componente (i, i) . De forma análoga se procede con b y ε , que son vectores de tamaño mq y N respectivamente, con N=Pm i=1 ni . La hipótesis distribucional también se puede reescribir como sigue, b∼Nmq(0, G ) y ε∼NN(0, R ) donde G y R son matrices diagonales por bloques vericando que en la diagonal i presentan a D y Σi con i= 1, . . . , m . Gracias a la hipótesis de independencia, V ar( Y ) = V ar( X β+ Z b+ε) = V ar( Z b) + V ar(ε) = ZGZ ′+ R . = V (4.3) De este modo V es también diagonal por bloques y V i= Z i DZ ′ i+Σi con i= 1 ...,m Observación 4.1. En (4.2), Y sigue siendo un vector aunque esté con notación matricial. Se ha elegido así para agilizar la notación. Observación 4.2. Z b modela la covarianza de Y y a través de ella, como covarían sus valores entre las unidades experimentales. Z es una matriz de diseño que denota la pertenencia de cada dato Y a un nivel de la variable modelada como efecto aleatorio.
Capítulo 4. Modelos Mixtos 57 Nota 4.1. En este capítulo volveremos a tomar β como un vector paramétrico p -dimensional. Se puede extender muy fácilmente esta teoría al caso (p+ 1) -dimensional. 4.2.1. Algunas cuestiones sobre las hipótesis del modelo Usando McCulloch and Searle (2004) vamos a ver el sentido de las hipótesis y cómo cuadran. Si Y = X β+ Z b+ε es obvio que E[ Y|b ] = X β+ Z b y V ar( Y |β) = V ar(ε) = R de donde se deduce que gracias a las propiedades de distribuciones condicionadas, proposición 1.2, E[ Y ] = E[E[ Y |b]] = X β+ Z E[b] = X β y V ar( Y ) = E[V ar( Y |b)] + V ar(E[ Y |b]) = R +V ar( X β+ Z b) = ZGZ ′+ R = V que es justamente lo que teníamos. Ahora bien, puede surgir la duda de por qué podemos considerar E[b] = 0 . Primeramente, b lo hemos introducido como un efecto aleatorio que nos modele la variabilidad de Y , luego tiene sentido que podamos asumir su esperanza nula al igual que ocurre con ε . Supongamos en cualquier caso que E[b] = τ y E[ Y |b] = X β+ Z b= X β+ Z τ+ Z (b−τ) Redenimos X ∗=h X Z i , β∗=hβ τi′ y b∗=b−τ obteniendo que E[ Y |b] = X ∗β∗+ Z b∗ donde ahora sí E[b∗] = 0 Antes de pasar a la siguiente sección, de ahora en adelante vamos a trabajar con tres posibles casos que pueden ocurrir en los Modelos Mixtos. El primero y más simple, que V ar(yi) = V i sea conocida. Después pasaremos a un segundo caso, que es el más parecido a una extensión de la Regresión Lineal que veremos y es que V ar(yi) = σ2 V i , con V i conocida. Por último, consideraremos el caso en que dicha covarianza es desconocida y veremos como podemos afrontar el problema. 4.3. Estimaciones de máxima verosimilitud En esta sección trabajaremos con el método de máxima verosimilitud. Aunque la máxima verosimilitud que solemos usar no es suciente para obtener estimadores para b debido a que Y ∼NN( X β, V ) , abordaremos el problema de otro modo. Cuando estamos en el tercer caso, además de necesitar métodos iterativos para aproximar la matriz de covarianzas, veremos otra forma de verosimilitud que se suele emplear, la verosimilitud restringida.
58 Sección 4.3. Estimaciones de máxima verosimilitud 4.3.1. Ecuaciones del Modelo Mixto de Henderson Fahrmeir et al. (2013) usan dichas ecuaciones para obtener los estimadores de β y b . En primer lugar, supondremos que tanto V como G y R son conocidas e invertibles (primer caso). Después veremos como abordar el caso en que no se conocen. Para llegar a ellas usaremos que Y |b∼NN( X β+ Z b, R ) y b∼Nmq(0, G ) . Además, gracias a que f( Y , b) = f( Y |b)f(b) , podremos obtener el término b que buscábamos. Nótese que como Y almacena la muestra de la variable a predecir, hallar la función conjunta es otra forma de abordar la máxima verosimilitud, de modo que L( Y ;β, b) = det( R )−1/2det( G )−1/2 (2π)N/2+mq/2exp −1 2(( Y − X β− Z b)′ R −1( Y − X β− Z b) + b′ G −1b) es la función que debemos maximizar. Como se suele hacer maximizaremos el logaritmo de la misma. Como la primera fracción no depende de lo que queremos estimar al tomar logaritmos lo podemos suprimir obteniendo que debemos maximizar ln(L( Y ;β, b)) . =l(β, b) = −1 2(( Y − X β− Z b)′ R −1( Y − X β− Z b) + b′ G −1b) . Para hallar el máximo derivamos e igualamos a 0 obteniendo las ecuaciones de Henderson ∂l(β, b) ∂β = 0 ⇐⇒ − X ′ R −1 X β+ X ′ R −1 Y − X ′ R −1 Z b= 0 ∂l(β, b) ∂b = 0 ⇐⇒ − Z ′ R −1 Z b+ Z ′ R −1 Y − Z ′ R −1 X β− G −1b= 0 En forma matricial " X ′ R −1 X X ′ R −1 Z Z ′ R −1 X Z ′ R −1 Z + G −1#"ˆ β ˆ b#=" X ′ R −1 Y Z ′ R −1 Y # que es equivalente a "ˆ β ˆ b#= ( C ′ R −1 C + B )−1 C ′ R −1 Y donde C =h X Z i y B ="0 0 0 G −1# El sistema anterior se resuelve de forma muy tediosa, aunque no es complicado usando un par de propiedades del álgebra de matrices, ver Correa Morales and Salazar Uribe (2016). Además, gracias a la expresión que se acaba de dar, es fácil ver que la matriz de segundas derivadas es denida negativa y por tanto la función l(β, b) cóncava, por tanto los estimadores que se obtendrán serán máximos globales. Finalmente, obtenemos que la expresión de los mismos es la siguiente:
Capítulo 4. Modelos Mixtos 59 ˆ β= ( X ′ V −1 X )−1 X ′ V −1 Y (4.4) ˆ b= GZ ′ V −1( Y − X ˆ β) (4.5) En el caso en que V y el resto de componentes de la varianza sean desconocidas, la solución es estimarlas y una vez hayamos obtenido dichas estimaciones, sustituir las componentes que aparezcan por estas. Se verá que no se pierden prácticamente propiedades si las aproximaciones son buenas (consistencia, que se vericará por ser EMV). Ahora bien de la expresión de ˆ β vamos a deducir algunas propiedades que nos serán de utilidad para los contrastes de los efectos jos. Proposición 4.1. Se verica que ˆ β∼Np(β, ( X ′ V −1 X )−1) Demostración: La normalidad está garantizada gracias a que es una transformación lineal de normales. Calculemos sus características. E[ˆ β] = E[( X ′ V −1 X )−1 X ′ V −1 Y ] = ( X ′ V −1 X )−1 X ′ V −1E[ Y ] = β V ar(β) = V ar(( X ′ V −1 X )−1 X ′ V −1 Y ) = ( X ′ V −1 X )−1 X ′ V −1V ar( Y ) V −1 X ( X ′ V −1 X )−1= ( X ′ V −1 X )−1 gracias a que V es simétrica por ser matriz de covarianzas. Esto concluye la prueba. Nótese que nuevamente gracias al lema 3.1, ∃( X ′ V −1 X )−1 , ya que V −1 es denida positiva y suponemos como en capítulos anteriores rg( X ) = p . 4.3.2. La estimación de máxima verosimilitud de σ2 Dado ahora el segundo caso, es decir V ar(yi) = σ2 V i , donde lo único desconocido del modelo es σ2 (el resto de componentes de la varianza también se conocen). Nótese que con unos cálculos análogos, la expresión de (4.4) y (4.5) no se ven afectadas. Ahora procedemos a calcular el estimador de σ2 por máxima verosimilitud y después veremos que se parecerá mucho al que obtuvimos en la Regresión Lineal. Para hallarlo vamos a usar la máxima verosimilitud
66 Sección 4.4. Inferencia en el vector paramétrico contraste individual en el que además podíamos dar intervalos de conanza sobre los coecientes se harían exactamente igual que en el capítulo pasado. Nos centraremos en resolver el más general de todos en cada uno de los casos que hemos venido trabajando. Sabemos que ˆ β∼Np(β, ( X ′ V −1 X )−1) . En el primer caso, es decir si V es conocida, directamente C ˆ β es normal k -dimensional por ser transformación lineal de una normal multivariante. Asumiendo H0 , como E[ C ˆ β] = C E[ˆ β] = C β=d y V ar( C ˆ β) = CV ar(ˆ β) C ′= C ( X ′ V −1 X )−1 C ′ , se tiene que C ˆ β∼H0Nk(d, C ( X ′ V −1 X )−1 C ′) o equivalentemente T= ( C ˆ β−d)′( C ( X ′ V −1 X )−1 C ′)−1( C ˆ β−d)∼H0χ2 k (4.11) debido a que podemos aplicar el lema 3.1 ( X ′ V −1 X es simétrica, denida positiva y rg( C ) = k ), la inversa tiene sentido. Dicho estadístico nos resuelve el contraste proporcionándonos la siguiente región crítica, Rc :T > χ2 k,1−α Ahora bien, si V ar(yi) = σ2 V i , como ˆ β∼Np(β, σ2( X ′ V −1 X )−1) , unos cálculos análogos nos llevan a que T σ2∼H0χ2 k donde T fue usado en (4.11). Utilizando ahora el teorema 4.1 que nos dice que SCE V σ2∼χ2 N−p independiente de ˆ β y se puede probar que por ello es además independiente de T , podemos hallar una cantidad pivotal, esta es F= T σ2k SCE V σ2(N−p) ∼H0Fk,N−p⇐⇒ ( C ˆ β−d)′( C ( X ′ V −1 X )−1 C ′)−1( C ˆ β−d) kˆ σ2 ∼H0Fk,N−p con la que hallamos la siguiente región crítica, Rc :F > Fk,N−p,1−α resolviendo el contraste. Ahora si V ar(yi) es desconocida, podremos resolver el contraste pero de forma asintótica. Nos basaremos en Verbeke et al. (1997) de aquí en adelante. Al no conocer V debemos estimarla ya sea con la máxima verosimilitud o la restringida, obteniendo ˆ θ estimador de θ , que nos da la estimación ˆ V = V (ˆ θ) . Si se da la consistencia de ˆ θ (en el caso de la máxima verosimilitud se dará gracias a las propiedades de este estimador), para un tamaño de muestra muy grande, entonces ˆ V = V (ˆ θ)≈ V (θ) = V por lo que V ar(ˆ β) = ( X ′ˆ V −1 X )−1 obteniendo gracias a los cálculos previos en (4.11), ˆ T= ( C ˆ β−d)′( C ( X ′ˆ V −1 X )−1 C ′)−1( C ˆ β−d)∼H0χ2 k una cantidad pivotal asintótica que nos resuelve el contraste gracias a la siguiente región crítica asintótica, Rc :ˆ T > χ2 k,1−α
Capítulo 4. Modelos Mixtos 67 Como notaron Dempster et al. (1981), el estadístico de Wald está basado en errores estándar estimados que subestiman la verdadera variabilidad de ˆ β debido a que no tienen en cuenta la variabilidad introducida al estimar θ . Este problema se puede resolver usando estadísticos t - y F - como ˆ F=( C ˆ β−d)′( C ( X ′ˆ V −1 X )−1 C ′)−1( C ˆ β−d) k en los que el número de grados de libertad del denominador tiene que ser estimado a partir de los datos (el del numerador es k ). Para más detalles sobre cómo estimar los grados de libertad del denominador, ver Verbeke et al. (1997). 4.4.2. Test de razón de verosimilitudes del vector paramétrico Este estadístico nos servirá para contrastar H0:β∈Θβ,0⊂Θβ , donde Θβ es el espacio paramétrico de β . Sea LML la función de verosimilitud de Y . El estadístico que se emplea es −2 ln λN=−2 ln "LML(ˆ βML,0) LML(ˆ βML)# (4.12) Donde ˆ βML,0 y ˆ βML son los máximos de LML en Θβ,0 y Θβ respectivamente. Se puede probar que −2 ln λN∼H0χ2 dim Θβ−dim Θβ,0 de forma asintótica. Este es un resultado clásico de la teoría de verosimilitud, véase Cox (1990). El resultado anterior no es del todo válido si usamos la verosimilitud restringida en lugar de la original, esto se debe a que como dijimos en secciones anteriores, la función a maximizar que obtenemos con la verosimilitud restringida no es una función de verosimilitud real. Sin embargo, a partir de un cierto n , sí sería válido, ya que es válido para la verosimilitud y esta acaba siendo equivalente a la restringida. Observación 4.5. Podríamos haber usado este resultado en el capítulo anterior para el estadístico desviación, pero no se hizo pues se tenían las herramientas sucientes para probarlo de forma más o menos sencilla. 4.5. Inferencia en los elementos de la varianza En muchas situaciones, la importancia de aproximar recae en la estructura de la esperanza en lugar de la varianza. Sin embargo, una modelización adecuada de la varianza puede ser útil para la interpretación de la variabilidad aleatoria en los datos y puede llegar a ser esencial para la estructura de la media. La sobreparametrización de la estructura de la varianza conduce a una
68 Sección 4.5. Inferencia en los elementos de la varianza estimación ineciente y posiblemente a una evaluación deciente de los errores estándar para las estimaciones de los perles de respuesta media. En esta parte, se discutirán los procedimientos inferenciales para componentes de la varianza en Modelos Mixtos. Lógicamente esto tiene sentido cuando la varianza no es conocida, ya que en los otros casos, al conocerse, no es necesario inferir en dichos elementos. 4.5.1. Test de Wald para componentes de la varianza Se sigue de la teoría clásica de verosimilitud, por ejemplo, Cox (1990), que bajo ciertas condiciones de regularidad, a partir de un cierto n , la dsitribución del estimador de máxima verosimilitud, y por equivalencia asintótica, el de máxima verosimilitud restringida, ˆ θ puede ser bien aproximado por una normal de media θ y matriz de covarianzas la inversa de la matriz de información de Fisher. Por ello, para estandarizar los errores aproximados en las estimaciones de las componentes de la varianza en θ , es factible calcularlos fácilmente mediante la inversión de la matriz de información de Fisher observada. Usando la normalidad asintótica de los parámetros estimados, podemos usar el test de Wald y obtener intervalos de conanza similares a los obtenidos anteriormente. Sin embargo, es necesario destacar que esto ocurre a nivel asintótico, es decir, necesitamos muestras muy grandes. Si el verdadero valor de θ es un valor de la frontera de Θθ , la aproximación falla. 4.5.2. Test de razón de verosimilitudes para componentes de la varianza Es muy similar al caso que vimos anteriormente para los efectos jos. Supongamos ahora que queremos contrastar H0:θ∈Θθ,0⊂Θθ y −2 ln λN denido como antes pero usando ˆ θML,0 y ˆ θML en lugar de ˆ βML,0 y ˆ βML respectivamente, que son los estadísticos de máxima verosimilitud en Θθ,0 y Θθ . Usando el mismo razonamiento, de forma asintótica, −2 ln λN∼H0χ2 dim Θθ−dim Θθ,0 Una de las condiciones de regularidad bajo las que la χ2 es válida es que el parámetro en H0 no esté en la frontera de Θθ . Por ello, este test sufre el mismo problema que el test de Wald. Gracias a la equivalencia asintótica, todo esto ocurre de igual forma si consideramos la máxima verosimilitud restringida y sus estimadores, pero el tamaño de la muestra debe ser lo sucientemente grande.
Capítulo 4. Modelos Mixtos 69 4.6. Inferencia en los efectos aleatorios Bien sabemos la utilidad de estimar β y las componentes de la varianza, sin embargo, también es útil calcular estimaciones de los bi , ya que reejan cuánto se desvían los perles especícos del sujeto del perl promedio general. Estimar estos efectos aleatorios nos será también de mucha ayuda para predicciones. Para esta sección nos apoyaremos en Verbeke et al. (1997). Como los efectos aleatorios son variables aleatorias, la forma más natural de estimarlos es con técnicas bayesianas,ver por ejemplo Liu et al. (1992). Como ya vimos, yi|bi∼Nni( X iβ+ Z ibi,Σi) y bi∼Nq(0, D ) . En teoría bayesiana, esta última es conocida como la distribución a priori del parámetro bi . Una vez los valores yi0 para yi han sido observados, podemos calcular la llamada distribución a posteriori de bi denida como la distribución bi|yi0 , basada en el teorema de Bayes, f(bi|yi0) = f(bi|yi=yi0) = f(yi0|bi)f(bi) Rf(yi0|bi)f(bi)dbi Usando la teoría de modelos lineales de Bayes, Lindley and Smith (1972), se puede probar que la expresión anterior es la densidad de una distribución normal multivariante. Con frecuencia, bi es estimado con la media de esa distribución, llamada la media a posteriori de bi . Está dada por ˆ bi=E[bi|yi] = DZ ′ i V −1 i(yi− X iˆ β) = DZ ′ i V −1 iei que tiene sentido pues es lo que obtuvimos estimando b en la expresión (4.5). En caso de que las componentes de la varianza no sean conocidas, basta con tomar ˆ D y ˆ V i . Además podemos calcular la varianza del estadístico anterior, pues V ar(e) = ( I N− H ∗ V −1) V ( I N− V −1 H ∗) , gracias a la demostración del teorema 4.1. Por tanto, V ar(ˆ b) = V ar( GZ ′ V −1e) = GZ ′ V −1V ar(e) V −1 ZG pues recordemos que G y V son simétricas por ser matrices de covarianzas. Una vez hallado ˆ b , podemos inferir en sus características de forma análoga a como hicimos en secciones anteriores, volviendo a tener en cuenta que si debemos aproximar la varianza, se subestima variabilidad introducida al estimar θ por ˆ θ . Por ello, en esos casos, la inferencia se basa en t -tests y F -tests aproximados en lugar de tests de Wald.
70 Sección 4.6. Inferencia en los efectos aleatorios 4.6.1. Mejor predicción lineal insesgada (BLUP) Suponemos ahora que queremos dar una estimación de la combinación u=x′β+z′bi donde x y z son vectores dados, p y q -dimensionales respectivamente. Lo más lógico sería tomar como estimador ˆu=x′ˆ β+z′ˆ bi Se puede probar que el estimador ˆu es la mejor predicción lineal insesgada, en el sentido de que es insesgado para u y que tiene la menor varianza de los insesgados que se escriben como combinación lineal, ver Harville (1976), por lo que si quisiéramos dar una predicción de u , una buena aproximación sería ˆu . Además ˆu=x′ˆ β+z′ˆ bi=x′ˆ β+z′ DZ ′ i V −1 i(yi− X iˆ β) = (x′−z′ DZ ′ i V −1 i X i)ˆ β+ (z′ DZ ′ i V −1 i)yi Y gracias a esta descomposición, podemos hallar fácil, pero tediosamente V ar(ˆu) , ya que los paréntesis se ven como constantes desde el punto de vista de la varianza y el problema se reduce a calcular la varianza de una combinación lineal de ˆ β e yi . Si las componentes de la varianza fueran desconocidas, bastaría con tomar sus aproximaciones. 4.6.2. Shrinkage Estos son estimadores conocidos como de contracción o encogimiento. Consideremos el problema de predicción de la evolución del sujeto i , ˆyi= X iˆ β+ Z iˆ bi= X iˆ β+ Z i DZ ′ i V −1 i(yi− X iˆ β) = ( I ni− Z i DZ ′ i V −1 i) X iˆ β+ Z i DZ ′ i V −1 iyi= Σi V −1 i X iˆ β+ ( I ni−Σi V −1 i)yi En la última igualdad se ha despejado Z i DZ ′ i de V = Z i DZ ′ i+Σi Hemos obtenido que ˆyi es un promedio ponderado del perl poblacional promedio X iˆ β y datos observados yi con pesos Σi V −1 i e I ni−Σi V −1 i respectivamente (que en caso de ser desconocidos se sustituye por sus estimadores). Observemos que X iˆ β recibe un peso mayor si la varianza residual es grande en comparación con la total. Este fenómeno es conocido como shrinkage (encogimiento, contracción), ya que los datos se contraen al perl promedio a priori X iˆ β .
Capítulo 4. Modelos Mixtos 71 4.7. Residuos Para esta parte nos basaremos en Correa Morales and Salazar Uribe (2016). Los residuos que denimos con anterioridad ( ei=yi− X iˆ β ), son correlados. Si queremos obtener unos que no lo sean, una solución la hallamos tras aplicar la factorización de Cholesky a la matriz V −1 i . Sea L i triangular superior vericando V −1 i= L i L ′ i podemos denir lo siguiente: Denición 4.3. Se dene el i -ésimo residuo de Cholesky como ˜ei= L ′ iei Estos residuos verican que su suma de cuadrados coincide con SCE V . En efecto, Pm i=1 ˜e′ i˜ei=Pm i=1 e′ i L i L ′ iei=Pm i=1 e′ i V −1 iei=SCE V Los residuos de Cholesky se usan para buscar outliers, correlación temporal, normalidad, etc. Si el modelo es correcto y la variabilidad de ˆ β despreciable, un cálculo sencillo nos dice que ˜ei∼Nni(0, I ni) (pues V ar(ei)≈V ar(yi) = V i ). Jacqmin-Gadda et al. (2007), usando resultados de simulación, sugieren que incluso para tamaños muestrales como n= 50 y 5 efectos jos, la variabilidad de ˆ β tenía un efecto despreciable. También existen otro tipo de residuos: Denición 4.4. Se dene el i -ésimo residuo condicionado como RCi=yi− X iˆ β− Z iˆ bi Los residuos condicionados pueden ser útiles para modelos homocedásticos, es decir, εij ∼N(0, σ2 ε) , por lo que para cada i , V i= Z i DZ ′ i+σ2 ε I ni . Usando esto y manipulando la expresión de los residuos condicionados, obtenemos que RCi=yi− X iˆ β− Z i DZ ′ i V −1 i(yi− X iˆ β) = σ2 ε V −1 i(yi− X iˆ β) Gracias a esto se demuestra fácilmente que si el modelo es correcto y V ar(ˆ β) despreciable, RCi∼Nni(0, σ4 ε V −1 i) . Antes de concluir esta sección daremos un coeciente para la bondad del ajuste que generaliza el R2 de la Regresión Lineal y por tanto tiene una interpretación análoga. Denición 4.5. Se dene el coeciente R2 en Modelos Mixtos como R2= 1 −Pm i=1 R′ CiRCi Pm i=1(yi−¯yi 1 ni)′(yi−¯yi 1 ni) , donde 1 ni es el vector ni -dimensional de unos e ¯yi es la media muestral asociada al individuo i , es decir, ¯yi=1 niPni j=1 yij
72 Sección 4.8. Modelos Lineales Mixtos Generalizados 4.8. Modelos Lineales Mixtos Generalizados En esta sección, basándonos en Fahrmeir et al. (2013), introduciremos un nuevo modelo que supone una mezcla de las características del Modelo Lineal Generalizado (visto en el capítulo anterior) y los Modelos Mixtos. Del primero utiliza que la variable a estudiar pertenece a la familia exponencial, así como la existencia de ciertas funciones link y respuesta que relacionan la esperanza con el predictor lineal. Dicho predictor lineal está sujeto a una muestra de la forma estudiada en este capítulo, es decir ηij =x′ ijβ+z′ ijbii= 1, . . . , m, j = 1, . . . , ni Las hipótesis del modelo son por tanto: -Hipótesis estructural: E[yij|bi] = µij =h(ηij) o ηij =g(µij)∀i= 1, . . . , m , j= 1, . . . , ni -Hipótesis distribucional: yij sigue una distribución de la familia exponencial al igual que en el Modelo Lineal Generalizado. Los efectos aleatorios, bi verican que son independientes e idénticamente distribuidos con bi∼Nq(0, D ) Como se dijo al comienzo del proyecto, el estudio de este modelo se omitirá, pues las herramientas teóricas fundamentales para el mismo se han visto ya en capítulos anteriores, es decir, estaríamos usando conceptos similares pero con operaciones mucho más tediosas, por tanto el modelo carece de interés teórico.
Conclusiones A lo largo del proyecto, tras introducir el modelo de regresión más simple, la Regresión Lineal, hemos estudiado en profundidad modelos más complejos y útiles como son el Modelo Lineal Generalizado y los modelos Mixtos. Finalmente, concluimos que fusionando estos dos modelos obtenemos el más general de todos, el Modelo Lineal Mixto Generalizado. Ahora veremos varias aplicaciones de toda la potencia teórica estudiada en casos de la vida cotidiana. En cuanto a la Regresión Lineal, es una primera opción muy buena si queremos estudiar datos de dependencia. Sin embargo, debido a la sencillez del modelo, no tiene un impacto directo muy grande en el día a día. Es decir la aplicación principal que tiene es que como modelo de inicio puede funcionar en algunas situaciones muy especicas aunque no es lo normal. Pasamos ahora al Modelo Lineal Generalizado que tiene mucho más impacto. Dicho modelo es mucho más completo a la hora de estudiar muestras independientes. El mayor impacto en el día a día lo tiene la Regresión Logística, que será un tipo de Modelo Lineal Generalizado cuando la variable a estudiar sea una Bernoulli. Tiene impacto en varios campos, aunque destaca en medicina, ya que este modelo sirve para detectar factores de riesgo de enfermedades gracias al conocido odds ratio (que está totalmente relacionado con el vector paramétrico). En este ámbito, también permite predecir la probabilidad que tiene un individuo con ciertas características de desarrollar una enfermedad (el padecer dicha enfermedad sería la variable a predecir). Gracias a la generalidad del modelo, también sirve para modelar todo tipo de situaciones en las que la variable que queremos estudiar se encuadre en la familia exponencial. Por otro lado, los Modelos Mixtos se pueden aplicar para estudiar variables normales en las que las muestras no tienen por qué ser independientes. Tiene impacto en la agricultura (si, por ejemplo, queremos estudiar el efecto de varios fertilizantes sobre un cultivo y tomar muestras a lo largo del tiempo para ver cuál nos interesa), biología (si queremos estudiar el efecto de un tratamiento para la producción de testosterona en ratas macho en su crecimiento, ver Verbeke et al. (1997)), etc. Por último, tenemos la fusión de los dos últimos modelos, los Modelos Lineales Mixtos Generalizados, que sirven para modelar casi cualquier situación en la que la variable a estudiar 73
presente aleatoriedad (basta con que la distribución de la variable a estudiar esté incluida en la familia exponencial). Por tanto, gracias a su generalidad, puede tener impacto en cualquier ámbito cotidiano y por ser la generalización de los modelos precedentes, puede modelar situaciones sencillas o situaciones muy complejas. Un buen resumen de todo se recoge en la siguiente gura. Regresión Lineal (RL) Modelo Lineal Generalizado (GLM) Modelo Mixto (MM) Modelo Lineal Mixto Generalizado (GLMM) Figura 4.1: Relaciones de contención entre modelos estadísticos. En conclusión, hemos estudiado cómo partiendo de un modelo tan simple cómo la Regresión Lineal, podemos llegar de forma natural a un modelo tan completo como el Modelo Lineal Mixto Generalizado. Esta es la evolución que han tenido los modelos lineales a lo largo del tiempo. 74
Bibliografía Agresti, A. (2015) Foundations of linear and generalized linear models . John Wiley & Sons. Correa Morales, J. C. and Salazar Uribe, J. C. (2016) Introducción a los modelos mixtos. Escuela de estadistica . Cox, C. (1990) Fieller's theorem, the likelihood and the delta method. Biometrics , pp. 709718. Cox, D. R. and Hinkley, D. V. (1979) Theoretical statistics . CRC Press. Cressie, N. and Lahiri, S. N. (1993) The asymptotic distribution of reml estimators. Journal of multivariate analysis , vol. 45(2): 217233. Demidenko, E. (2013) Mixed models: theory and applications with R . John Wiley & Sons. Dempster, A. P., Rubin, D. B. and Tsutakawa, R. K. (1981) Estimation in covariance components models. Journal of the American Statistical Association , vol. 76(374): 341353. Dobson, A. J. and Barnett, A. G. (2018) An introduction to generalized linear models . CRC press. Fahrmeir, L. and Kaufmann, H. (1985) Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models. The Annals of Statistics , vol. 13(1): 342368. Fahrmeir, L., Kneib, T., Lang, S., Marx, B., Fahrmeir, L., Kneib, T., Lang, S. and Marx, B. (2013) Regression models . Springer. Harville, D. (1976) Extension of the gauss-markov theorem to include the estimation of random eects. The Annals of Statistics , vol. 4(2): 384395. Harville, D. A. (1974) Bayesian inference for variance components using only error contrasts. Biometrika , vol. 61(2): 383385. Harville, D. A. (1977) Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American statistical association , vol. 72(358): 320338. Jacqmin-Gadda, H., Sibillot, S., Proust, C., Molina, J.-M. and Thiébaut, R. (2007) Robustness of the linear mixed model to misspecied error distribution. Computational Statistics & Data Analysis , vol. 51(10): 51425154. Laird, N. M. and Ware, J. H. (1982) Random-eects models for longitudinal data. Biometrics , pp. 963974. Lindley, D. V. and Smith, A. F. (1972) Bayes estimates for the linear model. Journal of the Royal Statistical Society Series B: Statistical Methodology , vol. 34(1): 118. Lindsey, J. K. (2000) Applying generalized linear models . Springer Science & Business Media. Liu, L.-M., Hudak, G. B., Box, G. E., Muller, M. E. and Tiao, G. C. (1992) Forecasting and time series analysis using the SCA statistical system , vol. 1. Scientic Computing Associates DeKalb, IL. 75