Modelos Estadísticos. Grado Biotecnología

Modelos ANOVA

En este tema se presentan los modelos ANOVA. Esto modelos surgen cuando la variable o variables predictoras son de tipo categórico, es decir, factores con diferentes niveles de clasificación (o categorías). La variable respuesta sigue siendo de tipo numérico y el objetivo principal es el estudio de la media para los diferentes niveles del factor. Se trata de una generalización de la comparaciones de medias estudiadas en la unidad de inferencia. No solo se está interesado en comparar los diferentes grupos sino además en cuantificar numéricamente esas diferencias.

La principal diferencia con los modelos de regresión es que en este caso no modelizamos toda la respuesta observada sino la media observada para cada nivel del factor. En los puntos siguientes veremos como escribir dicho modelo cuando tenemos uno o dos factores de clasificación. Introduciremos el concepto de interacción que nos permitirá plantear modelos más complejos.

Modelo ANOVA de una via

El modelo de una vía se plantea cuando tenemos una única variable predictora de tipo categórico y una respuesta de tipo numérico. La situación más general, para una muestra de tamaño \(n\), viene dada por:

  • Una variable respuesta, \(Y\), de tipo numérico con observaciones \(Y_1,...,Y_n\).
  • Una variable predictora, \(F\), de tipo categórico con \(k\) grupos o niveles distintos de tamaños muestrales \(n_1,n_2,...,n_k\), de forma que \(n = n_1 + n_2 + ... + n_k\).
  • \(\mu_i\) denota a la media de todas las observaciones de la respuesta asociadas con el nivel \(i\) del factor.

Con esta estructura el vector de respuestas se puede organizar separando las respuestas para cada uno de los niveles del factor: \[y_1,...,y_{n_1},y_{n_1+1},...,y_{n_1+n_2},...,y_{n_1+n_2+...+n_{k-1}+1}..,y_n\] Para construir un modelo para este tipo de datos debemos crear un conjunto de variables predictoras ficticias, de forma que, podamos representar que cada elemento de la respuesta pertenece a un grupo u a otro. Para ello se crean \(k-1\) variables ficticias, \(X_1,...X_{k-1}\) que toman los valores siguientes: \[X_i = \left\{\begin{array}{ll} 1 & \text{si la respuesta pertence al grupo i}\\ 0 & \text{en otro caso}\\ \end{array}\right.\]

El modelo se puede escribir entonces: \[ Y = \alpha_0 + \alpha_1 * X_1 + \alpha_2 * X_2 + ... + \alpha_{k-1} * X_{k-1} + \epsilon\] donde \(\alpha_0\) representa el efecto común a todos los grupos, y los \(\alpha_1,\alpha_2,...,\alpha_{k-1}\) representan el efecto asociado con cada uno de los grupos del factor. En este tipo de modelos siempre se añade la restricción de que el efecto asociado con el grupo \(k\) viene determinado únicamente por efecto común, es decir, \(\alpha_k = 0\). Esta estructura nos permite ver que los modelos ANOVA se construyen como cualquier otro modelo lineal, pero en la practica la forma de escribir este modelo viene dada por: \[\mu_i = \alpha_0 + \alpha_i + \epsilon_i, \text{ para } i = 1, 2,...,k, \text{ con } \alpha_k = 0\]

Las hipótesis de este modelo es que los errores se distribuyen de forma independiente mediante una distribución Normal de media cero y varianza constante \(\sigma^2\) para cada uno de los grupos que determina la variable predictora.

La especificación de este modelo se realiza a través de la ecuación: \[ Y \sim F\] En este tipo de modelos se plantean varias situaciones de modelos anidados que deberemos estudiar para la selección del modelo más adecuado: * Modelo con efecto de \(F\): \[ Y \sim F\]. Se detectan diferencias entre al menos dos medias para los niveles del factor \(F\). * Modelo sin efectos: \[ Y \sim 1\]. No se detectan diferencias entre las medias.

Modelo ANOVA de dos vías

El modelo ANOVA de dos vías aparece cuando consideramos dos variables predictoras de tipo categórico. Si tenemos dos factores \(F_1\) y \(F_2\) con \(k\) y \(l\) niveles respectivamente, el modelo viene dado por: \[\mu_{ij} = \alpha_0 + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ij}, \text{ para } i = 1,...,k; j = 1,...,l \text{ con } \alpha_k = \beta_l = \alpha\beta_{kl} = 0; \]

donde:

  • \(\mu_{ij}\) es la media de la respuesta para el nivel \(i\) de \(F_1\) y el nivel \(j\) de \(F_2\).
  • \(\alpha_0\) es el efecto común a todas las combinaciones de niveles para ambos factores.
  • \(\alpha_i\) son los incrementos asociados a cada uno de los niveles del factor \(F_1\)
  • \(\beta_j\) son los incrementos asociados a cada uno de los niveles del factor \(F_2\)
  • \(\alpha\beta_{ij}\) son los incrementos asociados a la combinación de los niveles \(i\) de \(F1\) y \(j\) de \(F_2\). Es lo que se denomina efecto de interacción que valora los cambios en las combinaciones de las medias para ambos factores. Si este efecto no es significativo la media para una combinación de niveles de los factores se construye sumando los incrementos de cada nivel de forma independiente.

La especificación de este modelo se realiza a través de la ecuación: \[ Y \sim F_1 + F_2 + F_1:F_2,\] donde \(F_1:F_2\) representa el efecto de interacción. En este tipo de modelos se plantean varias situaciones de modelos anidados que deberemos estudiar para la selección del modelo más adecuado:

  • Modelo saturado: \[ Y \sim F_1 + F_2 + F_1:F_2\]. Se detecta alguna diferencia entre las medias obtenidas para dos combinaciones diferentes de ambos factores.
  • Modelo con efectos de los factores independientes: \[ Y \sim F_1 + F_2\]. Se detecta alguna diferencia entre las medias obtenidas para un nivel de \(F_1\) y un nivel de \(F_2\).
  • Modelo con efecto de \(F_1\): \[ Y \sim F_1\]. Se detectan diferencias al menos entre dos medias para los niveles del factor \(F_1\) pero no para \(F_2\).
  • Modelo con efecto de \(F_2\): \[ Y \sim F_2\]. Se detectan diferencias al menos entre dos medias para los niveles del factor \(F_2\) pero no para \(F_1\).
  • Modelo sin efectos: \[ Y \sim 1\]. No se detectan diferencias entre las medias.

Ejemplos

Veamos los diferentes ejemplos con los que vamos a trabajar

# Cargamos las librerías
library(tidyverse)
library(forcats)
library(broom)
library(reshape2)
library(lmtest)
library(mgcv)
library(MASS)

Ejemplo 8. Los datos que se presentan a continuación (Dobson 2002) corresponden a un estudio en el que semillas genéticamente iguales son asignadas aleatoriamente, bien a un entorno enriquecido nutricionalmente (tratamiento), bien a condiciones estándar (control). Una vez han crecido todas las plantas, se recolectan, secan y pesan. El interés es investigar el efecto del tratamiento utilizado sobre le peso seco (en gramos) de las plantas en cuestión.

peso <- c(4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14, 4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69)
entorno <- c(rep("tratamiento",10),rep("control",10))
ejemplo08 <- data.frame(peso,entorno)

Antes de presentar el modelo, veamos como representar gráficamente los datos a través del diagrama de cajas.

# Diagrama de cajas
ggplot(ejemplo08,aes(x = entorno, y = peso)) + 
  geom_boxplot() + 
  theme_bw()

Visualmente se aprecia una media de peso superior en aquellas plantas que han crecido en un ambiente enriquecido (tratamiento). El modelo asociado con este experimento vendrá dado por:

\[\mu_i = \alpha_0 + \alpha_i + \epsilon_i, \text{ para } i = \text{ \{tratado, control\}} \text{ con } \alpha_{control} = 0\] De esta forma la media para cada nivel del factor tendrá la expresión: \[\left\{\begin{array}{ll} \widehat{\mu}_{tratado} & = \alpha_0 + \alpha_{tratado}\\ \widehat{\mu}_{control} & = \alpha_0 \\ \end{array}\right.\]

Estimando los \(\alpha\) nos resultará posible extraer todas las conclusiones sobre las medias. El modelo en forma reducida se expresa: \[peso \sim entorno\]

Ejemplo 9. Se desea estudiar la fabricación de cuatro tipos de máquinas automáticas en el cortado de piezas de embutido para envasado. Para ello se toman datos del número de envases sin defecto que produce cada máquina durante periodos de una hora.

maquina <- c(rep("M1",4),rep("M2",4),rep("M3",4),rep("M4",4))
produccion <- c(103,115,101,105,109,106,116,124,104,98,117,99,128,117,121,130)
ejemplo09 <- data.frame(maquina,produccion)

Representación gráfica

ggplot(ejemplo09,aes(x = maquina, y = produccion)) + 
  geom_boxplot() + 
  theme_bw()

Se puede apreciar como las máquinas 2 y 4 tienen una mayor producción. De hecho, se observa que la producción con la máquina 4 es muy superior al resto. Tenemos indicios para pensar que las medias pueden resultar distintas.

El modelo asociado con este experimento vendrá dado por:

\[\mu_i = \alpha_0 + \alpha_i + \epsilon_i, \text{ para } i = \text{ \{M1, M2, M3, M4\}} \text{ con } \alpha_{M1} = 0\] De esta forma la media para cada nivel del factor tendrá la expresión: \[\left\{\begin{array}{ll} \widehat{\mu}_{M1} & = \alpha_0 + \alpha_{M1} = \alpha_0\\ \widehat{\mu}_{M2} & = \alpha_0 + \alpha_{M2}\\ \widehat{\mu}_{M3} & = \alpha_0 + \alpha_{M3}\\ \widehat{\mu}_{M4} & = \alpha_0 + \alpha_{M4}\\ \end{array}\right.\]

Estimando los \(\alpha\) nos resultará posible extraer todas las conclusiones sobre las medias. El modelo en forma reducida se expresa: \[produccion \sim maquina\]

Ejemplo 10. Se ha realizado un experimento para comprobar la efectividad de diferentes antídotos (AA, AB, AC y AD) frente a diferentes venenos (VA, VB y VC). Para ello se recoge el tiempo que cada antídoto tarda en hacer efecto para cada veneno.

tiempo <- c(0.31,0.45,0.46,0.43,0.36,0.29,0.40,0.23,0.22,0.21,0.18,0.23,0.82,1.10,0.88,0.72,0.92,0.61,0.49,1.24,0.30,0.37,0.38,0.29,0.43,0.45,0.63,0.76,0.44,0.35,0.31,0.40,0.23,0.25,0.24,0.22,0.45,0.71,0.66,0.62,0.56,1.02,0.71,0.38,0.30,0.35,0.31,0.33)
antidoto <-c(rep("AA",12),rep("AB",12),rep("AC",12),rep("AD",12))
veneno <- c(rep("VA",4),rep("VB",4),rep("VC",4),rep("VA",4),rep("VB",4),rep("VC",4),rep("VA",4),rep("VB",4),rep("VC",4),rep("VA",4),rep("VB",4),rep("VC",4))
ejemplo10 <- data.frame(tiempo,antidoto,veneno)

En este caso nos enfrentamos a un problema donde tenemos dos variables predictoras de tipo categórico. Veamos como representar la información: realizamos el gráfico de cajas y el gráfico de interacción de medias (para ver como varaína las medias para cada combinación de los factores)

# Diagrama de cajas
ggplot(ejemplo10,aes(x = antidoto, y = tiempo, color = veneno)) + geom_boxplot() + theme_bw()

# Gráfico de interacción de medias
ggplot(ejemplo10, aes(x = antidoto, y = tiempo, group = veneno, color = veneno)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line") +
  theme_bw()

En el gráfico se puede apreciar que antídoto es más efectivo para cada tipo de veneno. demás se puede ver que las curvas del tiempo de reacción (líneas de colores) no tiene el mismo comportamiento lo que indica habitualmente que hay cierto efecto de interacción entre los factores considerados, es decir, existe dos combinaciones antídoto - veneno que presentan resultados significativos (medias distintas).

El modelo asociado con este experimento vendrá dado por:

\[\mu_{ij} = \alpha_0 + \alpha_i + \beta_j + \alpha\beta_{ij} +\epsilon_i,\] \[\text{ para } i = \text{ \{AA, AB, AC, AD\}}, j = \text{ \{VA, VB, VC\}} \text{ con } \alpha_{AA} = \beta_{VA} = \alpha\beta_{AA,VA} = 0\] De esta forma la media para cada combinación de los niveles del factor tendrá la expresión: \[\begin{array}{lll} \widehat{\mu}_{AA,VA} & = \alpha_0 + \alpha_{AA} + \beta_{VA} + \alpha\beta_{AA,VA} & = \alpha_0 \\ \widehat{\mu}_{AA,VB} & = \alpha_0 + \alpha_{AA} + \beta_{VB} + \alpha\beta_{AA,VB} & = \alpha_0 + \beta_{VB} + \alpha\beta_{AA,VB} \\ \widehat{\mu}_{AA,VC} & = \alpha_0 + \alpha_{AA} + \beta_{VC} + \alpha\beta_{AA,VC} & = \alpha_0 + \beta_{VC} + \alpha\beta_{AA,VC} \\ \widehat{\mu}_{AB,VA} & = \alpha_0 + \alpha_{AB} + \beta_{VA} + \alpha\beta_{AB,VA} & = \alpha_0 + \alpha_{AB} +\alpha\beta_{AB,VA} \\ \widehat{\mu}_{AB,VB} & = \alpha_0 + \alpha_{AB} + \beta_{VB} + \alpha\beta_{AB,VB} & = \alpha_0 + \alpha_{AB} + \beta_{VB} + \alpha\beta_{AB,VB} \\ \widehat{\mu}_{AB,VC} & = \alpha_0 + \alpha_{AB} + \beta_{VC} + \alpha\beta_{AB,VC} & = \alpha_0 + \alpha_{AB} + \beta_{VC} + \alpha\beta_{AB,VC} \\ \widehat{\mu}_{AC,VA} & = \alpha_0 + \alpha_{AC} + \beta_{VA} + \alpha\beta_{AC,VA} & = \alpha_0 + \alpha_{AC} + \alpha\beta_{AC,VA}\\ \widehat{\mu}_{AC,VB} & = \alpha_0 + \alpha_{AC} + \beta_{VB} + \alpha\beta_{AC,VB} & = \alpha_0 + \alpha_{AC} + \beta_{VB} + \alpha\beta_{AC,VB} \\ \widehat{\mu}_{AC,VC} & = \alpha_0 + \alpha_{AC} + \beta_{VC} + \alpha\beta_{AC,VC} & = \alpha_0 + \alpha_{AC} + \beta_{VC} + \alpha\beta_{AC,VC} \\ \widehat{\mu}_{AD,VA} & = \alpha_0 + \alpha_{AD} + \beta_{VA} + \alpha\beta_{AD,VA} & = \alpha_0 + \alpha_{AD} + \alpha\beta_{AD,VA} \\ \widehat{\mu}_{AD,VB} & = \alpha_0 + \alpha_{AD} + \beta_{VB} + \alpha\beta_{AD,VB} & = \alpha_0 + \alpha_{AD} + \beta_{VB} + \alpha\beta_{AD,VB} \\ \widehat{\mu}_{AD,VC} & = \alpha_0 + \alpha_{AD} + \beta_{VC} + \alpha\beta_{AD,VC} & = \alpha_0 + \alpha_{AD} + \beta_{VC} + \alpha\beta_{AD,VC} \\ \end{array}\]

Estimando los parámetros del modelo podremos realizar el correspondiente análisis de las medias. El modelo en forma reducida se expresa: \[tiempo \sim antidoto + veneno + antidoto:veneno\]

Estimación del modelo

La estimación del modelo se realiza de forma similar a los modelos de regresión. Utilizamos la función lm() especificando el modelo en su forma reducida.

Ejemplo 8

# Ajuste del modelo
ajuste <- lm(peso ~ entorno,data = ejemplo08)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef

El modelo ajustado viene dado por:

\[\left\{\begin{array}{ll} \widehat{\mu}_{tratamiento} & = 4.661 + 0.371 = 5.032\\ \widehat{\mu}_{control} & = 4.661 \\ \end{array}\right.\] En este caso no se interpretan las significatividades de cada coeficiente en el modelo, ya que al contrario de los que ocurría en los modelos de regresión, donde cada coeficiente representaba el efecto de una predictora, en los modelos ANOVA el efecto de la predictora va asociado con todos los coeficientes estimados, ya que cada uno hace referencia a un nivel del factor.

Podemos especificar la estimación de cada nivel del factor de forma directa si eliminamos de la forma reducida del modelo el efecto común a todos los niveles. En este caso si obtendremos los intervalos de confianza, ya que trabajamos directamente con las medias de cada grupo

# Ajuste del modelo eliminando efecto común
ajuste <- lm(peso ~ entorno - 1,data = ejemplo08)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef
# Intervalo de confianza
confint(ajuste)
##                      2.5 %  97.5 %
## entornocontrol     4.19834 5.12366
## entornotratamiento 4.56934 5.49466

Se puede comprobar que los resultados son los mismos que obtuvimos anteriormente. La principal ventaja de esta forma de proceder es que podemos obtener los intervalos de confianza de cada nivel, lo que nos permite compararlos de forma rápida. En este caso podemos ver que ambos intervalos comparten valores comunes (los intervalos no son disjuntos), lo que es un indicativo que las diferencias de ambos grupos no se pueden considerar distintas. Este resultado se podrá verificar en la bondad del ajuste del modelo, siempre y cuando se cumpla con las hipótesis.

Ejemplo 9

El ajuste para el modelo propuesto para estos datos viene dado por:

# Ajuste del modelo
ajuste <- lm(produccion ~ maquina, data = ejemplo09)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef

Las estimaciones de la medias de producción para cada máquina se obtienen como (recordemos que \(\alpha_{M1} = 0\):

\[\left\{\begin{array}{ll} \widehat{\mu}_{M1} & = \alpha_0 \\ \widehat{\mu}_{M2} & = \alpha_0 + \alpha_{M2}\\ \widehat{\mu}_{M3} & = \alpha_0 + \alpha_{M3}\\ \widehat{\mu}_{M4} & = \alpha_0 + \alpha_{M4}\\ \end{array}\right.\]

Siempre hay que comprobar que efectos son los que asimilamos a cero para no equivocarnos en la especificación del modelo. Para estimar las medias directamente ajustamos el modelo sin efecto común:

# Ajuste del modelo
ajuste <- lm(produccion ~ maquina -1, data = ejemplo09)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef
# Intervalo de confianza
confint(ajuste)
##               2.5 %   97.5 %
## maquinaM1  97.99607 114.0039
## maquinaM2 105.74607 121.7539
## maquinaM3  96.49607 112.5039
## maquinaM4 115.99607 132.0039

El hecho más relevante es que el intervalo de confianza para la producción de la máquina M4 es claramente distinto de las producciones de M1 y M3 (intervalos disjuntos), indicando que las medias se comportan de forma diferente. ¿Qué conclusión podemos extraer?

Ejemplo 10

El ajuste para el modelo propuesto para estos datos viene dado por:

# Ajuste del modelo
ajuste <- lm(tiempo ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef

Siguiendo la tabla tabla de estimación visto al inicio del ejemplo, la estimación de las medias se lleva a cabo con:

\[\begin{array}{lll} \widehat{\mu}_{AA,VA} & = \alpha_0 + \alpha_{AA} + \beta_{VA} + \alpha\beta_{AA,VA} & = \alpha_0 \\ \widehat{\mu}_{AA,VB} & = \alpha_0 + \alpha_{AA} + \beta_{VB} + \alpha\beta_{AA,VB} & = \alpha_0 + \beta_{VB} + \alpha\beta_{AA,VB} \\ \widehat{\mu}_{AA,VC} & = \alpha_0 + \alpha_{AA} + \beta_{VC} + \alpha\beta_{AA,VC} & = \alpha_0 + \beta_{VC} + \alpha\beta_{AA,VC} \\ \widehat{\mu}_{AB,VA} & = \alpha_0 + \alpha_{AB} + \beta_{VA} + \alpha\beta_{AB,VA} & = \alpha_0 + \alpha_{AB} +\alpha\beta_{AB,VA} \\ \widehat{\mu}_{AB,VB} & = \alpha_0 + \alpha_{AB} + \beta_{VB} + \alpha\beta_{AB,VB} & = \alpha_0 + \alpha_{AB} + \beta_{VB} + \alpha\beta_{AB,VB} \\ \widehat{\mu}_{AB,VC} & = \alpha_0 + \alpha_{AB} + \beta_{VC} + \alpha\beta_{AB,VC} & = \alpha_0 + \alpha_{AB} + \beta_{VC} + \alpha\beta_{AB,VC} \\ \widehat{\mu}_{AC,VA} & = \alpha_0 + \alpha_{AC} + \beta_{VA} + \alpha\beta_{AC,VA} & = \alpha_0 + \alpha_{AC} + \alpha\beta_{AC,VA}\\ \widehat{\mu}_{AC,VB} & = \alpha_0 + \alpha_{AC} + \beta_{VB} + \alpha\beta_{AC,VB} & = \alpha_0 + \alpha_{AC} + \beta_{VB} + \alpha\beta_{AC,VB} \\ \widehat{\mu}_{AC,VC} & = \alpha_0 + \alpha_{AC} + \beta_{VC} + \alpha\beta_{AC,VC} & = \alpha_0 + \alpha_{AC} + \beta_{VC} + \alpha\beta_{AC,VC} \\ \widehat{\mu}_{AD,VA} & = \alpha_0 + \alpha_{AD} + \beta_{VA} + \alpha\beta_{AD,VA} & = \alpha_0 + \alpha_{AD} + \alpha\beta_{AD,VA} \\ \widehat{\mu}_{AD,VB} & = \alpha_0 + \alpha_{AD} + \beta_{VB} + \alpha\beta_{AD,VB} & = \alpha_0 + \alpha_{AD} + \beta_{VB} + \alpha\beta_{AD,VB} \\ \widehat{\mu}_{AD,VC} & = \alpha_0 + \alpha_{AD} + \beta_{VC} + \alpha\beta_{AD,VC} & = \alpha_0 + \alpha_{AD} + \beta_{VC} + \alpha\beta_{AD,VC} \\ \end{array}\]

Ajustamos el modelo sin efecto común para obtener las medias de forma directa:

# Ajuste del modelo
ajuste <- lm(tiempo ~ antidoto + veneno + antidoto:veneno -1, data = ejemplo10)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(ajuste)
fitcoef
# Intervalo de confianza
confint(ajuste)
##                          2.5 %      97.5 %
## antidotoAA           0.2613254  0.56367465
## antidotoAB           0.7288254  1.03117465
## antidotoAC           0.4163254  0.71867465
## antidotoAD           0.4588254  0.76117465
## venenoVB            -0.3062932  0.12129324
## venenoVC            -0.4162932  0.01129324
## antidotoAB:venenoVB -0.2748493  0.32984930
## antidotoAC:venenoVB -0.4023493  0.20234930
## antidotoAD:venenoVB -0.1523493  0.45234930
## antidotoAB:venenoVC -0.6448493 -0.04015070
## antidotoAC:venenoVC -0.4323493  0.17234930
## antidotoAD:venenoVC -0.3873493  0.21734930

Al tener muchas combinaciones resulta más difícil comparar los intervalos de confianza. De hecho, en primer lugar deberíamos valorar si el efecto de interacción es significativo o no, es decir, si debe estar presente en el modelo. Para ello utilizaremos los procedimientos de bondad de ajuste y de comparación de modelos que veremos en los puntos siguientes.

Bondad de ajuste del modelo

La bondad de ajuste para este tipo se basa únicamente en el test \(F\) de la regresión. En este caso dicho test valora si existen al menos dos medias, de todas las posibles combinaciones de los niveles de los factores, que podamos considerar distintas.

Analizamos los resultados obtenidos para cada uno de los ejemplos.

Ejemplo 8

En este caso comparamos los resultados con los obtenidos en el tema anterior cuando ajustábamos un modelo de regresión lineal simple.

En este caso el pvalor resulta no significativo indicando que no podemos considerar como distintas las medias de peso de las plantas que han sido tratadas y aquellas que no lo han sido. El tratamiento utilizado no tiene efecto sobre el crecimiento de las plantas. Tan sólo nos restaría realizar la validación del modelo para establecer nuestras conclusiones finales.

Ejemplo 9

# Ajuste del modelo
ajuste <- lm(produccion ~ maquina, data = ejemplo09)
gof <- glance(ajuste)
gof

El pvalor resulta significativo indicando que existen al menos dos medias que son estadísticamente distintas. Tenemos por tanto que existen al menos dos máquinas cuyo número medio de envases con defectos pueden considerarse distintos. El problema de este test es que no nos proporciona que medias son las que resultan diferentes. Para resolver este problema debemos realizar las comparaciones múltiples entre todas las medias posibles. Este procedimiento se realizar con la función TukeyHSD() que os proporciona los pvalores asociados a todas las comparaciones dos a dos de todas las medias.

# Ajuste del modelo para que la función TukeyHSD lo pueda interpretar
ajuste <- aov(produccion ~ maquina, data = ejemplo09)
TukeyHSD(ajuste)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = produccion ~ maquina, data = ejemplo09)
## 
## $maquina
##        diff        lwr       upr     p adj
## M2-M1  7.75  -7.673887 23.173887 0.4716327
## M3-M1 -1.50 -16.923887 13.923887 0.9911761
## M4-M1 18.00   2.576113 33.423887 0.0210561
## M3-M2 -9.25 -24.673887  6.173887 0.3283778
## M4-M2 10.25  -5.173887 25.673887 0.2508228
## M4-M3 19.50   4.076113 34.923887 0.0126965

La columna p adj proporciona la significativa del contraste de comparación de medias que aparece en la primera columna. Sólo resultan significativos los contrastes entre M4 y M1, y entre M4 y M3. ¿Cómo interpretamos este resultado? ¿Podemos construir una tabla de comparación para visualizar las diferencias?

Ejemplo 10

Para este tipo de datos (ANOVA de dos vías) el test \(F\) solo nos indicara si hay diferencias entre las medias de al menos dos combinaciones de los niveles de ambos factores. Cuando el número de niveles es alto el número de combinaciones es muy alto, y puede resultar complicado realizar todas las comparaciones múltiples.

# Ajuste del modelo
ajuste <- lm(tiempo ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
gof <- glance(ajuste)
gof

¿Qué conclusiones podemos extraer de este análisis?

# Ajuste del modelo
ajuste <- aov(tiempo ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
TukeyHSD(ajuste)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = tiempo ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
## 
## $antidoto
##              diff         lwr        upr     p adj
## AB-AA  0.36250000  0.19858517  0.5264148 0.0000046
## AC-AA  0.07833333 -0.08558150  0.2422482 0.5769172
## AD-AA  0.21916667  0.05525184  0.3830815 0.0050213
## AC-AB -0.28416667 -0.44808150 -0.1202518 0.0002320
## AD-AB -0.14333333 -0.30724816  0.0205815 0.1045211
## AD-AC  0.14083333 -0.02308150  0.3047482 0.1136890
## 
## $veneno
##            diff        lwr         upr     p adj
## VB-VA -0.073125 -0.2019588  0.05570882 0.3580369
## VC-VA -0.341875 -0.4707088 -0.21304118 0.0000005
## VC-VB -0.268750 -0.3975838 -0.13991618 0.0000325
## 
## $`antidoto:veneno`
##                diff         lwr          upr     p adj
## AB:VA-AA:VA  0.4675  0.09956481  0.835435189 0.0041185
## AC:VA-AA:VA  0.1550 -0.21293519  0.522935189 0.9392346
## AD:VA-AA:VA  0.1975 -0.17043519  0.565435189 0.7669184
## AA:VB-AA:VA -0.0925 -0.46043519  0.275435189 0.9989703
## AB:VB-AA:VA  0.4025  0.03456481  0.770435189 0.0219844
## AC:VB-AA:VA -0.0375 -0.40543519  0.330435189 0.9999999
## AD:VB-AA:VA  0.2550 -0.11293519  0.622935189 0.4197520
## AA:VC-AA:VA -0.2025 -0.57043519  0.165435189 0.7390988
## AB:VC-AA:VA -0.0775 -0.44543519  0.290435189 0.9998030
## AC:VC-AA:VA -0.1775 -0.54543519  0.190435189 0.8638430
## AD:VC-AA:VA -0.0900 -0.45793519  0.277935189 0.9991978
## AC:VA-AB:VA -0.3125 -0.68043519  0.055435189 0.1614923
## AD:VA-AB:VA -0.2700 -0.63793519  0.097935189 0.3373794
## AA:VB-AB:VA -0.5600 -0.92793519 -0.192064811 0.0003192
## AB:VB-AB:VA -0.0650 -0.43293519  0.302935189 0.9999649
## AC:VB-AB:VA -0.5050 -0.87293519 -0.137064811 0.0014859
## AD:VB-AB:VA -0.2125 -0.58043519  0.155435189 0.6803714
## AA:VC-AB:VA -0.6700 -1.03793519 -0.302064811 0.0000138
## AB:VC-AB:VA -0.5450 -0.91293519 -0.177064811 0.0004874
## AC:VC-AB:VA -0.6450 -1.01293519 -0.277064811 0.0000282
## AD:VC-AB:VA -0.5575 -0.92543519 -0.189564811 0.0003426
## AD:VA-AC:VA  0.0425 -0.32543519  0.410435189 0.9999996
## AA:VB-AC:VA -0.2475 -0.61543519  0.120435189 0.4639955
## AB:VB-AC:VA  0.2475 -0.12043519  0.615435189 0.4639955
## AC:VB-AC:VA -0.1925 -0.56043519  0.175435189 0.7934572
## AD:VB-AC:VA  0.1000 -0.26793519  0.467935189 0.9979346
## AA:VC-AC:VA -0.3575 -0.72543519  0.010435189 0.0632967
## AB:VC-AC:VA -0.2325 -0.60043519  0.135435189 0.5563431
## AC:VC-AC:VA -0.3325 -0.70043519  0.035435189 0.1083714
## AD:VC-AC:VA -0.2450 -0.61293519  0.122935189 0.4790943
## AA:VB-AD:VA -0.2900 -0.65793519  0.077935189 0.2434424
## AB:VB-AD:VA  0.2050 -0.16293519  0.572935189 0.7247659
## AC:VB-AD:VA -0.2350 -0.60293519  0.132935189 0.5407312
## AD:VB-AD:VA  0.0575 -0.31043519  0.425435189 0.9999898
## AA:VC-AD:VA -0.4000 -0.76793519 -0.032064811 0.0233746
## AB:VC-AD:VA -0.2750 -0.64293519  0.092935189 0.3120799
## AC:VC-AD:VA -0.3750 -0.74293519 -0.007064811 0.0424758
## AD:VC-AD:VA -0.2875 -0.65543519  0.080435189 0.2541042
## AB:VB-AA:VB  0.4950  0.12706481  0.862935189 0.0019557
## AC:VB-AA:VB  0.0550 -0.31293519  0.422935189 0.9999935
## AD:VB-AA:VB  0.3475 -0.02043519  0.715435189 0.0788658
## AA:VC-AA:VB -0.1100 -0.47793519  0.257935189 0.9953187
## AB:VC-AA:VB  0.0150 -0.35293519  0.382935189 1.0000000
## AC:VC-AA:VB -0.0850 -0.45293519  0.282935189 0.9995273
## AD:VC-AA:VB  0.0025 -0.36543519  0.370435189 1.0000000
## AC:VB-AB:VB -0.4400 -0.80793519 -0.072064811 0.0085085
## AD:VB-AB:VB -0.1475 -0.51543519  0.220435189 0.9562021
## AA:VC-AB:VB -0.6050 -0.97293519 -0.237064811 0.0000887
## AB:VC-AB:VB -0.4800 -0.84793519 -0.112064811 0.0029419
## AC:VC-AB:VB -0.5800 -0.94793519 -0.212064811 0.0001810
## AD:VC-AB:VB -0.4925 -0.86043519 -0.124564811 0.0020941
## AD:VB-AC:VB  0.2925 -0.07543519  0.660435189 0.2330953
## AA:VC-AC:VB -0.1650 -0.53293519  0.202935189 0.9103689
## AB:VC-AC:VB -0.0400 -0.40793519  0.327935189 0.9999998
## AC:VC-AC:VB -0.1400 -0.50793519  0.227935189 0.9694983
## AD:VC-AC:VB -0.0525 -0.42043519  0.315435189 0.9999960
## AA:VC-AD:VB -0.4575 -0.82543519 -0.089564811 0.0053752
## AB:VC-AD:VB -0.3325 -0.70043519  0.035435189 0.1083714
## AC:VC-AD:VB -0.4325 -0.80043519 -0.064564811 0.0103293
## AD:VC-AD:VB -0.3450 -0.71293519  0.022935189 0.0832419
## AB:VC-AA:VC  0.1250 -0.24293519  0.492935189 0.9868596
## AC:VC-AA:VC  0.0250 -0.34293519  0.392935189 1.0000000
## AD:VC-AA:VC  0.1125 -0.25543519  0.480435189 0.9943551
## AC:VC-AB:VC -0.1000 -0.46793519  0.267935189 0.9979346
## AD:VC-AB:VC -0.0125 -0.38043519  0.355435189 1.0000000
## AD:VC-AC:VC  0.0875 -0.28043519  0.455435189 0.9993811

¿Cómo interpretamos los resultados de estas tablas? ¿Cuál de ellas nos interesa?

Comparación de modelos

La comparación de modelos se basa en determinar la capacidad explicativa de cada uno de los efectos presentes en el modelo. Mientras que en los modelos de regresión esos efectos van asociados con un único coeficiente del modelo (pendiente asociada a la predictora), en los modelos ANOVA el efecto de un factor va asociado con todos los incrementos establecidos. En realidad tendremos un efecto significativo cuando detectemos diferencias entre al menos dos medias de la respuesta para los diferentes niveles del factor. Utilizamos la tabla ANOVA y los procedimiento de selección automática para alcanza el mejor modelo en cada caso.

Análisis de los ejemplos

Pasamos al estudio de cada uno de los ejemplos

Ejemplo 8

La tabla ANOVA para este modelo se obtiene mediante

El test \(F\) muestra que no hay diferencias entre el peso de las plantas que han sido tratadas y las que no. Esto ya lo habíamos visto antes. Eso indica que el tratamiento no influye en el peso final de las plantas.

Ejemplo 9

La tabla ANOVA para este modelo se obtiene mediante

# Ajuste del modelo
ajuste <- lm(produccion ~ maquina, data = ejemplo09)
anova(ajuste)

La tabla ANOVA resultante muestra que existen al menos dos máquinas cuyo número medio de defectos se puede considerar como distinto. Hay al menos una máquina que funciona mejor que otra. Si recuperamos los test de comparación múltiples podemos ver que la máquina M4 funciona de forma distinta a como lo hacen la M1 y la M3.

No nos planteamos la selección de efectos o variables dado que sólo tenemos una predictora.

Ejemplo 10

La tabla ANOVA para este modelo se obtiene mediante

# Ajuste del modelo
ajuste <- lm(tiempo ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
anova(ajuste)

La tabla ANOVA muestra que el efecto asociado con la interacción no resulta significativo, lo que implica que no hay diferencias por esa combinación, y que podríamos ajustar un nuevo modelo sin dicho efecto. Construimos el nuevo modelo y comparamos el modelo saturado con el nuevo modelo sin interacción.

# Ajuste del modelo sin interacción
ajuste1 <- lm(tiempo ~ antidoto + veneno , data = ejemplo10)
# Comparamos los modelos anidados 
anova(ajuste, ajuste1)

La comparación de ambos modelos produce en pvalor (Pr(>F)) no significativo lo que indica que el efecto de interacción es prescindible y puede eliminarse del modelo, ya que el modelo sin interacción puede considerarse igual (hablando en términos de capacidad explicativa) al modelo que si tiene interacción. El modelo resultante (a falta de comprobar las hipótesis del modelo) es: \[tiempo \sim antidoto + veneno \] ¿Cómo interpretamos ese modelo? ¿Qué significa que no hay interacción?

Veamos el proceso de selección de variables

## Start:  AIC=-172.52
## tiempo ~ antidoto + veneno + antidoto:veneno
## 
##                   Df Sum of Sq    RSS     AIC
## <none>                         0.8001 -172.52
## - antidoto:veneno  6   0.25027 1.0504 -171.46
## 
## Call:
## lm(formula = tiempo ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
## 
## Coefficients:
##         (Intercept)           antidotoAB           antidotoAC  
##              0.4125               0.4675               0.1550  
##          antidotoAD             venenoVB             venenoVC  
##              0.1975              -0.0925              -0.2025  
## antidotoAB:venenoVB  antidotoAC:venenoVB  antidotoAD:venenoVB  
##              0.0275              -0.1000               0.1500  
## antidotoAB:venenoVC  antidotoAC:venenoVC  antidotoAD:venenoVC  
##             -0.3425              -0.1300              -0.0850

El proceso de selección no prescinde de la interacción (que es el efecto más complejo y por tanto e primer que debemos intentar eliminar) y nos quedaríamos con el modelo completo. Esta diferencia entre ambos procedimientos puede llegar a ser habitual. EN ese caso siempre debemos elegir el modelo más sencillo.

Diagnóstico del modelo

Los procedimientos de diagnóstico son los mismas que en los modelos de regresión salvo por el hecho de que la igualdad de varianzas se debe cumplir para los residuos de cada uno de los grupos determinados por los niveles de los factores, es decir, varianzas residuales iguales.

En el caso de la hipótesis de normalidad utilizaremos el test de shapiro para cada grupo, mientras que para la hipótesis de igualdad de varianzas utilizamos el test de Barlett. Obviaremos los procedimientos gráficos

Ejemplo 9.

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(produccion ~ maquina, data = ejemplo09)
fitinf <- fortify(ajuste)

Tests para verificar las hipótesis

# Varianza constante
bartlett.test(produccion ~ maquina, data = ejemplo09)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  produccion by maquina
## Bartlett's K-squared = 0.52592, df = 3, p-value = 0.9132
# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$maquina == "M1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$maquina == "M1"]
## W = 0.85422, p-value = 0.2401
shapiro.test(fitinf$.stdresid[fitinf$maquina == "M2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$maquina == "M2"]
## W = 0.95045, p-value = 0.7189
shapiro.test(fitinf$.stdresid[fitinf$maquina == "M3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$maquina == "M3"]
## W = 0.84246, p-value = 0.2027
shapiro.test(fitinf$.stdresid[fitinf$maquina == "M4"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$maquina == "M4"]
## W = 0.92714, p-value = 0.5777

Tanto el test de homogeneidad como el de normalidad proporcionan resultados no significativos indicando el cumplimiento de las hipótesis.

Ejemplo 10

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(tiempo ~ antidoto + veneno, data = ejemplo10)
fitinf <- fortify(ajuste)

Tests para verificar las hipótesis

En primer lugar creamos un factor para identificar todos los grupos posibles

# Añadimos un nuevo factor a los resultados con todas las combinaciones de los grupos
fitinf01 <- fitinf %>% mutate(grupo = paste(antidoto,veneno))
# Varianza constante
bartlett.test(tiempo ~ grupo, data = fitinf01)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  tiempo by grupo
## Bartlett's K-squared = 46.058, df = 11, p-value = 3.158e-06
# Normalidad: existen exesivas combinaciones para escribir el test.
# realizamos un análisis grañfico tenienco en cuenta el bajo número de observaciones
ggplot(fitinf01, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  facet_wrap(~ grupo, scales="free_x") +  
  theme_bw()

Se puede ver que no se verifican las hipótesis de homogeneidad de varianzas. En cuanto a la denormalidad al haber tan pocos datos por grupo resulta casi imposible asegurar dicha normalidad. Probamos a realizar un test global de los residuos:

shapiro.test(fitinf01$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf01$.stdresid
## W = 0.92213, p-value = 0.003538

El test global resulta significativo indicando falta de normalidad. Para tratar de corregir este defecto utilizamos las transformaciones de boxcox sobre la respuesta.

# Ajuste del modelo y obtención de los residuos
boxcox(ajuste)

Se puede ver que el intervalo de confianza para la transformación contiene al valor de -1. Esto nos da indicios de que podríamos utilizar la transformación inversa para corregir los problemas detectados de homogeneidad de varianzas. En ese caso habrá que tener en cuenta que todas nuestras conclusiones se basarán en dicha transformación.

Creamos la transformación y comenzamos de nuevo con la construcción del modelo.

# Creamos la variable transformada
ejemplo10 <- ejemplo10 %>% mutate(tiempoinv = 1/tiempo)
# Ajustamos el modelo de nuevo
ajuste <- lm(tiempoinv ~ antidoto + veneno + antidoto:veneno, data = ejemplo10)
# Valoramos el modelo con la tabla anova
anova(ajuste)
# Elinamos efecto de interacción y reajustamos
ajuste <- lm(tiempoinv ~ antidoto + veneno, data = ejemplo10)
anova(ajuste)
# Valoramos el cumplimiento de las hipótesis
fitinf <- fortify(ajuste)
# Añadimos un nuevo factor a los resultados con todas las combinaciones de los grupos
fitinf <- fitinf %>% mutate(grupo = paste(antidoto,veneno))
# Varianza constante
bartlett.test(tiempoinv ~ grupo, data = fitinf)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  tiempoinv by grupo
## Bartlett's K-squared = 10.511, df = 11, p-value = 0.4851
# Normalidad 
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.9789, p-value = 0.5339

Podemos ver como se verifican las hipótesis del modelo, de forma que el modelo resultante viene dado por: \[1/tiempo \sim antidoto + veneno\]

Predicción

El proceso de predicción en este tipo de modelos es muy simple ya que los posibles valores de las predictoras coincide con los posibles niveles o combinación de los niveles del factor o factores. En este caso solo tiene sentido predecir los valores de la media ya que nuestro objetivo es el estudio de dichas medias. Veremos además como representar gráficamente los resultados de la predicción. No predecimos resultados para el ejemplo 8 dado que el modelo no contiene variables predictoras.

Ejemplo 9

Queremos predecir la el número de nevases defectuosos en función de la máquina que produce los cortes:

Predecimos la media:

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(produccion ~ maquina, data = ejemplo09)
# Secuencia de valores de predicción
newdata <- data.frame(maquina = unique(ejemplo09$maquina))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
# Gráfico resultante (predicción y bandas de predicción). Ordenamos las clases de acuerdo al valor de fit
ggplot(newdata, aes(x = fct_reorder(maquina,fit), y= fit)) + 
    geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1) +
    geom_line() +
    geom_point() +
    labs(x = "Máquina", y = "Producción (número de envases sin defecto)") +
    coord_flip() +
    theme_bw()

Podemos ver como la predicción para la máquina 4 es diferente de las predicciones para las máquinas 1 y 3. No hay diferencia entre la M4 y la M2, ni entre la M2 y la M1 o M3. La máquina que podemos indicar como más efectiva es la M4 ya que es la que produce un mayor número de envases sin defectos.

Ejemplo 10

Predecimos la media:

# Ajuste del modelo y obtención de los residuos. Utilizamos el modelo obtenido el apartado de diagnóstico
ajuste <- lm(tiempoinv ~ antidoto + veneno, data = ejemplo10)
# Creamo la combinación de grupos con ambos factores
f1 <- unique(ejemplo10$antidoto)
f2 <- unique(ejemplo10$veneno)
# Generamos el grid. 
newdata <- data.frame(expand.grid(f1,f2))
colnames(newdata) <- c("antidoto","veneno")
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
# Creamos el factor con todas las combinaciones y ordenamos de acuerdo a la media, colocando el veneo en priemr lugar para poder localizar el mejor antídoto
newdata <- newdata %>% mutate(grupo = paste(veneno,antidoto)) 
# Gráfico resultante (predicción y bandas de predicción). Ordenamos las clases de acuerdo al valor de fit y coloreamos por veneno apra identificar el mejor antídoto
ggplot(newdata, aes(x = fct_reorder(grupo,fit), y= fit, colour = veneno)) + 
    geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1) +
    geom_line() +
    geom_point() +
    coord_flip() +
    labs(y = "1/ tiempo de reacción", x = "Combinación antídoto-veneno") +
    theme_bw()

para interpretar este gráfico hay que tener en cuenta que al transformar la respuesta mediante la función inversa, las mejores combinaciones son las que se sitúan en la parte superior del gráfico y no en la parte inferior Las conclusiones que podemos extraer son:

  • Para el veneno VA loa mejores antídotos son el AA y el AC con tiempos de reacción similares (intervalos de predicción no disjuntos), y distintos de los que proporcionan el antídoto AD y el AB.
  • Para el veneno VB loa mejores antídotos son el AA y el AC con tiempos de reacción similares, y distintos de los que proporcionan el antídoto AD y el AB.
  • Para el veneno VC loa mejores antídotos son el AA y el AC con tiempos de reacción similares, y distintos de los que proporcionan el antídoto AD y el AB.

A la hora de elegir deberíamos hacerlo entre AA y AC para cualquiera de los venenos.

Además podemos ver que el veneno VA es el más resistente a los antídotos ya que tiene mayores tiempo de reacción (sitúa más intervalos de predicción en la zona baja). El menos resistente sería el veneno VC.

Bibliografía

Dobson, A. 2002. An Introduction to Generalized Linear Models. CRC Press.

Copyright © 2018 Javier Morales. Universidad Miguel Hernández de Elche.