Modelos lineales mixtos en R

Luis Cayuela Mayo de 2011

Area de Biodiversidad y Conservacin, Universidad Rey Juan Carlos, o Departamental 1 DI. 231, c/ Tulipn s/n. E-28933 Mstoles (Madrid), a o Espaa. E-mail: [Link]@[Link]. n

128

Modelos lineales mixtos en R (versin 1.2) o

Se autoriza a cualquier persona a utilizar, copiar, distribuir y modicar esta obra con las siguientes condiciones: (1) que se reconozca la autor de la misma; a (2) que no se utilice con nes comerciales; y (3) que si se altera la obra original, el trabajo resultante sea distribuido bajo una licencia similar a sta. e

Para cualquier comentario o sugerencia por favor remitirse al autor de la obra.

129

 Indice
1. Introduccin o 132 1.1. Los paquetes nlme y lme4 . . . . . . . . . . . . . . . . . . . . . . 132 1.1.1. La funcin lme() . . . . . . . . . . . . . . . . . . . . . . . 133 o 1.1.2. La funcin lmer() . . . . . . . . . . . . . . . . . . . . . . . 134 o 1.2. Entendiendo los modelos mixtos . . . . . . . . . . . . . . . . . . 134 1.2.1. Un ejemplo con bentos marino . . . . . . . . . . . . . . . 134 1.2.2. Una aproximacin al anlisis de los datos con modelos o a lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 1.2.3. Re-analizando los datos con modelos lineales mixtos: El paquete nlme . . . . . . . . . . . . . . . . . . . . . . . . . 141 1.2.4. Re-analizando los datos con modelos lineales mixtos: El paquete lme4 . . . . . . . . . . . . . . . . . . . . . . . . . 143 1.2.5. Otros efectos aleatorios en modelos lineales mixtos . . . . 144 1.3. Seleccin de modelos . . . . . . . . . . . . . . . . . . . . . . . . . 146 o 1.4. Anlisis de los residuos . . . . . . . . . . . . . . . . . . . . . . . . 148 a 2. Dise o por bloques n 149

2.1. Paso 1: Aplicacin de un modelo lineal . . . . . . . . . . . . . . . 150 o 2.2. Paso 2: Ajuste de un modelo lineal mixto . . . . . . . . . . . . . 151 2.3. Paso 3: Eleccin de la estructura de los efectos aleatorios . . . . . 152 o 2.4. Paso 4: Eleccin de la estructura de los efectos jos . . . . . . . . 152 o 2.5. Paso 5: Presentacin del modelo nal con REML . . . . . . . . . 153 o 3. Dise o de tipo split-plot n 156

3.1. Paso 1: Aplicacin de un modelo lineal . . . . . . . . . . . . . . . 158 o 3.2. Paso 2: Ajuste de un modelo lineal mixto . . . . . . . . . . . . . 160 3.3. Paso 3: Eleccin de la estructura de los efectos aleatorios . . . . . 161 o 3.4. Paso 4: Eleccin de la estructura de los efectos jos . . . . . . . . 162 o 3.5. Paso 5: Presentacin del modelo nal con REML . . . . . . . . . 163 o 4. Dise os de medidas repetidas n 166

4.1. Paso 1: Aplicacin de un modelo lineal . . . . . . . . . . . . . . . 167 o 4.2. Paso 2: Ajuste de un modelo lineal mixto . . . . . . . . . . . . . 170 4.3. Paso 3: Eleccin de la estructura de los efectos aleatorios . . . . . 171 o 4.4. Paso 4: Eleccin de la estructura de los efectos jos . . . . . . . . 172 o 4.5. Paso 5: Presentacin del modelo nal con REML . . . . . . . . . 172 o 130

 INDICE

INDICE

5. Dise os factoriales anidados o jerarquizados n

175

5.1. Paso 1: Aplicacin de un modelo lineal . . . . . . . . . . . . . . . 177 o 5.2. Paso 2: Ajuste de un modelo lineal mixto . . . . . . . . . . . . . 178 5.3. Paso 3: Eleccin de la estructura de los efectos aleatorios . . . . . 180 o 5.4. Paso 4: Eleccin de la estructura de los efectos jos . . . . . . . . 181 o 5.5. Paso 5: Presentacin del modelo nal con REML . . . . . . . . . 181 o 6. Ms ejemplos a 7. Referencias 185 185

131

Luis Cayuela

Modelos lineales mixtos en R

1.

Introduccin o

Los modelos mixtos son usados cuando los datos tienen algn tipo de u estructura jerrquica o de agrupacin como los diseos de medidas repetidas, a o n las series temporales, los diseos anidados o por bloques aleatorizados. Los n modelos mixtos permiten tener coecientes jos (aquellos cuyos niveles son de inters para el experimentador) y aleatorios (aquellos cuyos niveles son slo e o una realizacin de todos los posibles niveles procedentes de una poblacin) y o o varios trminos de error. Pueden ser una herramienta muy util, pero son e potencialmente dif ciles de comprender y aplicar. En esta sesin intentaremos o dar una visin aplicada de los modelos lineales mixtos con especial referencia a o los distintos tipos de diseos vistos en la sesin anterior. n o Podis encontrar una buena introduccin al tema de los modelos mixtos en el e o cap tulo 8 de Zuur et al. (2007) y en el cap tulo 19 de Crawley (2002). El libro de Zuur et al. (2009) es tambin una buena referencia para profundizar ms en e a el tema de los modelos mixtos sin mucho aditamento matemtico. Otras a referencias con un fondo ms matemtico y, en mi opinin, bastante ms a a o a dif ciles de entender, son Pinheiro & Bates (2000) o Faraway (2006).

1.1.

Los paquetes nlme y lme4

En R hay dos paquetes que nos van a permitir ajustar modelos mixtos. El primero que apareci fue el paquete nlme(), escrito inicialmente por Jos o e Pinheiro (Bell Laboratories) y Douglas Bates (University of Wisconsin) y al que luego se han sumado otros autores. El cdigo de este paquete se escribi o o para que fuera compatible con S y S-Plus (las versiones comerciales de R). > citation("nlme") To cite package 'nlme' in publications use: Jose Pinheiro, Douglas Bates, Saikat DebRoy, Deepayan Sarkar and the R Development Core Team (2010). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-97. A BibTeX entry for LaTeX users is @Manual{, title = {nlme: Linear and Nonlinear Mixed Effects Models}, author = {Jose Pinheiro and Douglas Bates and Saikat DebRoy and Deepayan Sarkar and {R year = {2010}, note = {R package version 3.1-97}, } El otro paquete, que apareci posteriormente, se llama lme4 y ha sido escrito o originalmente por Douglas Bates. El cdigo de este paquete no es compatible o con S y S-Plus y la sintaxis, como veremos ms adelante, cambia ligeramente a con respecto al paquete nlme. 132

Luis Cayuela

Modelos lineales mixtos en R

> citation("lme4") To cite package lme4 in publications use: Douglas Bates <bates@[Link]>, Martin Maechler <maechler@[Link]> and Ben Bolker <bbolker@[Link]> (2011). lme4: Linear mixed-effects models using S4 classes. R package version 0.999375-39. [Link] A BibTeX entry for LaTeX users is @Manual{, title = {lme4: Linear mixed-effects models using S4 classes}, author = {Douglas Bates and Martin Maechler and Ben Bolker}, year = {2011}, note = {R package version 0.999375-39}, url = {[Link] } ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see help("citation") . Ambos paquetes contienen varias funciones relacionadas con modelos mixtos, no solamente modelos lineales mixtos. Las que ms nos van a interesar son las a funciones lme() y lmer() de los paquetes nlme y lme4 respectivamente. 1.1.1. La funcin lme() o

La especicacin de los efectos jos y aleatorios del modelo en la funcin lme() o o se realiza con dos argumentos distintos, xed y random respectivamente. En el caso de que no haya efectos jos en el modelo (todos los efectos son aleatorios) la componente ja del modelo estimar slo la constante o intercept: a o xed = y 1 El argumento xed no es totalmente necesario y se puede omitir en el caso de que no haya efectos jos. En este caso basta con especicar una frmula como o en el caso de la funcin lm() y R entiende que todos los factores explicativos o son aleatorios. La formulacin entonces ser as o a : ya+b+c dnde a, b y c ser factores aleatorios. En el caso de que haya factores jos y o an aleatorios ambos componentes deben ser denidos de manera individual. El componente aleatorio en este caso se especicar de la siguiente forma: a random = 1 |a/b/c Un detalle importante es que el nombre de la variable respuesta y no est a repetido en la frmula de los efectos aleatorios: en su lugar se deja un espacio o en blanco a la izquierda del s mbolo . En la mayor de los modelos mixtos se a 133

Luis Cayuela

Modelos lineales mixtos en R

asumen que los efectos aleatorios tienen una media de cero y que lo que interesa cuanticar es la variacin en la constante (esto es lo que signica el 1, o ver seccin 1.2) causada por las diferencias entre los niveles del factor de los o efectos aleatorios. Despus de la constante viene una barra vertical | que e signica dada la siguiente distribucin de las variables aleatorias. En este o ejemplo concreto hay tres efectos aleatorios con c anidado dentro de b, que a su vez est anidado dentro de a. Los factores estn separados por la barra a a oblicua / y las variables se enumeran de izquierda a derecha en orden decreciente siguiendo una jerarqu espacial o temporal. La formulacin a o completa del modelo utilizando la funcin lme() ser as o a : lme(xed = y 1, random = 1 | a/b/c) 1.1.2. La funcin lmer() o

En la funcin lmer() la frmula se especica en un unico argumento incluyendo o o ambos tipos de efectos. Los efectos jos se especican en primer lugar a la derecha del s mbolo en la forma habitual. A continuacin se especican los o efectos aleatorios precedidos por un + y en parntesis. R identica los efectos e aleatorios porque llevan la barra vertical |. La especicacin de factores o anidados se har con los : en vez de con la barra /. Sin embargo, al a especicar la anidacin de unos factores dentro de otros no se est asumiendo o a los efectos aleatorios de todos los factores incluidos en la frmula como en el o caso de la funcin lme(). Para el ejemplo anterior el modelo quedar o a especicado de la siguiente forma: lme4(y 1 + (1|a) + (1|a:b) + (1 | [Link])) A lo largo de las prximas secciones se vern formulaciones ms complejas o a a para la componente aleatoria del modelo mixto.

1.2.
1.2.1.

Entendiendo los modelos mixtos

Zuur et al. (2007) utilizan datos de bentos marino procedente de nueve zonas intermareales de la costa holandesa para presentar los modelos mixtos. Los datos fueron recogidos por el instituto holands RIKZ en el verano de 2002. En e cada zona intermareal (playa) se tomaron cinco muestras de la macro-fauna y variables abiticas siguiendo el siguiente esquema de muestreo. o

134

Luis Cayuela

Modelos lineales mixtos en R

En este ejemplo vamos a ver si existe alguna relacin entre la riqueza de o especies y el NAP, una variable que indica la altura de cada estacin de o muestreo con respecto al nivel medio de la marea. Como la riqueza de especies es un conteo, ser ms apropiado utilizar un modelo lineal generalizado a a (GLM) con una distribucin de errores de tipo Poisson. Sin embargo, para o simplicar las cosas utilizaremos un modelo de regresin lineal con una o distribucin de errores normal o Gausiana. o Vamos a leer los datos. Estos vienen en la forma de una matriz de abundancias de cada una de las especies observadas (columnas 2 a la 76). Para calcular la riqueza tenemos que convertir estas abundancias a presencia y sumar todas ellas para cada una de las las (muestras): > RIKZ <- [Link](file = "[Link] header = TRUE) > RIKZ$Richness <- apply(RIKZ[, 2:76] > 0, 1, sum) O alternativamente podr amos usar > RIKZ$Richness <- rowSums(RIKZ[, 2:76] > 0)

1.2.2.

Una aproximacin al anlisis de los datos con modelos o a lineales

Sin saber de la existencia de los modelos mixtos, tal vez nuestra primera aproximacin ser ajustar un modelo lineal asumiendo independencia de las o a muestras, de acuerdo al siguiente modelo: Riquezai = + N APi + i Dicho modelo originar la siguiente grca: a a i N (0, 2 )

135

Luis Cayuela

Modelos lineales mixtos en R

> tmp <- lm(Richness ~ NAP, data = RIKZ) > plot(RIKZ$NAP, RIKZ$Richness, xlab = "NAP", ylab = "Richness") > abline(tmp)

20

q q

15 Richness
q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q

10

1.0

0.5

0.0

0.5 NAP

1.0

1.5

2.0

Y contiene tres parmetros desconocidos: los dos parmetros de la regresin a a o (constante y pendiente) y la varianza residual ( 2 ). El modelo asume que la relacin entre riqueza y NAP es la misma en todas las playas. Sin embargo o podr ocurrir que dicho modelo no es el mismo en todas las playas, en cuyo a cayo tendr amos tres opciones: (1) que la pendiente fuera la misma pero la constante (Intercept) no lo fuera; (2) que la constante fuera la misma pero que la pendiente no lo fuera; y (3) que ni la constante ni la pendiente fueran las mismas para todas las playas. En el primer caso, el modelo vendr dado por la siguiente ecuacin: a o Riquezaij = j + N APij + ij ij N (0, 2 )

Y su grca correspondiente, en dnde tendr a o amos nueva rectas de regresin o (una por cada playa), todas con la misma pendiente pero distinta constante (Intercept):

136

Luis Cayuela

Modelos lineales mixtos en R

> > > > + + + + + +

tmp1 <- lm(Richness ~ NAP + factor(Beach), data = RIKZ) plot(RIKZ$NAP, RIKZ$Richness, xlab = "NAP", ylab = "Richness") abline(v = 0, lty = 3) for (i in 1:9) { J <- RIKZ$Beach == i x1 <- RIKZ$NAP[J] y1 <- tmp1$fitted[J] Ord <- order(x1) lines(x1[Ord], y1[Ord]) }

20

q q

15 Richness
q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q

10

1.0

0.5

0.0

0.5 NAP

1.0

1.5

2.0

En el segundo caso, el modelo vendr dado por la siguiente ecuacin: a o Riquezaij = + j N APij + ij ij N (0, 2 )

Y su grca correspondiente, en dnde tendr a o amos nueva rectas de regresin o (una por cada playa), todas con la misma constante pero distinta pendiente:

137

Luis Cayuela

Modelos lineales mixtos en R

> > > > + + + + + +

tmp2 <- lm(Richness ~ NAP + factor(Beach):NAP, data = RIKZ) plot(RIKZ$NAP, RIKZ$Richness, xlab = "NAP", ylab = "Richness") abline(v = 0, lty = 3) for (i in 1:9) { J <- RIKZ$Beach == i x1 <- RIKZ$NAP[J] y1 <- tmp2$fitted[J] Ord <- order(x1) lines(x1[Ord], y1[Ord]) }

20

q q

15 Richness
q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q

10

1.0

0.5

0.0

0.5 NAP

1.0

1.5

2.0

En el tercer caso, el modelo vendr dado por la siguiente ecuacin: a o Riquezaij = j + j N APij + ij ij N (0, 2 )

Y su grca correspondiente, en dnde tendr a o amos nueva rectas de regresin o (una por cada playa), todas con distinta constante y pendiente:

138

Luis Cayuela

Modelos lineales mixtos en R

> > > > + + + + + +

tmp3 <- lm(Richness ~ NAP * factor(Beach), data = RIKZ) plot(RIKZ$NAP, RIKZ$Richness, xlab = "NAP", ylab = "Richness") abline(v = 0, lty = 3) for (i in 1:9) { J <- RIKZ$Beach == i x1 <- RIKZ$NAP[J] y1 <- tmp3$fitted[J] Ord <- order(x1) lines(x1[Ord], y1[Ord]) }

20

q q

15 Richness
q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q

10

1.0

0.5

0.0

0.5 NAP

1.0

1.5

2.0

El modelo con una unica recta de regresin es sin duda el ms sencillo, y el o a modelo con nueve rectas de regresin con diferentes constantes y pendientes o para cada una de las nueve playas el ms complejo. El nmero de parmetros a u a en el primer modelo es de 3, en los dos modelos intermedios de 1 + 9 + 1 = 11, y en el ultimo modelo, el ms complejo, de 9 + 9 + 1 = 19. a Comparemos ahora los distintos modelos utilizando el estad stico F con la funcin anova(). o > anova(tmp, tmp1, test = "F") Analysis of Variance Table Model 1: Richness ~ NAP Model 2: Richness ~ NAP + factor(Beach) 139

Luis Cayuela

Modelos lineales mixtos en R

[Link] RSS Df Sum of Sq F Pr(>F) 1 43 744.12 2 35 327.74 8 416.37 5.5581 0.0001441 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 > anova(tmp, tmp2, test = "F") Analysis of Variance Table Model 1: Model 2: [Link] 1 43 2 35 Richness ~ NAP Richness ~ NAP + factor(Beach):NAP RSS Df Sum of Sq F Pr(>F) 744.12 699.28 8 44.831 0.2805 0.9681

> anova(tmp, tmp3, test = "F") Analysis of Variance Table Model 1: Richness ~ NAP Model 2: Richness ~ NAP * factor(Beach) [Link] RSS Df Sum of Sq F Pr(>F) 1 43 744.12 2 27 165.92 16 578.19 5.8804 2.993e-05 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 > anova(tmp1, tmp3, test = "F") Analysis of Variance Table Model 1: Richness ~ NAP + factor(Beach) Model 2: Richness ~ NAP * factor(Beach) [Link] RSS Df Sum of Sq F Pr(>F) 1 35 327.74 2 27 165.92 8 161.82 3.2915 0.009434 ** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Claramente el modelo ms complejo es el mejor modelo y el que explica mayor a cantidad de variabilidad (se minimiza la suma de cuadrados residual). Ahora bien, estamos interesados en la relacin general entre riqueza y NAP y no nos o importa tanto el efecto ocasionado por las caracter sticas particulares de cada playa. Pero si eliminamos la playa del anlisis entonces las diferencias en la a riqueza de especies que son atribuibles a las caracter sticas particulares de cada playa pasar a formar parte del trmino de la varianza residual. Adems, a e a los errores no ser totalmente independientes entre s No tener esto en an . 140

Luis Cayuela

Modelos lineales mixtos en R

consideracin podr afectar a la estimacin de los errores estndar y los o a o a p-valores de los efectos jos. Esto podr producir como resultado, por a ejemplo, la no deteccin de una relacin signicativa entre la riqueza y el NAP. o o Por tanto hace falta incorporar el efecto de la playa en el modelo. Pero esto implica incorporar 16 parmetros ms en el modelo, lo que conlleva la prdida a a e de 16 grados de libertad! Una alternativa posible para evitar sto es utilizar los e modelos mixtos. Pero el uso de los modelos mixtos ofrece, adems, otras a ventajas. Como vimos en la sesin anterior, si consideramos la playa como un o factor jo, nuestras conclusiones slo pueden referirse a esas playas concretas. o Sin embargo, si tratamos la playa como un factor aleatorio, entonces nuestras conclusiones se referirn a todas las playas, de las cules, estas nueve son slo a a o una muestra del total de la poblacin de playas. Por tanto los resultados del o anlisis al considerar la playa como un factor aleatorio permitir predecir la a a relacin entre la riqueza y el NAP en un sentido ms general, sin tener que o a limitarnos a estas nueve playas en concreto. 1.2.3. Re-analizando los datos con modelos lineales mixtos: El paquete nlme

Vamos a comenzar re-analizando el primero de los tres modelos anteriores, que consideraba la existencia de diferencias en la constante de la recta de regresin o entre la riqueza y el NAP para las distintas playas. Si consideramos la playa como un factor aleatorio, el modelo quedar formulado de la siguiente forma: a Riquezaij = + N APij + aj + ij
2 aj N (0, a )

ij N (0, 2 )

El ndice j (representando a las playas) toma valores de 1 a 9, e i (representando las muestras dentro de cada playa) toma valores de 1 a 5. En el modelo lineal anterior, acabbamos con nueve l a neas de regresin. La pendiente o estimada, los nueve puntos de corte y sus errores estndares nos dec cmo a an o era de diferente la relacin entre la riqueza y el NAP para las distintas playas o (sin que cambiase la pendiente). En este modelo, asumimos que slo hay una o l nea de regresin con una unica constante y una unica pendiente. La constante o y la pendiente son los parmetros jos del modelo. Adems, hay una a a constante aleatoria, aj , que aade cierta cantidad de variacin a la constante n o del modelo general en cada una de las playas. Se asume que esta constante 2 aleatoria sigue una distribucin normal de media 0 y varianza a . Por lo tanto, o los parmetros estimados en el modelo son cuatro, , , la varianza residual a 2 2 , y la varianza de la constante a . De esta forma tenemos un modelo equivalente al modelo lineal, pero en dnde slo es necesario estimar cuatro o o parmetros y no once, ya que en vez de tener nueve estimaciones de las a constantes en cada una de las rectas de regresin de las nueve playas, tenemos o nueve valores no estimados 1 , . . . , 9 que asumimos siguen una distribucin o normal. Lo que estimamos en este modelo es la varianza de esta distribucin. o Vamos a ajustar el modelo lineal mixto correspondiente con la funcin lme(): o

141

Luis Cayuela

Modelos lineales mixtos en R

> library(nlme) > lme5 <- lme(Richness ~ NAP, data = RIKZ, random = ~1 | factor(Beach)) > summary(lme5) Linear mixed-effects model fit by REML Data: RIKZ AIC BIC logLik 247.4802 254.5250 -119.7401 Random effects: Formula: ~1 | factor(Beach) (Intercept) Residual StdDev: 2.944065 3.05977 Fixed effects: Richness ~ NAP Value [Link] DF t-value p-value (Intercept) 6.581893 1.0957618 35 6.006682 0 NAP -2.568400 0.4947246 35 -5.191574 0 Correlation: (Intr) NAP -0.157 Standardized Within-Group Residuals: Min Q1 Med Q3 -1.4227495 -0.4848006 -0.1576462 0.2518966 Number of Observations: 45 Number of Groups: 9 > anova(lme5) numDF denDF F-value p-value 1 35 27.63494 <.0001 1 35 26.95244 <.0001

Max 3.9793918

(Intercept) NAP

La parte de los resultados que se reere a los efectos aleatorios nos muestra que la varianza residual es 2 = 3,062 y la varianza de la constante 2 a = 2,9442 . Para la parte de los efectos jos del modelo, + N APij , la constante se estima en = 6,582 y la pendiente en = 2,568. Ambos parmetros son signicativamente distintos de 0. La correlacin entre la a o constante y la pendiente es pequea. La tabla ANOVA tambin muestra la n e signicacin de la pendiente . o En resumen, el modelo estima una respuesta general de la forma 6,582 2,568 N APij . Estos parmetros estimados por el modelo son a signicativamente distintos de 0. Para cada playa, la constante aumenta o disminuye por un valor aleatorio. Este valor aleatorio sigue una distribucin o 2 normal con media esperada 0 y varianza a = 2,9442 . La varianza residual, es decir, el error que podemos aadir a cada una de nuestras predicciones, viene n dada por 2 = 3,06. 142

Luis Cayuela

Modelos lineales mixtos en R

1.2.4.

Re-analizando los datos con modelos lineales mixtos: El paquete lme4

El ajuste del modelo con la funcin lmer() del paquete lme4 es prcticamente o a idntico a cmo lo hemos hecho utilizando la funcin lme() del paquete nlme. e o o La principal diferencia es que no hace falta especicar el argumento random y los efectos aleatorios se especican como parte de la frmula, de la siguiente o forma: > library(lme4) > lmer5 <- lmer(Richness ~ NAP + (1 | Beach), data = RIKZ) > summary(lmer5) Linear mixed model fit by REML Formula: Richness ~ NAP + (1 | Beach) Data: RIKZ AIC BIC logLik deviance REMLdev 247.5 254.7 -119.7 241.9 239.5 Random effects: Groups Name Variance [Link]. Beach (Intercept) 8.6675 2.9441 Residual 9.3622 3.0598 Number of obs: 45, groups: Beach, 9 Fixed effects: Estimate Std. Error t value (Intercept) 6.5819 1.0957 6.007 NAP -2.5684 0.4947 -5.192 Correlation of Fixed Effects: (Intr) NAP -0.157 > anova(lmer5) Analysis of Variance Table Df Sum Sq Mean Sq F value NAP 1 252.33 252.33 26.952 Como vemos, los valores de los parmetros aleatorios y jos estimados por el a modelo son exactamente iguales utilizando cualquiera de las dos funciones. En el caso de la funcin lmer(), y al contrario de lo que ocurre con la funcin o o lme(), no se dan los valores de signicacin de los coecientes jos estimados o por el modelo. Buscando en la lista de distribucin de R, uno acabar o a encontrndose con un debate casi losco sobre la capacidad de estos a o modelos de realizar tests de signicacin utilizando el estad o stico F. En concreto, Douglas Bates, autor principal de este paquete, cree que no es posible estimar con abilidad los grados de libertad de los parmeros jos y, a 143

Luis Cayuela

Modelos lineales mixtos en R

por tanto, proponer una distribucin del estad o stico F lo sucientemente able como realizar un test de signicacin. En este contexto la signicacin de los o o parmetros jos en el modelo puede hacerse comparando modelos ms a a complejos frente a modelos menos complejos utilizando un criterio de informacin (AIC, BIC), siempre y cuando los efectos aleatorios no cambien. o En este ejemplo, el test de signicacin vendr dado de la siguiente forma: o a > lmer5.0 <- lmer(Richness ~ NAP + (1 | Beach), data = RIKZ, REML = F) > lmer5.1 <- lmer(Richness ~ 1 + (1 | Beach), data = RIKZ, REML = F) > anova(lmer5.0, lmer5.1) Data: RIKZ Models: lmer5.1: Richness lmer5.0: Richness Df AIC lmer5.1 3 269.30 lmer5.0 4 249.83 --Signif. codes: 0

~ 1 + (1 | Beach) ~ NAP + (1 | Beach) BIC logLik Chisq Chi Df Pr(>Chisq) 274.72 -131.65 257.06 -120.92 21.474 1 3.586e-06 *** *** 0.001 ** 0.01 * 0.05 . 0.1 1

La razn por la que se utiliza el argumento REML = F se explica en la seccin o o 1.3. Como vemos, el modelo completo (lmer5.0) produce un AIC que es inferior al AIC del modelo anidado (que no contempla la variacin en las pendientes), o y adems el test Chi-cuadrado realizado sobre el parmetro de verosimilitud a a (logLik) es signicativo, lo que indica que que el modelo completo es ms a parsimonioso. Todo ello apunta a que existe un efecto general y signicativo del NAP sobre la riqueza de especies, resultado coincidente con el obtenido anteriormente con la funcin lme(). o Y cmo comprobar la signicacin de las variables aleatorias? En el caso de o o las funciones lme() y lmer(), no se calculan p-valores para los trminos e aleatorios. La idea ser bsicamente la misma que para el caso de los efectos a a jos con la funcin lmer(): habr que comparar modelos alternativos con una o a estructura ja idntica pero que variasen en sus efectos aleatorios. Sin embargo e la estimacin de los parmetros de estos modelos ha de hacerse con el o a estimador de mxima verosimilitud restringida (REML) como se explica ms a a adelante. 1.2.5. Otros efectos aleatorios en modelos lineales mixtos

Siguiendo con el ejemplo anterior, es posible que contemplemos la posibilidad de que haya, adems de distintos puntos de corte, distintas pendientes para a cada una de las playas. En este caso el modelo quedar formulado de las a siguiente forma: Riquezaij = + aj + N APij + bj N APij + ij

144

Luis Cayuela

Modelos lineales mixtos en R

ij N (0, 2 )

2 aj N (0, a )

2 bj N (0, b )

El modelo es exactamente igual que el modelo anterior excepto por el trmino e bj N APij . Este nuevo trmino permite la variacin aleatoria de la pendiente e o en cada una de las playas. El modelo es muy parecido al modelo lineal completo con interaccin que ajustamos en la seccin 1.2.2, pero con muchos o o menos parmetros: dos para los efectos jos y , y tres para las varianzas de a 2 2 los trminos aleatorios 2 , a y b . En general, se permite que haya una e 2 2 correlacin alta entre las varianzas estimadas a y b . o Vamos a ajustar el modelo lineal mixto correspondiente (los dos modelos lme6 son equivalentes): > > > + > library(nlme) lme6 <- lme(Richness ~ NAP, data = RIKZ, random = ~NAP | factor(Beach)) lme6 <- lme(Richness ~ NAP, data = RIKZ, random = ~1 + NAP | factor(Beach)) summary(lme6)

Linear mixed-effects model fit by REML Data: RIKZ AIC BIC logLik 244.3839 254.9511 -116.1919 Random effects: Formula: ~1 + NAP | factor(Beach) Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 3.549068 (Intr) NAP 1.714957 -0.99 Residual 2.702822 Fixed effects: Richness ~ NAP Value [Link] DF t-value p-value (Intercept) 6.588704 1.2647621 35 5.209442 0e+00 NAP -2.830027 0.7229384 35 -3.914616 4e-04 Correlation: (Intr) NAP -0.819 Standardized Within-Group Residuals: Min Q1 Med Q3 -1.8213286 -0.3411047 -0.1674620 0.1921267 Number of Observations: 45 Number of Groups: 9 > anova(lme6)

Max 3.0397121

145

Luis Cayuela

Modelos lineales mixtos en R

(Intercept) NAP

numDF denDF F-value p-value 1 35 12.19512 0.0013 1 35 15.32422 0.0004

La parte de los resultados que se reere a los efectos aleatorios nos muestra que 2 la varianza residual es 2 = 2,702 , la varianza de la constante a = 3,552 y la 2 2 varianza de la pendiente b = 1,71 . Vemos que disminuye la varianza residual (lo que no explica el modelo) con respecto al modelo anterior. Para la parte de los efectos jos del modelo, + N APij , la constante estimada es = 6,588 y la pendiente = 2,83. Ambos parmetros son signicativamente distintos a de 0 y muy similares a los coecientes estimados en el modelo anterior.

1.3.

Seleccin de modelos o

Tenemos dos modelos mixtos: uno en dnde la constante var aleatoriamente o a (lme5) y otro en dnde tanto la constante y la pendiente var aleatoriamente o an (lme6). Qu modelo es mejor o ms apropiado? Hay que notar que ambos e a modelos son idnticos en cuanto a sus efectos jos y slo var en sus efectos e o an aleatorios. Para comparar los dos modelos con los mismos efectos jos pero con distintos efectos aleatorios, podemos usar un test de verosimilitud. > anova(lme5, lme6) Model df AIC BIC logLik Test [Link] p-value 1 4 247.4802 254.5250 -119.7401 2 6 244.3839 254.9511 -116.1919 1 vs 2 7.096378 0.0288

lme5 lme6

El AIC sugiere que el modelo lme6 es ms apropiado, pero el BIC sugiere que a el modelo lme5 es mejor. El p-valor del test de verosimilitud indica que el modelo ms complejo (lme6) es el ms parsimonioso y por tanto elegir a a amos ste. Sin embargo, hay un problema con el test de verosimilitud. En modelos e lineales se puede utilizar un test F para comparar modelos anidados, en dnde o la suma de cuadrados residual de los dos modelos obtenida por m nimos cuadrados es utilizada para calcular el estad stico de contraste F. En este caso, podemos comprobar hiptesis como H0 : 0 = 0 frente a H1 : 0 = 0. En los o modelos lineales mixtos el estimador utilizado no es m nimos cuadrados sino mxima verosimilitud (ML). El criterio de verosimilitud del modelo completo a L0 y el modelo anidado L1 son utilizados para comprobar hiptesis referentes o a la signicacin de los parmetros en el modelo. Tomando el logaritmo o, ms o a a comnmente, 2 log, se calcula un estad u stico de la forma L = 2(logL0 logL1 ). Se puede demostrar que bajo la hiptesis nula, este o estad stico sigue aproximadamente una distribucin Chi-cuadrado con o grados de libertad, en dnde es la diferencia entre el nmero de parmetros o u a de ambos modelos. Cul ser la hiptesis nula en este caso? La unica a a o diferencia entre el modelo lme6 y el modelo lme5 es el componente aleatorio 2 bj N APij , dnde bj N (0, b ). Si comparamos ambos modelos con un test de o 2 verosimilitud, estar amos comprobando la hiptesis nula de que H0 : b = 0 o 2 frente a H1 : b > 0. La hiptesis alternativa contiene un > porque se supone o 146

Luis Cayuela

Modelos lineales mixtos en R

que los componentes de la varianza no pueden tomar valores negativos. Esto es lo que se conoce como el problema del l mite (boundary problem); estamos comprobando si la varianza es cero, pero este valor est en el l a mite de todos los posibles valores que puede tomar la varianza. El problema es que la teor a subyacente al test de verosimilitud, que es la que genera un p-valor, asume que no hay un l mite en el espacio del parmetro estimado. En resumen, hay que a tener mucho cuidado cuando se interpreta el p-valor derivado de un test de verosimilitud cuando estamos comprobando la hiptesis nula en el l o mite de la distribucin del parmetro en cuestin. Si el p-valor es muy grande o muy o a o pequeo, no suele haber problema en interpretarlo, pero cuando encontramos n un p-valor en el borde de la signicacin hay que tener cuidado con cmo lo o o interpretamos. Pinheiro & Bates (2000) y Faraway (2006) argumentan que, en este sentido, las tcnicas de bootstraping son ms conables a la hora de e a estimar los p-valores si estamos comprobando una hiptesis en el l o mite de su distribucin. Otra alternativa es hacer la estimacin utilizando un estimador o o de mxima verosimilitud restringida (REML), que es lo que por defecto hacen a las funciones lme() y lmer(), que permite corregir en parte el problema de la estimacin de la varianza en el l o mite de su distribucin. o De lo visto hasta el momento concluimos que hay dos formas de estimar los parmetros de los modelos, mediante mxima verosimilitud (ML) y mediante a a mxima verosimilitud restringida (REML). Es importante mencionar que no es a posible comparar un modelo estimado mediante ML con otro modelo estimado mediante REML. En los modelos lineales mixtos, tenemos dos tipos de efectos, los efectos jos y los aleatorios. Aunque generalmente vamos a estar ms interesados en los a efectos jos, si tenemos una estructura de efectos aleatorios mal denida es posible que eso afecte a la estimacin de los efectos jos. Por ello es o importante seleccionar la mejor estructura para cada uno de estas dos componentes. Para ello se recomienda lo siguiente (Zuur et al. 2009): 1. Empieza ajustando un modelo en dnde la componente ja contenga o todas las variables explicativas e interacciones posibles. Esto es lo que se conoce como modelo ms all del ptimo (beyond optimal model ). Si no a a o es posible ajustar este modelo porque hay demasiados parmetros, hay a que intentar seleccionar las variables e interacciones que creamos inuyen ms sobre la variable respuesta. a 2. Utilizando el modelo ms all del ptimo, hay que encontrar la a a o estructura de la componente aleatoria ptima. Para ello propondremos o distintos modelos alternativos con la misma estructura en la componente ja pero que var en su componente aleatoria. Porque en la an componente aleatoria estamos estimando varianzas nos encontramos con el problema de la estimacin en el l o mite de la distribucin de los o parmetros estimados mencionado anteriormente. Por ello estos modelos a tienen que estar estimados utilizando REML. La comparacin entre o distintos modelos podemos hacerla utilizando la funcin anova(). o 3. Una vez que hemos denido la estructura ptima de la componente o aleatoria, tenemos que buscar la estructura ptima de la componente ja o

147

Luis Cayuela

Modelos lineales mixtos en R

del modelo. Para ello podemos utilizar el estad stico F o el estad stico t obtenido mediante el estimador REML con la funcin lme()1 o comparar o modelos anidados. Para comparar modelos que tienen la misma estructura en la componente aleatoria pero dieren en la componente ja se debe de utilizar un estimador LM y no un estimador de REML. 4. Cuando se ha seleccionado la estructura de la componente ja, se presenta el modelo nal utilizando un estimador REML.

1.4.

Anlisis de los residuos a

Para saber si el modelo es adecuado debemos mirar los residuos normalizados del modelo estimado a partir de REML. Al contrario que en el caso de los modelos lineales (lm) o modelos lineales generalizados (glm) no tenemos una funcin que nos genera todos los grcos de los residuos, as que tendremos que o a generarlos nosotros individualmente. Interesa dibujar los residuos frente a los valores predichos para ver si hay homocedasticidad. Si no hay homogeneidad de varianzas entonces se puede: (i) aplicar una transformacin; (ii) intentar o averiguar si el incremento de la varianza para determinados valores es debido a alguna covariable; (iii) utilizar otro tipo de modelos como los modelos lineales generalizados mixtos con una distribucin de errores diferente (por ejemplo o Poisson si la variable respuesta es un conteo)2 . Es interesante dibujar los grcos de los residuos frente a las variables explicativas. De nuevo, no tenemos a que observar ningn patrn en estas grcas. Si la variable explicativa es un u o a factor, la varianza ha de ser homognea entre los distintos niveles del factor. e Tambin se puede sacar un histograma de los residuos para ver si hay e normalidad. Dicha normalidad tambin se puede comprobar con tests e estad sticos como el test de Shapiro-Wilk. Finalmente podemos obtener el qqplot de los residuos con la funcin qqnorm(). o
1 Recordemos 2 Esto

que la funcin lmer() del paquete lme4 no hace estos contrastes de hiptesis. o o lo podemos hacer con la funcin glmer() del paquete lme4. o

148

Luis Cayuela

Modelos lineales mixtos en R

> > > > > > > > > >

Res <- residuals(lme6, type = "normalized") Fit <- fitted(lme6) par(mfrow = c(2, 2)) plot(Res ~ Fit, xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs. fitted" abline(h = 0) plot(Res ~ RIKZ$NAP, xlab = "NAP", ylab = "Residuals", main = "NAP") abline(h = 0) hist(Res, main = "Histogram of residuals", xlab = "Residuals") qqnorm(Res) qqline(Res)
Residuals vs. fitted
3
q

NAP
3
q q

Residuals

Residuals

q q

2 1

2 1

q q q q q q q q q qqq qq q q q qq q q q q q q q q q q q q

q q

q q q q q q qq q q q qq q q q q q q q qq qq q q q q q q q q q q q q q q

10 Fitted values

15

20

1.0

0.0 NAP

1.0

2.0

Histogram of residuals
20 3

Normal QQ Plot
q q

Sample Quantiles

Frequency

15

qq q q q q q qq qq qq q qq q qqq q q qq q qqq qq qq qq qq q q qqq

10

2 1

q q

0 2

Residuals

Theoretical Quantiles

2.

Dise o por bloques n

Los datos del RIKZ son un claro ejemplo de diseo por bloques. En este caso n no tenemos un diseo por bloques aleatorizados puesto que dentro de los n bloques no investigamos la respuesta de una variable (riqueza) frente a un factor, sino frente a una covariable (NAP). En este caso tendr amos por tanto un bloque que es consecuencia de la agrupacin de las unidades muestrales o dentro de un factor de agrupacin que son las playas (recordemos que se o tomaron cinco muestras por playa). Por tanto, no es del inters del e investigador ver si la playa tiene un efecto sobre la riqueza de especies, pero no se puede obviar este factor puesto que tiene una inuencia sobre la varianza de la variable respuesta y, por tanto, va a tener un efecto probable en la estimacin de los parmetros de inters en el modelo (NAP), adems de la o a e a

149

Luis Cayuela

Modelos lineales mixtos en R

clara violacin del supuesto hecho sobre el diseo experimental de que las o n muestras son independientes. Vamos a utilizar ahora otro ejemplo distinto. En este estudio un tcnico e agr cola investiga el efecto de seis tipos de dieta distintos (a,..., f) en la ganancia de peso de conejos domsticos. Para ello selecciona 3 conejos ms o e a menos uniformes en cuanto a tamao y peso procedentes de 10 camadas n distintas. En total hay por tanto 30 muestras. Como hay seis tratamientos y slo 3 conejos por camada, no se puede aplicar todos los tratamientos a todos o los conejos de cada camada. El diseo es, por tanto, un diseo factorial por n n bloques aleatorizados incompletos. Vamos a leer los datos: > library(faraway) > data(rabbit) > xtabs(gain ~ treat + block, data = rabbit) block treat b1 b10 b2 b3 b4 b5 b6 b7 b8 b9 a 0.0 37.3 40.1 0.0 44.9 0.0 0.0 45.2 44.0 0.0 b 32.6 0.0 38.1 0.0 0.0 0.0 37.3 40.6 0.0 30.6 c 35.2 0.0 40.9 34.6 43.9 40.9 0.0 0.0 0.0 0.0 d 0.0 42.3 0.0 37.5 0.0 37.3 0.0 37.9 0.0 27.5 e 0.0 0.0 0.0 0.0 40.8 32.0 40.5 0.0 38.5 20.6 f 42.2 41.7 0.0 34.3 0.0 0.0 42.8 0.0 51.9 0.0

2.1.

Paso 1: Aplicacin de un modelo lineal o

Como ya hemos visto, una forma de analizar el efecto de un bloque es considerando ste como un efecto jo por medio de un modelo lineal. Es e importante tener en cuenta que a la hora de formular el modelo el bloque tiene que ir en primer lugar. Esto es porque los modelos lineales en R se ajustan utilizando una suma de cuadrados de tipo I, que tiene en cuenta el orden de entrada de las variables en el modelo, y nos interesa analizar el efecto del factor (dieta) sobre la respuesta (ganancia de peso) una vez que se haya tenido en cuenta la variabilidad atribuible al bloque (camada). Tambin es e importante observar que el diseo no es cruzado, es decir que no todos los n niveles del factor estn representados dentro de cada uno de los niveles del a bloque. Por tanto, plantear una interaccin en este modelo no es posible. o > [Link] <- lm(gain ~ block + treat, data = rabbit) > anova([Link]) Analysis of Variance Table Response: gain Df Sum Sq Mean Sq F value 150

Pr(>F)

Luis Cayuela

Modelos lineales mixtos en R

block 9 730.39 81.154 8.0738 0.0002454 *** treat 5 158.73 31.745 3.1583 0.0381655 * Residuals 15 150.77 10.052 --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Observamos que tanto el bloque como el tratamiento son signicativos.

2.2.

Paso 2: Ajuste de un modelo lineal mixto

Aunque el modelo lineal no es necesariamente inadecuado, la inclusin del o bloque como un factor jo se hace a expensas de 9 grados de libertad. En este sentido el modelo mixto equivalente ser mucho ms parsimonioso en trmino a a e del nmero de parmetros. u a > library(nlme) > lme.rabbit1 <- lme(gain ~ treat, random = ~1 | block, data = rabbit) > anova(lme.rabbit1) numDF denDF F-value p-value 1 15 590.5297 <.0001 5 15 3.2818 0.0336

(Intercept) treat

> summary(lme.rabbit1) Linear mixed-effects model fit by REML Data: rabbit AIC BIC logLik 166.3569 175.7813 -75.17844 Random effects: Formula: ~1 | block (Intercept) Residual StdDev: 4.657853 3.175519 Fixed effects: gain ~ treat Value [Link] DF t-value p-value (Intercept) 39.53540 2.130338 15 18.558278 0.0000 treatb -2.50718 2.208700 15 -1.135139 0.2741 treatc -0.18407 2.208700 15 -0.083340 0.9347 treatd -0.88516 2.208700 15 -0.400759 0.6942 treate -5.64602 2.208700 15 -2.556263 0.0219 treatf 2.81003 2.208700 15 1.272254 0.2227 Correlation: (Intr) treatb treatc treatd treate treatb -0.518 treatc -0.518 0.500 treatd -0.518 0.500 0.500 151

Luis Cayuela

Modelos lineales mixtos en R

treate -0.518 treatf -0.518

0.500 0.500

0.500 0.500

0.500 0.500

0.500

Standardized Within-Group Residuals: Min Q1 Med Q3 -1.3794203 -0.5213453 -0.1701517 0.6456859 Number of Observations: 30 Number of Groups: 10

Max 1.4148990

2.3.

Paso 3: Eleccin de la estructura de los efectos o aleatorios

Recordemos que para elegir la estructura de los efectos aleatorios es necesario incluir todos los posibles trminos jos y sus interacciones en el modelo (ms e a all del ptimo). Luego se comparan distintos modelos que var en sus a o an efectos aleatorios pero que mantienen la misma estructura ja por medio de REML. Los modelos que podemos plantear en este caso son: (1) un modelo sin efecto aleatorio (equivalente a un modelo lineal pero utilizando un estimador REML). Este modelo lo calcularemos con la funcin gls(); (2) un modelo que o incluya la varianza aleatoria de la constante (gran media) dentro de cada bloque. > lme.rabbit0 <- gls(gain ~ treat, data = rabbit) > lme.rabbit1 <- lme(gain ~ treat, data = rabbit, random = ~1 | + block) > anova(lme.rabbit0, lme.rabbit1) Model df AIC BIC logLik Test [Link] p-value 1 7 174.2621 182.5085 -80.13107 2 8 166.3569 175.7813 -75.17844 1 vs 2 9.905255 0.0016

lme.rabbit0 lme.rabbit1

Vemos que el modelo con el efecto aleatorio bloque es ms parsimonioso que el a modelo sin el efecto aleatorio.

2.4.

Paso 4: Eleccin de la estructura de los efectos jos o

Ya tenemos la estructura de los efectos aleattorios. Podemos utilizar el test t o el test F para comprobar la signicacin de la(s) variable(s) ja(s). Pero tal o vez es mejor opcin comparar modelos anidados como en el caso de los efectos o aleatorios. Recordemos que, para hacer esto, debemos estimar los parmetros a de los modelos utilizando un estimador de ML. > lme.rabbit2 <- lme(gain ~ 1, data = rabbit, random = ~1 | block, + method = "ML") > lme.rabbit3 <- lme(gain ~ treat, data = rabbit, random = ~1 | + block, method = "ML") > anova(lme.rabbit2, lme.rabbit3) 152

Luis Cayuela

Modelos lineales mixtos en R

lme.rabbit2 lme.rabbit3

Model df AIC BIC logLik Test [Link] p-value 1 3 188.8307 193.0343 -91.41536 2 8 183.7730 194.9825 -83.88648 1 vs 2 15.05776 0.0101

Y, de nuevo, vemos que el modelo con el efecto jo del tratamiento es ms a parsimonioso que el modelo que contiene unicamente el efecto de la constante.

2.5.

Paso 5: Presentacin del modelo nal con REML o

Finalmente, presentamos el modelo nal utilizando REML. Es necesario que veamos cmo es de adecuado el modelo de acuerdo con los supuestos del o modelo, fundamentalmente normalidad y homocedasticidad. > [Link] <- lme(gain ~ treat, data = rabbit, random = ~1 | + block) > anova([Link]) numDF denDF F-value p-value 1 15 590.5297 <.0001 5 15 3.2818 0.0336

(Intercept) treat

> summary([Link]) Linear mixed-effects model fit by REML Data: rabbit AIC BIC logLik 166.3569 175.7813 -75.17844 Random effects: Formula: ~1 | block (Intercept) Residual StdDev: 4.657853 3.175519 Fixed effects: gain ~ treat Value [Link] DF t-value p-value (Intercept) 39.53540 2.130338 15 18.558278 0.0000 treatb -2.50718 2.208700 15 -1.135139 0.2741 treatc -0.18407 2.208700 15 -0.083340 0.9347 treatd -0.88516 2.208700 15 -0.400759 0.6942 treate -5.64602 2.208700 15 -2.556263 0.0219 treatf 2.81003 2.208700 15 1.272254 0.2227 Correlation: (Intr) treatb treatc treatd treate treatb -0.518 treatc -0.518 0.500 treatd -0.518 0.500 0.500 treate -0.518 0.500 0.500 0.500 treatf -0.518 0.500 0.500 0.500 0.500 153

Luis Cayuela

Modelos lineales mixtos en R

Standardized Within-Group Residuals: Min Q1 Med Q3 -1.3794203 -0.5213453 -0.1701517 0.6456859 Number of Observations: 30 Number of Groups: 10

Max 1.4148990

El modelo nal ser del tipo gain = 39,535 5,646 T rate (ya que los a coecientes para el resto de los tratamientos no son signicativamente distintos de cero) con una varianza residual de 2 = 3,1752 y una varianza para la 2 constante de a = 4,6582 (este es el efecto estimado del bloque en nuestro 3 modelo) .
3 Una forma ms elegante de comparar qu niveles del factor son signicativos ser juntar a e a los distintos niveles de dietas que hemos asumido que no tienen un efecto diferenciado entre s , volver a ajustar el modelo y comparar este nuevo modelo con el modelo anterior. Esto tambin e podr solucionar algunos de los problemas asociados con la violacin de los supuestos, como a o la heterocedasticidad.

154

Luis Cayuela

Modelos lineales mixtos en R

> > > > > > > > > >

Res <- residuals([Link], type = "normalized") Fit <- fitted([Link]) par(mfrow = c(2, 2)) plot(Res ~ Fit, xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs. fitted" abline(h = 0) boxplot(Res ~ rabbit$treat, ylab = "Residuals", main = "Treatment") abline(h = 0, lty = 3) hist(Res, main = "Histogram of residuals", xlab = "Residuals") qqnorm(Res) qqline(Res)
Residuals vs. fitted
q q q q q q q q q q q q q q q q q q q

Treatment
1.0

1.0

Residuals

0.0

qq

1.0

q q q

25

30

35

40

45

1.0

q q

0.0

q q

Residuals

Fitted values

Histogram of residuals
8 1.0

Normal QQ Plot
q q

Sample Quantiles

Frequency

qq q qq q q q q q q q q q qq qq qq qq q qq q q q

1.5

0.5

0.5

1.5

1.0 2

0.0

Residuals

Theoretical Quantiles

Parece que hay algo de heterocedasticidad debido a que no hay igualdad de varianzas entre los grupos. Adems, podr haber tambin problemas de a a e normalidad. No parece que haya problemas de linealidad. Soluciones al problema? Ninguna necesariamente ptima: probar otro tipo de modelos como o los modelos lineales generalizados con una distribucin de errores diferente o (por ejemplo, gamma), los modelos aditivos generalizados o modelos lineales mixtos no lineales. Tambin podr e amos utilizar alguna tcnica no paramtrica e e como los rboles de regresin o, simplemente, quedarnos con este modelo a a o pesar de sus limitaciones.

155

Luis Cayuela

Modelos lineales mixtos en R

3.

Dise o de tipo split-plot n

Vamos a usar otro ejemplo para ilustrar el diseo de tipo split-plot. n Recordemos que este tipo de diseos se basa en el diseo por bloques n n aleatorizados, pero existiendo un segundo factor cuyos niveles son aplicados entre los bloques y no dentro de los bloques. El objetivo de este experimento es investigar si distintos escenarios de cambio climtico (incluyendo ms y menos a a lluvia estival) afectan a distintos parmetros ecosiolgicos de varias especies a o leosas t n picas del bosque mediterrneo, y si este efecto var dependiendo del a a tipo de hbitat (pastizal, matorral, bosque)4 . En este ejercicio nos centraremos a en una unica especie (Quercus ilex ) y una unica variable respuesta (biomasa). El diseo viene ilustrado de la siguiente forma: n

Pero cuidado! La parcela 1 de pastizal no tiene absolutamente nada que ver con la parcela 1 de matorral o bosque. Preparamos los datos: > sn <- [Link]("[Link] header = T) > xtabs([Link] ~ Plot + Scenario + Habitat, data = sn) , , Habitat = Forest Scenario
L., Zamora, R., Castro, J. Sporadic rainy events are more critical than drought intensication for woody recruitment in Mediterranean mountains: a eld experiment simulating climate change. Ecology (submitted). Los datos han sido cedidos por Luis Mat as (Departamento de Ecolog Universidad de Granada). No se permite el uso de estos datos a, excepto para nes docentes sin el permiso del autor (lmatias@[Link]).
4 Mat as,

156

Luis Cayuela

Modelos lineales mixtos en R

Plot Control Dry Summer Wet Summer 1 1.389 1.084 0.985 2 1.422 1.667 1.793 3 1.189 1.356 1.636 4 1.953 1.371 1.490 5 0.849 1.369 1.143 6 0.950 1.351 1.951 7 1.410 2.913 2.303 8 1.203 1.538 1.853 , , Habitat = Grassland Scenario Plot Control Dry Summer Wet Summer 1 0.000 2.331 4.908 2 2.596 3.454 3.611 3 1.534 2.765 3.517 4 6.777 1.143 2.741 5 1.627 2.849 2.860 6 1.662 2.233 2.724 7 0.000 2.005 1.841 8 1.238 1.689 4.084 , , Habitat = Shrubland Scenario Plot Control Dry Summer Wet Summer 1 1.991 2.150 2.357 2 2.056 1.845 2.354 3 1.533 2.151 1.781 4 1.881 1.558 2.218 5 1.718 1.655 2.086 6 1.117 2.323 1.865 7 1.641 1.475 1.391 8 0.986 2.003 1.322 > sn$Plot <- rep(1:24, each = 3) > sn <- sn[!sn$[Link] == 0, ] Eliminamos los valores de biomasa que valen 0 porque las semillas de estas muestras no han germinado y el problema aqu ser otro. Si el sujeto muestral a es la parcela, tenemos un factor de medidas repetidas (escenario) con tres niveles y un factor inter-sujetos (habitat) tambin con tres niveles. Si e slamente tuviramos en cuenta el factor escenario tendr o e amos un diseo de n medidas repetidas, pero como tenemos los dos tipos de factores (intra- e inter-sujetos) entonces el diseo se denomina de tipo split-plot. n

157

Luis Cayuela

Modelos lineales mixtos en R

3.1.

Paso 1: Aplicacin de un modelo lineal o

Antes de ajustar un modelo lineal vamos a obtener algunos grcos a exploratorios. Como podemos ver, hay diferencias en la biomasa total dentro de cada parcela (pero recordemos que en cada parcela se experimentan los tres niveles del factor escenario, as que es ms fcil ver este efecto en el trmino a a e interaccin porque aqu estamos promediando los valores de biomasa de o plantas sometidas a sequ lluvia estival y al control). Tambin parece que hay a, e un efecto positivo de la lluvia estival en la produccin de biomasa leosa, y o n que se produce ms biomasa en el pastizal que en los otros dos tipos de a hbitat. Los grcos de interaccin parecen indicar una interaccin entre el a a o o hbitat y el escenario. La respuesta de la biomasa al escenario parece ser a homogneo dentro de cada parcela excepto por una parcela en dnde e o detectamos una respuesta muy distinta. Habr que ver qu pasa con estos a e datos. Tal vez haya algo anormal o algn error de medicin o de muestreo, en u o cuyo caso deber amos desestimar este dato en el anlisis. Por el momento lo a dejaremos. No tiene sentido plantear una interaccin entre la parcela y el o hbitat porque no son factores cruzados sino anidados (es decir, no todos los a niveles del factor parcela estn representados dentro del factor hbitat. De a a hecho las parcelas 1 a 8 se encuentran en un tipo de hbitat, las parcelas 9 a a 16 en otro y las parcelas 17 a 24 en otro, as que la interaccin no es posible). o

158

Luis Cayuela

Modelos lineales mixtos en R

> > > > > >

par(mfcol = c(3, 2)) plot([Link] ~ factor(Plot), data = sn) plot([Link] ~ Scenario, data = sn) plot([Link] ~ Habitat, data = sn) [Link](sn$Scenario, sn$Habitat, sn$[Link]) [Link](sn$Scenario, factor(sn$Plot), sn$[Link])

sn$Habitat mean of sn$[Link] 3.0 6 Grassland Shrubland Forest

[Link]

1 3 5 7 9

12

15

18

21

24

1.5 Control

2.0

2.5

Dry Summer sn$Scenario

factor(Plot)

factor(sn$Plot) mean of sn$[Link] 6 6 9 16 10 11 13 12 14 17 18 7 20 21 6 22 8 15 2 19 3 4 23 24 5 1

[Link]

Control

Dry Summer Scenario

1 Control

Dry Summer sn$Scenario

Wet Summer

[Link]

Forest

Grassland Habitat

Shrubland

Vamos a analizar los datos con un modelo lineal. Como slo hay un dato por o cada parcela y escenario no es posible calcular la signicacin del trmino o e interaccin entre el factor y el bloque (no hay rplicas). Por tanto el modelo o e que plantearemos ser el siguiente: a > [Link] <- lm([Link] ~ factor(Plot) + Scenario * Habitat, + data = sn) 159

Luis Cayuela

Modelos lineales mixtos en R

> anova([Link]) Analysis of Variance Table Response: [Link] Df Sum Sq Mean Sq F value factor(Plot) 23 29.5922 1.28662 1.8075 Scenario 2 2.7137 1.35683 1.9061 Scenario:Habitat 4 2.3363 0.58408 0.8205 Residuals 40 28.4733 0.71183 --Signif. codes: 0 *** 0.001 ** 0.01 *

Pr(>F) 0.04932 * 0.16192 0.51981

0.05 . 0.1 1

El bloque es signicativo pero ni el escenario ni el hbitat lo son. Con este a modelo estamos gastando 23 grados de libertad con el bloque (parcela) cuando podr amos utilizar slamente uno estimando el efecto de la parcela como la o varianza de todos los posibles coecientes de los niveles de la parcela en su efecto sobre la constante (Intercept).

3.2.

Paso 2: Ajuste de un modelo lineal mixto

Vamos ahora a ajustar el modelo mixto equivalente al modelo lineal anterior. > library(nlme) > [Link] <- lme([Link] ~ Scenario * Habitat, random = ~1 | + factor(Plot), data = sn) > anova([Link]) numDF denDF F-value p-value 1 40 431.7217 <.0001 2 40 2.5693 0.0892 2 21 13.9081 0.0001 4 40 0.9049 0.4703

(Intercept) Scenario Habitat Scenario:Habitat

> summary([Link]) Linear mixed-effects model fit by REML Data: sn AIC BIC logLik 186.9631 210.1827 -82.48156 Random effects: Formula: ~1 | factor(Plot) (Intercept) Residual StdDev: 2.953028e-05 0.804266 Fixed effects: [Link] ~ Scenario * Habitat 160

Luis Cayuela

Modelos lineales mixtos en R

Value [Link] DF t-value p-value (Intercept) 1.2956250 0.2843510 40 4.556429 0.0000 ScenarioDry Summer 0.2855000 0.4021330 40 0.709964 0.4818 ScenarioWet Summer 0.3486250 0.4021330 40 0.866940 0.3911 HabitatGrassland 1.2767083 0.4343533 21 2.939332 0.0078 HabitatShrubland 0.3197500 0.4021330 21 0.795135 0.4354 ScenarioDry Summer:HabitatGrassland -0.5492083 0.5919237 40 -0.927836 0.3591 ScenarioWet Summer:HabitatGrassland 0.3647917 0.5919237 40 0.616282 0.5412 ScenarioDry Summer:HabitatShrubland -0.0058750 0.5687019 40 -0.010331 0.9918 ScenarioWet Summer:HabitatShrubland -0.0422500 0.5687019 40 -0.074292 0.9411 Correlation: (Intr) ScnrDS ScnrWS HbttGr HbttSh SDS:HG ScenarioDry Summer -0.707 ScenarioWet Summer -0.707 0.500 HabitatGrassland -0.655 0.463 0.463 HabitatShrubland -0.707 0.500 0.500 0.463 ScenarioDry Summer:HabitatGrassland 0.480 -0.679 -0.340 -0.734 -0.340 ScenarioWet Summer:HabitatGrassland 0.480 -0.340 -0.679 -0.734 -0.340 0.538 ScenarioDry Summer:HabitatShrubland 0.500 -0.707 -0.354 -0.327 -0.707 0.480 ScenarioWet Summer:HabitatShrubland 0.500 -0.354 -0.707 -0.327 -0.707 0.240 SWS:HG SDS:HS ScenarioDry Summer ScenarioWet Summer HabitatGrassland HabitatShrubland ScenarioDry Summer:HabitatGrassland ScenarioWet Summer:HabitatGrassland ScenarioDry Summer:HabitatShrubland 0.240 ScenarioWet Summer:HabitatShrubland 0.480 0.500 Standardized Within-Group Residuals: Min Q1 Med Q3 -1.79635849 -0.52757734 -0.05789441 0.32727824 Number of Observations: 70 Number of Groups: 24 Aunque este no es todav el modelo denitivo es curioso observar que ahora a tanto el hbitat como el escenario (con un = 0,10) son signicativos (la a interaccin entre ambos no lo es). En cuanto a su componente ja, este ser el o a modelo ms all del ptimo (el modelo completo). Por lo tanto podemos a a o utilizar esta estructura de los efectos jos para elegir la estructura de los efectos aleatorios ms adecuada. a

Max 5.22795549

3.3.

Paso 3: Eleccin de la estructura de los efectos o aleatorios

Los modelos que podemos plantear en este caso son: (1) un modelo sin efecto aleatorio (equivalente a un modelo lineal pero utilizando un estimador REML). 161

Luis Cayuela

Modelos lineales mixtos en R

Este modelo lo calcularemos con la funcin gls(); (2) un modelo que incluya la o varianza aleatoria de la constante (gran media) dentro de cada bloque; (3) un modelo que incluya la varianza de la constante y las varianzas aleatorias para cada uno de los niveles del factor habitat menos uno. Este ultimo modelo asumir que la respuesta (biomasa) a los distintos niveles del factor escenario a no es igual en todas las parcelas, sino que var de manera aleatoria dentro de a cada una. > > + > + > lme.sn0 <- gls([Link] ~ Scenario * Habitat, data = sn) lme.sn1 <- lme([Link] ~ Scenario * Habitat, random = ~1 | factor(Plot), data = sn) lme.sn2 <- lme([Link] ~ Scenario * Habitat, random = ~Scenario | factor(Plot), data = sn) anova(lme.sn0, lme.sn1, lme.sn2) Model 1 2 3 df AIC BIC logLik Test [Link] p-value 10 184.9631 206.0719 -82.48156 11 186.9631 210.1827 -82.48156 1 vs 2 0.00000 0.9999 16 179.8083 213.5822 -73.90413 2 vs 3 17.15485 0.0042

lme.sn0 lme.sn1 lme.sn2

Vemos que el modelo lme.sn2, que incluye la varianza de la constante y las varianzas aleatorias para cada uno de los niveles del factor habitat es el ms a parsimonioso. Esto es como asumir que hay una interaccin entre la parcela y o el factor escenario. Con el modelo lineal no ten amos rplicas como para e calcular un coeciente del efecto de cada parcela en cada nivel del factor escenario (esto ser (24-1 * 3-1 coecientes = 69 coecientes ms en el an a modelo). Pero como aqu lo que estamos calculando son solamente dos coecientes (la varianza aleatoria para dos de los tres niveles del factor, que indicar como var estos coecientes de manera aleatoria dentro de cada a an parcela), entonces s que tenemos capacidad de estimacin. o

3.4.

Paso 4: Eleccin de la estructura de los efectos jos o

Ya tenemos la estructura de los efectos aleatorios. Podemos utilizar el test t o el test F para comprobar la signicacin de la(s) variable(s) ja(s). Pero tal o vez es mejor opcin comparar modelos anidados como en el caso de los efectos o aleatorios. Recordemos que, para hacer esto, debemos estimar los parmetros a de los modelos utilizando un estimador de ML. > + > + > + > + > lme.sn3 <- lme([Link] ~ 1, random = ~Scenario | factor(Plot), data = sn, method = "ML") lme.sn4 <- lme([Link] ~ Scenario, random = ~Scenario | factor(Plot), data = sn, method = "ML") lme.sn5 <- lme([Link] ~ Scenario + Habitat, random = ~Scenario | factor(Plot), data = sn, method = "ML") lme.sn6 <- lme([Link] ~ Scenario * Habitat, random = ~Scenario | factor(Plot), data = sn, method = "ML") anova(lme.sn3, lme.sn4, lme.sn5, lme.sn6) 162

Luis Cayuela

Modelos lineales mixtos en R

lme.sn3 lme.sn4 lme.sn5 lme.sn6

Model 1 2 3 4

df 8 10 12 16

AIC 188.0546 185.6914 171.5988 171.0547

BIC 206.0425 208.1763 198.5807 207.0307

logLik Test [Link] p-value -86.02728 -82.84568 1 vs 2 6.363196 0.0415 -73.79939 2 vs 3 18.092591 0.0001 -69.52736 3 vs 4 8.544054 0.0736

El modelo factorial sin interaccin ser el ms ptimo de acuerdo a los o a a o criterios de informacin AIC y BIC y al test de verosimilitud. Esto indicar o a un efecto signicativo de los dos factores sobre la biomasa, pero no de la interaccin entre ambos, si bien esta es marginalmente signicativa o (p-valor=0.0736).

3.5.

Paso 5: Presentacin del modelo nal con REML o

Finalmente, presentamos el modelo nal utilizando REML. > [Link] <- lme([Link] ~ Scenario + Habitat, random = ~Scenario | + factor(Plot), data = sn) > anova([Link]) numDF denDF F-value p-value 1 44 687.4992 <.0001 2 44 3.4570 0.0403 2 21 17.7136 <.0001

(Intercept) Scenario Habitat

> summary([Link]) Linear mixed-effects model fit by REML Data: sn AIC BIC logLik 181.0709 207.1635 -78.53544 Random effects: Formula: ~Scenario | factor(Plot) Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 1.0137234 (Intr) ScnrDS ScenarioDry Summer 1.3007906 -0.958 ScenarioWet Summer 1.2175833 -0.882 0.875 Residual 0.3670845 Fixed effects: [Link] Value (Intercept) 1.3760020 ScenarioDry Summer 0.0883947 ScenarioWet Summer 0.4440614 HabitatGrassland 1.0851344 HabitatShrubland 0.3064254 ~ Scenario + [Link] DF 0.2485840 44 0.2928944 44 0.2775897 44 0.1867898 21 0.1814308 21 163 Habitat t-value p-value 5.535361 0.0000 0.301797 0.7642 1.599704 0.1168 5.809389 0.0000 1.688938 0.1060

Luis Cayuela

Modelos lineales mixtos en R

Correlation: ScenarioDry Summer ScenarioWet Summer HabitatGrassland HabitatShrubland (Intr) ScnrDS ScnrWS HbttGr -0.843 -0.783 0.829 -0.308 -0.052 -0.055 -0.365 0.000 0.000 0.486

Standardized Within-Group Residuals: Min Q1 Med Q3 -1.64242419 -0.26997412 -0.04686347 0.18271625 Number of Observations: 70 Number of Groups: 24

Max 1.88194226

Tanto el escenario y el hbitat son signicativos y tambin hay un efecto claro a e de la parcela que afecta a todos los valores de biomasa dentro de cada parcela (constante) y a cada valor de biomasa dentro de cada nivel del factor escenario (coeciente de los niveles del factor escenario). Por ultimo, es necesario que veamos cmo es de adecuado el modelo de acuerdo con los supuestos del o modelo, fundamentalmente normalidad y homocedasticidad.

164

Luis Cayuela

Modelos lineales mixtos en R

> > > > > > > > > > > >

Res <- residuals([Link], type = "normalized") Fit <- fitted([Link]) par(mfrow = c(3, 2)) plot(Res ~ Fit, xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs. fitted" abline(h = 0) boxplot(Res ~ sn$Scenario, ylab = "Residuals", main = "Scenario") abline(h = 0, lty = 3) boxplot(Res ~ sn$Habitat, ylab = "Residuals", main = "Habitat") abline(h = 0, lty = 3) hist(Res, main = "Histogram of residuals", xlab = "Residuals") qqnorm(Res) qqline(Res)
Residuals vs. fitted
q q q q q

Scenario
q

1.5

1.5

q q

0.5

Residuals

q q q q qq q qq q q q q qq q q qq q q q q qq q q qq q q q qq q qq q q q qq q q q q q q q q qq qq q

q q

0.5

Residuals

1.5

1.5

0.5

0.5

qq q

Control

Dry Summer

Fitted values

Habitat
q

Histogram of residuals
30 Frequency

1.5

Residuals

0.5

0.5

1.5

Forest

Grassland

Shrubland

0 2

10

15

20

25

0 Residuals

Normal QQ Plot
q q q q q q q q qq q qq qq q qq q qqq qq qq qq qq qq qqq q qq qqq qqq qq qq qq q qq qq q qq q q qq qq q qq q

Sample Quantiles

0.5

0.5

1.5

1.5

Theoretical Quantiles

165

Luis Cayuela

Modelos lineales mixtos en R

Parece que hay algo de heterocedasticidad debido a que no hay igualdad de varianzas entre los grupos. Adems, podr haber tambin problemas de a a e normalidad y linealidad. Si incluimos la interaccin en el modelo (que era o marginalmente signicativa en nuestro test de verosimilitud) el modelo se ajusta mejor a los supuestos de homocedasticidad y normalidad. Esto ocurre a veces cuando hay alguna variable que guarda relacin con la estructura de los o residuos, tal podr ser el caso de la interaccin entre la parcela y el escenario. a o Por tanto, una opcin viable ser modicar el modelo e incluir dicha o a interaccin en los efectos jos. o

4.

Dise os de medidas repetidas n

En el ejemplo del diseo por bloques aleatorizados, tendr n amos un diseo de n medidas repetidas si en vez de tres conejos de cada camada, tomamos un unico conejo y le aplicamos distintos tratamientos a cada uno de ellos en distintos momentos. No obstante, el anlisis de los datos mediante el uso de modelos a lineales mixtos hubiera sido exactamente igual. En el diseo por bloques n asumimos que hay un efecto de la camada que afecta a la relacin entre la o ganancia de peso y el tratamiento. En el diseo de medidas repetidas n asumimos que hay un efecto del individuo que afecta a la relacin entre la o ganancia de peso y el tratamiento. Vamos a usar un ejemplo de fructicacin en cerezos. Se quiere ver si la o produccin de frutos (FRUIT) depende del tamao del rbol, medido como el o n a dimetro a la altura del pecho (DBH). Para ello se miden 179 rboles. Los a a mismos individuos (TAG) son remuestreados durante 3 aos (YEAR). El DBH n se mide una sola vez para los tres aos y, aunque este puede aumentar de un n ao para otro, es tan pequeo el cambio (son rboles adultos) que no se ha n n a tenido en cuenta. > fruit <- [Link]("[Link] header = T) > head(z <- xtabs(FRUIT ~ factor(TAG) + factor(YEAR), data = fruit), + n = 15L) factor(YEAR) factor(TAG) 1 2 3 86 45 662 0 101 6 7 0 155 700 225 96 160 67 0 43 192 61 164 12 203 631 306 156 352 3 0 0 415 17 0 0 486 118 3266 319 487 113 322 30 786 495 534 11 787 86 73 0 166

Luis Cayuela

Modelos lineales mixtos en R

836 889 906

30 80 10

301 0 111

15 0 0

4.1.

Paso 1: Aplicacin de un modelo lineal o

Antes de ajustar un modelo lineal vamos a obtener un grco exploratorio de a la variable respuesta frente a la variable explicativa (por ao). n > plot(fruit$FRUIT ~ fruit$DBH, xlab = "DBH", ylab = "Nmero de frutos", u + pch = fruit$YEAR)

Nmero de frutos

6000

8000

4000

2000

q q q q q qq q q q q q q qq q q q q q q qq qq q q q q q q q q q q q q qqqq q q q q q q qq q q qqqq qqqq q q q qqqqq qq qq q q q qq qq q q q q q q q q q qq qqqq qq q qq qqqqq qqq q q q q q q q q q q q q q q q q qq q q qq qq q q q q q q

200

300

400

500

600 DBH

700

800

900

Alternativamente podr amos haber usado la funcin xyplot() del paquete o lattice, que permite hacer grcas condicionadas. Aparentemente se observa a una relacin positiva entre el tamao del rbol y la produccin de frutos. No o n a o obstante, esta relacin no parece muy lineal. Los datos son conteos por lo que o podr amos utilizar un modelo lineal generalizado mixto con una distribucin o de errores de tipo Poisson o binomial negativo. Para este caso concreto, utilizaremos una distribucin de errores de tipo normal, pero transformaremos o la variable respuesta con logaritmos, para linealizar la relacin. o

167

Luis Cayuela

Modelos lineales mixtos en R

> fruit$LOGFRUIT <- log(fruit$FRUIT + 1) > plot(fruit$LOGFRUIT ~ fruit$DBH, xlab = "DBH", ylab = "Nmero de frutos", u + pch = fruit$YEAR)

q q

q q q q q q q q

q q q

q q q qq q q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq q q q q qq q q q q q q q q qq q q qq q q q q qq qq qq q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q

Nmero de frutos

q q q

2 200

300

400

500

600 DBH

700

800

900

Vamos a analizar los datos con un modelo de medidas repetidas. Por tanto el modelo que plantearemos ser el siguiente: a > > + + + > > > DBHtabs <- xtabs(DBH ~ factor(TAG) + factor(YEAR), data = fruit) for (i in 1:dim(DBHtabs)[1]) { DBHtabs[i, ] <- ifelse(DBHtabs[i, ] == 0, max(DBHtabs[i, ]), DBHtabs[i, ]) } DBHtabs <- apply(DBHtabs, 1, mean) mlmfruit <- lm(z ~ DBHtabs) summary(mlmfruit)

Response 1 : Call: lm(formula = `1` ~ DBHtabs) Residuals: Min 1Q Median -607.3 -259.1 -115.3

3Q Max 24.1 8926.6

168

Luis Cayuela

Modelos lineales mixtos en R

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -324.8865 203.8840 -1.593 0.11284 DBHtabs 1.2005 0.3925 3.058 0.00257 ** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 833.9 on 177 degrees of freedom Multiple R-squared: 0.0502, Adjusted R-squared: 0.04483 F-statistic: 9.354 on 1 and 177 DF, p-value: 0.002571

Response 2 : Call: lm(formula = `2` ~ DBHtabs) Residuals: Min 1Q Median -584.21 -230.86 -110.25

3Q Max 10.87 2826.69

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -238.999 134.021 -1.783 0.0763 . DBHtabs 1.050 0.258 4.069 7.1e-05 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 548.2 on 177 degrees of freedom Multiple R-squared: 0.08555, Adjusted R-squared: 0.08039 F-statistic: 16.56 on 1 and 177 DF, p-value: 7.093e-05

Response 3 : Call: lm(formula = `3` ~ DBHtabs) Residuals: Min 1Q Median -429.3 -162.7 -58.9

3Q Max 43.2 4359.8

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -229.2795 107.2247 -2.138 0.033863 * DBHtabs 0.7561 0.2064 3.663 0.000330 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 438.6 on 177 degrees of freedom 169

Luis Cayuela

Modelos lineales mixtos en R

Multiple R-squared: 0.07045, Adjusted R-squared: 0.0652 F-statistic: 13.41 on 1 and 177 DF, p-value: 0.0003296 Como podemos observar, el tamao del rbol condiciona la produccin de n a o frutos, pero esta relacin cambia de un ao a otro, como indican las diferencias o n entre las pendientes estimadas de los tres aos. El ao tambin tiene un efecto n n e (aunque realmente esta no es una hiptesis en la que estemos especialmente o interesados). Estas diferencias interanuales quedan reejadas en los distintos intercept estimados para cada ao. n

4.2.

Paso 2: Ajuste de un modelo lineal mixto

Vamos ahora a ajustar el modelo mixto equivalente al modelo de medidas repetidas anterior. > library(nlme) > [Link] <- lme(LOGFRUIT ~ factor(YEAR) * DBH, random = ~1 | + factor(TAG), data = fruit) > anova([Link]) numDF denDF F-value p-value 1 245 2218.4152 <.0001 2 245 26.1822 <.0001 1 177 48.5718 <.0001 2 245 0.6707 0.5123

(Intercept) factor(YEAR) DBH factor(YEAR):DBH

> summary([Link]) Linear mixed-effects model fit by REML Data: fruit AIC BIC logLik 1571.561 1603.922 -777.7807 Random effects: Formula: ~1 | factor(TAG) (Intercept) Residual StdDev: 0.9104415 1.190307 Fixed effects: LOGFRUIT ~ factor(YEAR) * DBH Value [Link] DF t-value p-value (Intercept) 2.1086711 0.4037734 245 5.222412 0.0000 factor(YEAR)2 0.4989693 0.5095819 245 0.979174 0.3285 factor(YEAR)3 -0.0666292 0.4849343 245 -0.137398 0.8908 DBH 0.0042968 0.0007590 177 5.660905 0.0000 factor(YEAR)2:DBH 0.0000580 0.0009390 245 0.061766 0.9508 factor(YEAR)3:DBH -0.0008920 0.0008967 245 -0.994741 0.3208 Correlation: (Intr) fc(YEAR)2 fc(YEAR)3 DBH f(YEAR)2: 170

Luis Cayuela

Modelos lineales mixtos en R

factor(YEAR)2 -0.535 factor(YEAR)3 -0.590 0.457 DBH -0.957 0.509 factor(YEAR)2:DBH 0.519 -0.959 factor(YEAR)3:DBH 0.567 -0.441

0.559 -0.443 -0.956

-0.538 -0.583

0.466

Standardized Within-Group Residuals: Min Q1 Med Q3 -2.75479404 -0.63189697 0.05649328 0.64647442 Number of Observations: 428 Number of Groups: 179

Max 2.38612044

Aunque este no es todav el modelo denitivo, este modelo parece indicar que a realmente no hay una interaccin entre DBH y ao. Es decir, que las diferentes o n pendientes estimadas para cada ao no son realmente diferentes, por lo que n podemos asumir una unica pendiente que resume la relacin entre la o produccin de frutos y el DBH. o

4.3.

Paso 3: Eleccin de la estructura de los efectos o aleatorios

Los modelos que podemos plantear en este caso son: (1) un modelo sin efecto aleatorio (equivalente a un modelo lineal pero utilizando un estimador REML). Este modelo lo calcularemos con la funcin gls(); (2) un modelo que incluya la o varianza aleatoria de la constante (gran media) dentro de cada bloque; (3) un modelo que incluya la varianza de la constante y las varianzas aleatorias para cada uno de los niveles del factor ao menos uno. Este ultimo modelo asumir n a que la respuesta (FRUIT) a los distintos niveles del factor YEAR no es igual en todos los rboles, sino que var de manera aleatoria en cada uno. a a > > + > + > lme.fruit0 <- gls(LOGFRUIT ~ factor(YEAR) * DBH, data = fruit) lme.fruit1 <- lme(LOGFRUIT ~ factor(YEAR) * DBH, random = ~1 | factor(TAG), data = fruit) lme.fruit2 <- lme(LOGFRUIT ~ factor(YEAR) * DBH, random = ~factor(YEAR) | factor(TAG), data = fruit) anova(lme.fruit0, lme.fruit1, lme.fruit2) Model df AIC BIC logLik Test [Link] p-value 1 7 1609.605 1637.920 -797.8026 2 8 1571.562 1603.921 -777.7807 1 vs 2 40.04369 <.0001 3 13 1575.523 1628.108 -774.7616 2 vs 3 6.03820 0.3025

lme.fruit0 lme.fruit1 lme.fruit2

Vemos que el modelo lme.fruit1 es el ms parsimonioso. Este modelo incluye la a varianza de la constante, pero no las varianzas aleatorias para cada uno de los niveles del factor YEAR. Esto indica que no hay una interaccin entre el rbol o a y el factor ao. n

171

Luis Cayuela

Modelos lineales mixtos en R

4.4.

Paso 4: Eleccin de la estructura de los efectos jos o

Ya tenemos la estructura de los efectos aleatorios. Podemos utilizar el test t o el test F para comprobar la signicacin de la(s) variable(s) ja(s). Pero tal o vez es mejor opcin comparar modelos anidados como en el caso de los efectos o aleatorios. Recordemos que, para hacer esto, debemos estimar los parmetros a de los modelos utilizando un estimador de ML. > + > + > + > + > lme.fruit3 <- lme(LOGFRUIT ~ 1, random = ~1 | factor(TAG), data = fruit, method = "ML") lme.fruit4 <- lme(LOGFRUIT ~ factor(YEAR), random = ~1 | factor(TAG), data = fruit, method = "ML") lme.fruit5 <- lme(LOGFRUIT ~ factor(YEAR) + DBH, random = ~1 | factor(TAG), data = fruit, method = "ML") lme.fruit6 <- lme(LOGFRUIT ~ factor(YEAR) * DBH, random = ~1 | factor(TAG), data = fruit, method = "ML") anova(lme.fruit3, lme.fruit4, lme.fruit5, lme.fruit6) Model df AIC BIC logLik 1 3 1608.961 1621.138 -801.4804 2 5 1566.002 1586.297 -778.0008 1 3 6 1523.951 1548.305 -755.9754 2 4 8 1526.595 1559.068 -755.2975 3 Test [Link] p-value <.0001 <.0001 0.5077

lme.fruit3 lme.fruit4 lme.fruit5 lme.fruit6

vs 2 46.95928 vs 3 44.05080 vs 4 1.35568

El modelo factorial sin interaccin (lme.fruit5) ser el ms ptimo de acuerdo o a a o a los criterios de informacin AIC y BIC y al test de verosimilitud. Esto o indicar un efecto signicativo del factor (ao) y de la covariable (DBH), pero a n no de la interaccin entre ambos. o

4.5.

Paso 5: Presentacin del modelo nal con REML o

Finalmente, presentamos el modelo nal utilizando REML. > [Link] <- lme(LOGFRUIT ~ factor(YEAR) + DBH, random = ~1 | + factor(TAG), data = fruit) > anova([Link]) numDF denDF F-value p-value 1 247 2225.7522 <.0001 2 247 26.2072 <.0001 1 177 48.7138 <.0001

(Intercept) factor(YEAR) DBH

> summary([Link]) Linear mixed-effects model fit by REML Data: fruit AIC BIC logLik 172

Luis Cayuela

Modelos lineales mixtos en R

1544.358 1568.657 -766.1792 Random effects: Formula: ~1 | factor(TAG) (Intercept) Residual StdDev: 0.9081265 1.190028 Fixed effects: LOGFRUIT ~ factor(YEAR) + DBH Value [Link] DF t-value p-value (Intercept) 2.2676767 0.31107380 247 7.289835 0e+00 factor(YEAR)2 0.5329282 0.14373756 247 3.707648 3e-04 factor(YEAR)3 -0.5283849 0.14240656 247 -3.710397 3e-04 DBH 0.0039911 0.00057184 177 6.979526 0e+00 Correlation: (Intr) f(YEAR)2 f(YEAR)3 factor(YEAR)2 -0.180 factor(YEAR)3 -0.214 0.446 DBH -0.927 -0.020 0.009 Standardized Within-Group Residuals: Min Q1 Med Q3 -2.78647922 -0.60881273 0.03793956 0.65215410 Number of Observations: 428 Number of Groups: 179 El DBH tiene un efecto signicativo sobre la produccin de frutos (en realidad, o sobre el logaritmo de la produccin de frutos). En el segundo ao hay en o n promedio ms frutos que en el primer ao, mientras que en el tercer ao hay a n n muchos menos. Tambin hay un efecto del rbol sobre la produccin de frutos. e a o Por ultimo, es necesario que veamos cmo es de adecuado el modelo de o acuerdo con los supuestos del modelo, fundamentalmente normalidad y homocedasticidad.

Max 2.38337641

173

Luis Cayuela

Modelos lineales mixtos en R

> > > > > > > > > > > >

Res <- residuals([Link], type = "normalized") Fit <- fitted([Link]) par(mfrow = c(3, 2)) plot(Res ~ Fit, xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs. fitted" abline(h = 0) plot(Res ~ fruit$DBH, ylab = "Residuals", main = "DBH") abline(h = 0, lty = 3) boxplot(Res ~ fruit$YEAR, ylab = "Residuals", main = "A~o") n abline(h = 0, lty = 3) hist(Res, main = "Histogram of residuals", xlab = "Residuals") qqnorm(Res) qqline(Res)
Residuals vs. fitted
q q q q

DBH

q q q q q q qq q q q q q q qq q q q qq q q q q q q qq q q q q q q qq q q qqq q q q q q q q q q q qqq q qq q qq q q qq q q q q q qq qq q q q q q qq qq q qq q q q q qq q q q q q q q q qqq q qq q q q q q qq q q qqqq qq q qq q q q q qq q q q qq q q q qq q q q q q qqq qq qq qqq qqq q q q q q q q q qq q q q qq q q q qq q q q q qq q q q q q q q q qqq qq qqq qq q q qq q q q q q qq q qq q q q q q q qq q q q qq qqq q q q q qq q q q qq qqqqqq q q q q q qq q q qq q q q qq q q q qq q qq q q q q q q q qq q qq q qq q qqq q q q q q qq q q q qq q q q q qq q q qq qq q q q q q qq q q q q q q q q q q qqq qq q q qq q q q q q q q q q q q qq q q q q q q qq qqq q q q q q q qq qq q q q q qq q qq q q q q q

q q q q

2
q

q q q q q q q q qq q q qq qq qq qq q q q q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q qq q qq qq qq q q q q qq q q q q qq q q q q qq q q q q q q q q qq q qq qq q q q q q q q q q q q q qq q q qq qq q q qq q q q q q q qq q qq q q qq q q qq q qq q q q qq qq q q qq q qqq q q q q q q q q qq q qqqq q qq q q q q qq q q q q qq q q qq q q q q qq q qq qq q q qq qq q q q q q q qq q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q qqqqqqq q q q qq q q q q q q q qq qq qq qq q q q qqq q qqq q q q q q q q qq q q q qq q q q q q q qq q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q qq q q q qq q qq q q q q q qq q q q q q q qq q qq q q q q q q q q q

Residuals

Residuals

q q

200 300 400 500 600 700 800 900 fruit$DBH

Fitted values

Ao
100

Histogram of residuals

Frequency
q

Residuals

0 3

20

40

60

80

Residuals

Normal QQ Plot
q q

q q

q qq qq qq q q qq qq q qq qq q qq q qq qq qq qq qq qq qq qq qq qq q q q q q q q qq qq qq qq qq q q q q q q qq qq qq q q q q qq q q qq qq q q q q q q q q q q qq qq qq q q q qq q qq q q q qq qq qq q qq q q q q q qq q q qq qq qq q q q q q q qq qq q qq qq qq qq qq qq qq q qq q qq

Sample Quantiles

1 3

Theoretical Quantiles

174

Luis Cayuela

Modelos lineales mixtos en R

El modelo es bastante adecuado. Se cumplen los supuestos de normalidad, homocedasticidad y linealidad.

5.

Dise os factoriales anidados o jerarquizados n

En este caso de estudio vamos a tomar un ejemplo recogido en el cap tulo 20 de Zuur et al. (2009) que hace referencia a la tesis doctoral de Patricia Lastra Luque (2008)5 . La estimacin de la edad en los odontocetos se basa en el o clculo de los grupos de capas de crecimiento que se depositan en estructuras a de registro como los dientes. Generalmente, las secciones dentales se obtienen utilizando un microtomo criostato. Sin embargo, algunos investigadores preeren obtener secciones nas utilizando un microtomo de parana tradicional. Hay escasa informacin disponible acerca de la aplicacin de esta o o tcnica sobre los dientes de los delnes. El objetivo de este estudio era e investigar si la tcnica de parana pod verse afectada por el mtodo de e a e tincin, controlando los efectos de la especie, el sexo y la regin de o o procedencia. Para ello se tomaron muestras de dientes de varios individuos de diferentes especies de odontocetos6 en Espaa y Escocia. El factor de inters es n e el mtodo de tincin (Mayer, Elrich, Toluidine) y la variable respuesta la edad e o estimada del cetceo. Otras variables explicativas son el sexo o la regin de a o procedencia del cetceo (Espaa, Escocia). El esquema de la estructura de los a n datos ser la siguiente: a

Leemos los datos. > cet <- [Link]("[Link] header = T) > cet$DolphinID <- factor(cet$DolphinID) > str(cet)

'[Link]': 180 obs. of 6 variables: $ DolphinID: Factor w/ 60 levels "1","2","3","4",..: 9 9 9 24 24 24 36 36 36 46 ... $ Species : Factor w/ 6 levels "Delphinusdelphis",..: 1 1 1 1 1 1 1 1 1 1 ...
5 Luque, P.L. (2008) Age determination and interpretation of mineralization anomalies in teeth of small cetaceans. Tesis doctoral, Universidad de Vigo, Espaa. n 6 Cetceos con dientes como los delnes, las marsopas, los cachalotes o las belugas. a

175

Luis Cayuela

Modelos lineales mixtos en R

$ $ $ $

Age Sex Stain Location

: : : :

num 11.5 11.5 11 1.5 1.5 1.5 2 2 3 12 ... int 1 1 1 1 1 1 1 1 1 2 ... Factor w/ 3 levels "Elrich","Mayer",..: 2 1 3 2 1 3 2 1 3 2 ... Factor w/ 2 levels "Scotland","Spain": 1 1 1 1 1 1 1 1 1 2 ...

Si exploramos ms de cerca los datos vemos que hay algunos 0 en la variable a Sex. Esto es, el sexo de algunos individuos no ha podido ser determinado correctamente. Tenemos que eliminar estos datos antes de continuar con el anlisis. a > I <- cet$Sex == 0 > cet <- cet[!I, ] > cet$Sex <- factor(cet$Sex) Antes de ajustar un modelo lineal vamos a obtener algunos grcos a exploratorios. Estas grcas muestran que hay mucha variacin en cuanto a la a o edad estimada para los distintos espec menes. Hay que resaltar que hay tres muestras por especimen. La misma grca para las especies muestra que hay a mucha menor variacin entre especies. Esto no es raro ya que cada animal va a o tener una unica edad estimada tres veces, pero es muy posible que dentro de cada especie tengamos muestras de espec menes que cubran gran parte del rango de edades presente en la poblacin. Incluso aunque alguna especie sea o ms longeva que otra, va a haber mucha superposicin en las edades presentes a o en las muestras de los distintos espec menes. Tambin parece haber mayor e variacin en las edades registradas en Escocia que en Espaa por lo que puede o n que sea necesario especicar distintas varianzas en las distintas regiones. Por otro lado, no parece, a simple vista, que haya una respuesta diferencial de la edad estimada frente al mtodo de tincin en funcin del nivel de ninguno de e o o los otros factores (especie, especimen, sexo, regin). o

176

Luis Cayuela

Modelos lineales mixtos en R

> > > > > > > > > >

par(mfcol = c(3, 3)) boxplot(Age ~ Species, data = cet) boxplot(Age ~ DolphinID, data = cet) boxplot(Age ~ Sex, data = cet) boxplot(Age ~ Location, data = cet) boxplot(Age ~ Stain, data = cet) [Link](cet$Stain, cet$Species, cet$Age) [Link](cet$Stain, cet$DolphinID, cet$Age) [Link](cet$Stain, cet$Location, cet$Age) [Link](cet$Stain, cet$Sex, cet$Age)

25

25

q q

25

cet$DolphinID 2 58 61 41 60 51 54 57 46 49 4 48 9 Mayer Toluidine 3 23 cet$Stain 27 31 47 28 7 22 29 cet$Location 35 Spain 14 Scotland 30 50 11 42 43 37 59 6 33 21 12 Mayer Toluidine 34 1 cet$Stain 5 8 25 32 36 38 39 cet$Sex 44 45 1 13 2 19 20 26 53 40 56 15 18 24 17 52 Mayer Toluidine 55 10 cet$Stain

20

20

mean of cet$Age Scotland Spain

15

15

10

10

Delphinusdelphis

Stenellafrontalis

0 Elrich
q q

25

25

20

20

mean of cet$Age Elrich Mayer Toluidine

15

15

10

10

1 6 12

20

27

34

41

48

55

4 Elrich mean of cet$Age Stenellacoeruleoalba Delphinusdelphis Stenellafrontalis Tursiopstruncatus Lagenorhynchusacutus Phocoenaphocoena 6.0 Elrich 6.5 Elrich Mayer Toluidine cet$Stain

25

12

20

mean of cet$Age

15

10

5.1.

Paso 1: Aplicacin de un modelo lineal o

Vamos a analizar los datos con un modelo lineal. En este modelo vamos a calcular el efecto del mtodo de tincin, la regin y el sexo, sin considerar la e o o especia ni el individuo. Tambin vamos a contemplar todas las posibles e interacciones entre los factores dos a dos y la interaccin entre los tres factores. o > f1 <- formula(Age ~ Sex * Stain * Location) > [Link] <- lm(f1, data = cet) > anova([Link]) Analysis of Variance Table Response: Age 177

10

7.0

q q

cet$Species

10

10

15

20

Luis Cayuela

Modelos lineales mixtos en R

Df Sum Sq Mean Sq F value Pr(>F) Sex 1 10.9 10.95 0.3354 0.5633 Stain 2 34.2 17.09 0.5234 0.5935 Location 1 1024.1 1024.15 31.3711 8.741e-08 *** Sex:Stain 2 1.2 0.58 0.0177 0.9825 Sex:Location 1 31.7 31.75 0.9724 0.3255 Stain:Location 2 17.7 8.85 0.2710 0.7629 Sex:Stain:Location 2 3.5 1.76 0.0539 0.9476 Residuals 165 5386.6 32.65 --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

5.2.

Paso 2: Ajuste de un modelo lineal mixto

Vamos ahora a ajustar el modelo mixto equivalente al modelo lineal anterior, pero considerando los efectos aleatorios, que en este caso vienen dados por el especimen y la especie. El especimen ser un factor de medidas repetidas ya a que dentro de cada especimen comprobamos todos los niveles del factor mtodo e de teido. Pero los espec n menes estn anidados dentro del factor especie. Esto a es porque no todos los niveles del factor especimen estn recogidos en todos los a niveles del factor especie. La formulacin del modelo ser as o a : > library(nlme) > [Link] <- lme(f1, random = ~1 | Species/DolphinID, data = cet) > anova([Link]) numDF denDF F-value p-value 1 110 60.19087 <.0001 1 50 0.05354 0.8180 2 110 24.86562 <.0001 1 50 8.18786 0.0061 2 110 0.83863 0.4350 1 50 0.28262 0.5973 2 110 12.87702 <.0001 2 110 2.56000 0.0819

(Intercept) Sex Stain Location Sex:Stain Sex:Location Stain:Location Sex:Stain:Location > summary([Link])

Linear mixed-effects model fit by REML Data: cet AIC BIC logLik 740.3277 786.9168 -355.1638 Random effects: Formula: ~1 | Species (Intercept) StdDev: 0.898073 178

Luis Cayuela

Modelos lineales mixtos en R

Formula: ~1 | DolphinID %in% Species (Intercept) Residual StdDev: 5.615181 0.828939 Fixed effects: list(f1) (Intercept) Sex2 StainMayer StainToluidine LocationSpain Sex2:StainMayer Sex2:StainToluidine Sex2:LocationSpain StainMayer:LocationSpain StainToluidine:LocationSpain Sex2:StainMayer:LocationSpain Sex2:StainToluidine:LocationSpain Correlation: Sex2 StainMayer StainToluidine LocationSpain Sex2:StainMayer Sex2:StainToluidine Sex2:LocationSpain StainMayer:LocationSpain StainToluidine:LocationSpain Sex2:StainMayer:LocationSpain Sex2:StainToluidine:LocationSpain Value 4.142943 -0.378897 0.375000 0.162500 4.480572 0.062500 0.170833 -0.902499 2.201923 1.029808 -1.407280 -0.738141 (Intr) -0.517 -0.093 -0.093 -0.642 0.057 0.057 0.356 0.058 0.058 -0.039 -0.039 Sx2:ST [Link] 1.4148661 2.0844159 0.2621335 0.2621335 2.1746964 0.4280622 0.4280622 3.0639805 0.4176455 0.4176455 0.6221848 0.6221848 DF 110 50 110 110 50 110 110 50 110 110 110 110 t-value p-value 2.928152 0.0041 -0.181776 0.8565 1.430569 0.1554 0.619913 0.5366 2.060321 0.0446 0.146007 0.8842 0.399085 0.6906 -0.294551 0.7696 5.272229 0.0000 2.465746 0.0152 -2.261836 0.0257 -1.186369 0.2380

Sex2 0.063 0.063 0.332 -0.103 -0.103 -0.682 -0.039 -0.039 0.071 0.071 Sx2:LS

StnMyr StnTld LctnSp Sx2:SM

0.500 0.060 -0.612 -0.306 -0.043 -0.628 -0.314 0.421 0.211 StM:LS

0.060 -0.306 -0.612 -0.043 -0.314 -0.628 0.211 0.421 StT:LS

-0.037 -0.037 0.500 -0.622 0.070 -0.096 0.384 -0.096 0.192 0.064 -0.688 0.064 -0.344 S2:SM:

Sex2 StainMayer StainToluidine LocationSpain Sex2:StainMayer Sex2:StainToluidine Sex2:LocationSpain 0.070 StainMayer:LocationSpain 0.192 0.068 StainToluidine:LocationSpain 0.384 0.068 0.500 Sex2:StainMayer:LocationSpain -0.344 -0.102 -0.671 -0.336 Sex2:StainToluidine:LocationSpain -0.688 -0.102 -0.336 -0.671 Standardized Within-Group Residuals: Min Q1 Med Q3 -2.86316216 -0.35443083 -0.00631148 0.29365611 Number of Observations: 177 179

0.500

Max 3.67902864

Luis Cayuela

Modelos lineales mixtos en R

Number of Groups: Species DolphinID %in% Species 6 59

5.3.

Paso 3: Eleccin de la estructura de los efectos o aleatorios

Los modelos que podemos plantear en este caso son: (1) un modelo sin efecto aleatorio (equivalente a un modelo lineal pero utilizando un estimador REML). Este modelo lo calcularemos con la funcin gls(); (2) un modelo que incluya la o varianza aleatoria de la constante (gran media) dentro de cada especimen (anidado dentro de especie) y de cada especie (por tanto ser dos coecientes an de varianza). En los grcos de perl no hemos visto aparentemente ninguna interaccin a o entre el efecto de la tincin y ninguno de los otros factores sobre la edad. En o un modelo lineal mixto la interaccin entre el mtodo de tincin y los factores o e o aleatorios se estimar por medio de las varianzas aleatorias estimadas para a cada uno de los niveles del factor (mtodo de teido) menos uno. Pero como a e n priori no parece haber una interaccin, podr o amos obviar este modelo ms a complejo o, simplemente incluirlo y ver si es mejor o peor que los otros dos modelos comparndolos mediante el test de verosimilitud. a > lme.cet0 <- gls(f1, data = cet) > lme.cet1 <- lme(f1, random = ~1 | Species/DolphinID, data = cet) Vemos que esto da un error para el modelo lme.cet2. Lo que ocurre es que hacen falta ms iteraciones para poder estimar los coecientes aleatorios por a medio del estimador REML. Para modicar el nmero de iteraciones en el u proceso de estimacin de los coecientes aleatorios del modelo podemos o utilizar la funcin lmeControl(). o > lmc <- lmeControl(niterEM = 5200, msMaxIter = 5200) > lme.cet2 <- lme(f1, random = ~Stain | Species/DolphinID, data = cet, + control = lmc) > anova(lme.cet0, lme.cet1, lme.cet2) Model 1 2 3 df AIC BIC logLik Test [Link] p-value 13 1101.4488 1141.8261 -537.7244 15 740.3277 786.9168 -355.1638 1 vs 2 365.1212 <.0001 25 704.9049 782.5536 -327.4525 2 vs 3 55.4227 <.0001

lme.cet0 lme.cet1 lme.cet2

Vemos que el modelo lme.cet1, que contiene dos trminos aleatorios para la e constante (uno para el especimen y otro para la especie) es mucho ms a adecuado que el modelo sin efectos aleatorios (el AIC disminuye de 1101.45 a 740.33). Cuando incorporamos adems el efecto de la varianza sobre los a coecientes de los niveles del factor principal (mtodo de teido), el modelo e n 180

Luis Cayuela

Modelos lineales mixtos en R

mejora un poco ms (lme.cet2), pero el cambio ms sustancial se produce al a a incluir la varianza atribuible al efecto de los dos factores aleatorios (uno anidado dentro del otro) sobre la constante del modelo. Como la estimacin de o los coecientes de este modelo es costosa y su interpretabilidad es menor, vamos a quedarnos con el modelo lme.cet1.

5.4.

Paso 4: Eleccin de la estructura de los efectos jos o

Ya tenemos la estructura de los efectos aleatorios. Vamos ahora a comparar modelos anidados como en el caso de los efectos aleatorios. Recordemos que, para hacer esto, debemos estimar los parmetros de los modelos utilizando un a estimador de ML. > + > + > + > + > + > + > + > + > + lme.cet3 <- lme(Age ~ 1, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet4 <- lme(Age ~ Sex, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet5 <- lme(Age ~ Stain, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet6 <- lme(Age ~ Stain + Location, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet7 <- lme(Age ~ Stain * Location, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet8 <- lme(Age ~ Stain * Location + Sex:Stain, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet9 <- lme(Age ~ Stain * Location + Sex:Location, random = ~1 | Species/DolphinID, data = cet, method = "ML") lme.cet10 <- lme(Age ~ Stain * Location + Sex:Location:Stain, random = ~1 | Species/DolphinID, data = cet, method = "ML") anova(lme.cet3, lme.cet4, lme.cet5, lme.cet6, lme.cet7, lme.cet8, lme.cet9, lme.cet10) Model 1 2 3 4 5 6 7 8 df 4 5 6 7 9 12 11 15 AIC 796.8521 798.8503 765.5389 760.8669 743.6632 744.9647 746.6293 745.2448 BIC 809.5567 814.7311 784.5958 783.0999 772.2486 783.0785 781.5669 792.8870 logLik -394.4260 -394.4252 -376.7695 -373.4334 -362.8316 -360.4824 -362.3146 -357.6224 Test 1 2 3 4 5 6 7 vs vs vs vs vs vs vs [Link] p-value 0.9665 <.0001 0.0098 <.0001 0.1953 0.0556 0.0522

lme.cet3 lme.cet4 lme.cet5 lme.cet6 lme.cet7 lme.cet8 lme.cet9 lme.cet10

2 0.00176 3 35.31138 4 6.67206 5 21.20366 6 4.69846 7 3.66451 8 9.38450

El modelo ms parsimonioso ser lme.cet7, con el mtodo de teido, la regin a a e n o y la interaccin entre ambos factores como efectos jos signicativos. o

5.5.

Paso 5: Presentacin del modelo nal con REML o

Finalmente, presentamos el modelo nal utilizando REML. 181

Luis Cayuela

Modelos lineales mixtos en R

> [Link] <- lme(Age ~ Stain * Location, random = ~1 | Species/DolphinID, + data = cet) > anova([Link]) numDF denDF F-value p-value 1 114 52.17758 <.0001 2 114 23.79974 <.0001 1 52 6.64360 0.0128 2 114 11.22038 <.0001

(Intercept) Stain Location Stain:Location

> summary([Link]) Linear mixed-effects model fit by REML Data: cet AIC BIC logLik 744.9385 773.2135 -363.4693 Random effects: Formula: ~1 | Species (Intercept) StdDev: 1.287742 Formula: ~1 | DolphinID %in% Species (Intercept) Residual StdDev: 5.515014 0.8472978 Fixed effects: Age ~ Stain * Location Value [Link] DF t-value p-value (Intercept) 4.049757 1.3567527 114 2.984889 0.0035 StainMayer 0.398438 0.2118245 114 1.880980 0.0625 StainToluidine 0.226563 0.2118245 114 1.069577 0.2871 LocationSpain 3.928176 1.8114631 52 2.168510 0.0347 StainMayer:LocationSpain 1.481192 0.3131268 114 4.730327 0.0000 StainToluidine:LocationSpain 0.671586 0.3131268 114 2.144772 0.0341 Correlation: (Intr) StnMyr StnTld LctnSp StM:LS StainMayer -0.078 StainToluidine -0.078 0.500 LocationSpain -0.720 0.058 0.058 StainMayer:LocationSpain 0.053 -0.676 -0.338 -0.086 StainToluidine:LocationSpain 0.053 -0.338 -0.676 -0.086 0.500 Standardized Within-Group Residuals: Min Q1 Med Q3 -3.17606813 -0.31200693 -0.04456998 0.24460500 Number of Observations: 177 Number of Groups: Species DolphinID %in% Species 6 59 182

Max 4.04733105

Luis Cayuela

Modelos lineales mixtos en R

El mtodo de teido puede afectar la estimacin que hagamos de la edad, con e n o el mtodo Mayer aumentando de manera (ligeramente) signicativa la e estimacin de la edad del cetceo con respecto a los otros dos mtodos. o a e Tambin la estimacin de la edad est condicionada por la regin geogrca, e o a o a con una edad estimada mayor en Espaa que en Escocia. Si vemos el grco n a de interaccin de estos dos factores podemos observar que el efecto ms o a destacado del mtodo de teido Mayer sobre la estimacin de la edad se e n o produce en las muestras de cetceos tomadas en Espaa mientras que en a n Escocia este efecto no es tan marcado y apenas se perciben diferencias entre los tres mtodos. Esto puede ser debido a que en Espaa los individuos e n muestreados son, en promedio, ms mayores y podr ser que la estimacin del a a o mtodo Mayer se vea afectada por la edad del especimen. e Debemos tambin comprobar la idoneidad del modelo desde el punto de vista e de los supuestos de normalidad, homocedasticidad y linealidad.

183

Luis Cayuela

Modelos lineales mixtos en R

> > > > > > > > > > > >

Res <- residuals([Link], type = "normalized") Fit <- fitted([Link]) par(mfrow = c(3, 2)) plot(Res ~ Fit, xlab = "Fitted values", ylab = "Residuals", main = "Residuals vs. fitted" abline(h = 0) boxplot(Res ~ cet$Stain, ylab = "Residuals", main = "Stain") abline(h = 0, lty = 3) boxplot(Res ~ cet$Location, ylab = "Residuals", main = "Location") abline(h = 0, lty = 3) hist(Res, main = "Histogram of residuals", xlab = "Residuals") qqnorm(Res) qqline(Res)
Residuals vs. fitted
4 4
q

Stain
q

q q

q q q q q q qq q q q q q qq q q q qq q qq q q q q q q q q qq q qq q q q qqqq q qq q qq q q q q q q q q q q qqq qqq qqq qq q q q qq q q qqq q q q qq qqq q q q q q q qq qq q q qq q qq q q qq qq q q q q q q qq qq q q q q q q q

q q

q q q q q q q q q q

Residuals

q
q q q q q

Residuals

q q q q

10

15

20

25

Elrich

Mayer

Toluidine

Fitted values

Location
4
q

Histogram of residuals
80 Frequency
q

Residuals

q q q q q

q q

q q q

Scotland

Spain

0 4

20

40

60

0 Residuals

Normal QQ Plot
4
q

Sample Quantiles

q q q qq qq qqq qq q q qq qq qq qq qq qq qq qq qq q q q q q q q qqq qqqq qqqq qqqq qq qqqq qqqq qqqq qqqq qqq qqq q qqq qqqq qq qq qq qq qq qq qq qq q q q q q q q q q qq qqq qqq q q q qq q q

Theoretical Quantiles

184

Luis Cayuela

Modelos lineales mixtos en R

Parece que hay heterocedasticidad en los residuos de las distintas regiones geogrcas y tambin en los mtodos de teido. En el conjunto de los residuos a e e n no se detecta este problema. Por otro lado, el histograma de los residuos tiene una forma un poco elongada (posibles problemas de leptokurtosis) y el qq-plot indica una falta de linealidad. Para corregir la falta de homocedasticidad se podr incluir una componente en el modelo que indicase diferencia de a varianzas entre diferentes regiones geogrcas con el argumento weights de la a funcin lme() (ver Zuur et al. 2009, pp. 464), si bien esto no lo vamos a ver o aqu Tambin mejorar el modelo si incorporamos la estructura de efectos . e a aleatorios ms compleja que vimos en el modelo lme.cet2. Ahora bien, si a decidimos quedarnos con esta estructura de efectos aleatorios deber amos repetir los pasos 4 y 5 para seleccionar de nuevo la estructura ms adecuada a para los efectos jos.

6.

Ms ejemplos a

Se pueden encontrar ms ejemplos resueltos en a [Link]

7.

Referencias
Zuur, A.F., Ieno, E.N., Walker, N.J., Saveliev, A.A. & Smith, G.M. 2009.

Mixed eects models and extensions in ecology with R. Springer, New York.
Zuur, A.F., Ieno, E.N. & Smith, G.M. 2007. Analysing ecological data.

Springer, New York.

Springer Verlag, New York.

linear, mixed eects and non parametric regression models. Chapman & Hall/CRC Press, Florida, USA.

185