Modelos lineales de efectos mixtos. Librerías R

UNIVERSIDAD DE SEVILLA Modelos lineales de efectos mixtos. Librer´ ıas R. Trabajo Fin de Grado - Grado en Matem´ aticas Departamento de Estad´ ıstica e Investigaci´ on Operativa Facultad de Matem´ aticas Junio 2025 TRABAJO REALIZADO POR: CRISTINA RUIZ MART´ IN TUTOR: RAFAEL PINO MEJ´ IAS Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 2 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Resumen Este Trabajo de Fin de Grado aborda el estudio te´ orico y pr´ actico de los Modelos Lineales de Efectos Mixtos. Se describe detalladamente el modelo cl´ asico y el modelo extendido, formulados rigurosamente desde el punto de vista matem´ atico. Se incluyen tambi´ en las estimaciones de sus par´ ametros por distintos m´ etodos, y tambi´ en m´ etodos de diagn´ ostico e inferencia. En cuanto a la parte pr´ actica del trabajo, usaremos R, con librer´ ıas como lme4 ynlme, para poner en pr´ actica el modelo lineal de efectos mixtos aplicado a un estudio real sobre la agudeza visual en pacientes con degeneraci´ on macular. Es muy importante tener en cuenta las diferencias entre los individuos dentro de un mismo grupo y entre grupos diferentes, para obtener resultados fiables. 3 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 4 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Abstract This Final Degree Project addresses the theoretical and practical study of Linear MixedEffects Models. Both the classical and extended models, which are detailed and rigorously described, are formulated from a mathematical perspective. The work also includes parameter estimation through various methods, as well as model diagnostic and inference techniques. Concerning about the practical part of the project, the Rprogramming language is used, using libraries such as lme4 and nlme to apply the linear mixed-effects model to a real study on visual acuity in patients with age-related macular degeneration. It is vital to pay attention to the differences between individuals both within and between groups in order to obtain reliable results. 5 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 6 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. ´ Indice general Cap´ ıtulo 1: Introducci´ on 11 Cap´ ıtulo 2: Modelo Lineal de Efectos Mixtos 15 2.1.-El modelo cl´ asico lineal de efectos Mixtos .................... 15 2.1.1.-Especificaci´ on a un nivel de un factor de agrupaci´ on ........... 15 2.1.2.-Especificaci´ on para todos los datos .................... 17 2.2.-El modelo lineal extendido de efectos mixtos ................... 18 2.3.-Distribuciones definidas por las variables aleatorias yyb............ 18 2.3.1.-Distribuci´ on incondicionada de efectos aleatorios ............. 18 2.3.2.-Distribuci´ on condicionada de ydados los efectos aleatorios ....... 19 2.3.3.-Distribuciones adicionales definidas por yyb.............. 19 2.3.3.1.-Distribuci´ on conjunta de yyb................. 20 2.3.3.2.-Distribuci´ on marginal de y................... 20 2.3.3.3.-Distribuci´ on de bdado y.................... 20 2.4.-Estimaci´ on ..................................... 21 2.4.1.-El modelo marginal impl´ ıcito en el LMM cl´ asico ............. 21 2.4.2.-Estimaci´ on por m´ axima verosimilitud ................... 22 2.4.3.-M´ ınimos cuadrados penalizados ...................... 23 2.4.4.-Incertidumbre en la estimaci´ on de par´ ametros ............... 26 2.4.5.-Enfoques de estimaci´ on alternativos .................... 26 2.5.-Diagn´ ostico del modelo .............................. 27 2.5.1.-Normalidad de efectos aleatorios ..................... 27 2.5.2.-Diagn´ ostico residual ............................ 28 2.5.3.-Diagn´ ostico influencia ........................... 29 2.6.-Inferencia y selecci´ on de modelos ......................... 29 2.6.1.-Contraste de hip´ otesis sobre los efectos fijos ............... 29 2.6.2.-Contraste de hip´ otesis sobre los par´ ametros de varianza-covarianza . . . 30 2.6.3.-Intervalos de confianza ........................... 31 Cap´ ıtulo 3: Implementaci´ on con R 35 3.1.-Modelo con interceptos aleatorios y varianza residual homog´ enea ........ 36 3.1.1.-Especificaci´ on del modelo ......................... 36 3.1.2.-Sintaxis y resultados en R......................... 38 3.2.-Modelo con interceptos aleatorios y la funci´ on de varianza residual varPower(·) . 43 3.2.1.-Especificaci´ on del modelo ......................... 43 3.2.2.-Sintaxis y resultados en R......................... 44 3.2.3.-Gr´ aficos de diagn´ ostico .......................... 47 7 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.3.-Modelos con interceptos y pendientes aleatorias y funci´ on de varianza residual varPower(·) .................................... 53 3.3.1.-Modelo con una matriz general D..................... 53 3.3.2.-Modelo con una matriz diagonal D.................... 57 3.3.3.-Modelo con una matriz diagonal Dy un efecto de tratamiento constante . 61 3.4.-Una funci´ on alternativa de varianza residual: varIdent(·) ............. 64 3.5.-Prueba de hip´ otesis sobre efectos aleatorios .................... 67 3.5.1.-Prueba de interceptos aleatorios ...................... 69 3.5.2.-Prueba de pendientes aleatorias ...................... 70 3.6.-An´ alisis utilizando la funci´ on lmer() ...................... 73 3.6.1.-Resultados b´ asicos ............................. 73 3.6.2.-P-valores basados en simulaci´ on: el m´ etodo simulate.mer() . . . . 78 3.6.3.-Prueba para los interceptos aleatorios ................... 82 3.6.4.-Prueba para las pendientes aleatorias ................... 83 3.7.-Conclusiones .................................... 86 Bibliograf´ ıa 89 8 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 9 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. donde yi,Xi,βyεison el vector de respuestas continuas, la matriz de dise˜ no, el vector de parametros de efectos fijos y el vector de errores residuales para el grupo i, respectivamente, mientras que Ziybison la matriz de covariables y el vector correspondiente de efectos aleatorios: Zi≡   z(1) i1z(2) i1. . . z(q) i1 . . .. . ..... . . z(1) iniz(2) ini. . . z(q) ini   =z(1) iz(2) i. . . z(q) i,bi≡   bi1 . . . biq   (2) De manera similiar a la matriz de dise˜ no Xi, la matriz Zicontiene valores conocidos de q covariables, con los correspondientes efectos no observados bi. Adem´ as, bi∼ Nq(0,D),εi∼ Nni(0,Ri)con bi⊥εi,(3) i.e., los errores residuales εipara el mismo grupo son independientes de los efectos aleatorios bi. Este supuesto juega un papel clave para distinguir un LMM cl´ asico de un LMM extendido. Adem´ as, suponemos que los vectores de efectos aleatorios y los errores residuales para diferentes grupos son independientes entre s´ ı, i.e., bies independiente de ε′ ipara i=i′. Tambi´ en especificamos que D=σ2DyRi=σ2Ri(4) donde σ2es un par´ ametro de escala desconocido. En general, asumiremos que DyRison definidas positivas, a menos que se indique lo contrario. La representaci´ on (4) no es ´ unica. Para que sea identificable, especificaremos la estructura de la matriz Rien t´ erminos de un conjunto de par´ ametros para una funci´ on de varianza y una matriz de correlaci´ on. Lo veremos en la secci´ on 2.3. La especificaci´ on implicar´ a reestricciones sobre Rique hagan que (4) sea identificable. Adem´ as de los par´ ametros de efectos fijos βpara las covariables utilizadas en la construcci´ on de la matriz de dise˜ no Xi, el modelo (1) incluye dos componentes aleatorias: los errores residuales εiy los efectos aleatorios bipara las covariables incluidas en la matriz Zi. La presencia de efectos fijos y aleatorios de variables conocidas da lugar al nombre del modelo. En muchos casos, los efectos (aleatorios) en bitienen sus correspondientes efectos (fijos), contenidos en β. En consecuencia, la matriz Zise crea a menudo seleccionando un subconjunto de columnas apropiado de la matriz Xi. En tal caso, se dice que los efectos aleatorios correspondientes est´ an “acoplados”. Denominaremos al LMM cl´ asico definido en (1)-(4), LMM de un solo nivel. Dicho modelo, se puede adaptar a datos agrupados en m´ ultiples niveles. Por ejemplo, un modelo para datos con dos niveles de agrupamiento, con observaciones agrupadas en Ngrupos de primer nivel (indexados por i= 1,...,N), cada uno con ni(sub)grupos de segundo nivel (indexados por j= 1,..., ni) que contienen nij observaciones, se puede escribir como yij =Xijβ+Z1,ijbi+Z2,ijbij +εij (5) con bi∼ Nq1(0,D1),bij ∼ Nq2(0,D2),yεij ∼ Nnij (0,Rij), 16 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. donde los vectores aleatorios bi,bij yεij son independientes entre s´ ı. En el modelo (5), bi son los efectos aleatorios asociados a los grupos de primer nivel, mientras que bij son los efectos aleatorios, independientes de los efectos aleatorios del primer nivel, asociados a los grupos de segundo nivel. Las matrices de dise˜ no Z1,ij yZ2,ij pueden ser id´ enticas, pero no tienen por qu´ e serlo. Denominaremos a este modelo LMM de dos niveles. 2.1.2 Especificaci´ on para todos los datos En esta secci´ on, describiremos la especificaci´ on del LMM de un solo nivel, dada por (1)-(4), para todos los datos. Sea y≡(y′ 1,y′ 2,...,y′ N)′el vector que contiene todos los n=PN i=1 nivalores observados de la variable dependiente. De la misma forma, sean b≡(b′ 1,b′ 2,...,b′ N)′yε≡(ε′ 1,ε′ 2,...,ε′ N)′ los vectores que contienen todos los Nq efectos aleatorios y nerrores residuales, respectivamente. Definamos las matrices X≡     X1 X2 . . . XN     ,Z≡     Z10. . . 0 0 Z2. . . 0 . . .. . ..... . . 0 0 . . . ZN     (6) donde 0denota una matriz con todos los elementos nulos. En general, Xes de dimensi´ on n×p, mientras que Zes de dimensi´ on n×Nq. El modelo dado por (1)-(4) puede entonces escribirse para todos los datos de la forma : y=Xβ+Zb +ε,(7) con b∼ NNq(0,σ2D)yε∼ Nn(0,σ2R),(8) donde D≡IN⊗D=     D0. . . 0 0D. . . 0 . . .. . ..... . . 0 0 . . . D     ,R≡     R10. . . 0 0R2. . . 0 . . .. . ..... . . 0 0 . . . RN     (9) donde ⊗denota el producto Kronecker. Es importante destacar que la estructura particular, en forma diagonal por bloques, de las matrices Z,DyR, presentadas en (6)y(9), respectivamente, se debe a la naturaleza del LMM de un solo nivel. Este modelo, definido por (1)-(4), implica una jerarqu´ ıa espec´ ıfica de datos y efectos aleatorios, la cual queda claramente reflejada en (3). En particular, el modelo asume que los efectos aleatorios para distintos grupos, definidos por los niveles de un factor espec´ ıfico, son independientes. Sin embargo, es posible formular modelos de efectos aleatorios utilizando la representaci´ on (7) con matrices Z,DyRque no sean diagonales por bloques. 17 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 2.2 El modelo lineal extendido de efectos mixtos En ciertos casos, el supuesto de que los errores residuales εisean independientes de los efectos aleatorios bipuede resultar demasiado restrictivo. Al relajar este supuesto, se obtiene el LMM extendido. Este modelo se especifica utilizando (1)-(2), y sustituyendo (3) por bi∼ Nq(0,D),εi|bi∼ Nni(0,Ri)(10) con DyRidescompuestas como en (4). Nos referimos a la especificaci´ on anterior como una especificaci´ on jer´ arquica. Observemos que si suponemos que εien (10) es independiente de los efectos aleatorios bi, obtenemos el LMM cl´ asico, especificado por (1)-(4). Por lo tanto, el LMM extendido ofrece un enfoque de modelado m´ as flexible y general en comparaci´ on con el LMM cl´ asico. Una especificaci´ on jer´ arquica de un LMM extendido de dos niveles, correspondiente a (5), ser´ ıa equivalente a suponer que bi∼ Nq1(0,D1),bij |bi∼ Nq2(0,D2),yεij |bi,bij ∼ Nnij (0,Rij). 2.3 Distribuciones definidas por las variables aleatorias yyb Tanto el LMM cl´ asico (secci´ on 2.1) como el LMM extendido (secci´ on 2.2) introducen dos variables aleatorias continuas: yyb. Ambos modelos se han descrito mediante dos funciones de densidad, las cuales son fundamentales para la definici´ on de los LMM. La primera es la distribuci´ on incondicionada de los efectos aleatorios b(no observados), definida en (8). La segunda es la distribuci´ on condicionada de la variable dependiente (aleatoria) y, bajo el supuesto de que los efectos aleatorios son conocidos. En las dos secciones siguientes, ofreceremos una descripci´ on m´ as detallada de ambas distribuciones, completando as´ ı la especificaci´ on del modelo para los LMM cl´ asicos y extendidos. 2.3.1 Distribuci´ on incondicionada de efectos aleatorios La distribuci´ on incondicionada fb(bi)de los efectos aleatorios bi, definida por (3), es una distribuci´ on normal multivariante de media cero y matriz de varianza-covarianza D. Teniendo en cuenta (4), escribimos D(σ2,θD) = σ2D(θD)(11) donde θDes un vector de par´ ametros que representa las varianzas y covarianzas de los elementos de bi. 18 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Obs´ ervese que, seg´ un (11), la matriz D, utilizada para definir la matriz de varianzacovarianza de los efectos aleatorios bi, esta parametrizada utilizando el vector de par´ ametros θD. En muchos casos, se supone que dos elementos cualesquiera del vector bipueden correlacionarse y no hay reestricciones impuestas a la matriz D, excepto que sea definida positiva y sim´ etrica. En ese caso, Dtiene una estructura general de una matriz definida positiva, con q(q+ 1)/2elementos distintos correspondientes a qvarianzas y q(q−1)/2covarianzas de los efectos aleatorios incluidos en bi. En consecuencia, θDcontiene q(q+ 1)/2par´ ametros distintos. Aunque qes t´ ıpicamente peque˜ no, estimar todos los par´ ametros puede ser dif´ ıcil si, por ejemplo, el tama˜ no de la muestra nes limitado. En tal situaci´ on, podemos elegir una estructura de Dsimplificada; por ejemplo, suponer una forma diagonal,i.e, suponer que todos los elementos del vector bisean independientes. As´ ı, θDcontiene s´ olo qpar´ ametros distintos. 2.3.2 Distribuci´ on condicionada de ydados los efectos aleatorios Notemos que, de (1)-(4), se deduce que , para los LMM cl´ asicos, la distribuci´ on condicionada, fy|b(yi|bi), de yidado bies una normal multivariante, con media y varianza definidas como: E(yi|bi)≡µi=Xiβ+Zibi(12) Var(yi|bi) = σ2Ri,(13) con µi≡(µi1,...,µini)′y E(yij |bi)≡µij =x′ ijβ+z′ ijbi,(14) donde xij ≡x(1) ij ,...,x(p) ij ′ yzij ≡z(1) ij ,...,z(q) ij ′ son vectores columna que contienen los valores de los predictores XyZpara la j-´ esima observaci´ on del i-´ esimo grupo. As´ ı, condicionalmente a los valores (desconocidos) de los efectos aleatorios bi, el valor medio del vector de la variable dependiente yise define como una combinaci´ on lineal de los vectores de las covariables XyZ, que est´ an incluidas como columnas en las matrices de dise˜ no espec´ ıficas del grupo XiyZi, correspondientes a los efectos fijos βy los efectos aleatorios bi, respectivamente. Adem´ as, la matriz de varianza-covarianza condicional de yies igual a la matriz de varianza-covarianza de los errores residuales εi. Generalmente, los LMM no son identificables debido a la no unicidad de la representaci´ on (4) y porque potencialmente contienen demasiados par´ ametros desconocidos. 2.3.3 Distribuciones adicionales definidas por yyb En esta secci´ on, presentamos distribuciones auxiliares adicionales relacionadas con los LMM. Se fundamentan en las distribuciones definidas en las secciones 2.3.1 y2.3.2, y son clave en diversos aspectos del ajuste del modelo y en la validaci´ on de sus resultados. 19 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 2.3.3.1 Distribuci´ on conjunta de yyb La distribuci´ on conjunta fy,b(yi,bi)de yybpara los LMM cl´ asicos se puede especificar tomando el producto de la distribuci´ on incondicionada de los efectos aleatorios b, dada en la secci´ on 2.3.1, y la distribuci´ on condicionada de y, definida en la secci´ on 2.3.2: fy,b(yi,bi) = fy|b(yi|bi)fb(bi). Dado que las distribuciones fb(b)yfy|b(y|b)fb(b)son normales multivariantes, la distribuci´ on conjunta tambi´ en sigue una normal. 2.3.3.2 Distribuci´ on marginal de y La funci´ on marginal fy(yi)de yise obtiene integrando los efectos aleatorios bide la distribuci´ on conjunta de yiybi. Es decir, calculamos la funci´ on de densidad de la distribuci´ on marginal de yicomo fy(yi) = Zfy,b(yi,bi)dbi=Zfy|b(yi|bi)fb(bi)dbi,(15) donde fy,b es la funci´ on de densidad de la distribuci´ on conjunta de yiybi,fy|bes la distribuci´ on condicionada de yidado bi, y fbes la funci´ on de densidad de la distribuci´ on incondicionada de bi. Dado que fy,b yfbson densidades de distribuciones normales multivariantes, la distribuci´ on marginal de ytambi´ en es normal multivariante y se puede derivar anal´ ıticamente. 2.3.3.3 Distribuci´ on de bdado y La distribuci´ on fb(bi)de efectos aleatorios bidefinida en (3) no depende de los valores observados yi. Luego, en el contexto bayesiano, es conocida como distribuci´ on a priori de bi. Suponiendo que los valores observados de yison iguales a y(obs) i, la denominada distribuci´ on a posteriori de bicondicionada a y(obs) ise puede calcular usando la siguiente f´ ormula general: fb|y(bi|yi)≡fb|y(bi|yi=yobs i) = fy|b(yi|bi)fb(bi) Rfy|b(yi|bi)fb(bi)dbi (16) Suponiendo que los par´ ametros βyθson conocidos, la distribuci´ on a posteriori fb|y(bi|yi) para los LMM cl´ asicos es normal multivariante. 20 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 2.4 Estimaci´ on En esta secci´ on, presentaremos m´ etodos para obtener un conjunto de estimaciones de los par´ ametros θ,σ2,θD, y θRpara el LMM cl´ asico, definido por (1)-(4). En la secci´ on 2.4.1, presentamos el modelo marginal impl´ ıcito en el LMM cl´ asico. En la secci´ on 2.4.2, describimos brevemente las modificaciones necesarias de los m´ etodos, con un enfoque particular en aquellos implementados en R. En la secci´ on 2.4.4 explicamos los m´ etodos para evaluar la incertidumbre en las estimaciones de los par´ ametros. Finalmente, para completar la descripci´ on de los enfoques de estimaci´ on, en la secci´ on 2.4.5 analizamos brevemente los enfoques alternativos a los presentados en la secci´ on 2.4.2. 2.4.1 El modelo marginal impl´ ıcito en el LMM cl´ asico Para el LMM cl´ asico, tenemos que la media marginal y la matriz de varianza-covarianza de yise dan de la siguiente manera: E(yi) = Xiβ,(17) Var(yi)≡ Vi(σ2,θ,vi) = σ2Vi(θ,vi) = σ2[ZiD(θD)Z′ i+Ri(θR;vi)],(18) donde θ′≡(θ′ D,θ′ R). Para simplificar la notaci´ on, de ahora en adelante, en general, vamos a suprimir el uso de θyven las f´ ormulas. De (17)y(18), se sigue que, marginalmente, yi∼ NniXiβ,σ2ZiDZ′ i+σ2Ri.(19) El valor medio marginal del vector de la variable dependiente yise define mediante una combinaci´ on lineal de los vectores de covariables incluidos, como columnas, en la matriz de dise˜ no Xi, con par´ ametros β. Adem´ as, la matriz de varianza-covarianza de yiconsta de dos componentes. El primero, σ2ZiDZ′ i, es aportado por los efectos aleatorios bi. El segundo, σ2Ri, esta relacionado con los errores residuales εi. Por lo tanto, estrictamente hablando, el modelo que emplea efectos aleatorios, especificado en (1)-(4), implica una distribuci´ on normal marginal, definida por (19). Es importante se˜ nalar que el modelo marginal, definido por (17)-(19), no incluye los efectos aleatorios bi. Por lo tanto, la matriz Dno necesita ser tratada como una matriz de varianzacovarianza. En consecuencia, no es necesario que sea definida positiva, siempre que la matriz Visea sim´ etrica. Esto nos lleva a concluir que, aunque cada modelo marginal de la forma especificada en (1)-(4) implica un modelo marginal definido por (19), no todos los modelos de la forma (19) pueden interpretarse como derivados de un LMM. Por lo tanto, ambos modelos no son equivalentes. 21 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 2.4.2 Estimaci´ on por m´ axima verosimilitud En t´ erminos generales, la estimaci´ on por m´ axima verosimilitud (MV) consiste en construir la funci´ on de verosimilitud a partir de una funci´ on de distribuci´ on de probabilidad adecuada para los datos observados. Las distribuciones incondicionada de biy condicionada de yidado bi, definidas en la secci´ on 2.3 para el LMM cl´ asico, no son apropiadas para construir la funci´ on de verosimilitud, ya que los efectos aleatorios bino son observables. Por una raz´ on similar, no es posible utilizar la distribuci´ on conjunta de yiybi. En cambio, la estimaci´ on del LMM se basa en la distribuci´ on marginal de yi, que coincide con la expresada en (19). Por ello, la estimaci´ on de los par´ ametros del LMM cl´ asico puede lograrse utilizando la estimaci´ on MV o la estimaci´ on por m´ axima verosimilitud restringida (REMV) para el modelo marginal impl´ ıcito. En particular, la estimaci´ on por MV se basa en la verosimilitud marginal resultante de (19). Podemos expresar la verosimilitud logar´ ıtmica como sigue: ℓ(β,σ2,θ)≡ −N 2log(σ2)−1 2 N X i=1 log[det(Vi)]− N X i=1 1 2σ2(yi−Xiβ)′V−1 i(yi−Xiβ),(20) donde Vi, definida en (18), depende de θ. Las estimaciones de los par´ ametros β,σ2yθse obtienen mediante el uso de la probabilidad de perfil logar´ ıtmico, una t´ ecnica que se emplea para estimar par´ ametros en modelos de efectos mixtos y otros modelos estad´ ısticos complejos. Dicha probabilidad resulta de introducir en (20) los estimadores de βyσ2, dados por ˆ βMV(θ)≡ N X i=1 X′ iV−1 iXi!−1 N X i=1 X′ iV−1 iyi!,(21) ˆσ2 MV(θ)≡1 n N X i=1 r′ iV−1 iri,(22) donde ri≡ri(θ) = yiˆ βMV(θ). Consideremos ahora la funci´ on de verosimilud restringida al logaritmo dada por ℓ∗ REMV(σ2,θR)≡lˆ βMV(θR),σ2,θR+p 2log(σ2)−1 2log "det N X i=1 X′ iR−1 iXi!#. A partir de dicha funci´ on, reemplazaremos el estimador ˆσ2 MV por d Var(b β)≡b σ2 N X i=1 X′ ib V−1 iXi!−1 (23) con ridefinido como en (22). Como resultado, se obtiene una versi´ on perfilada de la funci´ on de verosimilitud restringida, que est´ a expresada ´ unicamente en t´ erminos de θ: ℓ∗ REMV(θ)≡− n−p 2log N X i=1 r⊤ iri!−1 2 N X i=1 log [det(Vi)] −1 2log "det N X i=1 X′ iV−1 iXi!# (24) 22 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. La maximizaci´ on de la expresi´ on (24) proporciona un estimador para θ, que posteriormente se introduce en (21)y(23) para proporcionar estimadores de βyσ2, respectivamente. 2.4.3 M´ ınimos cuadrados penalizados En esta secci´ on presentamos un enfoque ligeramente distinto para la estimaci´ on de los par´ ametros β,σ2yθen el LMM cl´ asico definido por las ecuaciones (1)-(4) . Este m´ etodo se fundamenta en la funci´ on de verosimilitud restringida perfilada para θ, descrita en la secci´ on 2.4.2. No obstante, su implementaci´ on se optimiza mediante un algoritmo num´ erico basado en matrices dispersas, lo cual mejora notablemente la eficiencia computacional del enfoque basado en m´ ınimos cuadrados penalizados (PnLS). Consideramos un LMM de un solo nivel, especificado para el conjunto completo de datos (secci´ on 2.1.2). Adem´ as, asumimos independencia condicionada y homogeneidad en la varianza del error residual, es decir, R≡In. En el enfoque de estimaci´ on PnLS, el punto de partida es la funci´ on de densidad de la distribuci´ on conjunta de yy los efectos aleatorios b, los cuales son introducidos en t´ erminos generales en la secci´ on 2.3.3. El logaritmo de la funci´ on de densidad de la distribuci´ on conjunta de yy los efectos aleatorios bviene dado por h(y,b;β,σ2,θ)≡− n+Nq 2log(σ2)−1 2log [det(D)] −(y−Xβ−Zb)′(y−Xβ−Zb) + b′D−1b 2σ2 (25) donde XyZse definieron en (7), mientras que Dse especific´ o en (9). Cabe se˜ nalar que, bajo el supuesto de que R≡In, desde la ecuaci´ on (25) hasta el final de esta secci´ on, se tiene que θ=θD. Aplicando la siguiente repreresentaci´ on de Cholesky: D=TSST′(26) donde Tes una matriz triangular inferior con todos los elementos de la diagonal iguales a 1 y Ses una matriz diagonal con los elementos de la diagonal no negativos. De esta forma, expresamos bde la siguiente manera: b=TSu,con u∼ NNq0,σ2INq. Deducimos entonces que ycondicionada a usigue una normal con E(y|u) = Xβ+ZTSu ≡Xβ+A′u, Var(y|u) = σ2In(27) 23 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. mientras que la media y varianza marginales de yse expresan como E(y) = Xβ, Var(y) = σ2(A′A+In).(28) As´ ı, reescribimos (25) de la forma: hPnLS(y,b;β,σ2,θ)≡− n+Nq 2log(σ2) −(y−Xβ−A′u)′(y−Xβ−A′u) + u′u 2σ2 ≡− n+Nq 2log(σ2)−d(β,θ) 2σ2. (29) El t´ ermino d(β,θ)en (29) parece una suma de cuadrados penalizados, y podemos verlo como una suma de cuadrados residual en un modelo de regresi´ on lineal Ey 0=A′X INq 0u β≡X∗u β. La soluci´ on (˜ u′,˜ β′)′del problema de regresi´ on lineal satisface expl´ ıcitamente que AA′+INq AX X′A′X′X˜ u ˜ β≡Ay X′y.(30) Con el objetivo de optimizar el uso de memoria y reducir la complejidad num´ erica, se propone la utilizaci´ on de una matriz de Cholesky triangular inferior y dispersa L=LZ0 LZX LX, que satisface LL′=P(X∗)′X∗P′.(31) La matriz ortogonal Pcorresponde a una matriz de permutaci´ on dise˜ nada espec´ ıficamente para realizar una reducci´ on de relleno (fill-reducing permutation). Esta permutaci´ on se determina analizando el patr´ on de elementos distintos de cero en la matriz Z, con el objetivo de reorganizar sus filas y/o columnas de manera que se minimice la aparici´ on de nuevos elementos distintos de cero durante el proceso de descomposici´ on. Al aplicar esta permutaci´ on, se logra reducir significativamente el n´ umero de elementos no nulos en la matriz triangular inferior L, obtenida mediante la descomposici´ on de Cholesky. Como consecuencia, se disminuyen considerablemente los requisitos de almacenamiento asociados a L. Cabe destacar que la estructura de Ldepende de θ. Si suponemos que la matriz Pes de la forma P=PZ0 0 PX, 24 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. obtendremos que P′ ZLZ0 P′ XLZX P′ XLXL′ ZPZL′ ZX PX 0 L′ XPX=AA′+INq AX X′A′X′X.(32) Por tanto, podemos reescribir (29) de la forma: hPnLS(y,b;β,σ2,θ)≡− n+Nq 2log(σ2)−˜ d(θ) 2σ2 −1 2σ2PZ(u−˜ u) PX(β−˜ β)′ LL′PZ(u−˜ u) PX(β−˜ β), (33) donde ˜ d(θ)es el valor de la suma de cuadrados penalizada de d(β,θ), definida en (29), calculada en la soluci´ on (˜ u′,˜ β′)′del problema (30). Por lo tanto, ˜ d(θ)es el valor m´ ınimo de la suma de cuadrados penalizada, asumiendo que θes conocido. La funci´ on de verosimilitud marginal que corresponde a (33), viene dada por ℓMV(β,σ2,θ)≡− n 2log σ2−1 2log [det(LZ)]2−˜ d(θ) 2σ2 −1 2σ2hL′ XPX(β−˜ β)i′L′ XPX(β−˜ β). (34) Dado θ, el estimador de βes ˜ β, definido en (30), mientras que el estimador de σ2viene dado por ˜σ2 MV ≡˜ d(θ) n(35) Al sustituir ˜ βy˜σ2 MV en (34), se obtiene la funci´ on de verosimilitud del perfil logar´ ıtmico para el par´ ametro θ: ℓ∗ MV(θ)≡ −1 2log [det(LZ)]2−n 2log[ ˜ d(θ)].(36) La maximizaci´ on de la expresi´ on (36) con respecto a θda lugar al estimador de m´ axima verosimilitud del vector de par´ ametros.A continuaci´ on, este estimador se utiliza para obtener los estimadores MV de σ2yβa partir de las expresiones (35) y (30), respectivamente. Los estimadores obtenidos corresponden a los dados en (21)y(22). El estimador REMV de θse obtiene maximizando la funci´ on de verosimilitud restringida del perfil logar´ ıtmico: ˜ ℓ∗ REMV(θ)≡ −1 2log [det(LZ) det(LX)]2−n−p 2log[ ˜ d(θ)].(37) El estimador obtenido se emplea posteriormente para calcular el estimador de σ2, el cual coincide con el especificado en la ecuaci´ on (23): ˜σ2 REMV ≡˜ d(θ) n−p(38) La estimaci´ on de βse calcula a partir de (30). 25 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. En algunas circunstancias, puede ser de inter´ es construir un intervalo de confianza para σ. Este intervalo puede construirse usando la distribuci´ on chi-cuadrado (χ2). M´ as espec´ ıficamente, un intervalo de confianza del (1 −α)·100 % para σ, utilizando el estimador REMV (ecuaci´ on (23)), est´ a dado por: "ˆσ2 REMV(n−p) χ2 1−α/2, n−p ,ˆσ2 REMV(n−p) χ2 α/2, n−p# donde χ2 α/2, n−pes el percentil (α/2) ·100 de la distribuci´ on chi-cuadrado con n−pgrados de libertad. Si el intervalo de confianza se basa en el estimador de m´ axima verosimilitud (ecuaci´ on (22)), entonces se debe reemplazar n−ppor nen la f´ ormula anterior. Por otro lado, los intervalos de confianza para los par´ ametros θR, relacionados con la matriz Ri, y para σpueden obtenerse de forma similiar. Los intervalos de confianza para los par´ ametros que describen la estructura de la matriz D pueden obtenerse al representar dicha matriz en t´ erminos de varianzas (o desviaciones est´ andar) y correlaciones. 32 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 33 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 34 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Cap´ ıtulo 3 Implementaci´ on con R En este cap´ ıtulo implementaremos lo descrito anteriormente con R, utilizando la funci´ on lme() del paquete nlme(Jos´ e Pinheiro, Douglas Bates 2007) y la funci´ on lmer() del paquete lme4. Para ello, utilizaremos los datos del ensayo ARMD, del paquete nlmeU(Andrzej Galecki, Tomasz Burzykowski), para el modelado de la agudeza visual. Cuando hablamos de ARMD nos refereimos a la Degeneraci´ on Macular Relacionada con la Edad que es un deterioro de la regi´ on macular del ojo. Los datos de ARMD provienen de un ensayo cl´ ınico multic´ entrico aleatorizado que compara un tratamiento experimental (interfer´ on-α) con un placebo en pacientes diagnosticados con ARMD. Nos centraremos en la comparaci´ on entre el placebo y la dosis m´ as alta (6 millones de unidades diarias) de interfer´ on-α. Los pacientes con degeneraci´ on macular pierden progresivamente la visi´ on. En el ensayo, la agudeza visual de cada uno de los 240 pacientes fue evaluada al inicio (l´ ınea base) y en cuatro momentos posteriores a la aleatorizaci´ on: a las 4, 12, 24 y 52 semanas. La agudeza visual se evalu´ o seg´ un la capacidad del paciente para leer l´ ıneas de letras en carteles estandarizados de visi´ on. Estos carteles muestran l´ ıneas de cinco letras de tama˜ no decreciente, que el paciente debe leer de arriba (letras m´ as grandes) hacia abajo (letras m´ as peque˜ nas). Cada l´ ınea en la que se leen correctamente al menos cuatro letras se considera una “l´ ınea de visi´ on”. En nuestros an´ alisis, nos centraremos en la agudeza visual definida como el n´ umero total de letras le´ ıdas correctamente. Por lo tanto, para cada uno de los 240 pacientes, contamos con datos longitudinales en forma de hasta cinco mediciones de agudeza visual, recolectadas en distintos momentos que son comunes para todos los pacientes. Estos datos ser´ an ´ utiles para ilustrar el uso de Modelos Lineales Mixtos (LMM) para datos continuos y longitudinales. En este cap´ ıtulo se presentan distintos modelos lineales mixtos, explorando c´ omo var´ ıan las estructuras de los efectos aleatorios y la forma en que se modela la varianza de los errores residuales. En la secci´ on 3.1, se analiza un modelo con intercepto aleatorio bajo la suposici´ on de que los errores residuales son homoced´ asticos, es decir, que presentan varianza constante a lo largo del tiempo. Este modelo b´ asico permite capturar diferencias individuales en el nivel medio (intercepto), pero asume que la variabilidad residual no cambia con el tiempo. Luego, en la secci´ on 3.2, se introduce un modelo tambi´ en con intercepto aleatorio, pero en el que los errores residuales son heteroced´ asticos. En este caso, se permite que la varianza de los errores cambie con el tiempo, y espec´ ıficamente se modela como una funci´ on potencia del tiempo. Esto resulta ´ util cuando la variabilidad de las mediciones tiende a aumentar o disminuir a medida que avanza el estudio. Las secciones 3.3 y3.4 ampl´ ıan los modelos anteriores incorporando tanto interceptos como pendientes aleatorias asociadas al tiempo. Esto significa que no solo se permite que cada individuo tenga un punto de partida diferente (intercepto), sino tambi´ en una tasa 35 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. de cambio distinta a lo largo del tiempo (pendiente). Adem´ as, estos modelos consideran errores heteroced´ asticos, lo que proporciona una estructura m´ as flexible y realista para los datos longitudinales. En la secci´ on 3.5, se aborda el problema de la prueba de hip´ otesis sobre los efectos aleatorios. Este an´ alisis permite evaluar si ciertos componentes aleatorios son necesarios para describir adecuadamente la variabilidad entre sujetos, o si pueden ser eliminados del modelo sin perder capacidad explicativa. En la secci´ on 3.6, se repite el an´ alisis para algunos de los modelos seleccionados, esta vez utilizando la funci´ on lmer() del paquete lme4 en R, una herramienta muy utilizada para el ajuste de modelos lineales mixtos. Y por ´ ultimo, en la secci´ on 3.7 se exponen las conclusiones del cap´ ıtulo. 3.1 Modelo con interceptos aleatorios y varianza residual homog´ enea Los interceptos aleatorios se utilizan en modelos lineales mixtos para capturar la variabilidad entre individuos o grupos en el punto de partida (intercepto). En lugar de suponer que todos los individuos comienzan desde el mismo valor, los interceptos aleatorios permiten que cada sujeto o grupo tenga su propio valor inicial, lo que refleja mejor las diferencias individuales en la variable dependiente. Esto es ´ util cuando los datos est´ an agrupados o son repetidos, como en estudios longitudinales o jer´ arquicos. Comenzaremos con un modelo simple, que denominaremos M3.1. Este modelo contiene interceptos aleatorios espec´ ıficos para cada sujeto y errores residuales homoced´ asticos. En consecuencia, las observaciones de un mismo individuo, que comparten el intercepto aleatorio, est´ an correlacionadas con un coeficiente de correlaci´ on constante (positivo). Marginalmente, esto corresponde a una estructura de correlaci´ on de simetr´ ıa compuesta con un par´ ametro de correlaci´ on mayor que cero. La estructura de simetr´ ıa compuesta es demasiado simple para describir la estructura de varianza-covarianza de las mediciones de agudeza visual. Por lo tanto, presentamos el modelo M3.1 principalmente para ilustrar los pasos fundamentales en la especificaci´ on y ajuste de un modelo lineal mixto (LMM). 3.1.1 Especificaci´ on del modelo Definiremos el modelo M3.1 de la forma: VISUALit =β0+β1×VISUAL0i+β2×TIEMPOit +β3×TRATAMIENTOi +β4×TRATAMIENTOi×TIEMPOit +b0i+ϵit (44) El t´ ermino VISUALit en la ecuaci´ on (44) denota el valor de agudeza visual medido para el paciente i(i= 1,...,240) en el tiempo t(t= 1,2,3,4, correspondiente a las semanas 4, 12, 24 y 52, respectivamente). En la parte de efectos fijos del modelo, dada por las dos primeras l´ ıneas de la ecuaci´ on (44), VISUAL0irepresenta el valor de agudeza visual medido al inicio del estudio; TIEMPOit es el tiempo de la medici´ on, correspondiente al valor t; 36 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. TRATAMIENTOies el indicador de tratamiento, igual a 1 para el grupo activo y 0 en caso contrario; y TRATAMIENTOi×TIEMPOit es la interacci´ on entre ambos. El par´ ametro β0es el intercepto general del modelo; β1describe el cambio en la agudeza visual media debido a un aumento unitario en la agudeza visual al inicio; β2describe el cambio asociado a una variaci´ on de una semana en el tiempo; β3representa un efecto general del tratamiento, y β4describe el cambio adicional debido a una variaci´ on de una semana en el tiempo para los pacientes que recibieron el tratamiento activo. En la parte de efectos aleatorios del modelo, dada por la ´ ultima l´ ınea de la ecuaci´ on (44), b0irepresenta un intercepto aleatorio espec´ ıfico para cada paciente, que se asume distribuido normalmente con media cero y varianza d11. Por su parte, ϵit es un t´ ermino de error aleatorio residual, tambi´ en asumido como normalmente distribuido con media cero y varianza σ2. Cabe se˜ nalar que, desde un punto de vista formal, el intercepto aleatorio b0irepresenta una desviaci´ on individual respecto a la intercepto fijo β0. No obstante, es habitual referirse a b0icomo un intercepto aleatorio espec´ ıfico del sujeto, aunque en realidad tanto β0como b0iest´ an “acopladas” (v´ ease secci´ on 2.1.1), y ambas contribuyen conjuntamente al intercepto espec´ ıfico del paciente. La parte fija del modelo M3.1 asume que el perfil promedio es lineal en el tiempo, con diferentes interceptos y pendientes para los grupos de tratamiento con placebo y tratamiento activo. Los perfiles espec´ ıficos de cada sujeto tambi´ en se asumen lineales en el tiempo, con interceptos aleatorios (espec´ ıficos de cada sujeto) que desplazan los perfiles individuales con respecto a la tendencia lineal promedio. En notaci´ on matricial, el modelo para el sujeto icon un conjunto completo de cuatro mediciones de agudeza visual se expresa de la siguiente manera:     VISUALi1 VISUALi2 VISUALi3 VISUALi4    =    1VISUAL0i4TRATAMIENTOi4·TRATAMIENTOi 1VISUAL0i12 TRATAMIENTOi12 ·TRATAMIENTOi 1VISUAL0i24 TRATAMIENTOi24 ·TRATAMIENTOi 1VISUAL0i52 TRATAMIENTOi52 ·TRATAMIENTOi          β0 β1 β2 β3 β4       +    1 1 1 1    ×b0i+    ϵi1 ϵi2 ϵi3 ϵi4    (45) Podemos escribir (45) de la forma (1)-(3), si definimos: yi≡    VISUALi1 VISUALi2 VISUALi3 VISUALi4    ,(46) Xi≡    1VISUAL0i4TRATAMIENTOi4·TRATAMIENTOi 1VISUAL0i12 TRATAMIENTOi12 ·TRATAMIENTOi 1VISUAL0i24 TRATAMIENTOi24 ·TRATAMIENTOi 1VISUAL0i52 TRATAMIENTOi52 ·TRATAMIENTOi    ,Zi≡    1 1 1 1    ,(47) 37 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. ϵi≡    ϵi1 ϵi2 ϵi3 ϵi4    ,β≡      β0 β1 β2 β3 β4      ,bi≡b0i,(48) con D ≡ d11, y Ri≡σ2I4.(49) La parte aleatoria del modelo M3.1, especificada por la ecuaci´ on (49), conduce, de acuerdo con la ecuaci´ on (18), a la siguiente matriz de varianzas y covarianzas marginal para el sujeto i con cuatro observaciones: Vi≡ZiDZ′ i+σ2I4=    1 1 1 1    d11 1 1 1 1+    σ20 0 0 0σ20 0 0 0 σ20 0 0 0 σ2     =    σ2+d11 d11 d11 d11 d11 σ2+d11 d11 d11 d11 d11 σ2+d11 d11 d11 d11 d11 σ2+d11    .(50) 3.1.2 Sintaxis y resultados en R En el panel R1, usamos la funci´ on lme() para ajustar el modelo M3.1, el cual est´ a especificado en las ecuaciones (44)-(49). La f´ ormula formula1, utilizada en el panel R1, define la parte fija del modelo, tal como se describe en (44), incluyendo una interacci´ on entre el tiempo y el tratamiento. El factor treat.f est´ a parametrizado tomando “Placebo” como nivel de referencia. El argumento random = ˜1 | subject especifica interceptos aleatorios para cada sujeto. Por defecto, lme() asume que los errores residuales son independientes y tienen una varianza constante, denotada como σ2. Adem´ as, como en la llamada a lme() no se especifica ning´ un argumento method, se utiliza por defecto la estimaci´ on por REMV. Para cambiar a la estimaci´ on por MV, deber´ ıamos a˜ nadir el argumento method = "ML" en la llamada a la funci´ on. Adem´ as de especificar el modelo, en el panel R1 tambi´ en mostramos los resultados del ajuste del modelo M3.1. Los resultados presentados en el panel R1 indican que la desviaci´ on est´ andar √d11 de los interceptos aleatorios, tal como se especifica en (49), se estima en 8,9782, mientras que la desviaci´ on est´ andar de los residuos, σ, se estima en 8,6275. Cabe se˜ nalar que los p-valores que acompa˜ nan a los estad´ ısticos tpara los coeficientes de efectos fijos corresponden a pruebas basadas en el enfoque marginal. Un resumen de las estimaciones basadas en REMV para el modelo M3.1 tambi´ en se presenta en la tabla T1. 38 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Panel R1: Modelo M3.1 usando la funci´ on lme() > formula1 <- + formula(visual ˜ visual0 + time + treat.f + treat.f:time) > modelo1 <- lme(formula1,random = ˜1|subject, data = armd) > summary(modelo1) Linear mixed-effects model fit by REML Data: armd AIC BIC logLik 6591.971 6625.286 -3288.986 Random effects: Formula: ˜1 | subject (Intercept) Residual StdDev: 8.978212 8.627514 Fixed effects: list(formula1) Value Std.Error DF t-value p-value (Intercept) 9.288078 2.6818888 631 3.463260 0.0006 visual0 0.826440 0.0446670 231 18.502244 0.0000 time -0.212216 0.0229295 631 -9.255150 0.0000 treat.fActive -2.422000 1.4999667 231 -1.614703 0.1077 time:treat.fActive -0.049591 0.0335617 631 -1.477594 0.1400 Correlation: (Intr) visul0 time trt.fA visual0 -0.920 time -0.185 -0.003 treat.fActive -0.295 0.022 0.335 time:treat.fActive 0.126 0.002 -0.683 -0.476 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -4.18750513 -0.39692515 0.03204783 0.55138252 2.95132118 Number of Observations: 867 Number of Groups: 234 39 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Tabla T1: Estimaciones de par´ ametros basadas en REMV para los modelos M3.1 yM3.2 con interceptos aleatorios espec´ ıficos del sujeto Par´ ametro modelo1 modelo2 Modelo M3.1 M3.2 Valor de log-REMV -3288.986 -3260.563 Efectos fijos Intercepto β09.288 7.067 Agudeza visual en t=0 β10.826 0.866 Tiempo (en semanas) β2-0.212 -0.213 Tratamiento (Activo vs. Placebo) β3-2.422 -2.305 Tiempo × Tratamiento (Activo) β4-0.050 -0.051 reStruct(subject) SD(bi0)√d11 8.978 7.706 Funci´ on de varianza PowerTIEMPOδδ0.314 Escala σ8.627 3.607 Panel R2: Agrupaci´ on de datos/jerarqu´ ıa impl´ ıcita en el modelo M3.1 > getGroupsFormula(modelo1) #F´ ormula de agrupaci´ on ˜subject <environment: 0x000001d7c3fc9eb8> > str(grpF <- getGroups(modelo1)) #Factor de agrupaci´ on Factor w/ 234 levels "1","2","3","4",..: 1 1 2 2 2 2 3 3 3... - attr(*, "label")= chr "subject" > grpF[1:17] [1]11222233344446666 234 Levels: 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17 ... 240 > levels(grpF)[1:5] [1] "1" "2" "3" "4" "6" > range(xtabs(˜grpF)) #M´ ınimo y m´ aximo n´ umero de observaciones [1] 1 4 40 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. En el panel R2, mostramos c´ omo extraer informaci´ on acerca de la agrupaci´ on de los datos o, de manera equivalente, acerca de la jerarqu´ ıa de los datos impl´ ıcita en el modelo ajustado. Al utilizar la funci´ on getGroupsFormula(), obtenemos la expresi´ on de condicionamiento utilizada en la especificaci´ on del argumento aleatorio. Esta expresi´ on indica un solo nivel de agrupaci´ on, definido por los niveles del factor subjet. Al aplicar la funci´ on getGroups() al objeto de ajuste del modelo, extraemos el factor de agrupaci´ on y lo almacenamos en el objeto grpF. Con la ayuda de la funci´ on str(), mostramos la estructura del objeto. La salida implica que el factor de agrupaci´ on tiene 234 niveles (sujetos). Adem´ as, podemos concluir que, por ejemplo, para el primer sujeto tuvimos dos observaciones, para el segundo sujeto tuvimos cuatro observaciones, para el tercero 3 observaciones, etc. Informaci´ on similar se obtiene al listar un subconjunto de los elementos del factor de agrupaci´ on. El n´ umero m´ ınimo y m´ aximo de observaciones a trav´ es de todos los sujetos se obtiene aplicando la funci´ on range() al resultado de una tabla cruzada de los niveles del factor grpF, proporcionada por la funci´ on xtabs(). Para obtener una mejor comprensi´ on de la estructura estimada de varianza-covarianza del modelo M3.1, utilizamos las funciones getVarCov() yVarCorr(), como se muestra en el panel R3. Panel R3: Las matrices estimadas de varianza-covarianza para los efectos aleatorios (D) y para los errores residuales (Ri) del modelo M3.1 (a)Estimaci´ on de la matriz D > getVarCov(modelo1,individual = "2") Random effects variance covariance matrix (Intercept) (Intercept) 80.608 Standard Deviations: 8.9782 > VarCorr(modelo1) subject = pdLogChol(1) Variance StdDev (Intercept) 80.60829 8.978212 Residual 74.43401 8.627514 (b)Estimaci´ on de la matriz Ri > getVarCov(modelo1,type = "conditional",individual = "2") subject 2 Conditional variance covariance matrix 1234 1 74.434 0.000 0.000 0.000 2 0.000 74.434 0.000 0.000 3 0.000 0.000 74.434 0.000 4 0.000 0.000 0.000 74.434 Standard Deviations: 8.6275 8.6275 8.6275 8.6275 41 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Un gr´ afico modificado de los residuos por cada punto temporal y grupo de tratamiento puede resultar m´ as ´ util. Para este prop´ osito, utilizamos la forma del comando plot() mostrada en el panel R7 apartado (b). Es importante notar que, en la f´ ormula del gr´ afico, se utiliza el argumento type = "pearson" dentro de la funci´ on resid(), lo que indica que se emplear´ an los residuos de Pearson. Adem´ as, en la f´ ormula del gr´ afico se incluye el t´ ermino ˜ time |treat, lo que permite generar paneles separados del gr´ afico por grupo de tratamiento a lo largo del tiempo. Por ´ ultimo, mediante el argumento id = 0.05 en la instrucci´ on plot(), se etiquetan autom´ aticamente aquellos residuos cuyo valor absoluto supera el percentil 97.5 de la distribuci´ on normal est´ andar, utilizando el n´ umero de observaci´ on correspondiente del marco de datos armd. Obtenemos como resultado el siguiente gr´ afico: time Standardized residuals −4 −2 0 2 10 20 30 40 50 46 68 68 87 91 93 104 104 107 121 121 135 137 137 143 162 178 209 227 227 Placebo 10 20 30 40 50 40 51 56 56 70 73 73 73 75 77 90 112 112 120 151 165 191 200 Active En la figura F2 presentamos una versi´ on mejorada del gr´ afico, en la que se superponen diagramas de caja y bigotes (boxplots) sobre un gr´ afico de puntos (stripplot) de los residuos para cada punto temporal y grupo de tratamiento. Para ello, se utiliza la funci´ on bwplot() del paquete lattice (Deepayan Sarkar), como se muestra en el panel R7 apartado (b). En el primer argumento de bwplot() se emplea una f´ ormula que solicita un gr´ afico de los residuos de Pearson frente a los niveles del factor time.f, diferenciando los paneles por los niveles del factor treat.f. Los residuos se extraen del objeto ajustado del modelo modelo2 mediante la funci´ on resid(). La figura F2 permite evaluar la distribuci´ on de los residuos de Pearson condicionales para cada momento de medici´ on y grupo de tratamiento. A pesar de la estandarizaci´ on, la variabilidad de los residuos parece no ser constante. El gr´ afico tambi´ en revela la presencia de varios valores at´ ıpicos, es decir, residuos que superan, en valor absoluto, el percentil 97.5 de la distribuci´ on normal est´ andar. No obstante, dado el gran n´ umero de observaciones, es esperable encontrar algunos valores at´ ıpicos. Es importante destacar que estos outliers aparecen en todos los grupos de tratamiento y en todos los puntos temporales. 48 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. resid(modelo2, type = "p") −4 −2 0 2 4wks 12wks 24wks 52wks Placebo 4wks 12wks 24wks 52wks Active Fig. F2: Gr´ aficos de puntos (y diagramas de caja y bigotes) de los residuos de Pearson condicionales para cada punto temporal y grupo de tratamiento del modelo M3.2 Panel R8: Lista de los residuos de Pearson condicionales at´ ıpicos para el modelo M3.2 > id <- 0.05 > outliers.idx <- + within(armd, + { + resid.p <- resid(modelo2, type = "pearson") + idx <- abs(resid.p) > -qnorm(id/2) + }) > outliers <- subset(outliers.idx, idx) > nrow(outliers) [1] 38 > outliers$subject [1] 40 46 51 56 56 68 68 70 73 73 73 75 77 87 90 [16] 91 93 104 104 107 112 112 120 121 121 135 137 137 143 151 [31] 162 165 178 191 200 209 227 227 234 Levels: 1 2 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ... 240 49 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. El panel R8 muestra los sujetos para los cuales se etiquetaron residuos at´ ıpicos en la figura F2. Para ello, los residuos de Pearson condicionales se extraen del objeto ajustado del modelo modelo2 y se almacenan en el vector resid.p. Los ´ ındices correspondientes a los residuos cuyo valor absoluto es mayor que el percentil 97.5 de la distribuci´ on normal est´ andar se almacenan en el vector l´ ogico idx. El marco de datos outliers.idx contiene las variables seleccionadas del conjunto de datos armd, junto con los residuos y el vector de ´ ındices l´ ogicos. El marco de datos outliers es un subconjunto de outliers.idx y contiene las observaciones para las cuales el valor de la variable idx, dado como el segundo argumento de la funci´ on subset(), es igual a 1. Existen 38 observaciones de este tipo, para las cuales se imprime el n´ umero de sujeto correspondiente. Es importante destacar que, para varios sujetos, hay m´ as de un residuo at´ ıpico, ya que es posible que haya m´ as de una medici´ on de agudeza visual por sujeto. La figura F3 muestra el gr´ afico normal Q-Q de los residuos de Pearson condicionales para cada punto temporal. El gr´ afico se obtuvo utilizando la primera llamada a la funci´ on qqnorm() que aparece en el panel R7 apartado (c). Los patrones muestran algunas desviaciones de una tendencia lineal. Tambi´ en podemos observar el gr´ afico normal Q-Q de los efectos aleatorios predichos (interceptos aleatorios). Los efectos son estimados mediante EBLUPs (ver secci´ on 2.5.1). Estos efectos pueden ser extra´ ıdos del objeto ajustado del modelo modelo2 utilizando la funci´ on ranef(), como se muestra en la segunda llamada a la funci´ on qqnorm() en el panel R7 apartado (c). El gr´ afico Q-Q resultante se muestra en la figura F4 y tiene una ligera curvatura. Esto podr´ ıa interpretarse como una indicaci´ on de no normalidad de los efectos aleatorios. Sin embargo, como se menciona en la secci´ on 2.5.1, dicho gr´ afico no necesariamente refleja la verdadera distribuci´ on de los efectos aleatorios. Por lo tanto, debe interpretarse con cautela. Residuals Quantiles of standard normal −2 0 2 −40 −20 0 20 4wks 12wks 24wks −40 −20 0 20 −2 0 2 52wks Fig. F3: Gr´ aficos Q-Q normales de los residuos de Pearson condicionales para cada punto temporal del modelo M3.2 50 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Random effects Quantiles of standard normal −3 −2 −1 0 1 2 3 −20 −10 0 10 (Intercept) Fig. F4: El gr´ afico Q-Q normal de los interceptos aleatorios predichos para el modelo M3.2 Un gr´ afico diagn´ ostico importante se presenta en la figura F5. En ´ el se muestran los valores observados y predichos de las mediciones de agudeza visual para pacientes seleccionados. El panel R9 muestra c´ omo generar el objeto que contiene los datos necesarios para construir esta figura, utilizando la funci´ on augPred(). La funci´ on augPred() permite obtener valores predichos a partir de un objeto especificado como primer argumento. Si el objeto tiene una estructura de agrupamiento, los valores predichos se obtienen para cada grupo. De manera conveniente, la funci´ on tambi´ en a˜ nade las observaciones originales al objeto retornado, el cual es un data frame con cuatro columnas: los valores de la covariable principal, los grupos, los valores predichos u observados, y un indicador que se˜ nala si el valor en la tercera columna es observado o predicho. Entre los argumentos opcionales de la funci´ on augPred() se encuentran: primary,level,length.out, minimum ymaximum. El argumento primary es una f´ ormula unilateral que indica la covariable para la cual deben calcularse los valores predichos. En el ejemplo presentado en el panel R9, se especifica que los valores predichos deben calcularse en funci´ on de la variable time. Los argumentos minimum ymaximum permiten especificar los l´ ımites inferior y superior, respectivamente, para los valores de la covariable principal en los que deben calcularse los valores predichos. Por defecto, estos argumentos toman los valores m´ ınimo y m´ aximo observados de dicha covariable. En la llamada a la funci´ on presentada en el panel R9, se utilizan los valores por defecto, es decir, el m´ ınimo y el m´ aximo de la variable time, que corresponden a 4 y 52 semanas, respectivamente. El argumento level de la funci´ on augPred() es un vector de enteros que especifica los niveles de agrupamiento para los que se deben calcular las predicciones. En la llamada a augPred() mostrada en el panel R9, se emplea level = 0:1, lo que implica que las predicciones se obtienen tanto a nivel poblacional (level = 0) como a nivel de sujeto (level = 1). Por ´ ultimo, el argumento length.out es un n´ umero entero que indica cu´ antos valores de la covariable principal deben utilizarse para calcular las predicciones. El valor por defecto es 51. En el panel R9, se establece length.out = 2, es decir, las predicciones se calculan 51 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. en dos valores de time, el m´ ınimo (4 semanas) y el m´ aximo (52 semanas). Estos dos puntos son suficientes para describir la tendencia lineal (a nivel poblacional e individual) en el tiempo continuo, seg´ un la forma ajustada del modelo M3.2 Panel R9: Valores predichos de agudeza visual para el modelo M3.2 > aug.Pred <- augPred(modelo2,primary = ˜time,level = 0:1,length.out = 2) > plot(aug.Pred, layout = c(4, 4, 1), key = list(lines = list(lty = c(1,2)), text = list(c("Marginal", "Subject-specific")), columns = 2)) time visual 0 20 40 60 80 10 20 30 40 50 1 2 10 20 30 40 50 3 4 6 7 8 0 20 40 60 80 9 0 20 40 60 80 10 11 12 13 14 10 20 30 40 50 15 16 10 20 30 40 50 0 20 40 60 80 17 Marginal Subject−specific Fig. F5: Valores observados y predichos de agudeza visual para pacientes seleccionados seg´ un el modelo M3.2 (L´ ınea azul : Marginal, l´ ınea naranja : Subject-specific) Al aplicar, en el panel R9, la funci´ on plot() al objeto aug.Pred con el argumento level = 0:1, se obtiene un gr´ afico que representa los valores predichos tanto a nivel poblacional como a nivel individual (por sujeto). El argumento layout = c(4, 4, 1) solicita que se genere una ´ unica p´ agina de gr´ aficos, dispuestos en 4 filas y 4 columnas. Es decir, se visualizan 16 gr´ aficos en total, cada uno correspondiente a un sujeto diferente; por lo tanto, se representan las predicciones de los primeros 16 sujetos del estudio. Finalmente, el argumento key permite definir la leyenda del gr´ afico, la cual se ubica en la parte superior del mismo para facilitar la interpretaci´ on de los distintos tipos de predicciones representadas. 52 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. El gr´ afico resultante se presenta en la figura F5. Las medias poblacionales predichas, mostradas en el gr´ afico, disminuyen linealmente con el tiempo. Seg´ un la estructura del modelo asumido, las medias poblacionales se ajustan para cada individuo mediante interceptos aleatorios espec´ ıficos por sujeto. Como consecuencia de esto, las pendientes de los perfiles individuales son iguales para todos los sujetos. En otras palabras, todas las l´ ıneas espec´ ıficas por sujeto son paralelas a la l´ ınea que representa la media poblacional. No obstante, para algunos pacientes, los perfiles individuales predichos difieren notablemente de los valores observados. Por ejemplo, para los sujetos 4 y 15, las predicciones sugieren una disminuci´ on en la agudeza visual con el tiempo, mientras que los datos observados muestran un aumento. Una forma de mejorar estas predicciones individuales consiste en extender el modelo para incluir no solo interceptos aleatorios, sino tambi´ en pendientes aleatorias espec´ ıficas por paciente. Esta cuesti´ on ser´ a analizada en detalle en la siguiente secci´ on. 3.3 Modelos con interceptos y pendientes aleatorias y funci´ on de varianza residual varPower(·) En esta secci´ on se estudia un modelo que incluye dos efectos aleatorios espec´ ıficos por sujeto: un intercepto aleatorio y una pendiente aleatoria asociada al tiempo. Esto permite que cada individuo tenga su propio punto de partida y ritmo de cambio en el tiempo. Se consideran dos posibles estructuras para la matriz de varianza-covarianza Dde los efectos aleatorios: Una estructura general, que permite correlaci´ on entre el intercepto y la pendiente aleatoria. Una estructura diagonal, que asume independencia entre ambos efectos aleatorios. Adicionalmente, se utiliza la funci´ on de varianza varPower(·) para modelar la heterocedasticidad en los errores residuales, permitiendo que la varianza del error cambie con el tiempo. 3.3.1 Modelo con una matriz general D Para especificar el modelo M3.3 con una matriz de varianza-covarianza general D, se modifica la ecuaci´ on (44) del modelo M3.1 para incluir no solo un intercepto aleatorio, sino tambi´ en una pendiente aleatoria con posible correlaci´ on entre ambos: VISUALit =β0+β1×VISUAL0i+β2×TIEMPOit +β3×TRATAMIENTOi +β4×TRATAMIENTOi×TIEMPOit +b0i+b2i×TIEMPOit +ϵit (53) 53 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. La ecuaci´ on (53) se puede escribir en la forma (1)-(3), al definir yi,Xi,ϵiyβcomo en (46)-(48), pero con Zi=    1 4 1 12 1 24 1 52    ,bi=b0i b2i,(54) y con la estructura de varianza-covarianza de los efectos aleatorios dada por bi∼ N(0,D)yϵi∼ N(0,Ri)(55) donde D=d11 d12 d21 d22(56) yRidado por la ecuaci´ on (51). Obs´ ervese que la forma asumida de la matriz Dimplica que los interceptos aleatorios y las pendientes aleatorias est´ an correlacionados. Por ejemplo, una correlaci´ on positiva entre b0iyb2i significa que, para individuos con un mayor valor inicial de agudeza visual, las mediciones posteriores a la aleatorizaci´ on tender´ an a aumentar m´ as r´ apidamente o a disminuir m´ as lentamente que en pacientes con un valor inicial m´ as bajo. De acuerdo con este modelo, la covarianza marginal entre dos mediciones de agudeza visual del sujeto ien los tiempos t1yt2(t1, t2= 1,2,3,4) se puede expresar como: Cov(yit1, yit2) = 1TIEMPOit1D1 TIEMPOit2+I(t1=t2)σ2(TIEMPOit1)2δ =d11 +d12(TIEMPOit1+TIEMPOit2) + d22 TIEMPOit1TIEMPOit2 +I(t1=t2)σ2(TIEMPOit1)2δ(57) donde I(A)es la funci´ on indicadora que vale 1 si se cumple la condici´ on A, y 0 en caso contrario. Por tanto, la varianza marginal de las mediciones de agudeza visual para el sujeto ien el tiempo tes: Var(yit) = d11 + 2d12TIEMPOit +d22TIEMPO2 it +σ2(TIEMPOit)2δ(58) Esto implica que la varianza se comporta como una funci´ on potencial del tiempo de medici´ on, incluyendo un componente cuadr´ atico. En el panel R10 se ajusta el modelo M3.3, definido por las ecuaciones (53)-(56), actualizando el objeto modelo2. Se utiliza la sintaxis random = ˜1 + time | subject para especificar la estructura de efectos aleatorios, lo cual implica que para cada nivel de la variable subject se consideran un intercepto y una pendiente aleatoria para time. Esta formulaci´ on asume una matriz de varianzas-covarianzas general D, representada internamente mediante un objeto de clase pdLogChol, que garantiza su definici´ on positiva mediante una descomposici´ on de Cholesky logar´ ıtmica. 54 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Panel R10: La matriz ˆ Destimada e intervalos de confianza para los par´ ametros θDdel modelo M3.3 > (modelo3 <-update(modelo2,random = ˜1 + time | subject, + data = armd)) Linear mixed-effects model fit by REML Data: armd Log-restricted-likelihood: -3215.299 Fixed: formula1 (Intercept) visual0 time 4.74032694 0.90929518 -0.21506152 treat.fActive time:treat.fActive -2.26000666 -0.05746769 Random effects: Formula: ˜1 + time | subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 6.9789084 (Intr) time 0.2722504 0.138 Residual 5.1221953 Variance function: Structure: Power of variance covariate Formula: ˜time Parameter estimates: power 0.1074438 Number of Observations: 867 Number of Groups: 234 > getVarCov(modelo3, individual = "2") Random effects variance covariance matrix (Intercept) time (Intercept) 48.70500 0.26266 time 0.26266 0.07412 Standard Deviations: 6.9789 0.27225 > intervals(modelo3, which = "var-cov") Approximate 95% confidence intervals Random Effects: Level: subject lower est. upper sd((Intercept)) 5.9882014 6.9789084 8.1335212 sd(time) 0.2301014 0.2722504 0.3221201 cor((Intercept),time) -0.1272998 0.1382425 0.3852932 Variance function: lower est. upper power 0.01161225 0.1074438 0.2032754 55 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Within-group standard error: lower est. upper 3.958193 5.122195 6.628501 Los resultados b´ asicos del ajuste del modelo M3.3 se presentan en el panel R10, y se detallan m´ as en la tabla T2. En este panel tambi´ en se incluyen los intervalos de confianza (IC) del 95 % para los par´ ametros de la funci´ on de varianza y de la estructura de correlaci´ on, calculados seg´ un los m´ etodos descritos en la secci´ on 2.6.3. Los resultados indican un valor estimado bajo del coeficiente de correlaci´ on entre los efectos aleatorios b0iyb2i, igual a 0,138. El intervalo de confianza de dicho coeficiente sugiere que los dos efectos aleatorios podr´ ıan no estar correlacionados. Por lo tanto, en la siguiente secci´ on se considerar´ a una versi´ on simplificada de la matriz D. Tabla T2: Estimaciones basadas en REMV para modelos de efectos mixtos lineales con interceptos aleatorios y pendientes en el tiempo Par´ ametro modelo3 modelo4 modelo5 Modelo M3.3 M3.4 M3.5 Valor de log-REMV -3215.299 -3215.896 -3214.47 Efectos fijos Intercepto β04.740 5.262 5.441 Agudeza visual en t=0 β10.909 0.900 0.900 Tiempo (en semanas) β2-0.215 -0.215 0.242 Tratamiento (Activo vs. Placebo) β3-2.260 -2.279 -2.655 Tiempo × Tratamiento (Activo) β4-0.057 -0.056 reStruct(subject) SD(bi0)√d11 6.979 7.232 7.236 SD(bi1)√d22 0.272 0.281 0.281 cor((Intercerto),tiempo) ρ12 0.138 Funci´ on de varianza PowerTIEMPOδδ0.107 0.111 0.110 Escala σ15.122 5.031 5.039 56 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.3.2 Modelo con una matriz diagonal D En esta secci´ on, consideramos el modelo M3.4, el cual, de manera similar al modelo M3.3, est´ a definido por la ecuaci´ on (53), pero especificamos que D=d11 0 0d22.(59) Por lo tanto, asumimos que los interceptos aleatorios b0iy las pendientes aleatorias b1itienen varianzas diferentes y no est´ an correlacionados. Panel R11: Definici´ on del modelo M3.4 y definici´ on de los intervalos de confianza para los par´ ametros del mismo > (modelo4 <-update(modelo3, + random = list(subject = pdDiag(˜time)), data = armd)) Linear mixed-effects model fit by REML Data: armd Log-restricted-likelihood: -3215.896 Fixed: formula1 (Intercept) visual0 time 5.2622129 0.8999000 -0.2150309 treat.fActive time:treat.fActive -2.2787559 -0.0564510 Random effects: Formula: ˜time | subject Structure: Diagonal (Intercept) time Residual StdDev: 7.231948 0.2809649 5.031159 Variance function: Structure: Power of variance covariate Formula: ˜time Parameter estimates: power 0.1110751 Number of Observations: 867 Number of Groups: 234 > intervals(modelo4) Approximate 95% confidence intervals Fixed effects: lower est. upper (Intercept) 0.8127705 5.2622129 9.71165544 visual0 0.8246435 0.8999000 0.97515652 time -0.2795380 -0.2150309 -0.15052381 treat.fActive -4.5888197 -2.2787559 0.03130786 time:treat.fActive -0.1505476 -0.0564510 0.03764565 Random Effects: Level: subject 57 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. El panel R14 muestra las formas estimadas de las matrices D,RiyVipara el modelo M3.5. La matriz de varianzas-covarianzas marginal ˆ Viestimada indica una tendencia creciente de las varianzas de las mediciones de agudeza visual a lo largo del tiempo. Adem´ as, la matriz de correlaciones correspondiente sugiere una disminuci´ on en la correlaci´ on entre mediciones realizadas en puntos temporales m´ as distantes. 3.4 Una funci´ on alternativa de varianza residual: varIdent(·) Como se mencion´ o anteriormente, los modelos lineales mixtos presentados en las secciones 3.2 y 3.3 fueron especificados utilizando la funci´ on de varianza varPower(·) (v´ ease la definici´ on de la matriz Λien la ecuaci´ on (51)). Esta funci´ on puede resultar demasiado restrictiva, ya que asume que las varianzas de las mediciones de agudeza visual cambian como una funci´ on potencia del tiempo de medici´ on. En el contexto de los LMM, podr´ ıa ser m´ as adecuado utilizar una funci´ on de varianza m´ as general, con el objetivo de mejorar el ajuste del modelo en comparaci´ on con la funci´ on potencia. Para comprobar esta hip´ otesis, empleamos una prueba de raz´ on de verosimilitudes, basada en los modelos M3.3 yM3.6. Ambos modelos comparten la misma estructura de efectos fijos y aleatorios, dada por la ecuaci´ on (53), pero se diferencian en la especificaci´ on de la matriz Ri. En el modelo M3.3,Ri est´ a definida seg´ un una funci´ on potencia, mientras que en el modelo M3.6,Rise define como: Ri=σ2 1           1 0 0 0 0σ2 2 σ2 1 0 0 0 0 σ2 3 σ2 1 0 0 0 0 σ2 4 σ2 1           ≡σ2    δ2 10 0 0 0δ2 20 0 0 0 δ2 30 0 0 0 δ2 4    ,(61) donde δt≡σt/σ1para t= 1,...,4representa la raz´ on entre la desviaci´ on est´ andar en la ocasi´ on ty la de la primera ocasi´ on, y σ2≡σ2 1. Esta parametrizaci´ on corresponde a una funci´ on de varianza de clase varIdent y est´ a dise˜ nada para permitir la identificaci´ on expl´ ıcita de los par´ ametros δtque controlan la variaci´ on espec´ ıfica de cada ocasi´ on de medici´ on. Para ajustar el modelo M3.6, actualizamos el objeto modelo3 utilizando una forma adecuada del constructor varIdent() dentro del argumento weights de la funci´ on lme(). La sintaxis correspondiente y los resultados del ajuste se muestran en el panel R15 apartado (a), y resultados adicionales aparecen en la tabla T3. El panel R15 apartado (b) tambi´ en incluye el resultado de la prueba de raz´ on de verosimilitudes, calculada con la funci´ on anova() a partir de las verosimilitudes de los modelos M3.6 yM3.3. Obs´ ervese que el modelo M3.3 (nulo) est´ a anidado dentro del modelo M3.6 (alternativo), porque ambos comparten la misma estructura de efectos fijos y aleatorios, pero M3.6 permite mayor flexibilidad en la especificaci´ on de la varianza residual. El resultado de la prueba es estad´ ısticamente significativo al nivel del 5 %, lo cual indica que el uso de la funci´ on de varianza m´ as general varIdent(·) para definir la matriz Ri, seg´ un la ecuaci´ on (61), proporciona un mejor ajuste que la funci´ on varPower(·) usada en el modelo M3.3. Debemos ser cautelosos antes de aceptar esta conclusi´ on. Una inspecci´ on m´ as detallada de los resultados presentados en el panel R15 muestra que el valor estimado del par´ ametro δ4es extremadamente peque˜ no y difiere sustancialmente de los valores estimados de δ2yδ3. Esto resulta sorprendente, ya que todos los an´ alisis previos indicaban que la varianza de la ´ ultima medici´ on de agudeza visual (en la semana 52) era la m´ as alta. 64 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Tabla T3: Estimaciones basadas en REMV para modelos mixtos lineales con interceptos y pendientes aleatorias Par´ ametro modelo6 modelo7 Modelo M3.6 M3.7 Valor de log-REMV -3204.049 -3218.568 Efectos fijos Intercepto β05.103 5.349 Agudeza visual en t=0 β10.901 0.898 Tiempo (en semanas) β2-0.210 -0.215 Tratamiento (Activo vs. Placebo) β3-2.184 -2.314 Tiempo × Tratamiento (Activo) β4-0.059 -0.055 reStruct(subject) SD(bi0)√d11 7.353 SD(bi1)√d22 0.282 Escala σ6.683 Panel R15: Ajuste del modelo M3.6 y prueba de su funci´ on de varianza mediante un test de raz´ on de verosimilitud basado en REMV (a)Ajuste del modelo M3.6 > (modelo6 <- update(modelo3, weights = varIdent(form = ˜1 | time.f))) Linear mixed-effects model fit by REML Data: armd Log-restricted-likelihood: -3204.049 Fixed: formula1 (Intercept) visual0 time treat.fActive 5.10353829 0.90120343 -0.21040652 -2.18433949 treat.fActive time:treat.fActive -2.18433949 -0.05931033 Random effects: Formula: ˜1 + time | subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 7.3462142 (Intr) time 0.3110427 -0.132 Residual 4.6231094 Variance function: Structure: Different standard deviations per stratum Formula: ˜1 | time.f Parameter estimates: 4wks 12wks 24wks 52wks 1.0000000000 1.6252527188 1.7435761307 0.0005327589 Number of Observations: 867 Number of Groups: 234 65 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. (b)Prueba de raz´ on de verosimilitud basada en REMV para la funci´ on de varianza > anova(modelo3, modelo6) Model df AIC BIC logLik Test L.Ratio p-value modelo3 1 10 6450.598 6498.19 -3215.299 modelo6 2 12 6432.099 6489.21 -3204.049 1 vs 2 22.49859 <.0001 Una se˜ nal de los problemas con la estimaci´ on del modelo M3.6 tambi´ en se puede obtener, por ejemplo, intentando calcular los intervalos de confianza para los par´ ametros de la matriz de varianzacovarianza. En particular, al ejecutar el comando: > intervals(modelo6, which = "var-cov") se produce un mensaje de error que indica problemas con la estimaci´ on de la matriz de varianzacovarianza para las estimaciones de los par´ ametros. Finalmente, el problema con la convergencia del algoritmo de estimaci´ on para el modelo M3.6 tambi´ en se refleja claramente en el gr´ afico normal Q-Q de los residuos de Pearson condicionales, mostrado en la figura F10, y obtenido al ejecutar el comando: > qqnorm(modelo6, ˜resid(.)|time.f) Es importante notar que los residuos correspondientes a la semana 52 son todos iguales a 0. Residuals Quantiles of standard normal −2 0 2 −30 −20 −10 0 10 20 4wks 12wks 24wks −30 −20 −10 0 10 20 −2 0 2 52wks Fig. F10: El gr´ afico normal Q-Q de los residuos de Pearson condicionales para el modelo M3.6. El panel correspondiente a las 52 semanas indica el problema con el ajuste del modelo 66 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Para investigar la causa del problema detectado en el ajuste del modelo M3.6, nos fijaremos en las secciones transversales de la superficie de verosimilitud restringida (REMV) para los par´ ametros δ2, δ3, δ4yσ. Si mantenemos fijos los dem´ as par´ ametros en sus valores estimados mediante REMV, la secci´ on transversal correspondiente a δ4, que representa el cociente entre la desviaci´ on est´ andar residual de las mediciones de agudeza visual a las 52 semanas con respecto a la semana 4, mostrar´ a una l´ ınea horizontal pr´ acticamente plana cerca de cero. En concreto, la diferencia entre los valores de la log-verosimilitud restringida para los valores de δ4 dentro de un cierto intervalo y el valor de log-REML reportado (−3204,049) var´ ıa entre 2,4×10−7y −4,0×10−7. Esto sugiere que, bajo el modelo M3.6, los datos contienen muy poca informaci´ on acerca de este par´ ametro, ya que la superficie de verosimilitud es esencialmente plana en esa direcci´ on del espacio de par´ ametros. Adem´ as, la secci´ on transversal de δ4no sugiere un m´ aximo de la funci´ on de verosimilitud, lo cual implica que la estimaci´ on REMV reportada por la funci´ on lme() en el panel R15 apartado (a) no representa un valor ´ optimo. Dada la gran similitud estructural entre los modelos M3.6 yM3.3, surge la pregunta: ¿por qu´ e no hubo problemas evidentes al ajustar el segundo modelo? Aunque los modelos son similares, difieren en la forma de la estructura marginal de varianzacovarianza de las mediciones de agudeza visual. La forma de la covarianza para las mediciones del sujeto ien diferentes tiempos, implicada por el modelo M3.6, es la misma que la del modelo M3.3 y est´ a dada por: Var(yit) = d11 + 2d12TIEMPOit +d22TIEMPO2 it +σ2δ2 t.(62) Las ecuaciones (57)y(62) definen los diez elementos ´ unicos de la matriz de varianza-covarianza marginal Vipara el modelo M3.6 como funciones lineales de siete par´ ametros: d11, d12, d22,σ2, δ2, δ3 yδ4. Dado que el n´ umero de par´ ametros es cercano al n´ umero de ecuaciones, puede haber colinealidad entre los par´ ametros, lo que provocar´ ıa problemas de convergencia del algoritmo de estimaci´ on. Por otro lado, el lado derecho de la ecuaci´ on (58) en el modelo M3.3 involucra menos par´ ametros y utiliza una funci´ on potencia del tiempo, la cual es no lineal respecto al par´ ametro δ. Por lo tanto, es menos probable que ocurra colinealidad. En consecuencia, comparado con el modelo M3.6, el modelo M3.3 impone una restricci´ on adicional sobre la forma de la estructura de varianza-covarianza marginal. Esta restricci´ on limita el espacio param´ etrico, haciendo que los datos sean m´ as informativos para encontrar una soluci´ on ´ optima. Se deduce entonces que, para utilizar la funci´ on Ident(·), se necesitar´ ıan restricciones adicionales en el modelo M3.6. 3.5 Prueba de hip´ otesis sobre efectos aleatorios Como se mencion´ o en la secci´ on 2.6.2, se pueden realizar pruebas formales de hip´ otesis sobre la estructura de varianza-covarianza utilizando la prueba de raz´ on de verosimilitud, basada en la funci´ on de verosimilitud restringida. Un aspecto importante en este tipo de pruebas es la distribuci´ on nula del estad´ ıstico. En particular, cuando los valores de los par´ ametros de varianza-covarianza bajo la hip´ otesis nula se encuentran en el interior del espacio param´ etrico, la distribuci´ on nula sigue una distribuci´ on χ2con grados de libertad iguales a la diferencia en el n´ umero de par´ ametros independientes entre los modelos nulo y alternativo. Ejemplos de tales pruebas se presentan en las secciones 3.3.2 (panel R11)y3.4 (panel R15). No obstante, cuando los par´ ametros bajo la hip´ otesis nula se ubican en la frontera del espacio param´ etrico, la forma exacta de la distribuci´ on nula se vuelve dif´ ıcil de determinar. Como se indic´ o en la secci´ on 2.6.2, en ciertos casos la distribuci´ on del estad´ ıstico de prueba corresponde a una mezcla de distribuciones χ2. Este resultado se ha obtenido bajo la suposici´ on de errores residuales independientes y homoced´ asticos. 67 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. En otros escenarios, la ´ unica alternativa viable es simular la distribuci´ on nula. En el lenguaje R, esto puede lograrse usando la funci´ on simulate() del paquete nlme o la funci´ on exactRLRT() del paquete RLRsim(Fabian Scheipl). Ambas funciones mencionadas solo permiten errores residuales independientes y homoced´ asticos. Adem´ as, la funci´ on exactRLRT() solo admite efectos aleatorios independientes, mientras que simulate() no est´ a definida para objetos de clase gls. Estas limitaciones impiden, por ejemplo, probar si la inclusi´ on de interceptos aleatorios mejora el ajuste del modelo M3.2, en comparaci´ on con un modelo sin efectos aleatorios pero con varianzas residuales modeladas mediante la funci´ on de varianza Power(·). Por la misma raz´ on, tampoco es posible contrastar la significancia estad´ ıstica de extender el modelo M3.2 con pendientes aleatorias, lo cual da lugar al modelo M3.3. En estos casos, el grado de credibilidad de las modificaciones en la estructura de los efectos aleatorios debe evaluarse utilizando, por ejemplo, diagn´ osticos de residuos y/o aplicando criterios de informaci´ on. Panel R16: Valores del AIC para los modelos del M3.1 al M3.5 > AIC(modelo1, modelo2, modelo3, modelo4) df AIC modelo1 7 6591.971 modelo2 8 6537.126 modelo3 10 6450.598 modelo4 9 6449.792 > modelo4ml <- update(modelo4, method = "ML") > modelo5ml <- update(modelo5, method = "ML") > anova(modelo4ml, modelo5ml) Model df AIC BIC logLik Test L.Ratio p-value modelo4ml 1 9 6437.974 6480.859 -3209.987 modelo5ml 2 8 6437.371 6475.491 -3210.685 1 vs 2 1.39716 0.2372 El panel R16 presenta los valores del Criterio de Informaci´ on de Akaike (AIC) para los modelos del M3.1 al M3.5. Se puede observar, por ejemplo, que el AIC del modelo M3.2 es de 6537,126 , mientras que el del modelo M3.3 es considerablemente menor, con un valor de 6450,598. Esta diferencia sugiere que el modelo M3.3 se ajusta mejor a los datos. Adem´ as, como se muestra en la figura F9, los valores predichos por el modelo M3.3 siguen m´ as de cerca a los valores observados, en comparaci´ on con el modelo M3.2 (v´ ease figura F5). Es importante destacar que el menor valor de AIC se obtiene con el modelo M3.5, lo cual sugiere que este proporciona el mejor ajuste global. Este resultado refleja las decisiones adoptadas respecto a la estructura de efectos aleatorios durante el proceso de desarrollo del modelo. En el resto de esta secci´ on, se ilustra el uso de resultados anal´ ıticos y funciones de simulaci´ on en R para realizar pruebas de hip´ otesis sobre la estructura de efectos aleatorios, en situaciones donde los valores de los par´ ametros bajo la hip´ otesis nula se encuentran en la frontera del espacio param´ etrico. Para este fin, se consideran diversos modelos aplicados a los datos ARMD, bajo el supuesto de homocedasticidad de los errores residuales. 68 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.5.1 Prueba de interceptos aleatorios Comencemos considerando el modelo M3.1, el cual incluye interceptos aleatorios. Para evaluar si es necesario incorporar interceptos espec´ ıficos para cada sujeto, podemos utilizar una prueba de raz´ on de verosimilitud basada en el m´ etodo de m´ axima verosimilitud restringida (REMV). Esta prueba compara el modelo alternativo M3.1, que incluye interceptos aleatorios, con un modelo nulo que no considera efectos aleatorios y que asume errores residuales homoced´ asticos (es decir, con varianza constante). El objetivo es determinar si la inclusi´ on de interceptos aleatorios mejora significativamente el ajuste del modelo a los datos. Panel R17: La prueba de raz´ on de verosimilitud basada en REMV para ausencia de interceptos aleatorios en el modelo M3.1 (a)Usando 0,5χ2 0+ 0,5χ2 1como distribuci´ on nula > vis.gls1a <- gls(formula1, data = armd) > (anova.res <- anova(vis.gls1a, modelo1)) Model df AIC BIC logLik Test L.Ratio p-value vis.gls1a 1 6 6839.937 6868.492 -3413.968 modelo1 2 7 6591.971 6625.286 -3288.986 1 vs 2 249.9655 <.0001 > (anova.res[["p-value"]][2])/2 [1] 1.321121e-56 (b)Usando la funci´ on exactRLRT() para simular la distribuci´ on nula > library(RLRsim) > exactRLRT(modelo1) simulated finite sample distribution of RLRT. (p-value based on 10000 simulated values) data: RLRT = 249.97, p-value < 2.2e-16 En el panel R17, realizamos la prueba de raz´ on de verosimilitud (LR) basada en REMV, comparando el estad´ ıstico de la prueba con una distribuci´ on nula obtenida mediante dos enfoques: una mezcla de distribuciones χ2y una t´ ecnica de simulaci´ on. En el primer enfoque (panel R17 apartado (a)), ajustamos el modelo nulo sin interceptos aleatorios (objeto vis.gls1a), el cual comparte la misma estructura fija que el modelo alternativo M3.1 (modelo1). Posteriormente, usamos la funci´ on anova() para calcular el valor de la LR. Dado que la hip´ otesis nula sit´ ua la varianza del intercepto aleatorio en el l´ ımite del espacio param´ etrico, la p-valor reportada por anova() es incorrectamente basada en una distribuci´ on χ2 1. La distribuci´ on correcta es una mezcla 50-50 de χ2 0yχ2 1. Por tanto, ajustamos la p-valor dividi´ endola entre 2, concluyendo que el efecto es estad´ ısticamente significativo. Como alternativa (panel R17 apartado (b)), utilizamos la funci´ on exactRLRT() del paquete RLRsim, que estima la p-valor exacta mediante simulaci´ on Monte Carlo (10,000 r´ eplicas). Esta p-valor coincide con la de la mezcla te´ orica, demostrando la importancia de los interceptos aleatorios. Pese a que la funci´ on simulate() del paquete nlme podr´ ıa emplearse para simular la distribuci´ on nula, no est´ a implementada para objetos de clase gls. En la siguiente secci´ on, exploraremos su uso en la prueba de pendientes aleatorias. 69 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.5.2 Prueba de pendientes aleatorias Con fines ilustrativos, consideramos un modelo con interceptos y pendientes aleatorias espec´ ıficas por sujeto, no correlacionadas entre s´ ı, as´ ı como errores residuales independientes y homoced´ asticos. Es decir, consideramos un modelo definido por las ecuaciones (53)–(55), donde la matriz de varianzascovarianzas de los efectos aleatorios, D, se especifica seg´ un (59), y la matriz de errores residuales es: Ri=σ2×I4. Denominaremos a este modelo como M3.7. En esta secci´ on, aplicaremos una prueba de raz´ on de verosimilitud basada en REMV para determinar si es necesario incluir pendientes aleatorias en el modelo. La prueba consiste en comparar los siguientes modelos: Modelo nulo (M3.1): incluye ´ unicamente interceptos aleatorios. Modelo alternativo (M3.7): incluye tanto interceptos como pendientes aleatorias, sin correlaci´ on entre ellos. En el panel R18 se presentan tres enfoques distintos para realizar la prueba de raz´ on de verosimilitud basada en REMV, con el objetivo de determinar si es necesario incluir pendientes aleatorias en el modelo. En primer lugar, panel R18 apartado (a), se ajusta el modelo M3.7, que incluye pendientes aleatorias, partiendo del modelo M3.4 pero bajo el supuesto de varianza residual constante (homocedasticidad). Este modelo se almacena en el objeto modelo7, y sus resultados aparecen en la tabla T3. En el panel R18 apartado (b) uso de una mezcla de distribuciones χ2: En este enfoque, se realiza la prueba LR basada en REMV y se usa como distribuci´ on nula una mezcla del 50 % de χ2 1y 50 % de χ2 2. El estad´ ıstico de la prueba se extrae del objeto an.res (resultado de la funci´ on anova()) y se utiliza en la funci´ on pchisq() para calcular el p-valor. El p-valor ajustado indica que el efecto es estad´ ısticamente significativo, permitiendo rechazar la hip´ otesis nula de que la varianza de las pendientes aleatorias sea cero. A continuaci´ on, en el panel R18 apartado (c), se simula la distribuci´ on nula mediante la funci´ on exactRLRT(). Dado que esta funci´ on solo admite efectos aleatorios independientes, se considera una matriz de covarianzas Ddiagonal. Como el modelo contiene dos componentes de varianza (interceptos y pendientes aleatorios), se deben especificar los argumentos m,m0 ymA. El resultado de la simulaci´ on es un p-valor pr´ acticamente igual a cero, lo cual indica una clara significancia estad´ ıstica. Finalmente, en el panel R18 apartado (d), se utiliza la funci´ on simulate() para generar un gr´ afico que compara los p-valores emp´ ıricas (obtenidas por simulaci´ on del estad´ ıstico de LR) con los p-valores nominales. Este gr´ afico ayuda a decidir qu´ e distribuci´ on nula es m´ as adecuada, es decir, cu´ al refleja mejor lo que pasa en la realidad del modelo que est´ as usando. M´ as espec´ ıficamente, se aplica la funci´ on simulate() a los objetos modelo1 ymodelo7, donde el primero se especifica como el modelo nulo y el segundo, indicado mediante el argumento m2, como el modelo alternativo. El n´ umero de simulaciones del estad´ ıstico de prueba se fija en 100 mediante el argumento nsim. Posteriormente, la instrucci´ on plot() genera un gr´ afico que compara los p-valores emp´ ıricos del estad´ ıstico de la prueba de raz´ on de verosimilitud con los p-valores nominales, calculadas bajo tres distribuciones distintas: la distribuci´ on χ2 1, la distribuci´ on χ2 2, y una mezcla 50 %–50 % de las distribuciones χ2 1yχ2 2. Los grados de libertad requeridos se pasan al argumento df de la funci´ on plot() como un vector num´ erico. Para incluir, por ejemplo, una mezcla del 65 %–35 %, debe usarse expl´ ıcitamente el argumento weights = c(0.65, 0.35). 70 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Panel R18: Prueba de raz´ on de verosimilitud basada en REMV para pendientes aleatorias en el modelo M3.7 (a)Ajuste del modelo M3.7 > (modelo7 <- update(modelo4, weights = NULL,data = armd) ) Linear mixed-effects model fit by REML Data: armd Log-restricted-likelihood: -3218.568 Fixed: formula1 (Intercept) visual0 time 5.34880943 0.89846421 -0.21537026 treat.fActive time:treat.fActive -2.31374378 -0.05505963 Random effects: Formula: ˜time | subject Structure: Diagonal (Intercept) time Residual StdDev: 7.353234 0.2817028 6.683414 Number of Observations: 867 Number of Groups: 234 > intervals(modelo7, which = "var-cov") Approximate 95% confidence intervals Random Effects: Level: subject lower est. upper sd((Intercept)) 6.4144768 7.3532344 8.4293790 sd(time) 0.2434998 0.2817028 0.3258996 Within-group standard error: lower est. upper 6.254600 6.683414 7.141628 (b)Usando 0,5χ2 1+ 0,5χ2 2como distribuci´ on nula > (an.res <- anova(modelo1, modelo7)) Model df AIC BIC logLik Test L.Ratio p-value modelo1 1 7 6591.971 6625.286 -3288.986 modelo7 2 8 6453.137 6491.211 -3218.568 1 vs 2 140.8344 <.0001 > (RLRT <- an.res[["L.Ratio"]][2]) [1] 140.8344 > .5 *pchisq(RLRT, 1, lower.tail = FALSE) + .5 *pchisq(RLRT, 2, lower.tail = FALSE) [1] 1.397126e-31 71 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. (c)Usando la funci´ on exactRLRT() para simular la distribuci´ on nula > mAux <- update(modelo1, random = ˜0 + time|subject, data = armd) > exactRLRT(m = mAux, m0 = modelo1, mA = modelo7) simulated finite sample distribution of RLRT. (p-value based on 10000 simulated values) data: RLRT = 140.83, p-value < 2.2e-16 (d)Usando la funci´ on simulate() para simular la distribuci´ on nula > vis.lme2.sim <- simulate(modelo1, m2 = modelo7, nsim = 100) > plot(vis.lme2.sim, df = c(1, 2),abline = c(0,1, lty=2)) Empirical p−value Nominal p−value 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 df=1 ML Mix(1,2) ML 0.2 0.4 0.6 0.8 df=2 ML df=1 REML 0.2 0.4 0.6 0.8 Mix(1,2) REML 0.2 0.4 0.6 0.8 df=2 REML Fig. F11: P-valores emp´ ıricos y nominales para evaluar la necesidad de pendientes aleatorias en el modelo M3.7 El gr´ afico resultante se presenta en la figura F11. En ´ el se muestran dos filas, cada una compuesta por tres paneles: una fila corresponde a la estimaci´ on mediante REMV, y la otra a la estimaci´ on mediante m´ axima verosimilitud. El gr´ afico revela que los p-valores nominales, calculados utilizando distribuciones χ2 1,χ2 2o una mezcla 50 %–50 % de ambas, son consistentemente mayores que los p-valores obtenidos por simulaci´ on. Esto sugiere que el uso de cualquiera de estas distribuciones como referencia conduce a una prueba conservadora, es decir, con menor probabilidad de rechazar la hip´ otesis nula, incluso cuando hay evidencia que lo justifica. 72 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.6 An´ alisis utilizando la funci´ on lmer() En esta secci´ on, se vuelven a ajustar los modelos M3.1 yM3.7, descritos en las secciones 3.1.1 y3.5.2, respectivamente, utilizando la funci´ on lmer() del paquete lme4. La elecci´ on de estos modelos se basa en el hecho de que dicha funci´ on s´ olo permite considerar errores residuales independientes y homoced´ asticos. Es importante se˜ nalar que ambos modelos no proporcionan un ajuste adecuado a los datos del estudio ARMD, como se pudo concluir a partir de los an´ alisis previos realizados mediante la funci´ on lme(). Por consiguiente, los resultados presentados en esta secci´ on deben interpretarse principalmente como una ilustraci´ on del uso de la funci´ on lmer(), m´ as que como una evaluaci´ on definitiva de los modelos considerados. 3.6.1 Resultados b´ asicos En el panel R19, se muestra c´ omo ajustar el modelo M3.1 utilizando la funci´ on lmer(). Este modelo incluye interceptos aleatorios y asume que la varianza de los errores residuales es constante. Originalmente, este modelo fue ajustado en el panel R1 usando la funci´ on lme(), y su resultado se almacen´ o en el objeto modelo1. En el panel R19 apartado (a), se presenta la sintaxis empleada en la funci´ on lmer() para ajustar el modelo M3.1. Obs´ ervese que la estructura de efectos aleatorios se especifica directamente en el argumento Formula. Adem´ as, es importante destacar que el argumento data se pasa como un data frame simple, en lugar de un objeto de datos agrupados. A diferencia de la funci´ on lme(), el uso de objetos de datos agrupados no solo no es necesario, sino que tampoco se recomienda. El modelo se ajusta utilizando el m´ etodo REMV, que es la opci´ on por defecto. Los resultados del modelo ajustado se muestran mediante la funci´ on gen´ erica summary(). Es relevante notar que los valores de los estad´ ısticos tde los efectos fijos se presentan sin p-valores. Los m´ etodos para calcular p-valores ser´ an tratados m´ as adelante en esta secci´ on. Tambi´ en hemos hallado los valores de los criterios AIC y BIC, as´ ı como el valor del logar´ ıtmo de la verosmilitud. En el panel R19 apartado (b), se extrae directamente la matriz varianza-covarianza de los efectos fijos del objeto del modelo ajustado, utilizando la funci´ on vcov(). Luego, se convierte en matriz de correlaciones mediante la funci´ on cov2cor(). Los resultados mostrados en el panel R19 son coherentes con los presentados anteriormente en el panel R1 y en la tabla T1. 73 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. En el panel R23, se calcula el valor medio, la mediana y los percentiles 2.5 y 97.5 de las estimaciones de los coeficientes de efectos fijos y de los par´ ametros de varianza-covarianza obtenidos al reajustar el modelo M3.1 con los datos simulados. Para ello, en el panel R23 apartado (a), nos enfocamos primero en los coeficientes de los efectos fijos. Utilizamos la funci´ on apply() para calcular estos estad´ ısticos fila por fila en la matriz betaE, la cual contiene las estimaciones simuladas de los coeficientes. Adem´ as, empleamos esta funci´ on para calcular los p-valores emp´ ıricos. Esto se logra, para cada coeficiente (es decir, cada fila de betaE), computando la proporci´ on de estimaciones mayores que cero, y a partir de ello, el p-valor bilateral. Cabe destacar que este p-valor emp´ ırico no puede ser inferior a 1 nsim , lo que impone un l´ ımite inferior por razones de robustez num´ erica. La salida impresa con estos estad´ ısticos resumidos permite evaluar la distribuci´ on emp´ ırica de cada uno de los coeficientes de efectos fijos. En particular, la ´ ultima columna muestra los p-valores emp´ ıricos, que tienden a ser ligeramente mayores (y por tanto m´ as conservadores) que los calculados anal´ ıticamente en el panel R21 apartado (a). En el panel R23 apartado (b), se presenta el c´ odigo para calcular estad´ ısticos similares, media, mediana y percentiles 2.5 y 97.5, para las estimaciones simuladas de √d11 (la desviaci´ on est´ andar de los interceptos aleatorios) y de σ(la desviaci´ on est´ andar de los errores residuales). La estimaci´ on de √d11 se basa en la representaci´ on de la matriz de varianzas dada en la ecuaci´ on (26). Tanto las medias de √d11 como de σson muy cercanas a las estimaciones puntuales reportadas anteriormente en el panel R19. Los percentiles 2.5 y 97.5 pueden interpretarse como intervalos de confianza emp´ ıricos para evaluar la precisi´ on de las estimaciones, en una forma similar a la proporcionada por la funci´ on intervals() del paquete nlme. No obstante, es importante notar que el paquete lme4 no incluye una funci´ on an´ aloga a intervals(). Panel R24: Sintaxis para construir los gr´ aficos de densidad de las estimaciones basadas en simulaciones de los coeficientes de efectos fijos y los par´ ametros de varianza-covarianza para el modelo M3.1 > names(sigmaE) <- names(STe) <- NULL > parSimD1 <- rbind(betaE, ST1 = STe, sigma = sigmaE) > parSimD1t <- data.frame(t(parSimD1), check.names=FALSE) > parSimD1s <- subset(parSimD1t, select = -‘(Intercept)‘) > require(reshape) > densityplot(˜value | variable, + data = melt(parSimD1s), + scales = list(relation = "free"), + plot.points = FALSE) > detach(package:reshape) Podr´ ıa ser de inter´ es presentar la distribuci´ on de las estimaciones basadas en simulaciones de los par´ ametros del modelo M3.1. En el panel R24, se demuestra la sintaxis que puede utilizarse para crear gr´ aficos de las funciones de densidad correspondientes a las funciones de distribuci´ on emp´ ırica para los coeficientes de los efectos fijos y ST =qd11 σ2. 80 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Para ello, comenzamos creando la matriz parSimD1, que contiene, como columnas, las estimaciones basadas en simulaciones de los par´ ametros de inter´ es. Luego, la transponemos y la guardamos en el marco de datos parSimD1t. A continuaci´ on, utilizamos la funci´ on densityplot() del paquete lattice para crear los gr´ aficos. Este paquete se adjunta autom´ aticamente junto con lme4, por lo que no es necesario cargarlo por separado. Sin embargo, es importante se˜ nalar que, en la llamada a la funci´ on densityplot(), aplicamos la funci´ on melt() del paquete reshape. Por lo tanto, debemos cargar tambi´ en dicho paquete. La funci´ on melt() se utiliza para apilar las variables contenidas en el marco de datos parSimD1s en una sola columna, llamada (por defecto) value. Durante este proceso, se crea otra variable, llamada (por defecto) variable, que identifica las correspondientes variables originales de parSimD1s. La f´ ormula que se proporciona como primer argumento en la llamada a la funci´ on densityplot() solicita el gr´ afico de una estimaci´ on de la densidad de n´ ucleo Gaussiano (por defecto) para la distribuci´ on emp´ ırica de los valores de cada variable. Los gr´ aficos resultantes se presentan en la figura F12. Los gr´ aficos de densidad presentados en la figura F12 son relativamente sim´ etricos. Esto sugiere, por ejemplo, que los intervalos de confianza basados en la aproximaci´ on de la distribuci´ on normal de la distribuci´ on emp´ ırica podr´ ıan ser adecuados para la construcci´ on de las estimaciones de intervalo de los par´ ametros. value Density 0 2 4 6 8 0.7 0.8 0.9 1.0 visual0 0 5 10 15 −0.30 −0.25 −0.20 −0.15 time 0.00 0.10 0.20 −8 −6 −4 −2 0 2 treat.fActive 0 2 4 6 8 10 −0.15 −0.050.00 0.05 time:treat.fActive 0 1 2 3 4 5 0.8 0.9 1.0 1.1 1.2 1.3 ST1 0.0 0.5 1.0 1.5 8.0 8.5 9.0 9.5 sigma Fig. F12: Gr´ aficas de densidad para las estimaciones basadas en simulaci´ on del modelo M3.1 81 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.6.3 Prueba para los interceptos aleatorios Panel R25: Prueba de raz´ on de verosimilitud basada en REMV para la inexistencia de interceptos aleatorios en el modelo M3.1 (a)Usando 0,5χ2 0+ 0,5χ2 1como la distribuci´ on nula > formula3 <- visual ˜ visual0 + time + treat.f + treat.f:time > vis.lm2 <- lm(formula3, data = armd) > (RLRTstat <- + -2 *as.numeric(logLik(vis.lm2, REML=TRUE) + - logLik(modelo1mer))) [1] 249.9654 > 0.5 *pchisq(RLRTstat, 1, lower.tail = FALSE) [1] 1.321121e-56 (b)Usando la funci´ on exactRLRT() para simular la distribuci´ on nula > require(RLRsim) > exactRLRT(modelo1mer) simulated finite sample distribution of RLRT. (p-value based on 10000 simulated values) data: RLRT = 249.97, p-value < 2.2e-16 (c)Usando el m´ etodo simulate.mer() para obtener el p-valor emp´ ırico > lm2sim <- simulate(vis.lm2, nsim = 100) > RLRTstatSim <- apply(lm2sim, 2, + function(y){ + dfAux <- within(armd, visual <- y) + lm0 <- lm(formula(vis.lm2), data = dfAux) + llik0 <- as.numeric(logLik(lm0, REML=TRUE)) + llikA <- as.numeric(logLik(refit(modelo1mer, y))) + RLRTstat<- -2 *(llik0 - llikA) + }) > mean(RLRTstat <= RLRTstatSim) [1] 0 En el panel R25, se presentan diferentes enfoques para calcular el p-valor de la prueba de raz´ on de verosimilitudes basada en REMV, con el objetivo de evaluar la necesidad de incluir interceptos aleatorios en el modelo M3.1. 82 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Como se detalla en la secci´ on 2.6.2, la hip´ otesis nula plantea que la varianza de los efectos aleatorios es igual a cero. Dado que este valor est´ a en el l´ ımite del espacio param´ etrico, la distribuci´ on nula asint´ otica de la estad´ ıstica de la prueba corresponde a una mezcla de distribuciones χ2, en particular una combinaci´ on del 50 % de χ2 0y 50 % de χ2 1. En el panel R25 apartado (a), se muestra c´ omo calcular el p-valor utilizando esta distribuci´ on mixta. Se ajusta el modelo bajo la hip´ otesis nula (un modelo lineal cl´ asico homoced´ astico) y se guarda el ajuste en el objeto vis.lm2 de clase lm. Luego, se aplica la funci´ on logLik() para obtener los logaritmos de la verosimilitud REML del modelo nulo y del alternativo. Para obtener el p-valor, se divide por dos el p-valor asociado a una distribuci´ on χ2 1. Este resultado indica evidencia estad´ ısticamente significativa. En el panel R25 apartado (b), se eval´ ua el p-valor mediante simulaci´ on de la distribuci´ on del estad´ ıstico de LR en muestras de tama˜ no finito usando la funci´ on exactRLRT() del paquete RLRsim. Esta funci´ on se aplica al objeto modelo1mer, de clase mer, y el resultado coincide con el obtenido en el panel R17 apartado (b). Finalmente, en el panel R25 apartado (c) se estima un p-valor emp´ ırico simulando m´ ultiples muestras de la variable dependiente bajo el modelo nulo. Para ello, se generan nsim = 100 muestras con la funci´ on simulate() y se almacenan como columnas en el marco de datos lm2sim. Luego, mediante apply(), se calcula el estad´ ıstico de LR para cada muestra. Se crea un marco de datos auxiliar dfAux sustituyendo la variable visual en el conjunto de datos original por cada muestra simulada, se ajusta el modelo nulo y se obtiene su logverosimilitud con logLik(). Posteriormente, se vuelve a ajustar el modelo alternativo M3.1 usando refit() y se extrae su log-verosimilitud. Finalmente, se calcula el estad´ ıstico de LR y se estima el p-valor como la proporci´ on de valores simulados mayores o iguales al observado. En este caso, el p-valor resultante es 0. 3.6.4 Prueba para las pendientes aleatorias En esta secci´ on retomamos el modelo M3.7, definido en la secci´ on 3.5.2. Este modelo incluye interceptos aleatorios y pendientes aleatorias para la variable tiempo, bajo el supuesto de que la varianza residual es constante. Adem´ as, se asume que los interceptos y las pendientes aleatorias son independientes. En el panel R26, se ajusta el modelo M3.7 empleando la funci´ on lmer() y se realiza una prueba para evaluar el t´ ermino de interacci´ on. En particular, en el panel R26 apartado (a) se observa el uso de dos t´ erminos Z: (1|subject) y(0 + time|subject) en la f´ ormula de lmer(). Estos dos t´ erminos permiten emular una matriz diagonal 2×2, denotada como D. Tambi´ en se muestran algunos resultados seleccionados del modelo ajustado, extrayendo componentes adecuados del objeto modelo2mer. El panel R26 apartado (b) muestra la construcci´ on de la prueba de raz´ on de verosimilitudes para evaluar la interacci´ on treat.f:time. Para ello, se reajusta el modelo M3.7 omitiendo el t´ ermino de interacci´ on en la f´ ormula del modelo, y luego se aplica la funci´ on gen´ erica anova() para obtener el resultado de la prueba LR. Cabe destacar que los modelos ajustados, almacenados en los objetos lmer2Dd ylmer3Dd, fueron estimados usando el m´ etodo por defecto, es decir, REMV. Sin embargo, la prueba LR reportada por anova() se basa en m´ axima verosimilitud, dado que se trata de una hip´ otesis sobre un efecto fijo. El p-valor calculado indica que el resultado de la prueba no es estad´ ısticamente significativo al nivel del 5 %. 83 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Panel R26: Modelo M3.7 ajustado utilizando la funci´ on lmer() (a)Ajuste del modelo y extracci´ on de informaci´ on b´ asica > modelo2mer <- + lmer(visual ˜ visual0 + time + treat.f + treat.f:time + + (1|subject) + (0 + time|subject), + data = armd) > summ <- summary(modelo2mer) > coef(summ) Estimate Std. Error t value (Intercept) 5.35170677 2.32814866 2.298696 visual0 0.89841132 0.03924095 22.894740 time -0.21536129 0.03227106 -6.673512 treat.fActive -2.31384454 1.20777347 -1.915793 time:treat.fActive -0.05507348 0.04709814 -1.169334 > unlist(VarCorr(modelo2mer)) subject subject.1 53.73076997 0.07930878 > sigma(summ) [1] 6.690055 (b)Prueba de raz´ on de verosimilitudes para la interacci´ on treat.f:time > modelo2aux <- update(modelo2mer, . ˜ . - treat.f:time) > anova(modelo2aux, modelo2mer) refitting model(s) with ML (instead of REML) Data: armd Models: modelo2aux: visual ˜ visual0 + time + treat.f + (1 | subject) + (0 + time | subject) modelo2mer: visual ˜ visual0 + time + treat.f + treat.f:time + (1 | subject) + (0 + time | subject) npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) modelo2aux 7 6440.8 6474.2 -3213.4 6426.8 modelo2mer 8 6441.4 6479.5 -3212.7 6425.4 1.3764 1 0.2407 En el panel R27 se presentan diferentes enfoques para calcular el p-valor de la prueba de raz´ on de verosimilitudes basada en REMV, con el objetivo de evaluar la necesidad de incluir pendientes aleatorias en el modelo M3.7. En este caso, la distribuci´ on nula de la estad´ ıstica de prueba es, asint´ oticamente, una combinaci´ on al 50 %–50 % de las distribuciones χ2 1yχ2 2, tal como se describe en la secci´ on 2.6.2. En el panel R27 apartado (a) se ilustra c´ omo calcular el p-valor usando esta distribuci´ on mixta. Para ello, se emplea la funci´ on logLik() para extraer los logaritmos de la verosimilitud REMV de los modelos M3.1 yM3.7, representados por los objetos modelo1mer y modelo2mer, respectivamente. Cabe destacar que el valor del estad´ ıstico de LR coincide con el obtenido mediante la funci´ on anova() en el panel R18 apartado (b). El p-valor se obtiene sumando la mitad de los p-valores correspondientes a las distribuciones χ2 1yχ2 2. El resultado indica que la prueba LR es estad´ ısticamente significativa. 84 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. En el panel R27 apartado (b) se eval´ ua el p-valor mediante la simulaci´ on de la distribuci´ on del estad´ ıstico LR para un tama˜ no de muestra finito, utilizando la funci´ on exactRLRT() del paquete RLRsim. En este caso, es necesario ajustar un modelo auxiliar que contenga ´ unicamente pendientes aleatorias como efectos aleatorios, de acuerdo con los requisitos de la funci´ on exactRLRT(). El c´ odigo y los resultados obtenidos son comparables a los del Panel R18 apartado (c). Sin embargo, en el R27 apartado (b) la funci´ on exactRLRT() se aplica a objetos de clase mer. Panel R27: Prueba de raz´ on de verosimilitud basada en REMV para evaluar la ausencia de pendientes aleatorias en el modelo M3.7 (a)Usando 0,5χ2 0+ 0,5χ2 1como la distribuci´ on nula > RML0 <- logLik(modelo1mer) > RMLa <- logLik(modelo2mer) > (RLRTstat <- -2 *as.numeric(RML0 - RMLa)) [1] 140.8318 > .5 *pchisq(RLRTstat, 1, lower.tail = FALSE) + + .5 *pchisq(RLRTstat, 2, lower.tail = FALSE) [1] 1.398952e-31 (b)Usando la funci´ on exactRLRT() para simular la distribuci´ on nula > require(RLRsim) > mAux <- lmer(visual ˜ + visual0 + time + treat.f + treat.f:time + + (0 + time| subject), + data = armd) > exactRLRT(m = mAux, + m0= modelo1mer, + mA= modelo2mer) simulated finite sample distribution of RLRT. (p-value based on 10000 simulated values) data: RLRT = 140.83, p-value < 2.2e-16 85 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 3.7 Conclusiones En este cap´ ıtulo se han ajustado y comparado varios modelos lineales mixtos con el objetivo de analizar la evoluci´ on de la agudeza visual de los pacientes en funci´ on del tratamiento recibido y del tiempo. Se comenz´ o con modelos m´ as sencillos, como el modelo M3.1, que consideraba ´ unicamente interceptos aleatorios y asum´ ıa homogeneidad en la varianza residual. Este modelo, aunque simple, no capturaba adecuadamente la variabilidad entre sujetos ni la din´ amica temporal del proceso. Posteriormente, el modelo M3.2 introdujo heterocedasticidad en la varianza residual mediante la funci´ on varPower(), lo que mejor´ o significativamente el ajuste. A continuaci´ on, el modelo M3.3 incorpor´ o tanto interceptos como pendientes aleatorias con una matriz de varianza-covarianza no restringida para los efectos aleatorios. Este modelo ofreci´ o una descripci´ on muy detallada de la evoluci´ on individual, pero con un elevado n´ umero de par´ ametros, lo que lo hac´ ıa m´ as complejo y susceptible a sobreajuste. Para reducir dicha complejidad, el modelo M3.4 impuso una estructura diagonal en la matriz de efectos aleatorios, es decir, se asumi´ o independencia entre el intercepto y la pendiente aleatoria por sujeto. Esta simplificaci´ on mantuvo un buen nivel de ajuste y redujo la cantidad de par´ ametros a estimar, posicion´ andose como un candidato competitivo. Finalmente, el modelo M3.5 elimin´ o la interacci´ on entre tratamiento y tiempo en los efectos fijos, manteniendo interceptos y pendientes aleatorias con matriz diagonal y varianza heteroced´ astica. Este modelo logr´ o el menor valor de AIC entre todos los comparados, combinando un ajuste adecuado con una estructura m´ as parsimoniosa, i.e., logr´ o describir un fen´ omeno de la manera m´ as simple posible, sin perder precisi´ on ni poder explicativo. Por otro lado, el modelo M3.6 present´ o problemas de convergencia y el modelo M3.7 tuvo menor AIC y BIC que el modelo M3.1, pero no consigui´ o alcanzar al modelo M3.5. En conjunto, el proceso de comparaci´ on permiti´ o evaluar progresivamente la utilidad de cada componente del modelo. El modelo M3.5 fue seleccionado como el m´ as adecuado, al ofrecer el mejor equilibrio entre calidad de ajuste, simplicidad y relevancia interpretativa en el estudio de la agudeza visual en pacientes tratados a lo largo del tiempo. 86 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 87 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. 88 Modelos Lineales de Efectos Mixtos. Librer´ ıas R. Bibliograf´ ıa [1] Andrzej Gałecki • Tomasz Burzykowski. Linear Mixed-Effects Models Using R. Springer, 2013. [2] Juan Carlos Corra Morales, Juan Carlos Salazar Uribe. Introducci´ on a los modelos mixtos. Medell´ ın: Universidad Nacional de Colombia. Facultad de Ciencias. Escuela de Estad´ ıstica, 2016. [3] Douglas Bates,Martin M¨ achler, Benjamin M. Bolker, Steven C. Walker. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 2015. [4] Shonosuke Sugasawa • Tatsuya Kubokawa. Mixed-Effects Models and Small Area Estimation. Springer, 2023. [5] Bates D, Maechler M, Bolker B, Walker S (2014a). lme4: Linear Mixed-Effects Models Using Eigen and S4. R package version 1.1-37 Disponible en: http://CRAN. R-project.org/package=lme4. [6] Santos Nobre, J.,da Motta Singer, J. (2007). Residual analysis for linear mixed models. Biometrical Journal, 49(6), 863–875. [7] Gurka, M. (2006). Selecting the best linear mixed model under REML. The American Statistician, 60(1), 19–26. [8] Deepayan Sarkar. lattice: Trellis Graphics for R. R package version 0.22-7 Disponible en: https://doi.org/10.32614/CRAN.package.lattice. [9] Andrzej Galecki,Tomasz Burzykowski. nlmeU: Datasets and Utility Functions Enhancing Functionality of ’nlme’ Package. R package version 0.70-9 Disponible en: https: //doi.org/10.32614/CRAN.package.nlmeU. [10] Jos´ e Pinheiro,Douglas Bates (2007). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-168 Disponible en: https://doi.org/10.32614/CRAN. package.nlme. [11] Fabian Scheipl . RLRsim: Exact (Restricted) Likelihood Ratio Tests for Mixed and Additive Models. R package version 3.1-8 Disponible en: https://doi.org/10.32614/ CRAN.package.RLRsim. 89