1 Introducción

Anteriormente vimos la diferencia entre efectos fijos y aleatorios. Fijos son aquellos que, son parte de las condiciones que controlamos en el experimento y que afectan la media de cada tratamiento. En contraste los efectos aleatorios, son aquellos que no estamos controlando, generalmente relacionados con algún factor de agrupación como los bloques y que afectan la varianza del objeto de estudio. Entonces, modelos de efectos mixtos se refiere a aquellos que contienen tanto efectos aleatorios como fijos.

El programa original del curso contempla los modelos de efectos fijos y aleatorios por separado. Por motivos didpacticos yo decidí introducir los modelos de efectos mixtos antes porque considero que es más fácil aprender a distinguir sus diferencias de esta manera. Entonces ahora vamos a ver, sólo como un ejemplo, cómo sería un modelo únicamente con efectos aleatorios.

2 Modelo con efectos aleatorios

En los modelos con efectos mixtos vimos cómo se podía estimar la media de cada tratamiento y aislar los efectos de los bloques. En los modelos con efectos mixtos la meta es estimar la media global de la variable de respuesta y evitar que la variabilidad introducida por factores de agrupación alteren las estimaciones de la media. Supongamos que tenemos un diseño experimental para analizar con ANOVA de una vía, y antes de hacer el análisis queremos ver cuál es el promedio global de la variable de respuesta. La manera más correcta de hacer la estimación entonces sería con un modelo de efectos aleatorios. La manera más fácil de hacer esto en R con aov es:

modelo <- aov(formula = y ~ 1 + Error(Tratamientos), data = datos)

Es importante conocer este tipo de modelos porque en muchas áreas del conocimiento, incluyendo la ecología. Por ejemplo, la estimación de parámetros poblacionales como la densidad (cantidad de individuos por unidad de área) frecuente requiere de modelos de efectos aleatorios que se utilizan para tomar en cuenta los efectos del esfuerzo de muestreo o las características ambientales (densidad de vegetación) sobre las estimaciones de los parámetros poblacionales.

3 Modelos con efectos mixtos

Como ya hemos ajustado una serie de modelos con efectos mixtos, nos enfocaremos aquí en analizar diseños experimentales factoriales de dos vías con efectos mixtos utilizando métodos más robustos que aov.

3.1 Efectos mixtos con lme4

La principal diferencia entre aov y lmer radica en que lmer no utiliza las operaciones aritméticas de suma de cuadrados para estimar los efectos, sino el concepto de máxima verosimilitud. En breve, este método consiste en buscar por medio de una rutina computacional la combinación de parámetros del modelo que mejor se ajustan a los datos experimentales. Al igual que aov, lmer asume que la variable de respuesta tiene una distribución normal, de modo que se buscarán los parámetros \(\mu\) (media) y \(sigma\) (desviación estándar) que satisfacen la ecuación que describe una distribución normal:

\[ \frac{1}{\sigma\sqrt{2\pi}} \exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right) \]

Entonces, una vez estimados los parámetros se calculan las probabilidades de que los tratamientos afecten la media de los datos observados.

3.2 Ejemplos aplicados

Comenzaremos por re-analizar la base de datos npk de la última actividad. Si recordamos esta base de datos fue generada para su análisis con un ANOVA de tres vías. Aquí lo volveremos a analizar con ANOVA de dos vías pero tomando en cuenta los bloques. Para tener más grados de libertad, asumiremos que el factor N no existe, por lo tanto analizaremos la producción yield como función de la fertilización con nitrógeno N y potasio K en plantas que recibieron ambos tratamientos de P.

Comenzaremos al igual que antes por graficar los efectos fijos de ambos factores. En esta ocasión utilizaremos un método distinto llamado gráficas de violín utilizando un paquete más especializado y potenque que graphics, el paquete ggplot2.

npk.1 <- subset(npk, select = c("block", "N", "K", "yield")) #Quitando la columna N

library(ggplot2)

ggplot(npk.1) + geom_violin(aes(x = N, y = yield, colour = factor(K), fill = factor(K)), alpha = 0.3) +
      geom_boxplot(aes(x = N, y = yield, colour = factor(K), fill = factor(K)), alpha = 0.1) + 
      labs(colour= "K", fill = "K")
Gráfico de cajas y violín de `yield` en sustratamientos.

Figure 3.1: Gráfico de cajas y violín de yield en sustratamientos.

Estas gráficas de violín muestr la densidad de yield. Las zonas más anchas corresponden a zonas de mayor densidad de observaciones, y por lo tanto cercanas a las medidas de tendencia central (media, mediana). Para referencia en comparación con los gráficos de caja se muestra un boxplot de los mismos datos.

Como podemos ver la distribución de yield no es normal, por lo que e posible que haya problemas con ANOVA. Sin embargo, por se un diseño factorial con efectos mixtos, no podemos utilizar kruskal.test, pues esta sólo sirve para diseños de una vía. Entonces podemos ver si el modelo con efectos mixtos satiscafe nuestras necesidades.

3.2.1 Ajustando el modelo lmer

Comenzaremos por ajustar el modelo más complejo, es decir aquel que contiene los términos de \(1^{er}\) y \(2^o\) orden (N + K + N:K). Como vimos anteriormente, \(SS-II\) es un poco más fácil de interpretar que \(SS-I\) pues los resultados son los mismos independientemente de cómo especifiquemos el modelo:

library(lme4)

m1 <- lmer(yield ~ factor(N) * factor(K) + (1|block), npk.1) #Modelo de efectos mixtos
summary(m1) # Resumen del modelo
## Linear mixed model fit by REML ['lmerMod']
## Formula: yield ~ factor(N) * factor(K) + (1 | block)
##    Data: npk.1
## 
## REML criterion at convergence: 125
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4278 -0.5935  0.2340  0.5530  1.3855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  block    (Intercept) 13.57    3.684   
##  Residual             14.36    3.790   
## Number of obs: 24, groups:  block, 6
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)             52.883      2.158  24.508
## factor(N)1               7.967      2.188   3.641
## factor(K)1              -1.633      2.188  -0.746
## factor(N)1:factor(K)1   -4.700      3.094  -1.519
## 
## Correlation of Fixed Effects:
##             (Intr) fc(N)1 fc(K)1
## factor(N)1  -0.507              
## factor(K)1  -0.507  0.500       
## fc(N)1:(K)1  0.359 -0.707 -0.707
car::Anova(m1, type = 2) # Tabla de ANOVA
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: yield
##                       Chisq Df Pr(>Chisq)    
## factor(N)           13.1780  1  0.0002833 ***
## factor(K)            6.6281  1  0.0100385 *  
## factor(N):factor(K)  2.3069  1  0.1288009    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En el cómputo de la tabla de ANOVA podemos ver que la suma de cuadrados para la interacción (columna Chisq), tiene la menor suma de cuadrados, por lo que volveremos a ajustar el modelo sin la interacción, que además tiene la \(P = 0.13\) más grande.

m2 <- lmer(yield ~ factor(N) + factor(K) + (1|block), npk.1)
car::Anova(m2, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: yield
##             Chisq Df Pr(>Chisq)    
## factor(N) 12.1829  1  0.0004823 ***
## factor(K)  6.1275  1  0.0133091 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Y ahora sí, tenemos un modelo sólo con efectos significativos, con esto podemos empezar a sospechar que:

  1. Podemos rechazar lahipótesis nula de no diferencia entre tratamientos
  2. La fertilización con nitrógeno N explica gran parte de esa variabilidad
  3. Es poco probable que haya interacciones entre N y K.

Como de inicio sospechamos que la distribución de los datos dentro de los tratamientos podría se un problema, vamos a hacer sl diagnóstico visual del modelo (para lmer es un poco diferente que para aov y lm):

plot(m2)
Gráfico de diagnóstico de los residuales de `m2`.

Figure 3.2: Gráfico de diagnóstico de los residuales de m2.

Y vamos a ver si los residuales de m2 tienen una distribución normal centrada en cero:

resids <- residuals(m2) # Extrayendo los residuales de m2
mean(resids) # La media de los residuales
## [1] -8.585732e-15
shapiro.test(resids) # Probando la normalidad de los residuales
## 
##  Shapiro-Wilk normality test
## 
## data:  resids
## W = 0.97539, p-value = 0.7981

Con esto podemos ver que m2 parece robusto para las pruebas de hipótesis.

Regresar al índice del curso