Agro 6998 – Conferencia 3 Modelo mixto: estimación y prueba de hipótesis   Comenzaremos definiendo el modelo lineal de efectos fijos para luego extender dicha definición al caso del modelo lineal mixto. El modelo lineal es ampliamente utilizado en la experimentación para analizar la variabilidad de observaciones (respuestas) realizadas sobre características de importancia en función de una o más variables predictoras o factores. Los modelos de este tipo pueden ser especificados de la forma general:

yi = μ + β1 xi1 + β 2 xi 2 + ... + β p xip + ei Var(ei )=σ 2 . Los términos x de la forma general, asumen los valores 1 ó 0 y son usados para indicar a qué UE y a qué tratamiento corresponde la observación yi; por ejemplo si y3 fue observada sobre la unidad 1 bajo el tratamiento B, entonces los x correspondientes a la UE 1 y al tratamiento B serán 1 y los restantes cero. En notación matricial, el modelo lineal general tiene la forma:

y = Xβ + e donde y es un vector de observaciones, X es una matriz de diseño, β es el vector de parámetros (o efectos fijos) y

e

es el vector de errores, definido como

e = y − E (y ) = y − Xβ . El ejemplo anterior es un caso típico del modelo de ANOVA, donde los términos x representan a factores de clasificación (efectos categóricos) y por tanto la matriz X será una matriz de ceros y unos. Cuando los términos x representan covariables (variables medidas en una escala cuantitativa) en vez de factores, se tiene el modelo clásico de regresión lineal y en ese caso la matriz X contiene los valores de las variables regresoras para cada observación. Para modelar efectos de factores de clasificación se requieren varios parámetros mientras que el efecto de una covariable puede modelarse con uno o pocos parámetros. Los modelos que tienen ambos, factores y covariables, se denominan modelos de análisis de covarianza (ANCOVA).

  1   

Estimación de los parámetros β Utilizando el procedimiento de mínimos cuadrados ordinarios, se puede estimar el vector de parámetros β resolviendo las ecuaciones normales X´Xβ = X´y . La solución está dada por βˆ = ( X´X) − X´y , donde ( X´X) − es una inversa generalizada de

X´X (Searle, 1971). Para hallar una estimación del vector de parámetros, no hace falta hacer suposiciones distribucionales sobre el vector

e.

Si se asumen los

supuestos del modelo de muestreo ideal, i.e. términos de error independientes y normalmente distribuidos con media 0 y varianza σ 2 , entonces, la matriz de covarianzas de βˆ , utilizada para realizar inferencia estadística sobre β , es σ 2 ( X´X) − . Modelo Mixto Extendiendo el modelo lineal general presentado anteriormente a situaciones donde se incorporan efectos aleatorios se tiene el modelo lineal general mixto. La ecuación matricial para el modelo lineal mixto es:

y = Xβ + Zu + e donde y , X , β y e representan las mismas entidades del modelo de efectos fijos y los nuevos componentes son: 1) Z que representa una segunda matriz de diseño de dimensión nxq (matriz especificada exactamente en la misma forma que X , excepto que no incluye una columna para el término constante) y que asocia cada observación a los efectos aleatorios correspondientes y 2) el vector qx1 u de elementos aleatorios ( efectos o coeficientes) que usualmente se asume distribuido N ( 0 , G ). Sobre el vector e se supone distribución N ( 0 , R ), y este vector e es definido como:

e = y − E (y | u) = y − ( Xβ + Zu) Dado que la esperanza del vector aleatorio u es 0 , en el modelo lineal mixto, el valor esperado de una observación es la esperanza incondicional de la media de y (es decir promediada sobre todos los posibles valores de u ):

E (y ) = E ( Xβ + Zu) = Xβ Es decir, los niveles observados de un efecto aleatorio son una muestra aleatoria de la población de niveles y la esperanza incondicional es la media de y sobre toda esa población.   2   

Por otro lado, la esperanza condicional de y dado u es:

E (y | u) = Xβ + Zu esperanza que representa la media de y para el subconjunto específico de niveles del efecto aleatorio observados en el experimento. 2 La matriz R es modelada como R = σ I cuando se considera que los términos de

error (generalmente asociados a la UE) son independientes y tienen la misma varianza

σ 2 . Los términos aleatorios

u

se suponen independientes de los términos aleatorios

e. Resumiendo los supuestos usuales sobre la esperanza y la varianza de las componentes aleatorias, se tiene que:

⎡u ⎤ ⎡ 0 ⎤ E⎢ ⎥=⎢ ⎥ ⎣ e ⎦ ⎣0⎦

⎡u ⎤ ⎡G 0 ⎤ Var ⎢ ⎥ = ⎢ ⎥ ⎣e ⎦ ⎣ 0 R⎦

Cuando se asume distribución normal para el vector de observaciones, la función de densidad (verosimilitud) queda completamente determinada por el vector de valores esperados y la matriz de varianzas y covarianzas. La matriz de varianzas y covarianzas de y (marginal, o promedio para la población de efectos aleatorios) está dada por:

V (y ) = V = V ( Xβ + Zu + e) = ZV (u)Z´+V (e) = ZGZ´+ R Los supuestos clásicos de independencia y homogeneidad de varianzas para los términos aleatorios del modelo lineal general (muestreo ideal) se flexibilizan en el marco del modelo mixto general. La inclusión de efectos aleatorios produce observaciones correlacionadas. Tanto la estructura de correlaciones como la presencia de varianzas heterogéneas pueden ser especificadas a través de la modelación de las matrices de covarianza G y/o R . A través de G y R es posible modelar correlaciones entre efectos de tratamiento, entre parcelas experimentales ocasionadas por la distribución espacial y/o temporal de las mismas en el campo y/o considerar diferentes precisiones de ensayos cuando se combinan experimentos.

  3   

Modelos Marginales versus Modelos Jerárquicos El modelo mixto lineal general puede ser re-escrito como un modelo jerárquico (o modelo condicional):

y | u ~ N ( Xβ + Zu, R ) u ~ N (0, G ) Es decir existe un modelo para y dado u más un modelo para u . Esto sugiere que existen supuestos específicos sobre la dependencia de la media y la estructura de covarianza sobre las covariables en X y Z . La media marginal es Xβ y la estructura de covarianza es V = ZGZ´ + R . Es decir que el modelo implicado para la distribución marginal o incondicional de Y es N ( Xβ, ZGZ´ + R ) . Esta relación entre ambos modelos no se puede aplicar en general, y depende de propiedades de la distribución normal multivariada y de la linealidad del modelo.

Estimación   Las estimaciones por mínimos cuadrados generalizados pueden usarse para estimar los efectos fijos del modelo mixto. Estas estimaciones se obtienen minimizando

(y - Xβ)'V −1 (y - Xβ) ,

y

el

estimador

del

vector

de

efectos

fijos

β

es:

βˆ = ( X´V −1X) − X´V −1y . Si todas las componentes de varianza en V son conocidas este estimador es el mejor estimador lineal insesgado (BLUE) y se corresponde con el estimador máximo verosímil. En la práctica del análisis de datos experimentales V usualmente es desconocida y se reemplaza por su estimador

ˆ ´+ R ˆ = ZGZ ˆ . Si se V

puede asumir que u y e tienen distribución normal, la mejor aproximación para la estimación se logra con métodos basados en máxima verosimilitud. Los métodos de estimación más usados son máxima verosimilitud (ML) y máxima verosimilitud restringida (REML). La función de verosimilitud, L, puede pensarse como la probabilidad de observar los datos que tenemos si los parámetros del modelo fuesen los postulados. Se define usando la función de densidad de las observaciones, en este caso la función normal.

  4   



La estimación de los parámetros fijos será denotada como β ML y la de los parámetros



de la estructura de varianza como ξ ML

Estimador REML El simple ejemplo del estimador ML de la varianza σ 2 de una muestra aleatoria de variables normales, sugiere que cuando μ no es conocida y debe estimarse, dicha estimación introduce un sesgo en el estimador ML de la varianza. La pregunta entonces es, ¿cómo estimar las componentes de varianza sin tener que estimar los parámetros correspondientes a los efectos fijos? La respuesta conduce al estimador REML, sugerido por Patterson y Thompson (1971). En esta aproximación el vector de efectos fijos es eliminado de la función de verosimilitud, y por lo tanto le llamamos “verosimilitud restringida”, que nos sirve para estimar los parámetros de covarianza. Cuando los datos son balanceados, este método nos da estimadores insesgados iguales a los que nos daría un ANOVA. El estimador ML de ξ , basado en t se llama



 

estimador REML ( ξ REML ). La estimación resultante del vector de efectos fijos, β(ξ REML )



suele denotarse por β REML y se obtiene usando mínimos cuadrados generalizados. La idea del estimador REML es la siguiente: Primero se obtiene la verosimilitud basada en datos que en lugar de ser los observados son términos residuales, i.e.

 y - Xβ . Estos términos son conocidos como residuos completos ya que incluyen todas

las fuentes variación aleatoria; se demuestra que los mismos son independientes de



β.

Propiedades del estimador de efectos fijos El estimador de los efectos fijos se obtiene por mínimos cuadrados generalizados



usando ξ en lugar de ξ para construir V . Si E (y ) = Xβ , condicionando sobre las

 

componentes de varianza. Este estimador es insesgado, i.e. E (β(ξ)) = β . Luego, para obtener estimaciones insesgadas relacionadas a los efectos fijos es suficiente que la media de la respuesta sea correctamente especificada. Condicionando sobre ξ , el estimador del vector de efectos fijos tiene covarianza independiente de la Var( y ), si se asume que la matriz Var( y ) se modela correctamente como V = ZGZ´ + R . Por ello este estimador de covarianza suele   5   

llamarse “estimador naif o cándido”. La variabilidad incorporada por reemplazar las componentes de varianza por sus estimadores, no se tiene en cuenta en la construcción del estadístico de Wald que se presenta como candidato para contrastar hipótesis del tipo H 0 : Lβ = 0 , donde L es un arreglo de contrastes conocidos. El estadístico de Wald que se distribuye asintótica mente como una chi-cuadrado con grados de libertad iguales al rango de L , usa la siguiente expresión de varianza:

 −1 var(β) = ( X´V(ξ)X ) Luego, la prueba de Wald, solo proveerá de inferencia válida en caso de muestras grandes. Una alternativa práctica es reemplazar la distribución chi-cuadrado por una distribución F apropiada. El estadístico F para la hipótesis que contrasta efectos fijos mediante la matriz de contrastes L , es: -1  -1   β´L´ ⎡ L X´ V -1 (ξ)X L´ ⎤ Lβ ⎣⎢ ⎦⎥ F= rango(L)

(

)

Bajo la hipótesis nula, la distribución de F se aproxima a la distribución F con grados de libertad en el numerador igual al rango de L. Los grados de libertad del denominador se estiman desde los datos por diversos métodos: 1) método de containment (recomendado en modelos con efectos aleatorios y sin modelación de covarianza residual) , 2) aproximación de Sattherthwaite (casos donde existen efectos aleatorios y modelación de covarianza residual), 3) aproximación de Kenward-Roger (casos donde existen efectos aleatorios y modelación de covarianza residual), 4) Between-within (casos donde solo se modelación de covarianza residual; excepto que el tipo sea sin estructura donde se usa solo Between) y 5) Residual. Cuando existen varias observaciones por sujeto, los grados de libertad del denominador son en general muchos por lo que los tres métodos dan valores-p muy parecidos. Cuando la hipótesis es univariada, i.e. el rango de L es uno, la prueba F se reduce a la clásica prueba T.

Inferencia Si bien la inferencia respecto a la estructura de media es generalmente aquella donde se centra el interés, también la inferencia sobre componentes de varianza es importante ya que: 1) permite interpretar la variación aleatoria, 2) permite identificar estructuras de covarianza sobreparametrizadas que podrían conducir a ineficiencias   6   

en la inferencia para la estructura de medias o bien modelos muy restrictivos que también invalidan la inferencia sobre las medias y 3) en ocasiones, constituyen el objetivo de la investigación por ejemplo, en estudios genéticos para estimar heredabilidad y avance genético. Asintóticamente los estimadores ML y REML tienen distribución normal con la media correcta y matriz de covarianzas igual a la inversa de la matriz de información de Fisher. Luego es posible obtener errores estándares aproximados y realizar la prueba de Wald. En SAS, es posible obtener los estadísticos Z de esta prueba. En la interpretación de modelos jerárquicos algunas componentes de varianza deberían ser cero cuando otras de la distribución en la que se encuentran anidadas son iguales a cero. Por esto, estas pruebas deberían ser interpretadas completamente sólo en el caso de modelos marginales, es decir cuando no se supone una estructura de covarianza asociada a efectos aleatorios para representar la variación entre sujetos. En muchas circunstancias prácticas, los efectos (o coeficientes) de los niveles seleccionados del factor aleatorio no son de interés particular para el experimentador excepto por la información que ellos contienen acerca de la población de efectos. No obstante, existen situaciones donde el interés no se centra sólo en la distribución de los efectos aleatorios sino también en sus valores realizados, es decir en los niveles actualmente considerados en el estudio. En tales situaciones el análisis involucra el cálculo de predictores de esos efectos (o coeficientes) aleatorios. El mejor predictor lineal insesgado (BLUP) de efecto aleatorios constituye el predictor natural de efectos aleatorios en el contexto del modelo mixto lineal general, presentando propiedades óptimas en el sentido de minimizar el error cuadrático medio de la predicción dentro del conjunto de los predictores insesgados (Harville, 1990; Robinson, 1991). La idea subyacente en la predicción de efectos aleatorios implica determinar, para cada nivel del factor de efecto aleatorio, una predicción de su ubicación dentro de la distribución normal de la que proviene. El efecto aleatorio asociado al i-ésimo sujeto indica cómo éste se desvía del valor esperado. El BLUP produce un corrimiento de las predicciones hacia la media general de las observaciones que depende de la relación entre las componentes de varianza involucradas. En general, las predicciones son menos dispersas que las estimaciones. Mientras que la media asociada a factores de efectos fijos es un promedio realizado sobre todos los niveles del efecto en la población, el BLUP es una regresión hacia la media general basada en los componentes de varianza y covarianza   7   

asociados a los efectos aleatorios del modelo (shrinkage estimation). El predictor del vector de efectos aleatorios tiene la forma:

 = G Z'V ( ξ ) -1 (y - Xβ)   u El BLUP se distribuye normalmente con matriz de covarianza igual a: n

 ) = G Z´ ( V -1 − V -1 X( ∑ X´V -1 X )X´V -1 ) ZG var( u i =1

En la inferencia sobre el vector u se debe contemplar la variabilidad en u , por ello esta generalmente se basa en la

 − u ) = G − var( u ) var( u La raíz cuadrada de esta expresión es conocida como error de predicción (EP). Como los parámetros del modelo son generalmente desconocidos y deben ser estimados, en la práctica se calcula el estimador del BLUP( u ) o BLUP empírico como:

   ˆ u = G Z'V ( ξ ) -1 (y - Xβ) Usando la distribución normal antes mencionada, se pueden construir intervalos de predicción para los BLUP de efectos aleatorios, de manera análoga a la construcción de los intervalos de confianza para las medias, i.e. BLUP ± tν,

1-α

EP. También es

común en Agricultura, sobre todo si se desea ordenar material experimental en función de sus BLUP, usar los BLUPt o BLUP estandarizados, obtenidos dividiendo el BLUP empírico por su EP. Estos intervalos así como las pruebas T y F posibles son sólo aproximadas cuando se trabaja con los BLUP empíricos.

Comparando modelos diferentes Al ajustar distintos modelos a un mismo conjunto de datos, es necesario utilizar criterios para la comparación de los ajustes y por tanto para la selección de un modelo. Dos indicadores comúnmente usados son el criterio de información de Akaike (AIC) y el criterio de Schwarz (BIC). En las versiones más modernas de SAS MIXED (SAS Institute, 2001), los criterios AIC y BIC se definen como:

AIC = −2 L + 2d BIC = −2 L + d ln n   8   

donde L es el máximo valor de la función de verosimilitud (restringida), d=q+p es la dimensión del modelo, q es el número de parámetros de covarianza estimados y p es el rango de la matriz de diseño X . Bajo estas expresiones de AIC y BIC, el mejor modelo resulta ser aquel con menor valor para el indicador. Otra alternativa que puede usarse para comparar dos modelos anidados, por ejemplo con igual estructura de media pero diferente estructura de covarianza, o con diferente estructura de medias pero igual covarianza, es la prueba del cociente de verosimilitud, con base en la relación:

 ⎛ L(θ , reducido) ⎞ −2 ln λ = −2 ln ⎜  ⎟ ⎝ L(θ , completo) ⎠ Para la construcción de la prueba a la cantidad –2 ln(L) del modelo con más cantidad de parámetros (modelo completo) se le resta la cantidad –2 ln(L) del modelo reducido. La diferencia obtenida se compara con una distribución χ2 con grados de libertad igual a la diferencia entre el número de parámetros estimados por uno y otro modelo. Si la prueba resulta significativa, el modelo correcto es el más completo, en caso contrario, el modelo reducido es el adecuado. Si bien esta prueba se puede realizar tanto con lo estimadores de máxima verosimilitud como con los estimadores de máxima verosimilitud restringida, el uso de estimadores REML sólo es recomendable para comparar dos modelos que difieren en estructura de covarianza pero con igual media. Esta recomendación se basa en la naturaleza de los estimadores REML. Se discutió que éstos se obtienen maximizando la verosimilitud de un conjunto de contrastes de error derivados de una transformación de la variable respuesta. Si las estructuras de medias de los dos modelos a comparar son diferentes son diferentes, los vectores de contraste de error asociados a cada modelo también diferirán y por tanto si se aplica la prueba del cociente de verosimilitud basada en estimadores REML se estarán comparando verosimilitudes de variables diferentes. Otro problema relacionado con el uso de la prueba del cociente de verosimilitud se presenta cuando se usan para parámetros cuyo valor bajo la hipótesis nula está en la frontera del espacio de valores (por ejemplo, cuando probamos que una componente de varianza es cero).  

  9