Modelos Estadísticos. Grado Biotecnología

Modelos ANCOVA

En este tema se presentan los modelos ANCOVA. o modelos de análisis de la covarianza. Esto modelos surgen cuando entre las variables predictoras 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.

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_{11},...,X_{n1}\).
  • 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\).

El modelo teórico para esta situación viene dado por:

\[ Y = \alpha_0 + F + X + F:X +\epsilon\]

donde:

  • \(\alpha_0\) representa un efecto común a todo el modelo
  • \(F\) representa el efecto del factor, es decir, comparamos si la media 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\).

Tenemos las siguientes situaciones:

  • estimación de \(k\) curvas, con al menos dos de ellas no paralelas, para explicar el comportamiento de \(Y\) (efecto interacción significativo)
  • estimación de \(k\) curvas paralelas para explicar el comportamiento de \(Y\) (efecto interacción no significativo, pero efecto de \(X\) y \(F\) significativo)
  • estimación de una única curva para explicar el comportamiento de \(Y\) (efecto interacción y \(F\) no significativo, pero efecto de \(X\) significativo)
  • estimación únicamente de diferencia entre las medias para explicar el comportamiento de \(Y\) (efecto interacción y \(X\) no significativo, pero efecto de \(F\) significativo)
  • modelo sin capacidad explicativa, es decir, todos los efectos no significativos.

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, ya que dependen especialmente de si el factor está presente o no en el modelo. Otro aspecto fundamental es la comparación de modelos para determinar la situación experimental ante la que nos encontramos.

Veamos unos ejemplos antes de pasar a estudiar cada uno de los aspectos de este tipo de modelos.

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)
library(modelr)

Ejemplo 11. Se desea estudiar el tiempo de vida de una pieza cortadora de dos tipos, A y B, en función de la velocidad del torno 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.

vida <- c(18.73,    14.52,  17.43,  14.54,  13.44,  25.39,  13.34,  22.71,  12.68,  19.32,  30.16,  27.09,  25.40,  26.05,  33.49,  35.62,  26.07,  36.78,  34.95,  43.67)
velocidad <- c(610,950,720,840,980,530,680,540,980,730,670,770,880,1000,    760,590,910,650,810,500) 
tipo <- c(rep("A",10), rep("B",10))
ejemplo11 <- data.frame(vida,velocidad,tipo)

Tenemos dos variables predictoras (una de cada tipo) de forma que la secuencia de modelos (representados gráficamente) que nos podemos plantear son (obviamos el modelo sin capacidad explicativa:

  • Modelo con efecto del factor únicamente. Consideramos únicamente el modelo en que tendríamos el factor:

\[vida \sim tipo\]

# gráfico con error estándar
ggplot(ejemplo11, aes(x = tipo, y= vida)) + 
   geom_boxplot() +
    theme_bw()

Se aprecian diferencias en las medias del tiempo de vida en función del tipo de pieza.

  • Modelos con regresores. Obtenemos los tres modelos posibles de la combinación de la velocidad y el tipo de pieza

\[\begin{array}{ll} M0: & vida \sim velocidad\\ M1: & vida \sim velocidad + tipo\\ M2: & vida \sim velocidad + tipo + velocidad:tipo\\ \end{array}\]

¿Cómo se interpreta cada uno de estos modelos? ¿Cuantas tendencias ajustaremos para cada uno de ellos?

# Comenzamos desde le modelo más sencillo
# Modelo con una única recta
M0 <- lm(vida ~ velocidad,data = ejemplo11)
# M1: modelo con rectas paralelas
M1 <- lm(vida ~ tipo + velocidad,data = ejemplo11)
# M2: modelo con rectas no paralelas
M2 <- lm(vida ~ tipo + velocidad + tipo:velocidad,data = ejemplo11)
# grid de valores para construir los modelos
grid <- ejemplo11 %>% data_grid(tipo,velocidad) %>% gather_predictions(M0,M1,M2)
# Gráfico
ggplot(ejemplo11,aes(velocidad,vida,colour = tipo)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  facet_wrap(~ model) +
  labs(x = "Velocidad del torno", y = "Tiempo de vida") + 
  theme_bw()

Se puede apreciar que resulta necesario ajustar un modelo por cada nivel del factor (modelos M1 y M2), pero no está tan claro que la recta no tengan que ser paralelas. Se observa cierta pendiente diferente para el tipo B (modelo M2), pero habrá que comprobar si dicho efecto es significativo.

Ejemplo 12. Partridge y Farquhar realizan un experimento para relacionar la vida útil de las moscas de la fruta con su actividad sexual. 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 ya que se sospecha que afecta directamente a la longevidad de las moscas.

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,  44L)
actividad <- c(rep("G1",24),rep("G2",25),rep("G3",25),rep("G4",25),rep("G5",25))
ejemplo12 <- data.frame(thorax,longevidad,actividad)
  • Modelo con efecto del factor únicamente

\[longevidad \sim actividad\]

# gráfico con error estándar
ggplot(ejemplo12, aes(x = actividad, y= longevidad)) + 
   geom_boxplot() +
    theme_bw()

Se observa que no parece haber diferencias entre los diferentes grupos considerados, salvo tal vez para el grupo G5 con una longevidad inferior.

  • Modelos con regresores. Obtenemos los tres modelos posibles de la combinación de la velocidad y el tipo de pieza

\[\begin{array}{ll} M0: & longevidad \sim thorax\\ M1: & longevidad \sim thorax + actividad\\ M2: & longevidad \sim thorax + actividad + thorax:actividad\\ \end{array}\]

¿Cómo se interpreta cada uno de estos modelos? ¿Cuantas tendencias ajustaremos para cada uno de ellos?

# Comenzamos desde le modelo más sencillo
# Modelo con una única recta
M0 <- lm(longevidad ~ thorax,data = ejemplo12)
# M1: modelo con rectas paralelas
M1 <- lm(longevidad ~ actividad + thorax,data = ejemplo12)
# M2: modelo con rectas no paralelas
M2 <- lm(longevidad ~ actividad + thorax + actividad:thorax,data = ejemplo12)
# grid de valores para construir los modelos
grid <- ejemplo12 %>% data_grid(actividad,thorax) %>% gather_predictions(M0,M1,M2)
# Gráfico
ggplot(ejemplo12,aes(thorax,longevidad,colour = actividad)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  facet_wrap(~ model) +
  labs(x = "Longitud del tórax", y = "Longevidad") + 
  theme_bw()

Se observa como la nube de dispersión de los puntos, donde se ha identificado cada observación por su grupo de actividad, se encuentra muy mezclada. No se observan muchas diferencias entre los modelos M1 y M2, y de hecho, hasta parece posible considerar el modelo M0 como válido.

Ejemplo 13. Se conoce como infiltración el proceso por el cual el agua (riego o lluvia) se va introduciendo bajo la superficie de un terreno cultivado. Este proceso es vital para determinar las cantidades de agua de riego necesarias, para mantener el terreno en condiciones óptimas. Un parámetro habitual que sirve para estudiar dicho proceso es la carga hidráulica. Este depende tanto de la profundidad de la infiltración como del procedimiento de riego usado. Se diseña un experimento para estudiar la carga hidráulica de un terreno bajo diferentes condiciones de riego (denominados tratamientos).

tratamiento <- c(rep("A",15),rep("B",15),rep("C",15))
profundidad <-c(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150)
cargahid <- c(-406.90,-345.70,-335.50,-315.10,-304.90,-315.10,-323.26,-335.50,-345.70,-362.02,-374.26,-386.50,-421.18,-435.46,-447.70,-896.50,-737.38,-653.74,-470.14,-406.90,-388.54,-396.70,-396.70,-396.70,-406.90,-419.14,-437.50,-468.10,-466.06,-492.58,-896.50,-855.70,-818.98,-788.38,-678.22,-590.50,-545.62,-515.02,-498.70,-496.66,-517.06,-555.82,-619.06,-623.14,-623.14)
ejemplo13 <- data.frame(tratamiento,profundidad,cargahid)
  • Modelo con efecto del factor únicamente
# gráfico con error estándar
ggplot(ejemplo13, aes(x = tratamiento, y= cargahid)) + 
   geom_boxplot() +
    theme_bw()

Se aprecian diferencias entre los tres tratamientos utilizados, aunque la variabilidad de cada grupo pueden parecer diferentes.

  • Modelos con regresores. Obtenemos los tres modelos posibles de la combinación de la velocidad y el tipo de pieza

\[\begin{array}{ll} M0: & cargahid \sim profundidad\\ M1: & cargahid \sim profundidad + tratamiento\\ M2: & cargahid \sim profundidad + tratamiento + profundidad\\ \end{array}\]

¿Cómo se interpreta cada uno de estos modelos? ¿Cuantas tendencias ajustaremos para cada uno de ellos?

# Comenzamos desde le modelo más sencillo
# Modelo con una única recta
M0 <- lm(cargahid ~ profundidad,data = ejemplo13)
# M1: modelo con rectas paralelas
M1 <- lm(cargahid ~ tratamiento + profundidad,data = ejemplo13)
# M2: modelo con rectas no paralelas
M2 <- lm(cargahid ~ tratamiento + profundidad + tratamiento:profundidad,data = ejemplo13)
# grid de valores para construir los modelos
grid <- ejemplo13 %>% data_grid(tratamiento,profundidad) %>% gather_predictions(M0,M1,M2)
# Gráfico
ggplot(ejemplo13,aes(profundidad,cargahid,colour = tratamiento)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  facet_wrap(~ model) +
  labs(x = "Profundidad", y = "Carga hidráulica") + 
  theme_bw()

Se observa de forma muy clara que la relación entre la carga hidráulica y la profundidad no se puede analizar mediante una recta. Resulta más factible considerar un polinomio de grado 3 (se observan dos cambios en la curvatura) para ajustar las tendencias observadas en cualquiera de los tres modelos. Proponemos un nuevo conjunto de modelos basados en un modelo polinómico de grado 3.

\[\begin{array}{ll} M0: & cargahid \sim profundidad + profundidad^2 + profundidad^3\\ M1: & cargahid \sim profundidad + profundidad^2 + profundidad^3 + tratamiento\\ M2: & cargahid \sim (profundidad + profundidad^2 + profundidad^3)*tratamiento \\ \end{array}\]

Utilizamos el símbolo \(*\) para construir un modelo con todos los efectos posibles, es decir,

\[ X + F +X:F = X*F\] Tenemos un modelos con las dos predictoras y la interacción correspondiente

¿Cómo se interpreta cada uno de estos modelos? ¿Cuantas tendencias ajustaremos para cada uno de ellos?

# Comenzamos desde le modelo más sencillo
# Modelo con una única recta
M0 <- lm(cargahid ~ profundidad + I(profundidad^2) + I(profundidad^3),data = ejemplo13)
# M1: modelo con rectas paralelas
M1 <- lm(cargahid ~ tratamiento + profundidad + I(profundidad^2) + I(profundidad^3), data = ejemplo13)
# M2: modelo con rectas no paralelas
M2 <- lm(cargahid ~ tratamiento*(profundidad + I(profundidad^2) + I(profundidad^3)), data = ejemplo13)
# grid de valores para construir los modelos
grid <- ejemplo13 %>% data_grid(tratamiento,profundidad) %>% gather_predictions(M0,M1,M2)
# Gráfico
ggplot(ejemplo13,aes(profundidad,cargahid,colour = tratamiento)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  facet_wrap(~ model) +
  labs(x = "Profundidad", y = "Carga hidráulica") + 
  theme_bw()

Se aprecia claramente que el modelo M2 es el que parece adaptarse al comportamiento de la respuesta.

En este tipo de modelos procederemos de forma un poco distinta a lo hecho en el resto de modelos lineales. Antes de realizar la estimación del modelo, obtendremos el mejor modelo utilizando la tabla ANOVA o el criterio AIC para comparar los diferentes modelos anidados. Pasaremos después al diagnóstico del modelo.

Selección del modelo

Obtenemos el mejor modelo comparando el modelo saturado con modelos más simples. En este caso el modelo saturado es aquel que contiene todas las predictoras junto con los efectos de interacción. El primer paso siempre consiste en testar si podemos prescindir de las interacciones, para poder plantear un modelo donde todas las curvas o tendencias sean paralelas.

Ejemplo 11

Consideramos los modelos: \[\begin{array}{ll} M1: & vida \sim velocidad + tipo\\ M2: & vida \sim velocidad + tipo + velocidad:tipo\\ \end{array}\]

# M1: modelo con rectas paralelas
M1 <- lm(vida ~ tipo + velocidad,data = ejemplo11)
# M2: modelo con rectas no paralelas
M2 <- lm(vida ~ tipo + velocidad + tipo:velocidad,data = ejemplo11)
## Comparación de ambos modelos mediante la tabla ANOVA
## Siempre debemos colocar en primer luagr el modelo más simple
# Comparación con ANOVA
anova(M1,M2)
# Comparación con AIC
AIC(M1,M2)

El pvalor asociado con el test F de ANOVA resulta no significativo, lo que implica que ambos modelos pueden considerarse con la misma capacidad explicativa. En este tipo de situaciones siempre nos quedaremos con el modelo más sencillo que correspondería la M1. Sin embargo, si utilizamos el criterio AIC podemos ver que el valor más bajo se obtiene para el modelo 2, indicando que este es el mejor. ¿Cuál elegir? Esta pregunta no tiene una contestación nada fácil y debe ser el investigador el que decida finalmente. El procedimiento habitual pasa por hacer el diagnóstico de ambos modelos y si los dos cumplen las hipótesis nos quedaríamos con el más simple.

El test ANOVA esta diseñado para testar modelos anidados de forma que no resulta posible realizar la comparación de todos ellos de una vez. Por el contrario, el AIC nos permite valorar todos los modelos y seleccionar aquel con un valor más bajo. Procedamos con todos los modelos:

# M0: modelo sin predictoras
M0 <- lm(vida ~ 1,data = ejemplo11)
# M1: modelo con factor únicamente
M1 <- lm(vida ~ tipo,data = ejemplo11)
# M2: modelo con predictora numérica
M2 <- lm(vida ~ velocidad,data = ejemplo11)
# M3: modelo con rectas paralelas (efectos aditivos)
M3 <- lm(vida ~ tipo + velocidad,data = ejemplo11)
# M4: modelo con rectas no paralelas (efecto interacción)
M4 <- lm(vida ~ tipo + velocidad + tipo:velocidad,data = ejemplo11)
# Comparación con AIC
AIC(M0,M1,M2,M3,M4)

El AIC más pequeño se encuentra para el modelo con interacción. Obtenemos ahora las curvas o tendencias estimadas para este modelo:

tidy(M4)

que nos proporciona una tendencia para cada uno de los niveles del factor y con diferentes pendientes:

\[\begin{array}{ll} \text{Vida piezas tipo A:} & \widehat{y} = 32.87 - 0.02 * velocidad\\ \text{Vida piezas tipo B:} & \widehat{y} = 32.87 + 23.87 - 0.02 * velocidad - 0.01 * velocidad = 56.75 - 0.03 * velocidad\\ \end{array}\]

El gráfico del ajuste es:

# grid de valores para construir los modelos
grid <- ejemplo11 %>% data_grid(tipo,velocidad) %>% gather_predictions(M4)
# Gráfico
ggplot(ejemplo11,aes(velocidad,vida,colour = tipo)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  labs(x = "Velocidad del torno", y = "Tiempo de vida") + 
  theme_bw()

Ejemplo 12

Construimos la secuencia de modelos y los valoramos con el AIC

# Modelo sin efectos
M0 <- lm(longevidad ~ 1, data = ejemplo12)
# M1: modelo con factor
M1 <- lm(longevidad ~ actividad, data = ejemplo12)
# M2: modelo con predictora numérica
M2 <- lm(longevidad ~ thorax, data = ejemplo12)
# M3 : modelo aditivo
M3 <- lm(longevidad ~ actividad + thorax, data = ejemplo12)
# M4 : modelo con interacción
M4 <- lm(longevidad ~ actividad + thorax + actividad:thorax,data = ejemplo12)
# Resultado
AIC(M0,M1,M2,M3,M4)

El AIC más pequeño se obtiene para el modelo aditivo donde ajustamos tendencias paralelas para cada nivel del factor. Veamos el modelo obtenido:

tidy(M3)

\[\begin{array}{ll} \text{longevidad actividad G1:} & \widehat{y} = -44.61 + 134.34 * torax\\ \text{longevidad actividad G2:} & \widehat{y} = -44.61 - 4.14 + 134.34 * torax = -48.75 + 134.34 * torax\\ \text{longevidad actividad G3:} & \widehat{y} = -44.61 - 1.50 + 134.34 * torax = -46.11 + 134.34 * torax\\ \text{longevidad actividad G4:} & \widehat{y} = -44.61 - 11.15 + 134.34 * torax= -55.76 + 134.34 * torax\\ \text{longevidad actividad G5:} & \widehat{y} = -44.61 - 24.14 + 134.34 * torax= -68.75 + 134.34 * torax\\ \end{array}\]

Todas las rectas tienen la misma pendiente (coeficiente asociado a tórax) pero con distinto punto de origen debido al tratamiento al que son sometidos los sujetos.

El gráfico del ajuste es:

# grid de valores para construir los modelos
grid <- ejemplo12 %>% data_grid(actividad,thorax) %>% gather_predictions(M3)
# Gráfico
ggplot(ejemplo12,aes(thorax,longevidad,colour = actividad)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  labs(x = "Longitud del tórax", y = "Longevidad") + 
  theme_bw()

Ejemplo 13

En este caso tenemos muchos posibles modelos teniendo en cuenta la tendencia polinómica observada. Aceleramos el proceso de selección (esto siempre es posible si tenemos cierta certeza sobre el modelo) para encontrar el mejor modelo. Consideramos diferentes modelos polinómicos (grados distintos) para valorar el mejor modelo:

\[\begin{array}{ll} \text{Modelo grado 1} & cargahid \sim tratamiento * profundidad\\ \text{Modelo grado 2:} & cargahid \sim tratamiento * (profundidad + profundidad^2)\\ \text{Modelo grado 3:} & cargahid \sim tratamiento * (profundidad + profundidad^2 + profundidad^3)\\ \end{array}\]

# M1: grado 1
M1 <- lm(cargahid ~ tratamiento * profundidad, data = ejemplo13)
# M2: grado 2
M2 <- lm(cargahid ~ tratamiento * (profundidad + I(profundidad^2)), data = ejemplo13)
# M3 : modelo aditivo
M3 <- lm(cargahid ~ tratamiento * (profundidad + I(profundidad^2) + I(profundidad^3)), data = ejemplo13)
# Resultado
AIC(M1,M2,M3)

El mejor modelo es el de grado 3 (menos valor de AIC). Estudiamos ahora si podemos eliminar la interacción y dejar un modelo aditivo:

# Modelo sin interacción
M3.1 <- lm(cargahid ~ tratamiento + profundidad + I(profundidad^2) + I(profundidad^3), data = ejemplo13)
# Modelo con interacción
M3 <- lm(cargahid ~ tratamiento * (profundidad + I(profundidad^2) + I(profundidad^3)), data = ejemplo13)
# Resultado
AIC(M3.1,M3)

No podemos eliminar la interacción como ya habíamos podido observar en la representación gráfica del modelo. Veamos las curvas obtenidas:

tidy(M3)

¿Cuáles son las ecuaciones estimadas (una por tratamiento)?

Gráficamente:

# grid de valores para construir los modelos
grid <- ejemplo13 %>% data_grid(tratamiento,profundidad) %>% gather_predictions(M3)
# Gráfico
ggplot(ejemplo13,aes(profundidad,cargahid,colour = tratamiento)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  labs(x = "Profundidad", y = "Carga hidráulica") + 
  theme_bw()

¿Se podrían considerar modelos polinómicos de grado mayor? ¿Como procederíamos en ese caso?

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 esta presente en el modelo). Las hipótesis son linealidad, normalidad y varianza constante. Para verificar la hipótesis de linealidad utilizamos el gráfico de residuos vs ajustados y residuos vs predictora numérica, mientras que usamos los tests de normalidad y homogeneidad para el resto de hipótesis.

Ejemplo 11.

Utilizamos el modelo obtenido en el punto anterior.

ajuste <- lm(vida ~ tipo + velocidad + tipo:velocidad,data = ejemplo11)
fitinf <- fortify(ajuste)

Gráfico de residuos

Residuos versus ajustados

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(x = .fitted, y =.stdresid)) + 
  geom_point() +
  facet_wrap(~ tipo, scales="free_x") +  
  theme_bw()

No se observan problemas con los residuos, aunque si se puede ver que en el gráfico para el tipo A que la variabilidad de los residuos aumenta cuando aumenta el valor ajustado, indicando posibles problemas con la homogeneidad de varianzas.

Residuos versus predictora numérica

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(x = velocidad, y =.stdresid)) + 
  geom_point() +
  facet_wrap(~ tipo, scales="free_x") +  
  theme_bw()

No se observa ningún tipo de tendencia que pueda indicar falta de linealidad.

Tests para verificar las hipótesis

# Varianza constante
bartlett.test(.stdresid ~ tipo, data = fitinf)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  .stdresid by tipo
## Bartlett's K-squared = 0.7861, df = 1, p-value = 0.3753
# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$tipo == "A"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$tipo == "A"]
## W = 0.92232, p-value = 0.3767
shapiro.test(fitinf$.stdresid[fitinf$tipo == "B"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$tipo == "B"]
## W = 0.94236, p-value = 0.5796

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

Ejemplo 12

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(longevidad ~ actividad + thorax, data = ejemplo12)
fitinf <- fortify(ajuste)

Gráfico de residuos

Residuos versus ajustados

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(x = .fitted, y =.stdresid)) + 
  geom_point() +
  facet_wrap(~ actividad, scales="free_x") +  
  theme_bw()

Para alguna combinación se observa como los residuos aumentan o disminuyen en función del valor ajustado indicando posibles problemas con la homogeneidad de varianzas.

Residuos versus predictora numérica

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(x = thorax, y =.stdresid)) + 
  geom_point() +
  facet_wrap(~ actividad, scales="free_x") +  
  theme_bw()

No se observa ningún tipo de tendencia que pueda indicar falta de linealidad, aunque que si el aumento de os residuos en función de la predictora, indicando falta de homogeneidad.

Tests para verificar las hipótesis

# Varianza constante
bartlett.test(.stdresid ~ actividad, data = fitinf)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  .stdresid by actividad
## Bartlett's K-squared = 10.551, df = 4, p-value = 0.0321
# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G1"]
## W = 0.98698, p-value = 0.9836
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G2"]
## W = 0.91164, p-value = 0.03317
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G3"]
## W = 0.97543, p-value = 0.7824
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G4"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G4"]
## W = 0.95104, p-value = 0.2646
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G5"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G5"]
## W = 0.96176, p-value = 0.4506

Se puede ver que no se verifica la hipótesis de homogeneidad de varianzas y tampoco la de normalidad para algún grupo. Para tratar de corregir este defecto utilizamos las transformaciones de boxcox sobre la respuesta.

# Transformación
boxcox(ajuste)

Se puede ver que el intervalo de confianza para la transformación contiene al valor de 0.5. Esto nos da indicios de que podríamos utilizar la transformación raíz cuadrada 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
ejemplo12 <- ejemplo12 %>% mutate(rlonge = sqrt(longevidad))
# Ajustamos los modelos con interacción y sin interacción para elegir
M0 <- lm(rlonge ~ actividad + thorax, data = ejemplo12)
M1 <- lm(rlonge ~ actividad + thorax + actividad:thorax, data = ejemplo12)
AIC(M0,M1)

Nos quedamos con el modelo sin interacción y valoramos las hipótesis del modelo.

ajuste <- lm(rlonge ~ actividad + thorax, data = ejemplo12)
fitinf <- fortify(ajuste)
# Verificamos las hipótesis
# Añadimos un nuevo factor a los resultados con todas las combinaciones de los grupos
# Varianza constante
bartlett.test(.stdresid ~ actividad, data = fitinf)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  .stdresid by actividad
## Bartlett's K-squared = 5.061, df = 4, p-value = 0.2811
# Normalidad 
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G1"]
## W = 0.97533, p-value = 0.7969
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G2"]
## W = 0.93471, p-value = 0.1117
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G3"]
## W = 0.97797, p-value = 0.8422
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G4"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G4"]
## W = 0.93645, p-value = 0.1225
shapiro.test(fitinf$.stdresid[fitinf$actividad == "G5"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$actividad == "G5"]
## W = 0.94806, p-value = 0.2266

Podemos ver como se verifican las hipótesis del modelo, de forma que el modelo resultante viene dado por: \[\sqrt{longevidad} \sim actividad + thorax\]

El gráfico del modelo ajustado viene dado por:

# Ajuste
ajuste <- lm(rlonge ~ actividad + thorax, data = ejemplo12)
# grid de valores para construir los modelos
grid <- ejemplo12 %>% data_grid(actividad,thorax) %>% gather_predictions(ajuste)
# Gráfico
ggplot(ejemplo12,aes(thorax,rlonge,colour = actividad)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) +
  labs(x = "Longitud del tórax", y = "Raíz cuadrada de Longevidad") + 
  theme_bw()

Ejemplo 13

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(cargahid ~ tratamiento * (profundidad + I(profundidad^2) + I(profundidad^3)), data = ejemplo13)
fitinf <- fortify(ajuste)

Gráfico de residuos

Residuos versus ajustados

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(x = .fitted, y =.stdresid)) + 
  geom_point() +
  facet_wrap(~ tratamiento, scales="free_x") +  
  theme_bw()

Para alguna combinación se observa como los residuos aumentan o disminuyen en función del valor ajustado indicando posibles problemas con la homogeneidad de varianzas.

Residuos versus predictora numérica

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(x = profundidad, y =.stdresid)) + 
  geom_point() +
  facet_wrap(~ tratamiento, scales="free_x") +  
  theme_bw()

Todavía se observan ciertas curvaturas en el comportamiento de los residuos o que puede considerarse como indicativo de que necesitamos aumentar el grado del polinomio.

Tests para verificar las hipótesis

# Varianza constante
bartlett.test(.stdresid ~ tratamiento, data = fitinf)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  .stdresid by tratamiento
## Bartlett's K-squared = 26.19, df = 2, p-value = 2.056e-06
# Normalidad (Obtenemos los pvalores del test de normalidad)
shapiro.test(fitinf$.stdresid[fitinf$tratamiento == "A"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$tratamiento == "A"]
## W = 0.94196, p-value = 0.4077
shapiro.test(fitinf$.stdresid[fitinf$tratamiento == "B"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$tratamiento == "B"]
## W = 0.96646, p-value = 0.8026
shapiro.test(fitinf$.stdresid[fitinf$tratamiento == "C"])
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid[fitinf$tratamiento == "C"]
## W = 0.97885, p-value = 0.961

No se verifica la igualdad de varianzas. Dada la complejidad del modelo no probamos la transformación de boxcox. En temas posteriores veremos como analizar este tipo de datos sin necesidad de verificar las hipótesis de los modelos lineales.

Predicción

El proceso de predicción en este tipo de modelos es muy simple a partir de las ecuaciones de los modelos obtenidas. De hecho en los puntos anteriores ya hemos visto gráficamente la predicción para todos estos modelos. Si queremos 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 intervalo de predicción.

Vamos a construir la predicción y las bandas de confianza de predicción para cada uno de los ejemplos.

Ejemplo 11.

Utilizamos el modelo obtenido tras la fase de diagnóstico

ajuste <- lm(vida ~ tipo + velocidad + tipo:velocidad,data = ejemplo11)
# Creamos grid de predicción
# Utlizamos la opción each para indicar el número de repeticiones (tantas como niveles del factor)
newdata <- data.frame(velocidad = rep(seq(min(ejemplo11$velocidad), max(ejemplo11$velocidad), .01),each=2), 
                      tipo = factor(c("A", "B")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence"))
# Gráfico de la predicción
ggplot(newdata, aes(x = velocidad, y = fit, color = tipo)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = tipo), alpha = 1/5) +
  labs(x = "Velocidad", y = "Tiempo de vida", title = "Bandas de predicción") +
  theme_bw()

¿Cómo interpretamos las bandas de predicción obtenidas?

Ejemplo 12.

Utilizamos el modelo obtenido tras la fase de diagnóstico

ajuste <- lm(rlonge ~ actividad + thorax, data = ejemplo12)
# Creamos grid de predicción
newdata <- data.frame(thorax = rep(seq(min(ejemplo12$thorax), max(ejemplo12$thorax), .01),each=5), 
                      actividad = factor(c("G1", "G2", "G3", "G4", "G5")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence"))
# Gráfico de la predicción
ggplot(newdata, aes(x = thorax, y = fit, color = actividad)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = actividad), alpha = 1/5) +
  labs(x = "Longitud del thorax", y = "Raíz cuadrada de la longevidad", title = "Bandas de predicción") +
  theme_bw()

Obtenemos ahora el gráfico de predicción en la escala original. Deshacemos el cambio de ráiz cuadrada y representamos de nuevo las bandas de confianza

# Eliminamos la raíz cuadrado de las predicciones 
newdata$fit <- newdata$fit^2
newdata$lwr <- newdata$lwr^2
newdata$upr <- newdata$upr^2
# Gráfico de la predicción
ggplot(newdata, aes(x = thorax, y = fit, color = actividad)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = actividad), alpha = 1/5) +
  labs(x = "Longitud del thorax", y = "Longevidad", title = "Bandas de predicción") +
  theme_bw()

¿Cómo interpretamos las bandas de predicción obtenidas?

Ejemplo 13

Utilizamos el modelo obtenido tras la fase de diagnóstico

ajuste <- lm(cargahid ~ tratamiento * (profundidad + I(profundidad^2) + I(profundidad^3)), data = ejemplo13)
# Creamos grid de predicción
newdata <- data.frame(profundidad = rep(seq(min(ejemplo13$profundidad), max(ejemplo13$profundidad), .01),each=3), 
                      tratamiento = factor(c("A", "B", "C")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence"))
# Gráfico de la predicción
ggplot(newdata, aes(x = profundidad, y = fit, color = tratamiento)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = tratamiento), alpha = 1/5) +
  labs(x = "Profundidad", y = "Carga hidráulica", title = "Bandas de predicción") +
  theme_bw()

¿Cómo interpretamos las bandas de predicción obtenidas?

Bibliografía

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