1 Modelos de análisis de la varianza

Los modelos ANOVA surgen cuando la variable o variables predictoras son de tipo categórico, es decir, factores con diferentes niveles de clasificación, de forma que cada sujeto es medido (respecto de la respuesta) para una combinación específica de los factores considerados. 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 o combinaciones de los factores. 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 del conjunto de observaciones para cada nivel del factor. Además, veremos como estos modelos pueden ser descritos en términos de un modelo lineal de regresión lo que nos permite utilizar parte de los procesos de descripción, análisis y diagnóstico utilizamos en los capítulos anteriores. En cada punto de este capítulo estudiaremos las similitudes y diferencias con los modelos de regresión.

1.1 Bancos de datos

Veamos los diferentes bancos de datos que iremos analizando a lo largo de la unidad.

1.1.1 Ejemplo 1

Datos de Insecticidas. Este ejemplo contiene los datos de un experimento agronómico para conocer el efecto de diferentes insecticidas (spray) sobre el número de insectos vivos (count) tras un periodo de tratamiento.

Cargamos los datos y realizamos el gráfico descriptivo:

Se aprecia como el número de insectos vivos es superior cuando usamos los sprays A, B, y F (con medias muy similares). El resto de sprays muestran (medias) una supervivencia inferior. Se observan como dos grupos de sprays, uno con medias altas y otro con medias bajas.

1.1.2 Ejemplo 2

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

Cargamos los datos y realizamos el gráfico descriptivo:

Cada máquina muestra un nivel de envases sin defecto diferente. La máquina M4 es la que muestra mejores resultados y la M3 los peores.

1.1.3 Ejemplo 3

Datos de venenos. 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 de reacción (tiempo) que cada antídoto tarda en hacer efecto para cada veneno.

Cargamos los datos y realizamos el gráfico descriptivo:

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).

1.2 Modelos ANOVA 1 vía

El modelo ANOVA de una vía se plantea cuando tenemos una única variable predictora de tipo categórico y una respuesta de tipo numérico. Dicho modelo se puede describir de forma general, para una muestra de tamaño \(n\), como:

  • Una variable respuesta, \(Y\), de tipo numérico con observaciones \(y_1, \ldots ,y_n\).
  • Una variable predictora, \(F\), de tipo categórico con \(I\) grupos o niveles distintos de tamaños muestrales \(n_1, n_2, \ldots , n_I\), de forma que \(n = n_1 + n_2 + \ldots + n_I\), y el vector de observaciones de la respuesta se puede escribir como:

\[ y_{11},\ldots,y_{1n_1},y_{21},\ldots,y_{2n_2},\ldots,y_{I1},\ldots,y_{In_I} \] donde el primer subíndice indica el nivel del factor y el segundo la posición dentro del conjunto de datos de dicho nivel del factor.

  • Conjunto \(\mu_i\) de medias de todas las observaciones de la respuesta asociadas con el nivel \(i\) del factor, es decir:

\[ \mu_i = \frac{\sum_{j = 1}^{n_j} y_{ij}}{n_i}; \text{ i = 1, 2,..., I} \]

  • Media global de la respuesta, \(\mu\), que se puede obtener como:

\[ \mu = \frac{\sum_{j = 1}^{I} \mu_{j}}{I} \]

  • Incrementos, \(\alpha_i\), de cada una de las medias de cada grupo con respecto a la media global, es decir:

\[ \alpha_i= \mu - \mu_i; \text{ i = 1, 2,..., I} \]

El objetivo básico en este tipo de modelos es la comparación de las medias \(\mu_i\) para detectar posibles diferencias entre os niveles del factor, es decir, plateamos el contraste de comparación de medias: \[ \left\{ \begin{array}{ll} H_0: & \mu_1 = \mu_2 = \ldots = \mu_I\\ H_a: & \mbox{Existen dos grupos al menos con medias distintas} \end{array} \right. \]

Este contraste es equivalente al que se puede escribir en términos de los incrementos de las medias: \[ \left\{ \begin{array}{ll} H_0: & \alpha_1 = \alpha_2 = \ldots = \alpha_I = 0\\ H_a: & \mbox{Al menos hay un incremento distinto de cero} \end{array} \right. \]

que tiene una estructura similar al test \(F\) de la regresión, donde queremos contrastar si existe algún coeficiente distinto de cero. Esto implica que el modelo ANOVA se puede escribir en términos de un modelo lineal sin más que considerar un conjunto de variables ficticias \(X_1,...,X_{I}\) que toman los valores siguientes:

\[ X_i= \left\{ \begin{array}{ll} 1 & \mbox{si la respuesta pertenece al grupo i}\\ 0 & \mbox{en otro caso} \end{array} \right. \]

y un efecto común para todos los niveles del factor que viene dado por \(\mu\) que denotamos por \(\alpha_0\). La única diferencia entre ambas estructuras de comparación es que la referida a las medias \(\mu_i\) tiene \(I\) parámetros, mientras que la de los incrementos tiene \(I+1\) parámetros. Para hacer equivalentes ambas estructuras añadimos la conocida como restricción de identificabilidad que establece de partida que uno de los incrementos debe ser igual a cero, es decir, que una media de un grupo coincide con la media global. Este es el nivel de referencia que utilizamos como base de comparación con el resto de niveles o grupos del factor.

En esta situación si consideramos \(Y_i = \{y_{i1},...,y_{in_i}\}\), \(1_{n_i} = \{1,...,1\}\) un vector de \(n_i\) unos, \(0_{n_i} = \{0,...,0\}\) un vector de \(n_i\) ceros, para cada uno de los niveles \(i\) del factor, \(\mu = \alpha_0\), y la restricción de identificabilidad \(\alpha_I = 0\) (podríamos elegir cualquier otro \(\alpha_j =0\) y obtendríamos el mismo modelo), podemos escribir el modelo ANOVA como:

\[ Y = \left( \begin{array}{c} Y_1 \\ Y_2 \\ \ldots \\ Y_I \\ \end{array} \right) = \left( \begin{array}{ccccc} 1_{n_1} & 1_{n_1} & 0_{n_1} & \ldots & 0_{n_1} \\ 1_{n_2} & 0_{n_2} & 1_{n_2} & \ldots & 0_{n_2} \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1_{n_I} & 0_{n_I} & 0_{n_I} &\ldots & 1_{n_I} \\ \end{array} \right) \left( \begin{array}{c} \alpha_0 \\ \alpha_1 \\ \alpha_2 \\ \ldots \\ \alpha_{I_1} \\ \end{array} \right) + \left( \begin{array}{c} e_1 \\ e_2 \\ \ldots \\ e_n \\ \end{array} \right) = X\beta + \epsilon \]

o equivalentemente en términos de las medias como:

\[\begin{equation} \mu_i = \alpha_0 + \alpha_i + \epsilon_i, \text{ para } i = 1, 2,...,k, \text{ con } \alpha_k = 0 \end{equation}\]

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. En esta situación el contraste sobre los incrementos es equivalente al test \(F\) de regresión de los modelos RLM.

1.3 Modelos ANOVA dos vías

El modelo ANOVA de dos vías se presenta cuando consideramos dos variables predictoras de tipo categórico. Si tenemos dos factores \(F_1\) y \(F_2\) con \(I\) y \(J\) niveles respectivamente (y utilizando lo visto en el punto anterior), el modelo expresado en términos de las medias de las combinaciones de los diferentes niveles de los factores considerados viene dado por:

\[\begin{equation} \mu_{ij} = \alpha_0 + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ij}, \ \mbox{para} \ i = 1,\ldots,I; \ j = 1,\ldots,J \ \mbox{con} \ \alpha_1 = \beta_1 = \alpha\beta_{11} = 0 \end{equation}\]

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 \(F_1\) 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.

Dado que tenemos dos factores se hace necesario añadir más restricciones de identificabildiad. Añadimos una por cada factor y la que aparece de forma inmediata en la interacción de los dos efectos que estamos tomando como cero.

En este caso no escribimos la matriz de diseño del modelo pero se puede escribir considerado la matriz de coeficientes de dimensiones \(IJ \times IJ\) que surge para el modelo ANOVA de dos vías. Bastará con sustituir los unos y ceros por vectores de unos y ceros de dimensiones adecuadas en función de los tamaños de cada uno de los niveles de los dos factores considerados. Dicha matriz de coeficientes se puede extraer a partir de las ecuaciones completas del modelo ANOVA de dos vías (teniendo en cuenta las restricciones de identificabilidad) dadas por:

\[\begin{array}{lll} \mu_{11} & = \alpha_0 + \alpha_{1} + \beta_{1} + \alpha\beta_{11} + \epsilon_{11}& = \alpha_0 + \epsilon_{11}\\ \mu_{12} & = \alpha_0 + \alpha_{1} + \beta_{2} + \alpha\beta_{12} + \epsilon_{12}& = \alpha_0 + \beta_{2} + \alpha\beta_{12} + \epsilon_{12}\\ \ldots & \ldots &\ldots\\ \mu_{1J} & = \alpha_0 + \alpha_{1} + \beta_{J} + \alpha\beta_{1J} + \epsilon_{1J}& = \alpha_0 + \beta_{J} + \alpha\beta_{1J} + \epsilon_{1J}\\ \mu_{21} & = \alpha_0 + \alpha_{2} + \beta_{1} + \alpha\beta_{21} + \epsilon_{21}& = \alpha_0 + \alpha_{2} + \alpha\beta_{21} +\epsilon_{11}\\ \mu_{22} & = \alpha_0 + \alpha_{2} + \beta_{2} + \alpha\beta_{22} + \epsilon_{22}& = \alpha_0 + \alpha_{2} + \beta_{2} + \alpha\beta_{22} + \epsilon_{22}\\ \ldots & \ldots &\ldots\\ \mu_{2J} & = \alpha_0 + \alpha_{2} + \beta_{J} + \alpha\beta_{2J} + \epsilon_{2J}& = \alpha_0 + \alpha_{2} + \beta_{J} + \alpha\beta_{2J} + \epsilon_{2J}\\ \ldots & \ldots &\ldots\\ \ldots & \ldots &\ldots\\ \ldots & \ldots &\ldots\\ \mu_{I1} & = \alpha_0 + \alpha_{I} + \beta_{1} + \alpha\beta_{I1} + \epsilon_{I1}& = \alpha_0 + \alpha_{I} + \alpha\beta_{I1} + \epsilon_{I1}\\ \mu_{I2} & = \alpha_0 + \alpha_{I} + \beta_{2} + \alpha\beta_{I2} + \epsilon_{I2}& = \alpha_0 + \alpha_{I} + \beta_{2} + \alpha\beta_{I2} + \epsilon_{I2}\\ \ldots & \ldots &\ldots\\ \mu_{IJ} & = \alpha_0 + \alpha_{I} + \beta_{J} + \alpha\beta_{IJ} + \epsilon_{IJ}& = \alpha_0 + \alpha_{I} + \beta_{J} + \alpha\beta_{IJ} + \epsilon_{IJ}\\ \end{array}\]

En este caso el objetivo vuelve a ser la comparación de medias pero el contraste a resolver depende de la especificación del modelo considerado. Hay que tener en cuenta que en este caso no tenemos una única variable predictora, y además hay que valorar la posibilidad de la interacción entre ambos factores.

1.4 Especificación en R

A continuación se presenta la especificación en R de los modelos anteriores y de todos los anidados que pueden surgir a partir de ellos, así como su interpretación en términos de contraste de hipótesis o coeficientes asociados a cada efecto presente en el modelo.

1.4.1 ANOVA de una vía

La especificación en R del modelo ANOVA de una vía para su estimación se realiza a través de la expresión reducida del modelo dada por: \[ Y \sim F \]

En este tipo de modelos se plantean varias situaciones de modelos anidados que deberemos estudiar:

  • Modelo con efecto de \(F\):

\[\begin{equation} Y \sim F, \end{equation}\]

Si \(F\) resulta significativo, tenemos que se detectan diferencias entre al menos dos medias para los niveles del factor \(F\).

  • Modelo sin efectos:

\[\begin{equation} Y \sim 1 \end{equation}\]

En este caso \(F\) no resulta significativo y podmeos conluir que las medias de todos los niveles del factor se pueden considerar iguales.

En la práctica el ajuste de este tipo de modelos se puede plantear como una comparación entre los modelos anteriores, dado que están anidados, de forma similar al proceso de comparación de modelos planteado en los modelos RLM y MP.

1.4.2 ANOVA de dos vías

La especificación en R del modelo ANOVA de dos vías se realiza a través de la expresión reducida del modelo dada por:

\[\begin{equation} Y \sim F_1 + F_2 + F_1 : F_2, \end{equation}\]

donde \(F_1 : F_2\) representa el efecto de interacción, y \(F_1\) y \(F_2\) son los denominados efectos principales asociados con los factores considerados. En este tipo de modelos se plantean varias situaciones de modelos anidados que deberemos estudiar para establecer cual es el tipo de contraste sobre las medias que estamos resolviendo. Alternativamente ese modelo se puede escribir de forma más reducida como \(Y \sim F_1*F_2\), que es un modelo con los efectos principales de cada factor y con la interacción.

  • Modelo saturado o modelo con todos los factores y la interacción dado por:

\[\begin{equation} Y \sim F_1 + F_2 + F_1 : F_2 \end{equation}\]

Si la interacción está presente en el modelo el contraste planteado es: \[ \left\{ \begin{array}{ll} H_0: & \mbox{No hay diferencias entre las medias conjuntas de ambos niveles del factor}\\ H_a: & \mbox{Al menos hay dos combinaciones de ambos factores con medias distintas} \end{array} \right. \]

que nos permite concluir sobre si el efecto de interacción es necesario para explicar el comportamiento de la respuesta. En caso de que dicho contraste resulte no significativo podemos plantear el modelo con ambos efectos principales.

  • Modelo con ambos efectos principales:

\[ Y \sim F_1 + F_2 \]

donde se parte de que la interacción es no significativa y se plantean dos contrastes (uno por cada factor)

\[ \left\{ \begin{array}{ll} H_0: & \mbox{No hay diferencias entre las medias de } F_1\\ H_a: & \mbox{Al menos hay dos medias de } F_1 \mbox{ distintas} \end{array} \right. \] Si no rechazamos este contraste deberemos considerar que no hay diferencias en las medias de los niveles establecidos en \(F_1\) y debemos considerar un modelo con efecto principal de \(F_2\).

\[ \left\{ \begin{array}{ll} H_0: & \mbox{No hay diferencias entre las medias de } F_2\\ H_a: & \mbox{Al menos hay dos medias de } F_2 \mbox{ distintas} \end{array} \right. \] Si no rechazamos este contraste deberemos considerar que no hay diferencias en las medias de los niveles establecidos en \(F_2\) y debemos considerar un modelo con efecto principal de \(F_1\).

Estos contrastes nos permiten concluir sobre cada factor de forma independiente (no hay interacción entre ellos).

  • Modelo con efecto principal de \(F_1\):

\[ Y \sim F_1 \]

donde se parte de que la interacción y el efecto principal de \(F_2\) son no significativos y se plantea el contraste

\[ \left\{ \begin{array}{ll} H_0: & \mbox{No hay diferencias entre las medias de } F_1\\ H_a: & \mbox{Al menos hay dos medias de } F_1 \mbox{ distintas} \end{array} \right. \]

Si no rechazamos este contraste deberemos considerar que no hay diferencias en las medias de los niveles establecidos en \(F_1\) y nuestro modelo final no contiene ningún tipo de efecto predictivo.

  • Modelo con efecto principal de \(F_2\):

\[ Y \sim F_2 \]

donde se parte de que la interacción y el efecto principal de \(F_1\) son no significativos y se plantea el contraste

\[ \left\{ \begin{array}{ll} H_0: & \mbox{No hay diferencias entre las medias de } F_2\\ H_a: & \mbox{Al menos hay dos medias de } F_2 \mbox{ distintas} \end{array} \right. \] Si no rechazamos este contraste deberemos considerar que no hay diferencias en las medias de los niveles establecidos en \(F_2\) y nuestro modelo final no contiene ningún tipo de efecto predictivo.

  • Modelo sin efectos:

\[ Y \sim 1 \]

donde concluimos que todas las medias consideradas son iguales, es decir, los niveles de los factores considerados no tienen efecto sobre la respuesta.

Esta combinación de modelos nos da una estructura secuencial de modelización (modelos alternativos) que deberemos estudiar para elegir el modelo con una mayor capacidad predictiva. Dicha estructura viene dada por:

\[\begin{array}{lll} M1: & Y \sim F_1 + F_2 + F_1:F_2\\ M2: & Y \sim F_1 + F_2 \\ M3.1: & Y \sim F_1 \\ M3.2: & Y \sim F_2 \\ M4: & Y \sim 1\\ \end{array}\]

Más adelante veremos como llevar a cabo esta comparación secuencial y como interpretarla para determinar los efectos que finalmente están presentes en el modelo y los que no.

1.5 Estimación y bondad de ajuste del modelo

Dado que los modelos ANOVA se pueden escribir como modelos lineales de regresión con una matriz de diseño que solo contiene unos y ceros, se pueden utilizar los procedimientos de estimación vistos en el capitulo anterior para estimar los parámetros del modelo. Podemos utilizar las ecuaciones normales para obtener los coeficientes del modelo que en este caso representarán los incrementos de cada nivel con respecto al de referencia (que hemos tomado como cero en la restricción de identificabilidad), que nos permiten obtener las estimaciones de las medias asociadas a cada nivel o niveles del factor o factores. Además, podemos obtener los intervalos de confianza de dichos coeficientes.

Sin embargo, en este caso no se pueden interpretar de forma directa 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 con el factor o interacción de factores. Por eso es necesario proporcionar una medida de bondad de ajuste con dichas estimaciones que en este caso es el test F ANOVA que se obtiene como el test F de la regresión de los modelos de unidades anteriores. Este test nos permite extraer conclusiones directas sobre los efectos presentes en el modelo y no sobre los coeficientes estimados. En el punto siguiente veremos como utilizar el test F ANOVA para obtener el mejor modelo. En este caso no utilizaremos el \(R^2\) como medida de bondad de ajusta ya que se encuentra diseñada para modelos de regresión donde las predictoras son de tipo numérico.

Para resolver la estimación y la bondad de ajuste en este tipo de modelos utilizaremos las funciones utilizadas en la unidad anterior. En algunos casos modificaremos la sintaxis para adaptarnos a las variables predictoras de tipo factor.

1.5.1 Ejemplos

Pasamos a estudiar cada uno de los ejemplos ajustando el modelo saturado, es decir, aquel que contiene todos los efectos posibles en el modelo. Por defecto el programa R utiliza como restricción de identificabilidad que el incremento asociado con la primera categoría (por orden numérico o texto de los niveles del factor) es igual a cero, es decir, que la estimación de dicho incremento corresponde con la interceptación o efecto común del modelo.

En las tablas de estimaciones eliminamos el p-valor individual asociado a cada incremento, dado que en este caso esos efectos no deben interpretarse de forma individual.

1.5.1.1 Datos de Insecticidas

Ajustamos un modelo ANOVA de una vía para el conjunto de datos de insecticida. En este caso la restricción de identificabilidad es fijar el incremento del spray A igual a cero, es decir, \(\alpha_{spray,A} = 0\).

  count
Predictors Estimates CI
(Intercept) 14.50 12.24 – 16.76
spray [B] 0.83 -2.36 – 4.03
spray [C] -12.42 -15.61 – -9.22
spray [D] -9.58 -12.78 – -6.39
spray [E] -11.00 -14.20 – -7.80
spray [F] 2.17 -1.03 – 5.36
Observations 72

Las ecuaciones de estimación de las medias del número de insectos vivos se pueden extraer a partir de la tabla de coeficientes estimados. Sin embargo, podemos obtener las estimaciones y representar las medias de todos los sprays haciendo uso de la función update() quitando el término de interceptación del modelo. De esta forma se prescinde de la restricción de identificabilidad y se obtienen las medias directamente.

  count
Predictors Estimates CI
spray [A] 14.50 12.24 – 16.76
spray [B] 15.33 13.07 – 17.59
spray [C] 2.08 -0.18 – 4.34
spray [D] 4.92 2.66 – 7.18
spray [E] 3.50 1.24 – 5.76
spray [F] 16.67 14.41 – 18.93
Observations 72

Podemos ver que los sprays A, B, y F son los menos efectivos (medias de insectos vivos más altas), mientras que los otros tres tienen medias significativamente más bajas. El spray C es el que se muestra más efectivo con una media de insectos vivos de 2 (redondeando a entero).

Con respecto a la bondad del ajuste podemos ver que el test \(F\) resulta significativo indicando que al menos alguna de las medias es distinta del resto (algún incremento distinto de cero). Por tanto, existe un efecto del tipo de spray en el número de insectos que sobreviven.

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.724         0.704  3.92      34.7 3.18e-17     6  -197.  409.  425.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

1.5.1.2 Datos de Envasado

Ajustamos un modelo ANOVA de una vía para el conjunto de datos de envasado. En este caso la restricción de identificabilidad es fijar el incremento de Máquina M1 igual a cero, es decir, \(\alpha_{maquina,M1} = 0\).

  produccion
Predictors Estimates CI
(Intercept) 106.00 98.00 – 114.00
maquina [M2] 7.75 -3.57 – 19.07
maquina [M3] -1.50 -12.82 – 9.82
maquina [M4] 18.00 6.68 – 29.32
Observations 16

A la vista de los coeficientes la máquina el orden de las máquinas (de mejor a peor, es decir, de mayor número de envases sin defecto a menor) es \(M4 > M2 > M1 > M3\), con media estimadas dadas por:

  produccion
Predictors Estimates CI
maquina [M1] 106.00 98.00 – 114.00
maquina [M2] 113.75 105.75 – 121.75
maquina [M3] 104.50 96.50 – 112.50
maquina [M4] 124.00 116.00 – 132.00
Observations 16

A la vista del gráfico de estimación no se aprecian comportamientos muy diferenciados entre las diferentes máquinas.

Con respecto a la bondad del ajuste podemos ver que el test \(F\) resulta significativo indicando que al menos alguna de las medias es distinta del resto (algún incremento distinto de cero). Por tanto, existe un efecto del tipo de máquina en el número de envases sin defecto.

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.596         0.496  7.35      5.91  0.0102     4  -52.3  115.  118.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

1.5.1.3 Datos de Venenos

Ajustamos un modelo ANOVA de dos vías con interacción entre veneno y antidoto para el conjunto de datos de venenos. En este caso las restricciones de identificabilidad son:

  • \(\alpha_{antidoto,AA} = 0\),
  • \(\beta_{veneno,VA} = 0\),
  • \(\alpha\beta_{Antidoto,AA; Veneno,VA} = 0\),

de forma que el ajuste para este modelo es

  tiempo
Predictors Estimates CI
(Intercept) 0.41 0.26 – 0.56
antidoto [AB] 0.47 0.25 – 0.68
antidoto [AC] 0.15 -0.06 – 0.37
antidoto [AD] 0.20 -0.02 – 0.41
veneno [VB] -0.09 -0.31 – 0.12
veneno [VC] -0.20 -0.42 – 0.01
antidoto [AB] * veneno
[VB] 		0.03 -0.27 – 0.33
antidoto [AC] * veneno
[VB] 		-0.10 -0.40 – 0.20
antidoto [AD] * veneno
[VB] 		0.15 -0.15 – 0.45
antidoto [AB] * veneno
[VC] 		-0.34 -0.64 – -0.04
antidoto [AC] * veneno
[VC] 		-0.13 -0.43 – 0.17
antidoto [AD] * veneno
[VC] 		-0.08 -0.39 – 0.22
Observations 48

Desde luego interpretar la tabla de coeficientes en este caso es más complicado, dada la gran cantidad de parámetros involucrados, por lo que resulta más fácil obtener las medias y el gráfico asociado. En este caso no podemos actualizar el modelo eliminando la interceptación dado que tenemos tres restricciones y no una. Modificamos la tabla de estimaciones y el gráfico para obtener todo el conjunto de medias. Además, como no hemos descartado el efecto de interacción, en el gráfico sólo representaremos las medias asociadas con dicho efecto. Para representar el efecto de interacción añadimos el parámetro "int" a la función plot_model.

## 
## # Predicted values of tiempo
## # x = antidoto
## 
## # veneno = VA
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      0.41 | 0.07 |        VA | [0.27, 0.56]
## 2 |      0.88 | 0.07 |        VA | [0.73, 1.03]
## 3 |      0.57 | 0.07 |        VA | [0.42, 0.71]
## 4 |      0.61 | 0.07 |        VA | [0.46, 0.76]
## 
## # veneno = VB
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      0.32 | 0.07 |        VB | [0.17, 0.47]
## 2 |      0.82 | 0.07 |        VB | [0.67, 0.96]
## 3 |      0.38 | 0.07 |        VB | [0.23, 0.52]
## 4 |      0.67 | 0.07 |        VB | [0.52, 0.81]
## 
## # veneno = VC
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      0.21 | 0.07 |        VC | [0.06, 0.36]
## 2 |      0.33 | 0.07 |        VC | [0.19, 0.48]
## 3 |      0.24 | 0.07 |        VC | [0.09, 0.38]
## 4 |      0.32 | 0.07 |        VC | [0.18, 0.47]

Con respecto a la bondad del ajuste podemos ver que el test \(F\) resulta significativo indicando que al menos alguna de las medias es distinta del resto (algún incremento distinto de cero). Por tanto, existe un efecto de interacción entre venenos y antídoto que proporciona un tiempo de reacción distinto.

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.734         0.653 0.149      9.03 1.93e-7    12   30.2 -34.3 -9.98
## # … with 2 more variables: deviance <dbl>, df.residual <int>

1.6 Comparación y selecció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. Dado que el número de modelos que podemos construir en estas situaciones es limitado, podríamos construirlos todos ellos y compararlos todos de golpe para decidir el mejor modelo. Esto se puede hacer con el test F de comparación de modelos que es una versión del test F para un contraste especifico o bien utilizar el criterio del AIC o BIC y quedarnos con aquel modelo que de un valor más pequeño. En modelos más complejos optaremos por un procedimiento automático por pasos empezando por el modelo saturado. A continuación, presentamos las diferentes opciones con los bancos de datos que venimos trabajando en esta unidad.

Se recomienda leer el apartado de Selección de variables de la unidad anterior para completar los aspectos teóricos de este apartado.

1.6.1 Ejemplos

Para el desarrollo de los ejemplos utilizamos las funciones presentadas en la unidad anterior.

1.6.1.1 Datos de insecticidas

Dado que en el modelo asociado a este banco de datos es un ANOVA de una vía sólo existen dos posibles modelos:

  • \(count \sim spray\): Existe al memnos un spray que se proporciona un nivel medio de supervivientes diferente.
  • \(count \sim 1\): Todos los sprays se comportan de la misma forma.

Para comparar ambos modelos utilizamos el test \(F\) parcial.

## Analysis of Variance Table
## 
## Model 1: count ~ 1
## Model 2: count ~ spray
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     71 3684.0                                  
## 2     66 1015.2  5    2668.8 34.702 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test \(F\) parcial resulta significativo (p-valor < 0.05) indicando que el modelo con el factor spray es preferible al modelo sin spray. Dado que el factor considerado es significativo resulta necesario conocer que medias son realmente distintas y cuales podemos considerar iguales. Para realizar esta tarea utilizaremos los denominados tests de comparaciones múltiples que nos permiten comparar dos a dos todas las combinaciones del factor para determinar las que son significativamente distintas o no.

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = x)
## 
## $spray
##            diff        lwr       upr     p adj
## B-A   0.8333333  -3.866075  5.532742 0.9951810
## C-A -12.4166667 -17.116075 -7.717258 0.0000000
## D-A  -9.5833333 -14.282742 -4.883925 0.0000014
## E-A -11.0000000 -15.699409 -6.300591 0.0000000
## F-A   2.1666667  -2.532742  6.866075 0.7542147
## C-B -13.2500000 -17.949409 -8.550591 0.0000000
## D-B -10.4166667 -15.116075 -5.717258 0.0000002
## E-B -11.8333333 -16.532742 -7.133925 0.0000000
## F-B   1.3333333  -3.366075  6.032742 0.9603075
## D-C   2.8333333  -1.866075  7.532742 0.4920707
## E-C   1.4166667  -3.282742  6.116075 0.9488669
## F-C  14.5833333   9.883925 19.282742 0.0000000
## E-D  -1.4166667  -6.116075  3.282742 0.9488669
## F-D  11.7500000   7.050591 16.449409 0.0000000
## F-E  13.1666667   8.467258 17.866075 0.0000000

Para interpretar la tabla resultante nos debemos fijar en la columna p adj. Todos aquellos inferiores a 0.05 indican que las medias son distintas. Analizando todos los p-valores podemos establecer dos grupos de sprays, el formado por los tipos A, B y F y el formado por los tipos C, D, y E. Dentro de cada grupo no hay diferencias pero sin entres los sprays de grupos distintos.

También podemos representar gráficamente las comparaciones realizadas:

Los intervalos de confianza que contienen al cero muestran aquellos pares de sprays que pueden considerarse iguales. Los intervalos de confianza a la izquierda de cero muestran las combinaciones de sprays donde el primer spray tiene un efecto menor que el segundo (media negativa), mientras que en la parte derecha tenemos las combinaciones donde el primer spray es superior al segundo. Por ejemplo, podemos concluir:

  • los sprays A y B tienen comportamiento similar,
  • los sprays A y C tiene comportamientos distintos con la media del número de supervivientes en A superior que en con el spray C.
  • los sprays F y E tiene comportamientos distintos con la media del número de supervivientes en F superior que en con el spray E.

1.6.1.2 Datos de envasado

Consideramos los modelos:

  • \(produccion \sim maquina\): Existe al menos un spray que se proporciona un nivel medio de supervivientes diferente.
  • \(produccion \sim 1\): Todos los sprays se comportan de la misma forma.

Para comparar ambos modelos utilizamos el test \(F\) parcial.

## Analysis of Variance Table
## 
## Model 1: produccion ~ 1
## Model 2: produccion ~ maquina
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
## 1     15 1604.94                              
## 2     12  647.75  3    957.19 5.9108 0.01024 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test resulta significativo indicando que ambos son distintos. Esto implica directamente que hay al menos una máquina con un comportamiento diferente del resto. Realizamos las comparaciones múltiples para comparar todas las máquinas dos a dos.

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = x)
## 
## $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

Tan sólo las combinaciones M4-M1 y M4-M3 muestran resultados significativos indicando que la máquina M4 es comparable con la M2, es decir, la media de envases sin defecto es igual para dichas máquinas, pero distintas de la M1 y M3 que también puede considerarse iguales. Veamos los resultados gráficamente:

1.6.1.3 Datos de venenos

En este caso la especificación de todos los posibles modelos así como su comparación resulta más costosa ya que están involucrados dos factores de clasificación. Los posibles modelos son:

  • \(tiempo \sim veneno * antidoto\): Modelo con interacción, donde hay al menos una combinación de veneno y antídoto con tiempo de reacción medio distinto de otras combinaciones.
  • \(tiempo \sim veneno + antidoto\): Modelo sin interacción con efectos principales, donde el efecto del veneno y el antídoto tiene un efecto aditivo en e tiempo medio de reacción.
  • \(tiempo \sim veneno\): Modelo con efecto principal marginal del veneno, donde hay almenos un tipo de veneno con un tiempo medio de reacción distinto pero los antidotos no tienen efecto.
  • \(tiempo \sim antidoto\): Modelo con efecto principal marginal del antídoto, donde hay al menos un tipo de antídoto con un tiempo medio de reacción distinto pero los venenos no tienen efecto.
  • \(tiempo \sim 1\): Modelo sin efectos, donde el tiempo de reacción no cambia con el antídoto y el tipo de veneno.

Vamos a realizar el contraste \(F\) parcial para determinar si la interacción es necesaria en este modelo.

## Analysis of Variance Table
## 
## Model 1: tiempo ~ veneno + antidoto
## Model 2: tiempo ~ antidoto * veneno
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     42 1.0504                           
## 2     36 0.8001  6   0.25027 1.8768 0.1118

Dado que el p-valor no es significativo (p-valor > 0.05) podemos concluir que ambos modelos pueden ser considerarse como iguales, es decir, el efecto de interacción no resulta significativo. Vamos a comparar las estimaciones de ambos modelos y representamos gráficamente el modelo sin interacción. En este caso obtendremos un gráfico por cada factor dado que no actúan de forma conjunta sino de forma aditiva. Utilizamos el parámetro "pred" para obtener las estimaciones de las medias asociadas a cada factor.

  tiempo tiempo
Predictors Estimates CI p Estimates CI p
(Intercept) 0.41 0.26 – 0.56 <0.001 0.45 0.34 – 0.57 <0.001
antidoto [AB] 0.47 0.25 – 0.68 <0.001 0.36 0.23 – 0.49 <0.001
antidoto [AC] 0.15 -0.06 – 0.37 0.150 0.08 -0.05 – 0.21 0.232
antidoto [AD] 0.20 -0.02 – 0.41 0.069 0.22 0.09 – 0.35 0.002
veneno [VB] -0.09 -0.31 – 0.12 0.386 -0.07 -0.19 – 0.04 0.198
veneno [VC] -0.20 -0.42 – 0.01 0.063 -0.34 -0.45 – -0.23 <0.001
antidoto [AB] * veneno
[VB] 		0.03 -0.27 – 0.33 0.855
antidoto [AC] * veneno
[VB] 		-0.10 -0.40 – 0.20 0.507
antidoto [AD] * veneno
[VB] 		0.15 -0.15 – 0.45 0.321
antidoto [AB] * veneno
[VC] 		-0.34 -0.64 – -0.04 0.028
antidoto [AC] * veneno
[VC] 		-0.13 -0.43 – 0.17 0.389
antidoto [AD] * veneno
[VC] 		-0.08 -0.39 – 0.22 0.572
Observations 48 48 
## $veneno

## 
## $antidoto

Como era de esperar, dado que la interacción no es significativa, los coeficientes estimados (incrementos respecto de las medias de referencia dadas por el antídoto AA y el veneno VA) no varían sustancialmente entre ambos modelos. En el gráfico podemos ver que el veneno con un tiempo de reacción más pequeño es el asociado con el veneno VC, mientras que si nos fijamos en el antídoto, el menor tiempo de reacción se da para el AA. Dado que no hay interacción el modelo sugiere que independientemente del veneno utilizado, el antídoto más efectivo es el AA, seguido por el AC, AD, y por último el AB. Para confirmar este hecho deberíamos obtener las ecuaciones de estimación asociadas que vienen dadas por:

\[\begin{array}{lll} \mu_{AA,VA} & = \alpha_0 + \alpha_{AA} + \beta_{VA} & = 0.45 + 0 + 0 = 0.45\\ \mu_{AA,VB} & = \alpha_0 + \alpha_{AA} + \beta_{VB} & = 0.45 + 0 - 0.07 = 0.38\\ \mu_{AA,VC} & = \alpha_0 + \alpha_{AA} + \beta_{VC} & = 0.45 + 0 - 0.34 = 0.11\\ \mu_{AB,VA} & = \alpha_0 + \alpha_{AB} + \beta_{VA} & = 0.45 + 0.36 + 0 = 0.81\\ \mu_{AB,VB} & = \alpha_0 + \alpha_{AB} + \beta_{VB} & = 0.45 + 0.36 - 0.07 = 0.74\\ \mu_{AB,VC} & = \alpha_0 + \alpha_{AB} + \beta_{VC} & = 0.45 + 0.36 - 0.34 = 0.47\\ \mu_{AC,VA} & = \alpha_0 + \alpha_{AC} + \beta_{VA} & = 0.45 + 0.08 + 0 = 0.53\\ \mu_{AC,VB} & = \alpha_0 + \alpha_{AC} + \beta_{VB} & = 0.45 + 0.08 - 0.07 = 0.46\\ \mu_{AC,VC} & = \alpha_0 + \alpha_{AC} + \beta_{VC} & = 0.45 + 0.08 - 0.34 = 0.19\\ \mu_{AD,VA} & = \alpha_0 + \alpha_{AD} + \beta_{VA} & = 0.45 + 0.22 + 0 = 0.67\\ \mu_{AD,VB} & = \alpha_0 + \alpha_{AD} + \beta_{VB} & = 0.45 + 0.22 - 0.07 = 0.60\\ \mu_{AD,VC} & = \alpha_0 + \alpha_{AD} + \beta_{VC} & = 0.45 + 0.22 - 0.34 = 0.33\\ \end{array}\]

Para determinar el modelo final sin necesidad de comparar todos los posibles modelos utilizamos los procedimientos secuenciales. Fijamos el valor de la significatividad del test \(F\) para descartar efectos no significativos mediante el parámetro prem en la función ols_step_backward_p.

## 
## 
##                               Elimination Summary                               
## -------------------------------------------------------------------------------
##         Variable                         Adj.                                      
## Step        Removed        R-Square    R-Square     C(p)       AIC        RMSE     
## -------------------------------------------------------------------------------
##    1    antidoto:veneno      0.6508      0.6092    5.2608    -33.2407    0.1581    
## -------------------------------------------------------------------------------

Para utilizar el criterio AIC utilizamos la función step() porque la función ols_step_backward_aic no esta diseñada especificamente para trabajar con interacciones de factores y puede lleavr a soluciones erróneas. La única opción pasa por evaluar todos los modelos con la función ols_step_all_possible y elegir el mejor de ellos teniendo en cuenta la construcción del modelo. Mostraremos como hacerlo para este ejemplo tan sencillo pero no será práctico en modelos más complejos.

## Start:  AIC=-172.52
## tiempo ~ 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, data = venenos)
## 
## 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 indica que no podemos descartar el efecto de interacción contradiciendo los resultados obtenidos con el test \(F\). En este tipo de modelos siempre es preferible utilizar el test \(F\) en lugar del \(AIC\). Veamos la solución con todos los modelos:

##   mindex n                      predictors        aic
## 3      1 1                 antidoto:veneno -32.304440
## 2      2 1                          veneno  -9.042063
## 1      3 1                        antidoto  -4.274275
## 5      4 2        antidoto antidoto:veneno -34.304440
## 6      5 2          veneno antidoto:veneno -34.304440
## 4      6 2                 antidoto veneno -33.240672
## 7      7 3 antidoto veneno antidoto:veneno -34.304440

Lo primero que debemos hacer es ignorar aquellos modelos que no resultan posibles desde el punto de vista de la modelización estadística. Estos modelos son todos aquellos que consideran la interacción pero no los dos factores principales asociados. Estos son los modelos (mindex) 1, 4, y 5. Podemos ver en ese caso que los modelos con menor AIC son el 7 (modelo con interacción) y el 6 (modelo sin interacción). Concretamente el seleccionado sería el de interacción contradiciendo en cierta forma los resultados obtenidos con el test \(F\) como hemos visto anteriormente.

1.7 Diagnóstico

Los procedimientos de diagnóstico son los mismos que en los modelos de regresión salvo por el hecho de que la igualdad de varianzas de los residuos se debe cumplir para 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 los tests vistos en la unidad anterior, mientras que para la hipótesis de igualdad de varianzas tenemos diferentes alternativas:

  • bartlett.test(): Test de Bartlett para homogeneidad de varianzas utilizando el test \(F\). Este test puedo funcionar incorrectamente cuando los residuos no son normales.
  • leveneTest(): Test de Levene para igualdad de varianzas. Funciona bien en cualquier situación.

También se pueden realizar los diagnósticos gráficos de residuos para cada hipótesis, pero en para este tipo de modelos nos centraremos principalmente en la realización de los tests. El único gráfico que realizaremos será el gráfico de cajas de residuos versus niveles del factor para validar la homogeneidad de varianzas y la linealidad de las medias estimadas. La amplitud de las cajas debería ser similar para verificar homogeneidad y las medias deberían situarse en torno al cero para validar linealidad.

En cuanto al estudio de influencia, nos centraremos principalmente en la distancia de Cook dado que no tiene sentido considerar la influencia sobre los coeficientes del modelo ni sobre los valores ajustados dada la estructura agrupada del banco de datos a analizar.

1.7.1 Ejemplos

Para realizar el diagnóstico partiremos del modelo ajustado en el apartado anterior para cada banco de datos. En caso de detectar algún problema con las hipótesis del modelo deberemos ajustar un nuevo modelo y verificar que cumple las hipótesis.

1.7.1.1 Datos de insecticidas

Obtenemos los valores de diagnóstico y realizamos los correspondientes tests de hipótesis y análisis de influencia. En este caso utilizamos la versión del test de barlett más precisa (función barlett.test), sin la aproximación que realiza la función ols_test_bartlett.

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9601         0.0223 
## Kolmogorov-Smirnov        0.1301         0.1747 
## Cramer-von Mises          6.6076         0.0000 
## Anderson-Darling          1.2015         0.0037 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  5  3.8214 0.004223 **
##       66                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Las medias presentan pocas desviaciones respecto del valor cero y la amplitud de las cajas es bastante similar. Haciendo uso de los test se verifican las hipótesis de normalidad e igualdad de varianzas. Se detectan observaciones potencialmente influyentes, pero ningún valor de la distancia de Cook es superior a 1 y no clasificamos ninguna observación como influyente.

El modelo final para este conjunto de datos es el analizado en la sección de estimación de esta unidad.

1.7.1.2 Datos de envasado

Obtenemos los valores de diagnóstico y realizamos los correspondientes tests de hipótesis y análisis de influencia.

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9074         0.1058 
## Kolmogorov-Smirnov        0.176          0.6420 
## Cramer-von Mises          1.5833         0.0000 
## Anderson-Darling          0.5552         0.1269 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.2091 0.8882
##       12

Los tests de hipótesis resultan no significativos indicando que el modelo propuesto verifica dichas hipótesis. No se detectan además observaciones influyentes, con lo que el modelo propuesto y que fue analizado en la sección de estimación es el modelo final para este banco de datos.

1.7.1.3 Datos de venenos

Utilizamos el modelo sin interacción que hemos obtenido tras el proceso de selección de efectos tratada en la sección anterior.

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9221         0.0035 
## Kolmogorov-Smirnov        0.113          0.5353 
## Cramer-von Mises         12.0813         0.0000 
## Anderson-Darling          0.8952         0.0205 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value    Pr(>F)    
## group 11  4.1582 0.0005539 ***
##       36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se verifica la hipótesis de normalidad pero no así la de igualdad de varianzas. No se presentan oes influyentes. Se plantea utilizar Box-Cox para obtener una posible transformación de la respuesta que permita verificar las hipótesis del modelo.

A la vista del intervalo de confianza obtenido para el parámetro de Box-Cox, la transformación que deberíamos usar es la inversa. Transformamos la variable respuesta y procedemos de nuevo con la selección del modelo.

## 
## 
##                              Elimination Summary                               
## ------------------------------------------------------------------------------
##         Variable                         Adj.                                     
## Step        Removed        R-Square    R-Square     C(p)       AIC       RMSE     
## ------------------------------------------------------------------------------
##    1    antidoto:veneno      0.8454       0.827    0.4181    75.5483    0.4911    
## ------------------------------------------------------------------------------

El modelo resultante es aquel que no contiene la interacción. Obtenemos los coeficientes de ese modelo, así como las medias estimadas para la combinación de veneno-antídoto. Expresamos las medias con la transformación propuesta.

  tiempoinv
Predictors Estimates CI p
(Intercept) 2.70 2.35 – 3.05 <0.001
antidoto [AB] -1.66 -2.06 – -1.25 <0.001
antidoto [AC] -0.57 -0.98 – -0.17 0.007
antidoto [AD] -1.35 -1.76 – -0.95 <0.001
veneno [VB] 0.47 0.12 – 0.82 0.010
veneno [VC] 2.00 1.65 – 2.35 <0.001
Observations 48

Las ecuaciones de estimación son:

\[\begin{array}{lll} 1/\mu_{AA,VA} & = \alpha_0 + \alpha_{AA} + \beta_{VA} & = 2.70 + 0 + 0 = 2.70\\ 1/\mu_{AA,VB} & = \alpha_0 + \alpha_{AA} + \beta_{VB} & = 2.70 + 0 + 0.47 = 3.17\\ 1/\mu_{AA,VC} & = \alpha_0 + \alpha_{AA} + \beta_{VC} & = 2.70 + 0 + 2.00 = 4.70\\ 1/\mu_{AB,VA} & = \alpha_0 + \alpha_{AB} + \beta_{VA} & = 2.70 - 1.66 + 0 = 1.04\\ 1/\mu_{AB,VB} & = \alpha_0 + \alpha_{AB} + \beta_{VB} & = 2.70 - 1.66 + 0.47 = 1.51\\ 1/\mu_{AB,VC} & = \alpha_0 + \alpha_{AB} + \beta_{VC} & = 2.70 - 1.66 + 2.00 = 3.04\\ 1/\mu_{AC,VA} & = \alpha_0 + \alpha_{AC} + \beta_{VA} & = 2.70 - 0.57 + 0 = 2.13\\ 1/\mu_{AC,VB} & = \alpha_0 + \alpha_{AC} + \beta_{VB} & = 2.70 - 0.57 + 0.47 = 2.60\\ 1/\mu_{AC,VC} & = \alpha_0 + \alpha_{AC} + \beta_{VC} & = 2.70 - 0.57 + 2.00 = 4.13\\ 1/\mu_{AD,VA} & = \alpha_0 + \alpha_{AD} + \beta_{VA} & = 2.70 - 1.35 + 0 = 1.35\\ 1/\mu_{AD,VB} & = \alpha_0 + \alpha_{AD} + \beta_{VB} & = 2.70 - 1.35 + 0.47 = 1.82\\ 1/\mu_{AD,VC} & = \alpha_0 + \alpha_{AD} + \beta_{VC} & = 2.70 - 1.35 + 2.00 = 3.35\\ \end{array}\]

Veamos si se verifican las hipótesis del modelo:

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9789         0.5339 
## Kolmogorov-Smirnov        0.1057         0.6192 
## Cramer-von Mises          6.0617         0.0000 
## Anderson-Darling          0.2948         0.5834 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group 11  1.1672 0.3429
##       36

Se verifican las hipótesis del modelo. La transformación planteada permite obtener un modelo válido para la fase de predicción, y nos permite extraer conclusiones sobre la relación entre veneno y antídoto con el tiempo de reacción.

1.8 Predicción

El proceso de predicción en este tipo de modelos es muy simple ya que los posibles valores de las predictoras coinciden con las medias de los posibles niveles o combinación de los niveles del factor o factores que aparezcan en el modelo final. En este caso solo tiene sentido predecir los valores de la media ya que nuestro objetivo es el estudio de dichas medias. Dichas medias se obtienen a partir de los incrementos estimados para cada modelo.

1.8.1 Ejemplos

Obtenemos las tablas de medias estimadas para cada uno de los ejemplos tratados en esta unidad.

1.8.1.1 Datos de insecticidas

Podemos obtener la tabla de predicción y el gráfico asociado con el código siguiente:

## 
## # Predicted values of count
## # x = spray
## 
## x | Predicted |   SE | group_col |         95% CI
## -------------------------------------------------
## 1 |     14.50 | 1.13 |         1 | [12.28, 16.72]
## 2 |     15.33 | 1.13 |         1 | [13.11, 17.55]
## 3 |      2.08 | 1.13 |         1 | [-0.14,  4.30]
## 4 |      4.92 | 1.13 |         1 | [ 2.70,  7.14]
## 5 |      3.50 | 1.13 |         1 | [ 1.28,  5.72]
## 6 |     16.67 | 1.13 |         1 | [14.45, 18.89]

1.8.1.2 Datos de envasado

Podemos obtener la tabla de predicción y el gráfico asociado con el código siguiente:

## 
## # Predicted values of produccion
## # x = maquina
## 
## x | Predicted |   SE | group_col |           95% CI
## ---------------------------------------------------
## 1 |    106.00 | 3.67 |         1 | [ 98.80, 113.20]
## 2 |    113.75 | 3.67 |         1 | [106.55, 120.95]
## 3 |    104.50 | 3.67 |         1 | [ 97.30, 111.70]
## 4 |    124.00 | 3.67 |         1 | [116.80, 131.20]

1.8.1.3 Datos de venenos

En este caso debemos tener en cuenta que los valores de predicción vendrán expresados en términos de la inversa de la respuesta. Habrá que deshacer la transformación para obtener los valores de predicción en la escala real de la predictora.

## 
## # Predicted values of tiempoinv
## # x = veneno
## 
## # antidoto = AA
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      2.70 | 0.17 |        AA | [2.36, 3.04]
## 2 |      3.16 | 0.17 |        AA | [2.82, 3.50]
## 3 |      4.70 | 0.17 |        AA | [4.36, 5.04]
## 
## # antidoto = AB
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      1.04 | 0.17 |        AB | [0.70, 1.38]
## 2 |      1.51 | 0.17 |        AB | [1.17, 1.85]
## 3 |      3.04 | 0.17 |        AB | [2.70, 3.38]
## 
## # antidoto = AC
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      2.12 | 0.17 |        AC | [1.78, 2.46]
## 2 |      2.59 | 0.17 |        AC | [2.25, 2.93]
## 3 |      4.13 | 0.17 |        AC | [3.78, 4.47]
## 
## # antidoto = AD
## 
## x | Predicted |   SE | group_col |       95% CI
## -----------------------------------------------
## 1 |      1.34 | 0.17 |        AD | [1.00, 1.68]
## 2 |      1.81 | 0.17 |        AD | [1.47, 2.15]
## 3 |      3.35 | 0.17 |        AD | [3.01, 3.69]

Tomando los valores de las tablas anteriores y calcula la inversa tendríamos la estimación en la escala de la variable original (tiempo).

2 Modelos de análisis de la covarianza

En este punto se presentan los modelos ANCOVA o modelos de análisis de la covarianza. Esto modelos surgen cuando entre las posibles variables predictoras de la respuesta (de tipo numérico) consideramos tanto variables numéricas como factores. El objetivo principal de este tipo de modelos es estudiar si la relación entre las predictoras numéricas y la respuesta viene condicionada por el factor o factores de clasificación considerados, es decir, si debemos construir:

  • un único modelo entre la respuesta y las predictoras de tipo numérico,
  • un único modelo entre la respuesta y las predictoras de tipo categórico (factores),
  • un modelo diferente entre la respuesta y las predictoras numéricas para cada nivel o combinaciones de niveles de los factores.

Este tipo de modelos permiten una versatilidad que nos posibilita el estudio de situaciones experimentales más complejas. Sin embargo, no están exentos de dificultades sobre todo en lo que tiene que ver con el cumplimiento de las hipótesis del modelo. En función del modelo final las hipótesis de normalidad y homogeneidad varían en su aplicación. Además, el proceso de selección del mejor modelo requiere de un proceso de análisis más profundo debido a la inclusión de diferentes tipos de variables en el conjunto de posibles predictoras.

Para introducir los conceptos básicos de este tipo de modelos y mostrar todas sus posibilidades de análisis comenzaremos con el modelo ANCOVA más sencillo donde únicamente consideramos dos variables predictoras, una de tipo numérico y la otra un factor. La formulación presentada se puede generalizar rápidamente a situaciones más complejas donde el número de predictoras sea mayor.

2.1 Bancos de datos

Veamos los diferentes bancos de datos que iremos analizando a lo largo de la unidad.

Ejemplo 1. Datos de tiempo de vida. Se desea estudiar el tiempo de vida de una pieza (vida) cortadora de dos tipos, A y B (herramienta), en función de la velocidad del torno (velocidad) en el que está integrada (en revoluciones por segundo). El objetivo del análisis es describir la relación entre el tiempo de vida de la pieza y la velocidad del torno, teniendo en cuenta de qué tipo es la pieza.

Cargamos los datos y realizamos el gráfico descriptivo:

A la vista de la Figura anterior se puede apreciar que el comportamiento del tiempo de vida con respecto a la velocidad disminuye al aumentar esta última, y además ese descenso es distinto en función de la herramienta utilizada. A la misma velocidad la herramienta A proporciona un tiempo de vida inferior que la herramienta B. Para ambas máquinas se aprecia una tendencia lineal bastante clara que deberá ser objeto de estudio. Esto implicaría que nuestro modelo se desdoblaría en dos rectas de regresión lineal simple en función de la herramienta utilizada.

Ejemplo 2. Datos de longevidad. Partridge y Farquhar realizan un experimento para relacionar la vida útil (longevidad) de las moscas de la fruta con su actividad sexual (actividad). La información recogida es la longevidad en días de 125 moscas macho, divididas en cinco grupos bajo diferentes condiciones ambientales para medir su actividad sexual. Asimismo, se recoge la longitud del tórax (thorax) ya que se sospecha que afecta directamente a la longevidad de las moscas.

Cargamos los datos y realizamos el gráfico descriptivo:

# Carga de datos
thorax <- c(0.68, 0.68, 0.72, 0.72, 0.76, 0.76, 
            0.76, 0.76, 0.76, 0.8, 0.8, 0.8, 0.84, 0.84, 
            0.84, 0.84, 0.84, 0.84, 0.88, 0.88, 0.92, 0.92, 
            0.92, 0.94, 0.64, 0.7, 0.72, 0.72, 0.72, 0.76, 
            0.78, 0.8, 0.84, 0.84, 0.84, 0.84, 0.84, 0.88, 
            0.88, 0.88, 0.88, 0.88, 0.92, 0.92, 0.92, 0.92, 
            0.92, 0.92, 0.94, 0.64, 0.68, 0.72, 0.76, 0.76, 
            0.8, 0.8, 0.8, 0.82, 0.82, 0.84, 0.84, 0.84, 0.84, 
            0.84, 0.84, 0.88, 0.88, 0.88, 0.88, 0.88, 0.88, 
            0.88, 0.92, 0.92, 0.68, 0.68, 0.72, 0.76, 0.78, 
            0.8, 0.8, 0.8, 0.84,  0.84, 0.84, 0.84, 0.84, 
            0.84, 0.88, 0.88, 0.88, 0.9, 0.9, 0.9, 0.9, 0.9, 
            0.9, 0.92, 0.92, 0.64, 0.64, 0.68, 0.72, 0.72, 
            0.74, 0.76, 0.76, 0.76, 0.78, 0.8, 0.8, 0.82, 
            0.82, 0.84, 0.84, 0.84, 0.84, 0.88, 0.88, 0.88, 
            0.88, 0.88, 0.88, 0.92)
longevidad <- c(37, 49,  46,  63,  39,  46,  56,  63,  65,  
                56,  65,  70,  63,  65, 70,  77,  81,  86,  
                70,  70,  77,  77,  81,  77,  40,  37,  44, 
                47,  47,  47,  68,  47,  54,  61,  71,  75,  
                89,  58,  59,  62, 79,  96,  58,  62,  70,  
                72,  75,  96,  75,  46,  42,  65,  46, 58,  
                42,  48,  58,  50,  80,  63,  65,  70,  70,  
                72,  97,  46, 56,  70,  70,  72,  76,  90,  
                76,  92,  21,  40,  44,  54,  36, 40,  56,  
                60,  48,  53,  60,  60,  65,  68,  60,  81,  
                81,  48, 48,  56,  68,  75,  81,  48,  68,  
                16,  19,  19,  32,  33,  33, 30,  42,  42,  
                33,  26,  30,  40,  54,  34,  34,  47,  47, 
                42, 47,  54,  54,  56,  60,  44)
actividad <- c(rep("G1", 24), rep("G2", 25), rep("G3", 25), rep("G4", 25), rep("G5", 25))
longevidad <- data.frame(thorax, longevidad, actividad)
# Gráfico
ggplot(longevidad, aes(x = thorax, y = longevidad, color = actividad)) + 
  geom_point()

Se puede ver como la longevidad aumenta cuando aumenta la longitud del thorax pero ese crecimiento no parece distinto según actividad, dado que las nubes de puntos están bastante mezcladas. En este caso no parece adecuado un modelo lineal para cada grupo de actividad.

2.2 Modelo ANCOVA

Consideramos el modelo ANCOVA más sencillo donde consideramos dos variables predictoras: una numérica y otra un factor. Consideramos una muestra de tamaño \(n\) donde tenemos:

  • Una variable respuesta, \(Y\), de tipo numérico con observaciones \(y_1,...,y_n\).
  • Una variable predictora, \(X\), de tipo numérico con observaciones \(x_1,...,x_n\).
  • Una variable predictora, \(F\), de tipo categórico con \(I\) grupos o niveles distintos de tamaños muestrales \(n_1,n_2,...,n_I\), de forma que \(n = n_1 + n_2 + ... + n_I\), de forma que el vector de observaciones de la respuesta y de la predictora numérica se pueden escribir como:

\[(Y_1, Y_2,...,Y_I) = y_{11},\ldots,y_{1n_1},y_{21},\ldots,y_{2n_2},\ldots, y_{I1},\ldots,y_{In_I}\]

\[ (X_1, X_2,...,X_I) = x_{11},\ldots,x_{1n_1},x_{21},\ldots,x_{2n_2},\ldots,x_{I1},\ldots,x_{In_I} \]

donde el primer subíndice indica el nivel del factor y el segundo la posición dentro del conjunto de datos de dicho nivel del factor. * Conjunto \(\mu_i\) de medias de todas las observaciones de la respuesta asociadas con el nivel \(i\) del factor, es decir:

\[ \mu_i = \frac{\sum_{j = 1}^{n_j} y_{ij}}{n_i}; \text{ i = 1, 2,..., I} \]

  • Media global de la respuesta, \(\mu\), que se puede obtener como:

\[ \mu = \frac{\sum_{j = 1}^{I} \mu_{j}}{I} \]

  • Incrementos, \(\alpha_i\), de cada una de las medias de cada grupo con respecto a la media global, es decir:

\[ \alpha_i= \mu - \mu_i; \text{ i = 1, 2,..., I} \]

  • Pendiente común, \(\beta\), que representa la posible relación entre las variables de tipo numérico.

  • Pendientes diferentes entre las predictoras numéricas asociadas a cada nivel del factor, \(\gamma_i\) con \(i = 1,...,I\).

En esta situación el modelo que describe la posible relación entre respuesta y predictoras se puede escribir como: \[\begin{equation} y_{ij} = \alpha_0 + \alpha_i + \beta x_{ij} + \gamma_i x_{ij} + \epsilon_{ij};\quad i=1,...,I \quad \text{con} \quad \alpha_I = 0 \end{equation}\]

Que en forma matricial se puede escribir fácilmente (de forma análoga al ANOVA de un vía) sin más que considerar \(1_(n_i )={1,...,1}\) un vector de \(n_i\) unos, \(0_(n_i )={0,...,0}\) un vector de \(n_i\) ceros, para cada uno de los niveles i del factor la matriz de diseño viene dada por: \[ Y = \left( \begin{array}{c} Y_1 \\ Y_2 \\ \ldots \\ Y_I \\ \end{array} \right) = \left( \begin{array}{ccccccccc} 1_{n_1} & 1_{n_1} & 0_{n_1} & \ldots & 0_{n_1} & X_1 & X_1 & \ldots & 0_{n_1}\\ 1_{n_2} & 0_{n_2} & 1_{n_2} & \ldots & 0_{n_2} & X_2 & 0_{n_2} & \ldots & 0_{n_2}\\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & 0_{n_3}\\ 1_{n_I} & 0_{n_I} & 0_{n_I} &\ldots & 1_{n_I} & X_I & 0_{n_I} & \ldots & X_I\\ \end{array} \right) \left( \begin{array}{c} \alpha_0 \\ \alpha_1 \\ \alpha_2 \\ \ldots \\ \alpha_{I_1} \\ \beta \\ \gamma_1\\ \ldots\\ \gamma_I\\ \end{array} \right) + \left( \begin{array}{c} e_1 \\ e_2 \\ \ldots \\ e_n \\ \end{array} \right) = X\beta + \epsilon \]

donde las primeras \(I\) columnas representan el efecto del factor (ANOVA de una vía), la siguiente columna representa el efecto común entre predictoras numéricas (Regresión lineal simple), y las últimas \(I\) columnas representan el efecto distinto de la predictora numérica para cada nivel del factor. Estas columnas se obtienen fácilmente multiplicando la columna de \(X\) por cada de las columnas desde la \(2\) a la \(I\).

Las posibles modelos anidados que se pueden obtener a partir de la ecuación del modelo ANCOVA, así como sus interpretaciones se presentan a continuación. Además, se ofrece la interpretación de dichos modelos en términos de los coeficientes que resulten significativos. La estructura secuencial de contrastes pasa por estudiar si debemos considerar una pendiente distinta para la predictora numérica en función de los niveles del factor, que es equivalente a plantear el contraste:

\[\begin{equation} \left\{ \begin{array}{ll} H_0: & \gamma_1 = \gamma_2 = \ldots = \gamma_I = 0\\ H_a: & \mbox{Al menos hay una pendiente distinto de cero}\\ \end{array} \right. \end{equation}\]

  • Si rechazamos el contraste anterior tendremos un modelo de regresión entre la respuesta y la predictora numérica para cada nivel del factor, es decir, I rectas de regresión distintas con ecuaciones:

\[\begin{equation} \begin{array}{ll} y_{1j} &= (\alpha_0 + \alpha_1) + (\beta + \gamma_1) x_{1j} + \epsilon_{1j}\\ y_{2j} &= (\alpha_0 + \alpha_2) + (\beta + \gamma_2) x_{2j} + \epsilon_{2j}\\ \ldots &= \ldots \\ y_{Ij} &= (\alpha_0 + \alpha_I) + (\beta + \gamma_I) x_{Ij} + \epsilon_{Ij}\\ \end{array} \end{equation}\]

con interceptaciones \(\alpha_0 + \alpha_i\) y pendientes \(\beta + \gamma_i\) para cada uno de los \(I\) niveles del factor. Tenemos I modelos de regresión distintos (uno por cada nivel del factor), es decir, tenemos I rectas perpendiculares o que se cortan.

  • Si no rechazamos el contraste sobre los efectos \(\gamma_i\) tendremos un único modelo de regresión entre la respuesta y la predictora pero con diferentes interceptaciones (la pendiente es la misma), es decir, I rectas paralelas cuyas ecuaciones vienen dadas por:

\[\begin{equation} \begin{array}{ll} y_{1j} &= (\alpha_0 + \alpha_1) + \beta x_{1j} + \epsilon_{1j}\\ y_{2j} &= (\alpha_0 + \alpha_2) + \beta x_{2j} + \epsilon_{2j}\\ \ldots &= \ldots \\ y_{Ij} &= (\alpha_0 + \alpha_I) + \beta x_{Ij} + \epsilon_{Ij}\\ \end{array} \end{equation}\]

En la situación donde el contraste sobre los efectos \(\gamma_i\) es no significativo podemos definir diferentes modelos anidados en función de los incrementos del factor (`s) y la pendiente \(\beta\). Los contrastes son:

  • Contraste sobre el efecto del factor

\[\begin{equation} \left\{ \begin{array}{ll} H_0: & \alpha_1 = \alpha_2 = \ldots = \alpha_I = 0\\ H_a: & \mbox{Al menos hay un incremento distinto de cero}\\ \end{array} \right. \end{equation}\]

  • Contraste sobre el efecto de la predictora numérica

\[\begin{equation} \left\{ \begin{array}{ll} H_0: & \beta = 0\\ H_a: & \beta \neq 0 \end{array} \right. \end{equation}\]

  • Si rechazamos la hipótesis nula sobre los \(\alpha_i\) pero no rechazamos la hipótesis nula sobre \(\beta\) diríamos que no hay efecto de la predictora numérica pero que si podemos establecer diferencias entre las medias de la respuesta dadas por los diferentes niveles del factor (Modelo ANOVA de una vía).
  • Si no rechazamos la hipótesis nula sobre los \(\alpha_i\) pero si rechazamos la hipótesis nula sobre \(\beta\) diríamos que no hay efecto del factor pero si que podemos establecer un modelo de regresión entre la respuesta y la predictora numérica (Modelo de Regresión Lineal Simple).
  • Si no rechazamos la hipótesis nula sobre los \(\alpha_i\) y no rechazamos la hipótesis nula sobre \(\beta\) estaríamos ante un modelo nulo donde el comportamiento de la respuesta no viene explicado por las predictoras consideradas.

Las ecuaciones de los modelos resultantes son:

  • Modelo ANOVA:

\[\begin{equation} \begin{array}{ll} y_{1j} &= (\alpha_0 + \alpha_1) + \epsilon_{1j}\\ y_{2j} &= (\alpha_0 + \alpha_2) + \epsilon_{2j}\\ \ldots &= \ldots \\ y_{Ij} &= (\alpha_0 + \alpha_I) + \epsilon_{Ij}\\ \end{array} \end{equation}\]

  • Modelo de Regresión Lineal Simple:

\[\begin{equation} \begin{array}{ll} y_{ij} &= \alpha_0 + \beta x_{ij} + \epsilon_{ij}\\ \end{array} \end{equation}\]

Atendiendo a las ecuaciones obtenidas las posibles modelizaciones resultantes del modelo saturado son:

  • I rectas que se cortan (modelo con interacción factor - numérica ).
  • I rectas paralelas (modelo sin interacción pero con factor y numérica).
  • Una recta de regresión (modelo con la predictora numérica únicamente).
  • Diferencias entre las medias (modelo con el factor únicamente).

Vemos gráficamente los posibles modelos que podemos construir sobre los dos ejemplos presentados anteriormente:

  • Datos de tiempo de vida: ¿Cómo interpretamos los gráficos obtenidos?

  • Datos de longevidad: ¿Cómo interpretamos los gráficos obtenidos?

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. Estas hipótesis se adaptarán en función del tipo de modelo que finalmente alcanzemos en el proceso de contraste. Todos ellos se pueden resolver mediante un test \(F\), y más adelante veremos como construir la secuencia de modelos para decidir sobre el modelo final.

2.3 Especificación del modelo en R

El modelo ANCOVA planteado en el punto anterior para un factor F y una variable predictora numérica X, se puede escribir en R en su formato reducido como:

\[Y \sim F + X + F:X\]

donde:

  • \(F\) representa el efecto del factor, es decir, comparamos si las medias de la respuesta para cada grupo pueden considerarse iguales.
  • \(X\) representa el efecto de regresión asociado con la variable numérica, es decir, la respuesta y X están relacionadas mediante una única pendiente que deberemos estimar.
  • \(F:X\) representa el efecto de interacción entre predictoras, es decir, que la respuesta se relaciona con la predictora numérica mediante tantas curvas (generalmente líneas) como niveles tenga el factor \(F\).

A continuación, se presentan los modelos reducidos para diferentes situaciones experimentales en el número y tipo de variables predictoras:

  • Modelo para dos factores (\(F_1\) y \(F_2\)) y una numérica (\(X\))

\[Y \sim F_1 + F_2 + F_1:F_2 + X + F_1:X + F_2:X + F_1:F_2:X\] o en forma más simplificada \(Y \sim F_1*F_2*X\)

  • Modelo para un factor (\(F\)) y dos numéricas (\(X_1\) y \(X_2\))

\[Y \sim F + X_1 + X_2 + F:X_1 + F:X_2\] o en forma más simplificada \(Y \sim F*(X_1 + X_2)\)

  • Modelo para dos factores (\(F_1\) y \(F_2\)) y dos numéricas (\(X_1\) y \(X_2\))

\[Y \sim F_1*F_2*(X_1 + X_2)\]

Como se puede ver la complejidad del modelo aumenta sustancialmente con la consideración de más variables predictoras. La forma de expresar el modelo saturado debe contemplar tanto os efectos principales asociados a cada predictora, como los efectos de interacción entre factores y entre factores y numéricas.

2.4 Estimación y Selección del modelo

Dado que hemos expresado el modelo ANCOVA como un modelo de tipo lineal con una ecuación similar a los modelos de regresión múltiple, la estimación de los parámetros del modelo se puede realizar utilizando las ecuaciones normales. En el proceso de selección del mejor modelo actuaremos como en los modelos ANOVA, es decir partiremos del modelo saturado y veremos que efectos puede ser considerados como irrelevantes, y por tanto deben desaparecer del modelo. Esta selección nos permitirá elegir el modelo final resultante. En este caso más sencillo podemos escribir todos los modelos posibles y elegir el mejor de ellos, bien mediante la comparación con el test F o con el AIC, pero veremos como utilizar los procedimientos secuenciales automáticos para problemas más complejos.

2.4.1 Ejemplos

Realizamos la selección y estimación del mejor modelo para cada uno de los conjuntos de datos considerados.

2.4.1.1 Datos de tiempo de vida

Construimos el modelo saturado y seleccionamos mediante el test \(F\).

## 
## 
##                                  Elimination Summary                                  
## -------------------------------------------------------------------------------------
##         Variable                               Adj.                                      
## Step           Removed           R-Square    R-Square     C(p)       AIC        RMSE     
## -------------------------------------------------------------------------------------
##    1    velocidad:herramienta      0.8969      0.8847    3.9652    106.6591    3.0919    
## -------------------------------------------------------------------------------------

El proceso de selección identifica el efecto de interacción entre velocidad y herramienta como no significativo, de forma que el modelo final viene dado por: \[vida \sim velocidad + herramienta\] Ajustamos el modelo y estudiamos los parámetros obtenidos:

  vida
Predictors Estimates CI
(Intercept) 36.93 29.56 – 44.30
velocidad -0.03 -0.04 – -0.02
herramienta [B] 14.67 11.75 – 17.58
Observations 20

En este caso tenemos un modelo con dos rectas paralelas (una por cada tipo de herramienta) cuyas ecuaciones de estimación vienen dadas por:

\[ \left\{ \begin{array}{lll} \mbox{Herramienta A}:& \widehat{Vida_{A}} = 36.93 + 0 - 0.03*velocidad &= 36.93 -0.03*velocidad\\ \mbox{Herramienta B}:& \widehat{Vida_{B}} = 36.93 + 14.67 - 0.03*velocidad &= 51.60 - 0.03*velocidad \\ \end{array} \right. \]

Tenemos una interceptación mayor para la herramienta A indicando que la recta asociada con dicha herramienta está por encima de la de la herramienta B. Además, la pendiente negativa asociada con la velocidad indica que conforme aumenta esta disminuye el tiempo de vida. Ese efecto es mayor si usamos la herramienta de tipo B dado que la recta ajustada está por debajo de la de la herramienta A.

2.4.1.2 Datos de longevidad

Construimos el modelo saturado y seleccionamos mediante el test \(F\).

## 
## 
##                                Elimination Summary                                 
## ----------------------------------------------------------------------------------
##         Variable                          Adj.                                        
## Step        Removed         R-Square    R-Square     C(p)        AIC        RMSE      
## ----------------------------------------------------------------------------------
##    1    thorax:actividad      0.6527       0.638    -3.7881    943.8165    10.5394    
## ----------------------------------------------------------------------------------

El proceso de selección identifica el efecto de interacción entre thorax y actividad como no significativo, de forma que el modelo final viene dado por: \[longevidad \sim thorax + actividad\] Ajustamos el modelo y estudiamos los parámetros obtenidos:

  longevidad
Predictors Estimates CI
(Intercept) -44.61 -65.53 – -23.69
thorax 134.34 109.13 – 159.55
actividad [G2] -4.14 -10.13 – 1.86
actividad [G3] -1.50 -7.48 – 4.47
actividad [G4] -11.15 -17.15 – -5.16
actividad [G5] -24.14 -30.12 – -18.17
Observations 124

En este caso tenemos un modelo con dos rectas paralelas (una por cada tipo de actividad) cuyas ecuaciones de estimación vienen dadas por:

\[ \left\{ \begin{array}{ll} \mbox{G1}:& \widehat{logevidad_{G1}} = - 44.61 + 134.34*thorax\\ \mbox{G2}:& \widehat{logevidad_{G2}} = - 48.75 + 134.34*thorax \\ \mbox{G3}:& \widehat{logevidad_{G3}} = - 46.11 + 134.34*thorax \\ \mbox{G4}:& \widehat{logevidad_{G4}} = - 55.76 + 134.34*thorax \\ \mbox{G5}:& \widehat{logevidad_{G5}} = - 68.75 + 134.34*thorax \\ \end{array} \right. \]

Se observa un efecto de aumento de la longevidad conforme aumenta la longitud del thorax. El grupo con mayor longevidad es el G1, ya que tiene la interceptación más grande, mientras que el que tiene menor longevidad es G5. El orden vendría dado por \(G1 > G3 > G2 > G4 > G5\).

2.5 Diagnóstico del modelo

En este caso el diagnóstico es similar al de los modelos de regresión pero teniendo en cuenta que las hipótesis se deben verificar para los residuos asociados a cada nivel del factor (si este está presente en el modelo). Las hipótesis son linealidad, normalidad y varianza constante. Para verificar las hipótesis utilizamos los procedimientos gráficos y tests tratados en las unidades anteriores.

De nuevo utilizaremos la distancia der Cook para establecer posibles observaciones influyentes.

2.5.1 Ejemplos

A continuación, realizamos el diagnóstico de los ejemplos que venimos trabajando en esta unidad. Para obtener el diagnóstico partimos del modelo obtenido en la sección anterior y realizamos el gráfico de residuos versus ajustados y los tests de diagnóstico.

Como en los dos ejemplos el factor considerado se encuentra en el modelo final, utilizaremos el test de Levene para verificar la igualdad de varianzas entre los niveles del factor.

2.5.1.1 Datos de tiempo de vida

Obtenemos los valores de diagnóstico y realizamos los correspondientes tests de hipótesis y análisis de influencia.

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9715         0.7858 
## Kolmogorov-Smirnov        0.1232         0.8859 
## Cramer-von Mises          1.4412          2e-04 
## Anderson-Darling          0.2652         0.6555 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.3888  0.254
##       18

Los residuos tienen un comportamiento aleatorio, se verifican las hipótesis del modelo y no se detectan observaciones influyentes. El modelo obtenido parece adecuado para estudiar el tiempo de vida en función de la velocidad y la herramienta utilizada.

2.5.1.2 Datos de longevidad

Obtenemos los valores de diagnóstico y realizamos los correspondientes tests de hipótesis y análisis de influencia.

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9916         0.6607 
## Kolmogorov-Smirnov        0.0538         0.8654 
## Cramer-von Mises         10.2413         0.0000 
## Anderson-Darling          0.3224         0.5241 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   4  1.2925 0.2769
##       119

Aunque se verifican las hipótesis del modelo y no se detectan observaciones influyentes, si que es cierto que los gráficos de residuos muestran cierto comportamiento de embudo con variabilidades más pequeñas en valores más pequeños de thorax, y mayor dispersión al aumentar la longitud del thorax. Sería recomendable probar Box-Cox para tratar de obtener una transformación de la respuesta que nos permita obtener gráficos sin esos efectos indeseables.

La transformación raíz cuadrada parece adecuada en esta situación. Obtenemos la nueva variable y ajustamos de nuevo el modelo.

## 
## 
##                                Elimination Summary                                
## ---------------------------------------------------------------------------------
##         Variable                          Adj.                                       
## Step        Removed         R-Square    R-Square     C(p)        AIC        RMSE     
## ---------------------------------------------------------------------------------
##    1    thorax:actividad      0.6868      0.6736    -3.0694    267.2943    0.6888    
## ---------------------------------------------------------------------------------

De nuevo el modelo seleccionado prescinde del efecto de interacción. Ajustamos y estudiamos el nuevo modelo.

  rlongevidad
Predictors Estimates CI
(Intercept) 0.34 -1.03 – 1.71
thorax 9.41 7.77 – 11.06
actividad [G2] -0.30 -0.69 – 0.09
actividad [G3] -0.12 -0.51 – 0.27
actividad [G4] -0.76 -1.15 – -0.37
actividad [G5] -1.73 -2.12 – -1.34
Observations 124

¿Cuáles son las ecuaciones de estimación en este caso?

El proceso de diagnóstico para el nuevo modelo permite verificar el cumplimiento de las hipótesis y la leve mejora de los gráficos de residuos.

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.9951         0.9465 
## Kolmogorov-Smirnov        0.0474         0.9432 
## Cramer-von Mises         10.7431         0.0000 
## Anderson-Darling          0.2137         0.8484 
## -----------------------------------------------
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   4  0.6856 0.6033
##       119

2.6 Predicción

El proceso de predicción en este tipo de modelos es muy simple a partir de las ecuaciones de los modelos obtenidos. De hecho, en las secciones anteriores ya hemos visto gráficamente la predicción para todos estos modelos en los ejemplos que hemos ido trabajando. Básicamente, si queremos obtener una predicción especifica deberemos dar un valor del factor y otro de la predictora numérica para calcular el valor de predicción y su correspondiente intervalo. En este caso nos imitamos a representar las bandas de predicción que podemos obtener para cada modelo.

2.6.1 Ejemplos

Construimos la predicción y bandas de confianza para cada uno de los ejemplos a partir de los modelos determinados tras la fase de diagnóstico.

2.6.1.1 Datos de tiempo de vida

A continuación, se presentan las rectas de predicción para el modelo ajustado a este banco de datos.

2.6.1.2 Datos de longevidad

A continuación, se presentan las rectas de predicción para el modelo ajustado a este banco de datos.

3 Ejercicios

A continuación, se presenta una colección de ejercicios referidos a los modelos ANOVA y ANCOVA. Los pasos a seguir para la obtención del modelo son los que hemos ido desarrollando: representación gráfica y propuesta de modelo, ajuste, bondad de ajuste, diagnóstico y predicción. En caso de encontrar problemas con el diagnóstico se deberá proponer un nuevo modelo alternativo.

No olvides cargar las librerías para realizar los ejercicios.

Ejercicio 1. Se realiza una investigación para conocer los niveles de fosfato inorgánico en plasma (mg / dl) una hora después de una prueba de tolerancia a la glucosa estándar para sujetos obesos, con o sin hiperinsulinemia, y controles. Los datos corresponden con la tabla 6.18 de Dobson (2002).

Ejercicio 2. Se lleva a cabo un estudio sobre el contenido promedio de grasa (en porcentaje) en la leche del ganado de cinco razas distintas canadienses. Para ello se consideran veinte ejemplares de pura raza (diez de dos años y diez maduras de más de cuatro años de cada una de cinco razas.

Ejercicio 3. Los datos que se presentan a continuación [@Dobson02] 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.

Ejercicio 4. Se realiza un estudio experimental para estudiar el rendimiento de un proceso químico en función de la concentración del compuesto base (A, B, C y D), el catalizador usado (C1, C2, y C3), y la temperatura usada en el proceso (T1, T2). Se quiere estudiar como afectan estos factores en el rendimiento del proceso.

Ejercicio 5. Se quiere evaluar la eficacia de distintas dosis de un fármaco contra la hipertensión arterial, comparándola con la de una dieta sin sal. Para ello se seleccionan al azar 25 hipertensos y se distribuyen aleatoriamente en 5 grupos. Al primero de ellos no se le suministra ningún tratamiento (T1), al segundo una dieta con un contenido pobre en sal (T2), al tercero una dieta sin sal (T3), al cuarto el fármaco a una dosis determinada (T4) y al quinto el mismo fármaco a otra dosis (T5). Las presiones arteriales (Presion) sistólicas de los 25 sujetos al finalizar los tratamientos son:

Ejercicio 6. La convección es una forma de transferencia de calor por los fluidos debido a sus variaciones de densidad por la temperatura; las partes calientes ascienden y las frías descienden formando las corrientes de convección que hacen uniforme la temperatura del fluido. Se ha realizado un experimento para determinar las modificaciones de la densidad de fluido al elevar la temperatura en una determinada zona. Los resultados obtenidos han sido los siguientes:

Ejercicio 7. Un laboratorio de reciclaje controla la calidad de los plásticos utilizados en bolsas. Se desea contrastar si existe variabilidad en la calidad de los plásticos que hay en el mercado. Para ello, se eligen al azar cuatro plásticos y se les somete a una prueba para medir el grado de resistencia a la degradación ambiental. De cada plástico elegido se han seleccionado ocho muestras y los resultados de la variable que mide la resistencia son los de la tabla adjunta.

Ejercicio 8. Se realiza un estudio sobre el efecto que produce la descarga de aguas residuales de un planta sobre la ecología del agua natural de un río. En el estudio se utilizaron dos lugares de muestreo. Un lugar está aguas arriba del punto en el que la planta introduce aguas residuales en la corriente; el otro está aguas abajo. Se tomaron muestras durante un periodo de cuatro semanas y se obtuvieron los datos sobre el número de diatomeas halladas. Los datos se muestran en la tabla adjunta:

Ejercicio 9. Disponemos de los datos de peso de 24 niños recién nacidos (peso), su sexo (sexo; “H” = Hombres y “M” = Mujeres) y la edad de sus madres (edad). Nos gustaría ser capaces de determinar un modelo que explique el peso de los niños recién nacidos en función de su sexo y de la edad de sus madres.

Ejercicio 10. Se lleva a cabo una investigación sobre diversas malformaciones del sistema nervioso central registradas en nacidos vivos en Gales del Sur, Reino Unido. El estudio fue diseñado para determinar el efecto de la dureza del agua sobre la incidencia de tales malformaciones. La información registrada son: NoCNS = recuento de nacimientos sin problema CNS; An = conteo de nacimientos de Anencephalus; Sp = conteo de nacimientos de espina bífida; Otro = recuento de otros nacimientos del SNC; Agua = endurecimiento del agua; Trabajo = un factor con niveles Manual no manual en función del tipo de trabajo realizado por los padres Se está interesado en predecir el número total de malformaciones en función de la calidad del agua y el trabajo realizado por los padres.

Ejercicio 11. Se ha realizado un estudio para establecer la calidad de los vinos de la variedad Pino Noir en función de un conjunto de características analizadas. Las características analizadas son claridad, aroma, cuerpo, olor y matiz. Para medir la calidad se organiza una cata ciega a un conjunto de expertos y se calcula la puntuación final de cada vino a partir de la información de todos ellos. Además se registra la región (region) de procedencia del vino por si puede influir en la calidad del vino.

Ejercicio 12. Un fabricante de ropa que suministra uniformes militares debe cortar chaquetas, camisas, pantalones (variable Prenda) y otros complementos (en muchas tallas diferentes), de rollos de tela. La tela es cara, de modo que el desperdicio (Desperdicio) tiene un efecto muy grande en los beneficios. El fabricante tiene que elegir entre tres máquinas (Maquina) cortadoras asistidas por computadora: A, B y C. El fabricante decide experimentar haciendo que cada máquina corte varios lotes de chaquetas, varios más de camisas otros más de pantalones y complementos para determinar que máquina es más eficiente en cada caso, es decir, tratamos de conocer el desperdicio que se producirá para cada prenda y máquina.

Ejercicio 13. Una empresa dedicada a la fabricación de aislantes térmicos y acústicos establece un experimento que mide la pérdida de calor (Calor) a través de cuatro tipos diferentes de cristal para ventanas (Cristal) utilizando cinco graduaciones diferentes de temperatura exterior (TempExt). Se prueban tres hojas de cristal en cada graduación de temperatura, y se registra la pérdida de calor para cada hoja.

Ejercicio 14. El grupo de asesores LearnStatistics ha realizado un estudio para comprobar si las empresas destinan parte de los beneficios de sus ventas en la formación de sus empleados para mejorar su competitividad. Para ellos se recoge la información sobre ventas (Ventas) en miles de euros, capital invertido en formación (Capital) en miles de euros, y el nivel de productividad de la empresa establecido por un asesor externo (Productividad).

Javier Morales, email: j.morales@.umh.es.

Mª Asunción Martínez, email: am.mayoral@umh.es.