Modelos Estadísticos. Grado Biotecnología

Modelo de Regresión Lineal Simple

En este tema se presenta el primer tipo de modelos con los que vamos a trabajar. Para poder ir introduciendo toda la terminología, así como las funciones para la construcción y su análisis nos centraremos en los modelos más simples: Modelos de Regresión Lineal Simple. Este tipo de modelos tratan de explicar el comportamiento de una variable respuesta de tipo numérico a partir de una única variable predictora de tipo numérico, a través de una función \(f\) se expresa de tipo lineal: \[Y = \beta_0 + \beta_1 X + \epsilon\]

donde \(\beta_0\) y \(\beta_1\) son los parámetros desconocidos del modelo y \(\epsilon\) es el error aleatorio con distribución \(F\). 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\).

Los parámetros desconocidos de este modelo son \((\beta_0, \beta_1, \sigma^2)\) donde:

  • \(\beta_0\) se conoce como interceptación y representa el valor de la respuesta cuando la variable predictora toma el valor cero, interpretándose como un efecto común en la relación entre la predictora y la respuesta.
  • \(\beta_1\) es la pendiente de la recta y representa el aumento o disminución del valor de la respuesta cuando aumentamos en una unidad el valor de la predictora. En este tipo de modelos dicho parámetro se conoce también como el efecto de la predictora sobre la respuesta.
  • \(\sigma^2\) es la varianza residual del modelo, es decir, la variabilidad asociada con los errores aleatorios. Cuanto menor sea está variabilidad mejor será el modelo obtenido. Hemos de tener en cuenta que la variabilidad de la respuesta (VT) es la suma de la variabilidad aportada por la parte sistemática del modelo (VM) y la variabilidad de los errores aleatorios (VE), es decir: \[VT = VM + VE,\] de forma que si \(VE\) es muy pequeña en comparación con \(VM\) eso implicará directamente que la variable predictora es capaz de explicar el comportamiento de la respuesta a través del modelo propuesto.

Dada un muestra de \(n\) sujetos de la variable respuesta \((y_1,...,y_n)\) y de la variable predictora \((x_{11},...,x_{n1})\), el modelo de regresión lineal simple se puede escribir como:

\[Y = \left(\begin{array}{c} y_1 \\ y_2 \\ ...\\ y_n\\ \end{array} \right) = \left(\begin{array}{cc} 1 & x_{11}\\ 1 & x_{21}\\ ...& ... \\ 1 & x_{n1} \\ \end{array} \right) \left(\begin{array}{c} \beta_0 \\ \beta_1 \\ \end{array} \right) + \left(\begin{array}{c} e_1 \\ e_2 \\ ...\\ e_n \end{array} \right) = X \beta + \epsilon\]

donde \(X\) se denomina matriz del diseño, representando el efecto común (columna de 1’s) y el efecto de la predictora (columna con los valores de la variable), y los \(e_i\) representan los errores aleatorias para cada uno de los sujetos de la muestra.

Ejemplos

Antes de comenzar con el desarrollo del modelo se presentan los ejemplos con los que iremos trabajando. Los dos primeros ya los hemos visto en el tema anterior, pero los recuperamos aquí para realizar todo el análisis:

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

Ejemplo 1. Treinta aleaciones del tipo 90/10 Cu-Ni, cada una con un contenido específico de hierro son estudiadas bajo un proceso de corrosión. Tras un período de 60 días se obtiene la pérdida de peso (en miligramos al cuadrado por decímetro y día) de cada una de las aleaciones debido al proceso de corrosión. Los resultados obtenidos aparecen a continuación:

# carga de datos
hierro <- c(0.01, 0.48, 0.71, 0.95, 1.19, 0.01, 0.48, 1.44, 0.71, 1.96, 0.01, 1.44, 1.96)
peso <- c(127.6, 124, 110.8, 103.9, 101.5, 130.1, 122, 92.3, 113.1, 83.7, 128, 91.4, 86.2)
ejemplo01 <- data.frame(hierro,peso)

Como ya pudimos en el tema anterior parece existir una relación de tipo lineal entre el contenido en hierro y la pérdida de peso. Planteamos el modelo:

\[ peso = \beta_0 + \beta_1 hierro + \epsilon\]

Ejemplo 2. Se propone a una empresa que fabrica vasos de cristal un nuevo proceso de control de calidad. Hasta ahora la empresa seleccionaba una caja de vasos al final de la fabricación y observaba si había alguno roto. Esto provocaba un gran gasto ya que en caso de encontrar algún defecto la caja se desembala y los vasos vuelven a la cadena de embalaje. Ahora se propone seleccionar vasos antes de embalar y determinar así el porcentaje de defectos. Se desea saber si ambos porcentajes están relacionados. Los datos aparecen a continuación:

# carga de datos
cajas <- c(3.0, 3.1, 3.0, 3.6, 3.8, 2.7, 3.1, 2.7, 2.7, 3.3, 3.2, 2.1, 3.0, 2.6)
vasos <- c(3.1, 3.9, 3.4, 4.0, 3.6, 3.6, 3.1, 3.6, 2.9, 3.6, 4.1, 2.6, 3.1, 2.8)
ejemplo02 <- data.frame(cajas,vasos)

Como ya pudimos en el tema anterior parece existir una relación de tipo lineal entre el porcentaje de vasos defectuosos y el número de cajas con vasos defectuosos. Planteamos el modelo:

\[ cajas = \beta_0 + \beta_1 vasos + \epsilon\]

Ejemplo 3. Queremos estudiar la relación existente entre la concentración de madera contenida en la pulpa a partir de la que se elabora papel, y la resistencia (en términos de tensión que soporta) del papel resultante. El objetivo del análisis es describir la tendencia observada.

madera <- c(1.0, 1.5, 2.0, 3.0, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0)
resistencia <- c(6.3, 11.1, 20.0, 24.0, 26.1, 30.0, 33.8, 34.0, 38.1, 39.9, 42.0, 46.1, 53.1, 52.0, 52.5, 48, 42.8, 27.8, 21.9) 
ejemplo03 <- data.frame(madera,resistencia)

El objetivo es tratar de explicar la resistencia en función de la concentración de madera. En primer lugar representamos los datos obtenidos:

# gráfico
ggplot(ejemplo03,aes(madera,resistencia)) + 
  geom_point() + 
  labs(x = "Concentración de madera", y = "Resistencia del papel") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) # Por el momento añadimos la opción se = FALSE

Aunque se ha añadido el posible modelo lineal, se puede observar que los datos no presentan ese comportamiento sino más bien cierto grado de curvatura, que nos debe hacer sospechar que le modelo lineal no resultará válido. A pesar de ello planteamos el modelo lineal para mostrar como funciona la fase de diagnóstico del modelo cuando este es inválido. El modelo propuesto es:

\[ resistencia = \beta_0 + \beta_1 madera + \epsilon\]

En ocasiones disponemos de modelos teóricos que surgen en física, química, biología o medicina que se pueden transformar en un modelo lineal de forma muy sencilla. A continuación se presentan dos ejemplos de este tipo. Más adelante veremos como construir un modelo lineal a partir de las ecuaciones teóricas e interpretaremos los coeficientes del modelo.

Ejemplo 4. En 1929 Edwin Hubble investiga la relación entre la distancia (en megaparsecs. 1 parsec=3.26 años luz) de una galaxia a la tierra y la velocidad (en Km/sg) con la cual parece retroceder. Observo que las galaxias aparecen alejarse de nosotros no importa cual sea la dirección en que miremos. Se piensa de hecho que esto es el resultado del ”Big Bang”. Hubble esperaba aportar conocimiento sobre la formación del universo y sobre cual podría ser su evolución en el futuro. Los datos que recogió incluyen las distancias de la tierra a 24 galaxias, así como sus correspondientes velocidades de retroceso. El objetivo principal del estudio que Hubble llevo a cabo era determinar la edad del universo, y para ello establece lo que se conoce hasta ahora como la ley de Hubble: \[V = H_0*D\] donde \(V\) es la velocidad, \(D\) es la distancia, y \(H_0\) es la constante de Hubble.

distancia <- c(.032,.034,.214,.263,.275,.275,.45,.5,.5,.63,.8,.9,.9,.9,.9,1.0,1.1,1.1,1.4,1.7,2.0,2.0,2.0,2.0)
velocidad <- c(170,290,-130,-70,-185,-220,200,290,270,200,300,-30,650,150,500,920,450,500,500,960,500,850,800,1090) 
ejemplo04 <- data.frame(distancia,velocidad)

El modelo teórico se parece aun modelo lineal donde \(\beta_0 = 0\), y \(H_0\) corresponderá con el valor estimado de \(\beta_1\). Veamos el gráfico para estudiar la posibilidad de plantear un modelo de regresión lineal simple a estos datos.

# gráfico
ggplot(ejemplo04,aes(distancia,velocidad)) + 
  geom_point() + 
  labs(x = "Distancia", y = "Velocidad") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) # Por el momento añadimos la opción se = FALSE

Si que parece observarse un modelo de tipo lineal.

Ejemplo 5. Un grupo de ingenieros agrónomos esta estudiando la evolución de un proceso infeccioso sobre un conjunto de plantas. En concreto durante un período de días se calcula de forma aproximada la fracción de infección de dichas plantas. Por este motivo, dos días consecutivos pueden tener fracciones de infección inferiores. Para evitar esto se proponen dos modelos teóricos: \[(I) \text{ } y = \beta_0 x^{\beta_1}\] \[(II) \text{ } y = \beta_0 + \beta_1 log(x)\]

donde y es el grado de infección y x el número de días transcurridos

dia <- 1:200
infecc <- c(0.02063762, 0.05637206, 0.11889346, 0.11376972, 0.13772089, 0.19055130, 0.18509772, 0.16887353, 0.19753055, 0.21835949, 0.26677238, 0.26360653, 0.27772497, 0.28447172, 0.28300082, 0.34108177, 0.32594214, 0.28675461, 0.34971616, 0.33537176, 0.33217888, 0.35748261, 0.34925483, 0.36278067, 0.37211460, 0.35783240, 0.41498583, 0.40769174, 0.38800213, 0.44174297, 0.43087261, 0.42190606, 0.45097338, 0.45570700, 0.45946960, 0.46153400, 0.46340109, 0.45549254, 0.45487365, 0.45750686, 0.45521341, 0.46881463, 0.45141048, 0.52372846, 0.50803021, 0.46482596, 0.48254773, 0.48449405, 0.51255671, 0.49833262, 0.50802495, 0.50526564, 0.50777983, 0.53873568, 0.50950327, 0.54693458, 0.48815063, 0.53327501, 0.52645577, 0.53063462, 0.53618898, 0.52077545, 0.52633078, 0.51474555, 0.51575426, 0.54528626, 0.55015967, 0.54419107, 0.56346905, 0.58787669, 0.53886562, 0.50427534, 0.57230842, 0.53970820, 0.54179538, 0.57769618, 0.55308538, 0.53593047, 0.56550374, 0.56060245, 0.56496884, 0.57400395, 0.56030226, 0.58199322, 0.56606007, 0.57844415, 0.59505931, 0.58311616, 0.56916054, 0.59989923, 0.59801493, 0.59031303, 0.58529898, 0.56912505, 0.61003513, 0.57193641, 0.62878888, 0.61677661, 0.58247460, 0.56770688, 0.57505675, 0.59541545, 0.58634056, 0.58530427, 0.57418797, 0.59326985, 0.57940758, 0.56266765, 0.58932866, 0.61620602, 0.58719856, 0.61173102, 0.56806743, 0.60015458, 0.61248238, 0.60893367, 0.60582869, 0.59169408, 0.58829584, 0.58557803, 0.60917339, 0.58862023, 0.59849746, 0.60391548, 0.64663334, 0.59742612, 0.61587231, 0.61341390, 0.59329848, 0.61178139, 0.64276168, 0.62355522, 0.61599580, 0.60735889, 0.57537341, 0.63968664, 0.58846078, 0.63307852, 0.65706008, 0.59059116, 0.63408846, 0.61538542, 0.58975607, 0.59146830, 0.59028687, 0.61224876, 0.59417456, 0.63770437, 0.66647829, 0.59925939, 0.64127259, 0.64141050, 0.63317968, 0.60686830, 0.62514131, 0.62241142, 0.63976259, 0.62153212, 0.64899315, 0.62242964, 0.65143794, 0.60985758, 0.60609047, 0.69656194, 0.62384676, 0.63858650, 0.64578674, 0.62380855, 0.64424572, 0.64170765, 0.63043627, 0.63646096, 0.63488074, 0.67853395, 0.62153691, 0.61483840, 0.63790480, 0.64374543, 0.64664922, 0.62913057, 0.61740673, 0.66430865, 0.63241999, 0.62246721, 0.63541282, 0.63655235, 0.66304830, 0.64289529, 0.65662894, 0.63190605, 0.64652159, 0.63607656, 0.64479640, 0.62532881, 0.61734833, 0.68383389, 0.65622608, 0.61950582, 0.63262438, 0.62145169) 
ejemplo05 <- data.frame(dia,infecc)

Los modelos teóricos anteriores se pueden transformar a modelos de regresión lineal simple sin muchos problemas. Para el modelo (I) consideramos la transformación logarítmica:

\[log(y) = log(\beta_0) + \beta_1 log(x) \longrightarrow y' = \beta^`_0 + \beta_1 x`\]

donde el \('\) indica la nueva variable o parámetro. Podemos estudiar el comportamiento lineal obteniendo dichas variables, y volver al modelo original utilizando las transformaciones.

Para el modelo (II) ya viene expresado como un modelo lineal donde únicamente tenemos que obtener las variable predictora transformada mediante la función logaritmo.

Calculamos las nuevas variables:

ejemplo05 <- ejemplo05 %>% mutate(ynew = log(infecc), xnew = log(dia))

Para cada situación representamos el posible modelo lineal con las variables transformadas:

# gráfico
ggplot(ejemplo05,aes(xnew,ynew)) + 
  geom_point() + 
  labs(x = "log(dia)", y = "log(infección)") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) # Por el momento añadimos la opción se = FALSE

Se puede ver que la tendencia lineal no es apropiada para estos datos, ya que se observa cierta curvatura. Más adelante veremos que podríamos plantear un modelo de tipo no lineal para este tipo de tendencia. La conclusión puede ser que el modelo teórico no se corresponde con la información muestral recogida.

# gráfico
ggplot(ejemplo05,aes(xnew,infecc)) + 
  geom_point() + 
  labs(x = "log(dia)", y = "Infección") + 
  theme_bw() + 
  geom_smooth(method = "lm", se = FALSE) # Por el momento añadimos la opción se = FALSE

En este caso si se aprecia la tendencia lineal. Por tanto, este modelo teórico si parece adecuado para representar la información muestral recogida, relacionando los días transcurridos (en escala logarítmica) con la fracción de infección.

Una vez establecidos los posibles modelos y las correspondientes hipótesis, pasamos a la fase de análisis inferencial de los parámetros del modelo o fase de estimación del modelo.

Estimación del modelo

En esta fase procedemos con la estimación del modelo, es decir, obtener estimaciones puntuales,intervalos de confianza, y contrastes de hipótesis sobre cada uno de los parámetros del modelo. Todo el desarrollo teórico se puede consultar en los libros recomendados en la bibliografía. Aquí nos ocupamos únicamente de como obtener y analizar dicho proceso de estimación.

Para el establecimiento de un modelo utilizamos la nomenclatura habitual para indicar la relación entre respuesta y predictora: \[y \sim x\]

Veamos la expresión para cada uno de los ejemplos propuestos.

\[\left\{ \begin{array}{lc} \text{Ejemplo 1} & peso \sim hierro\\ \text{Ejemplo 2} & cajas \sim vasos\\ \text{Ejemplo 3} & resistencia \ sim madera\\ \text{Ejemplo 4} & velocidad \sim distancia - 1\\ \text{Ejemplo 5. (I)} & log(infecc) \sim log(dia)\\ \text{Ejemplo 5}. (II) & infecc \sim log(dia)\\ \end{array} \right. \]

El \(-1\) que aparece en la expresión del modelo 4 es para indicar que dicho modelo se debe obtener teniendo en cuenta que la interceptación debe ser igual 0, es decir, \(\beta_0 = 0\).

Dada una muestra de la variable predictora y respuesta, debemos obtener las estimaciones \(\widehat{\beta}_0\) y \(\widehat{\beta}_1\) de cada uno de los coeficientes del modelo. Una vez obtenidos tendremos la denominada ecuación del modelo ajustada que nos permite obtener el valor ajustado de la respuesta \(\widehat{y}_i\) para cada uno de los sujetos de la muestra: \[ \widehat{y}_i = \widehat{\beta}_0 + \widehat{\beta}_1 x_{i1}, \mbox{ } i=1,...,n\] A partir de esos valores ajustados será posible obtener los residuos del modelo sin más que calcular: \[r_i = y_i - \widehat{y}_i, \mbox{ } i=1,...,n\]

Para el ajuste de un modelo de tipo lineal utilizamos la funciónlm(), mientras que las funciones tidy(modelo), confint(modelo) y fortify(modelo nos permiten realizar todo el análisis inferencial sobre los coeficientes del modelo. Más concretamente:

  • tidy(modelo) nos proporciona la estimación puntual del parámetro junto con su error estándar. Las últimas dos columnas son utilizadas para resolver el contraste: \[\left\{ \begin{array}{lc} H_0: & \beta_i = 0\\ H_1: & \beta_i \neq 0\\ \end{array} \right. \] para cada coeficientes el modelo. Nos fijaremos en el pvalor asociado para concluir si el coeficiente analizado puede ser considerado distinto de cero. En el caso de \(\beta_1\) esto implicará directamente que la variable respuesta y la predictora están relacionadas.

  • confint(modelo) nos proporciona los intervalos de confianza para cada uno de los coeficientes del modelo.

  • fortify(modelo) nos proporciona diferentes columnas de datos entre las que por el momento utilizaremos la columna de valores ajustados (.fitted), los residuos (.resid), y los residuos estandarizados (.stdresid) obtenidos al dividir cada residuo por la desviación típica de todos ellos.

Análisis del ejemplo 1

Para los datos del ejemplo 1 la estimación del modelo propuesto nos proporciona los resultados siguientes:

# Ajuste del modelo
fit <- lm(peso ~ hierro,data = ejemplo01)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(fit)
fitcoef
# Intervalo de confianza
confint(fit)
##                 2.5 %    97.5 %
## (Intercept) 126.69920 132.87400
## hierro      -26.83664 -21.20314
# Valores ajustados y residuos
fitinf <- fortify(fit)
head(fitinf)

El modelo ajustado viene dado por: \[\widehat{peso} = 129.79 - 24.02* hierro\] donde tenemos que la pérdida de peso es igual a 129.79 cuando el contenido de hierro es igual a cero, mientras que un aumento de una unidad en el contenido de hierro implica una reducción de peso de 24.02. Fijando una significatividad de 0.05 podemos concluir que ambos coeficientes son distintos de cero ya que el pvalor obtenido es inferior a dicho valor. Por tanto, tenemos evidencias para establecer que es posible explicar la pérdida de peso por la oxidación en función del contenido de hierro inicial. A continuación representamos gráficamente el modelo obtenido:

# gráfico
ggplot(fitinf,aes(hierro,peso)) + 
  geom_point() + 
  geom_line(aes(hierro,.fitted),colour = "red") +
  theme_bw()

Análisis del ejemplo 2

Para los datos del ejemplo 2 la estimación del modelo propuesto nos proporciona los resultados siguientes:

# Ajuste del modelo
fit <- lm(cajas ~ vasos,data = ejemplo02)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(fit)
fitcoef
# Intervalo de confianza
confint(fit)
##                  2.5 %   97.5 %
## (Intercept) -0.5915737 2.362561
## vasos        0.1898898 1.054966
# Valores ajustados y residuos
fitinf <- fortify(fit)
head(fitinf)

El modelo ajustado viene dado por: \[\widehat{cajas} = 0.89 + 0.62* vasos\]

En vista de los resultados obtenidos (pvalores) se puede concluir que la interceptación puede ser cero. Es un resultado lógico ya que si no tenemos vasos con defectos no deberíamos tener cajas con defectos. Por tanto, resulta necesario proponer un segundo modelo donde consideramos la eliminación de ese parámetro:

# Ajuste del modelo
fit <- lm(cajas ~ vasos - 1,data = ejemplo02)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(fit)
fitcoef
# Intervalo de confianza
confint(fit)
##           2.5 %    97.5 %
## vasos 0.8220822 0.9369524
# Valores ajustados y residuos
fitinf <- fortify(fit)
head(fitinf)

El modelo resultante viene dado por: \[\widehat{cajas} = 0.88* vasos\]

El coeficiente sigue siendo significativo (pvalor inferior a 0.05), con lo que podemos concluir que el porcentaje de cajas con defectos aumenta un 0.88 cada vez que encontramos un nuevo vaso con defectos. En el apartado de bondad de ajuste y comparación de modelos veremos cuál de los modelos es el más adecuado para la información muestral recogida. A continuación representamos gráficamente el último modelo obtenido:

# gráfico
ggplot(fitinf,aes(vasos,cajas)) + 
  geom_point() + 
  geom_line(aes(vasos,.fitted),colour = "red") +
  theme_bw()

Análisis del ejemplo 3

Como ya pudimos ver en el análisis gráfico de estos datos, un modelo lineal no parece el más adecuado. Sin embargo, lo ajustamos para poder realizar un estudio más completo de las hipótesis del modelo en el apartado de diagnóstico del modelo.

# Ajuste del modelo
fit <- lm(resistencia ~ madera,data = ejemplo03)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(fit)
fitcoef
# Intervalo de confianza
confint(fit)
##                 2.5 %    97.5 %
## (Intercept) 9.8645891 32.777936
## madera      0.4042179  3.137753
# Valores ajustados y residuos
fitinf <- fortify(fit)
head(fitinf)

El modelo ajustado viene dado por: \[\widehat{resistencia} = 21.32 + 1.77* madera\] No interpretamos los coeficientes ya que como hemos dicho este modelo no es adecuado. Pasamos a representar gráficamente el modelo obtenido:

# gráfico
ggplot(fitinf,aes(madera,resistencia)) + 
  geom_point() + 
  geom_line(aes(madera,.fitted),colour = "red") +
  theme_bw()

Análisis del ejemplo 4

Para los datos del ejemplo 4 la estimación del modelo propuesto nos proporciona los resultados siguientes:

# Ajuste del modelo
fit <- lm(velocidad ~ distancia - 1,data = ejemplo04)
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(fit)
fitcoef
# Intervalo de confianza
confint(fit)
##              2.5 %   97.5 %
## distancia 336.7348 511.1398
# Valores ajustados y residuos
fitinf <- fortify(fit)
head(fitinf)

El modelo ajustado viene dado por: \[\widehat{velocidad} = 423.94* distancia\]

Dado que el coeficiente puede considerarse distinto de cero (pvalor inferior a significatividad), podemos concluir que la estimación de la constante de Hubble en base a la información muestral recogida es 423.94. Veamos gráficamente el modelo ajustado:

# gráfico
ggplot(fitinf,aes(distancia,velocidad)) + 
  geom_point() + 
  geom_line(aes(distancia,.fitted),colour = "red") +
  theme_bw()

Análisis del ejemplo 5

Se presentan los resultados para el modelo II propuesto inicialmente.

# Ajuste del modelo
fit <- lm(infecc ~ xnew,data = ejemplo05) # xnew = log(dia) 
# Análsis inferencial sobre los coeficientes del modelo
fitcoef <- tidy(fit)
fitcoef
# Intervalo de confianza
confint(fit)
##                   2.5 %      97.5 %
## (Intercept) -0.05723695 -0.02443335
## xnew         0.13016297  0.13758473
# Valores ajustados y residuos
fitinf <- fortify(fit)
head(fitinf)

El modelo ajustado viene dado por: \[\widehat{infecc} = -0.04 + 0.13 * log(dia)\]
¿Cómo interpretamos los coeficientes del modelo? ¿Cuál es la gráfica del modelo ajustado? ¿Qué modelo obtenemos is ajustamos el modelo I a estos datos?

Bondad de ajuste del modelo

La bondad de ajuste trata de cuantificar la capacidad explicativa del modelo obtenido. La capacidad explicativa nos indica lo buena que es la variable predictora para explicar el comportamiento de la respuesta. Sin embargo, hay que tener cuidado con la interpretación de estos procedimientos porque en ocasiones nos pueden llevar a conclusiones erróneas, lo que provoca el establecimiento de modelos que son incorrectos desde el punto de vista de lo que se está midiendo en el diseño experimental. Las dos medidas más habituales para el estudio de la bondad del ajuste son el \(R^2\) ajustado y el test \(F\) de la regresión.

El \(R^2\) cuadrado ajustado es una medida numérica que varía entre 0 y 1, de forma que valores próximos a 1 indican alta capacidad explicativa de la predictora sobre la respuesta. La dificultad estriba en que no tenemos una medida estándar para determinar lo cerca o lejos que nos encontramos del valor de 1, con lo que en muchas ocasiones, resulta en una medida arbitraría sujeta a la interpretación del investigador. Esta medida nunca se presenta de forma individual.

El test \(F\) de la regresión es un procedimiento de contraste que generaliza los contrastes individuales sobre cada uno de los coeficientes del modelo. El contraste considerado viene dado por:

\[\left\{ \begin{array}{lc} H_0: & \text{Todos los coeficientes del modelo son iguales a } 0\\ H_1: & \text{Algún coeficiente del modelo es distinto de } 0\\ \end{array} \right. \]

El test \(F\) se basa en la comparación de la variabilidad explicada por el modelo frente a la variabilidad residual (\(\sigma^2\)). Para el modelo de regresión lineal simple esto implica que si no rechazamos \(H_0\) los coeficientes podrían ser cero, y la capacidad explicativa de la predictora sería cero. La respuesta y la predictora no estarían relacionadas. Si rechazamos \(H_0\) alguno de los coeficientes será cero, pero sin saber cual de ellos exactamente. Para ello deberíamos acudir a los test individuales de cada coeficiente que hemos visto en el punto anterior. Veremos el funcionamiento de este test para más de un predictor en el tema siguiente.

La función glance(modelo) nos proporciona tanto las medidas de bondad de ajuste del modelo como las necesarias para la comparación de modelos que veremos . Los resultados que proporciona esta función son:

  • r.square: Valor de \(R^2\)
  • adj.r.squared: Valor de \(R^2\) ajustado
  • sigma: Varianza residual del modelo
  • statistic: Estimación del estadístico del test \(F\) de regresión, es decir, el cociente entre la variabilidad explicada por el modelo y la residual corregidas por sus respectivos grados de libertas. Cuanto mayor es el valor del estadístico mayor será la capacidad explicativa del modelo construido.
  • p.value: pvalor para el contraste \(F\) de la regresión. Si su valor es inferior a 0.05 el modelo tiene capacidad explicativa, es decir, existe algún coeficiente distinto de cero. Si es superior a 0.05, el modelo establecido carece de capacidad explicativa.
  • df: grados de libertad de la varianza explicada por el modelo. Coincide con el número de parámetros estimados.
  • logLik: logaritmo de la verosimilitud evaluada en los datos muestrales.
  • AIC, BIC, y deviance: Criterios para la valoración de la capacidad explicativa del modelo y la selección del mejor modelo.
  • df.residual: grados de libertad de la varianza residual.

En este punto nos ocupamos únicamente de adj.r.squared y pvalue que son los relacionados con las medidas de bondad del ajuste.

Ejemplo 1

# Ajuste del modelo
fit <- lm(peso ~ hierro,data = ejemplo01)
gof <- glance(fit)
gof

El \(R^2\) ajustado de 0.97 nos indica que el modelo tiene una alta capacidad explicativa, es decir, la variable predictora considerada explica el comportamiento de la respuesta. Habitualmente, cuanto mayor es el \(R^2\) ajustado menores son los residuos, y más cerca están los puntos del diagrama de dispersión del modelo ajustado.

En cuanto al test \(F\) de regresión podemos ver que el pvalor es inferior a 0.05, de forma que, podemos concluir que existe algún coeficiente del modelo que es distinto de cero. Si volvemos sobre el análisis inferencial de los coeficientes podemos ver que ambos eran considerados distintos de cero.

Ejemplo 2

# Ajuste del modelo
fit <- lm(cajas ~ vasos,data = ejemplo02)
gof <- glance(fit)
gof

En este caso el \(R^2\) ajustado toma un valor de 0.40 (poca capacidad explicativa), pero sin embargo el test \(F\) de la regresión resulta significativo (pvalor inferior a 0.05), lo que implica que alguno de los coeficientes del modelo es distinto de cero. Dada la poca capacidad explicativa, aunque el modelo resulte válido los valores ajustados tendrán residuos altos. Por tanto, la componente aleatoria del modelo todavía será muy relevante. Para mejorar la capacidad explicativa se pueden considerar otros tipos de modelos con funciones no lineales o aumentar el número de variables predictoras consideradas.

Ejemplo 3

# Ajuste del modelo
fit <- lm(resistencia ~ madera,data = ejemplo03)
gof <- glance(fit)
gof

De nuevo el valor del \(R^2\) ajustado es pequeño (0.26) pero el test \(F\) resulta significativo. Estos datos resultan confusos, ya que si parecen indicar que el modelo es adecuado cuando pudimos ver con el gráfico asociado que la tendencia de este modelo no era de tipo lineal. Por este motivo es necesario realizar un análisis de los residuos del modelo, que veremos más adelante, ya que de lo contrario podríamos plantear modelos claramente incorrectos.

Ejemplo 4

¿Qué conclusiones podemos extraer al obtener las medidas de bondad de ajuste para el modelo de estos datos?

Ejemplo 5

En este caso vamos a obtener las medidas de bondad de ajuste para los dos posibles modelos considerados. No se trata de un procedimiento de comparación de modelos, pero nos da una primera aproximación sobre ese punto que desarrollaremos en el apartado siguiente.

# Ajuste del modelo
fit1 <- lm(ynew ~ xnew,data = ejemplo05) # ynew = log(infecc); xnew = log(dia) (I)
fit2 <- lm(infecc ~ xnew,data = ejemplo05) # xnew = log(dia) (II)
gof1 <-glance(fit1) # Bondad de ajuste modelo I
gof2 <-glance(fit2) # Bondad de ajuste modelo II
gof <- rbind(gof1,gof2)  # Combinación de ambos resultados
gof

Podemos ver que para ambos modelos el test \(F\) de la regresión resulta significativo, pero sin embargo el modelo II presenta un \(R^2\) ajustado más alto, lo que indica que dicho modelo tiene mayor capacidad explicativa, y es más adecuado para estos datos.

Comparación de modelos

Aunque en los modelos de regresión lineal simple no tiene sentido realizar una comparativa de modelos, ya que sólo disponemos de una variable predictora, mostramos como utilizar las herramientas que proporciona la función glance() a la hora de elegir entre uno y otro. En el tema siguiente generalizamos este apartado a los modelos de regresión lineal múltiple (múltiples predictoras) e introducimos los procedimientos automáticos de comparación de modelos.

Las tres herramientas que proporciona la función son \(logLik\), \(AIC\) y \(BIC\). Si tenemos dos modelos \(M1\) y \(M2\) en competencia (debemos elegir entre dos modelos sobre los mismos datos) los valores de estos estadísticos nos indican que:

  • Si \(logLik\) de \(M1\) es mayor que el de \(M2\) tomamos como mejor modelo el \(M1\).
  • Si \(AIC\) de \(M1\) es menor que el de \(M2\) tomamos como mejor modelo el \(M1\).
  • Si \(BIC\) de \(M1\) es menor que el de \(M2\) tomamos como mejor modelo el \(M1\).

Ejemplo 5

Si recuperamos los resultados del ejemplo 5 podemos comprobar que los tres criterios seleccionan el modelo II frente al modelo I.

Diagnóstico del modelo

El objetivo de las técnicas que estudiamos en este tema es verificar el cumplimiento de las hipótesis básicas del modelo que justificaban los procedimientos de estimación e inferencia planteados en temas anteriores. En concreto: hipótesis de linealidad, incorrelación o independencia, homocedasticidad (variabilidad constante) y normalidad de los errores. Todos estos procedimientos están basados en el estudio del vector de residuos obtenidos tras el ajuste del modelo propuesto. Más concretamente nos centraremos en el estudio de los residuos estandarizados.

Para la verificación de las hipótesis asumidas utilizamos dos tipos de técnicas:

  • Análisis mediante procedimientos gráficos, que nos indicarán a primera vista si hay un incumplimiento de la hipótesis o nos permitirán detectar lo que denominaremos observaciones anómalas.
  • Análisis mediante tests estadísticos, que plantean diferentes contrastes de hipótesis para verificar cada una de las hipótesis del modelo.

Procedimientos gráficos

Consideramos diferentes tipos de gráficos que son utilizados para la comprobación de las diferentes hipótesis:

  • Gráficos de residuos versus valores ajustados, para detectar problemas de no linealidad o varianza no constante. Cuando representamos los residuos frente a los valores ajustados esperamos que dichos valores se distribuyan aleatoriamente sobre todo el gráfico. Si dicho comportamiento no es aleatorio o se aprecia algún tipo de tendencia podemos tener problemas relacionados con la homogeneidad o linealidad. En ocasiones resulta más fácil evaluar este gráfico en términos del valor absoluto de los residuos o de la raíz cuadrada del valor absoluto de los residuos.

  • Gráficos de residuos versus variables explicativas. Para variables predictoras presentes en el modelo podemos detectar problemas de no linealidad o varianza no constante, mientras que para variables predictoras no presentes, la necesidad de que sean incluidas en el modelo. Cuando representamos los residuos frente a variables explicativas en el modelo esperamos que el gráfico no presente ningún tipo de tendencia o comportamiento no aleatorio, ya que esto indicaría falta de homogeneidad y linealidad. Cuando representamos los residuos frente a variables explicativas que no están en el modelo esperamos que el gráfico no presente ningún tipo de tendencia, ya que esto indicaría falta de ajuste y la necesidad de incluir dicha variable en el modelo.

  • Gráficos de normalidad de los residuos, para verificar la hipótesis de normalidad de los errores. Cuando realizamos un gráfico qq de normalidad esperamos que los residuos se ajusten a la recta que aparece en la diagonal. La presencia de otros comportamientos indica la falta de normalidad.

  • Secuencia temporal de los residuos (si ésta es conocida), para identificar autocorrelación, y por lo tanto problemas con la independencia de los errores. Cuando representamos los residuos frente a la secuencia temporal esperamos que no haya ningún tipo de tendencia, ya que esto indicaría falta de independencia, y por tanto que os errores están correlados.

Tests estadísticos

Aunque los métodos gráficos para el análisis de los residuos son muy usados, solo resultan prácticos cuando el incumplimiento de las hipótesis esta más o menos clara. En general, apreciar cierto comportamiento en los residuos resulta una tarea muy subjetiva. Para evitar eso se plantean en la literatura diferentes contrastes para comprobar cada una de las hipótesis.

  • Homogeneidad de varianzas. Hay numerosos tests en la literatura para reconocer heterocedasticidad. Dichos tests se plantean como:

\[\left\{ \begin{array}{lc} H_0: & \text{Varianza constante } \\ H_1: & \text{Varianza no constante} \\ \end{array} \right. \] Si rechazamos la hipótesis nula estaremos concluyendo que nuestro modelo incumple la hipótesis de varianza constante. Los tests utilizados son el de Breusch-Pagan para variables predictoras de tipo numérico, y el de Bartlett para variables predictoras de tipo categórico. Para realizar el test de Breusch-Pagan utilizamos la función bptest(modelo), mientras que para realizar el test de Bartlett utilizamos la función bartlett.test(modelo)

  • Normalidad de los errores. La hipótesis de normalidad de los errores en el modelo lineal justifica la utilización de los procedimientos de inferencia presentados. En muestras pequeñas, la no normalidad de los errores es muy difícil de diagnosticar a través del análisis de los residuos, pues éstos pueden diferir notablemente de los errores aleatorios. La forma habitual de diagnosticar la hipótesis de normalidad es a través del test como el de Shapiro-Wilks. El test utilizado es: \[\left\{ \begin{array}{lc} H_0: & \text{residuos normales } \\ H_1: & \text{residuos no normales} \\ \end{array} \right. \] Rechazar la hipótesis nula implica el rechazo de la hipótesis de normalidad. Para realizar dicho contraste utilizamos la función shapiro.test(residuos).

  • Independencia de los errores. Se trata aquí de detectar la posible dependencia de los errores entre si. La desviación de la incorrelación más estudiada se refiere a la posible influencia de los errores con los que le preceden o suceden en la secuencia con que se han obtenido las observaciones. Es lo que se denomina correlación serial. Un contraste de correlación serial rutinario es el test de Durbin-Watson para contrastar: \[\left\{ \begin{array}{lc} H_0: & \text{residuos no correlados } \\ H_1: & \text{residuos correlados} \\ \end{array} \right. \] Rechazar la hipótesis nula implica el rechazo de la hipótesis de incorrelación. Para realizar dicho contraste utilizamos la función dwtest(modelo).

Incumplimiento de hipótesis

Una vez identificado el incumplimiento de alguna de las hipótesis del modelo, hay que tratar de identificar porque se produce dicho incumplimiento. Se estudia si el incumplimiento es debido a:

  • Subconjunto de los datos que influye desproporcionadamente en el ajuste del modelo propuesto, con lo cual las estimaciones y predicciones dependen mucho de él. En primer lugar, el objetivo es identificar dichas observaciones. Una vez detectadas la forma de proceder es la siguiente:
    • Comprobar si la influencia se debe a un error en la toma de observaciones, si es así se corrigen los defectos encontrados y se comienza de nuevo.
    • Si los datos son correctos y el subconjunto de influyentes es pequeño se opta casi siempre por su eliminación del banco de datos. En otras ocasiones se puede optar por estudiar de forma separada a dichas observaciones.
  • Comportamiento sistemático del modelo. Este caso es más complicado y requiere de procedimientos más sofisticados para corregir los defectos que aparezcan en el modelo.

De la primera parte se encarga de analizarla los diagnósticos de influencia, mientras que en el segundo caso se trata principalmente de realizar transformaciones de las variables involucradas en el modelo.

Análisis de influencia

El análisis de influencia pretende detectar aquellas observaciones cuya inclusión/exclusión en el ajuste altera sustancialmente los resultados. Es interesante siempre, localizar este tipo de datos, si existen, y evaluar su impacto en el modelo. Si estos datos influyentes son “malos” (provienen de errores en la medición, o de condiciones de experimentación diferentes, etc.) habrían de ser excluidos del ajuste; si son “buenos”, contendrán información sobre ciertas características relevantes a considerar en el ajuste.

A primera vista, observaciones que dan lugar a un residuo grande, pueden influir notablemente en el ajuste. Las denominaremos OBSERVACIONES ALEJADAS. Su existencia puede indicar también la inadecuación del modelo asumido a la realidad experimental. Si dicha observación tiene un residuo exageradamente grande la denominamos ANÓMALA (outlier en inglés). Por otra parte, observaciones que adoptan valores extremos de alguna o varias variables explicativas pueden tener más influencia que las usuales. Las denominaremos OBSERVACIONES ATÍPICAS.

Sin embargo, las dos características no siempre suponen que las observaciones que las manifiestan sean también influyentes. Generalmente se dice que una observación es alejada si el valor absoluto del residuo es mayor que 2. Se considera anómala si el valor absoluto del residuo es mayor que 3.

Un criterio para valora la influencia de una observación sobre los coeficientes del modelo es el cálculo de la distancia de CooK. Se consideran como obsrvaciones influyentes todas aquellas cuyo valor de la distancia sea mayor que 1. Dicha distancia se obtiene directamente con la función fortify(modelo) en al columna denomianda .cooksd. Existen otro tipo de medidas de influencia (se pueden obtener con la función influence.measures(ajuste)) pero no la estudiaremos aquí.

Si el incumplimiento de las hipótesis no es debido a la presencia de observaciones influyentes, sino más bien a un comportamiento sistemático del modelo, los remedios para corregir estas deficiencias pasan principalmente por:

  • Propuesta de otros modelos adecuados a la distribución de la respuesta y su relación con los predictores. Este punto o trataremos ampliamente más adelante
  • Transformar la variable respuesta (si es de tipo continuo), o las variables predictoras (si son de tipo continuo).

Transformaciones

Dividimos este punto en tres apartados:

  • Transformaciones de variable respuesta y predictora. Existen situaciones experimentales (ya hemos visto algún ejemplo de ello) donde ya partimos de un tipo de modelo de carácter no lineal pero que se podría convertir en lineal con una sencilla transformación de la respuesta o la predictora:
Modelo Teórico Modelo transformado (lineal)
\(Y = \beta_0 X^{\beta_1}\) \(log(Y) \sim log(X)\)
\(Y = \beta_0 exp^{\beta_1 X}\) \(log(Y) \sim X\)
\(Y = \beta_0 + \beta_1 log(X)\) \(Y \sim log(X)\)
\(Y = \frac{X}{\beta_0 + \beta_1 X}\) \(1/Y \sim 1/X\)
\(Y = \frac{1}{\beta_0 + \beta_1 X}\) \(1/Y \sim X\)
\(Y = \beta_0 + \beta_1 \frac{1}{X}\) \(Y \sim 1/X\)

Este tipo de modelos teóricos suelen aparecer en diseños experimentales sobre concentraciones de compuestos.

  • Transformaciones sobre la predictora. Se utilizan principalmente ante la falta de linealidad, y se basan principalmente en la construcción de modelos de predicción polinómicos (ejemplo 3). Estos modelos los veremos en el tema siguiente.

  • Transformaciones sobre la respuesta. Esta suele ser la opción más habitual. Obtener una transformación adecuada de la respuesta sin alterar las variables predictoras. Se suelen utilizar ante el incumplimiento de las hipótesis de normalidad o varianza constante. Como buscara una transformación adecuada es un tema que puede resulta costoso, vamos a utilizar un procedimiento automático que nos da una transformación rápida. Dicho procedimiento se conoce con el nombre de transformaciones de Box-Cox. Dicho procedimiento consiste en obtener un intervalo de confianza para un parámetro (\(\lambda\)) que refleja la transformación de la respuesta a utilizar. Las transformaciones más habituales son:

\(\lambda\) Transformación
-2 \(1/Y^2\)
-1 \(1/Y\)
-1/2 \(1/\sqrt{Y}\)
0 \(log(Y)\)
1/2 \(\sqrt{Y}\)
1 \(Y\)
2 \(Y^2\)

Una vez realizado el estudio de influencia la forma de proceder consiste en eliminar las observaciones influyentes y obtener un nuevo modelo sin ellas, o bien realizar alguna de las transformaciones planteadas y ajustar el nuevo modelo. Una vez construido deberemos realizar un nuevo diagnóstico para verificar que se cumple las hipótesis en ese nuevo modelo. Se trata pues de un proceso circular donde a cada modificación debemos obtener un nuevo modelo y analizarlo completamente hasta llegar a un modelo que cumpla con todas las especificaciones.

Diagnóstico de los ejemplos

Pasamos a realizar el proceso de diagnóstico de cada uno de los ejemplos.

Ejemplo 1

Ajustamos el modelo y obtenemos los residuos

# Ajuste del modelo y obtención de los residuos
fit <- lm(peso ~ hierro,data = ejemplo01)
fitinf <- fortify(fit)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","hierro",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method= mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

No se aprecia ninguna tendencia (recta horizontal) entre residuos y valores ajustado o variables predictoras. Además, los residuos presentan un comportamiento aleatorio. Ambos gráficos parecen verificar las hipótesis de varianza constante y linealidad.

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

A pesar de que hay pocos puntos para realizar un estudio más preciso, si parece apreciarse que los puntos de distribuyen a lo largo de la línea que marca la hipótesis de normalidad. Podemos asumir por tanto que se verifica la hipótesis de normalidad.

Pasamos ahora a la verificación de las hipótesis a través del los diferentes tests estadísticos asociados con las hipótesis del modelo.

# Varianza constante
bptest(fit)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 0.024539, df = 1, p-value = 0.8755
# Normalidad
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.92905, p-value = 0.3312
# Independencia
dwtest(fit)
## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 2.5348, p-value = 0.8524
## alternative hypothesis: true autocorrelation is greater than 0

Todos los tests resultan no significativos (pvalores mayores que 0.05) y por tanto tenemos evidencias para concluir que el modelo propuesto verifica todas las hipótesis necesarias.

Por último realizamos un análisis de influencia para conocer si hay alguna observación que se pueda considerar anómala:

# Calculamos los valores absolutos
abs(fitinf$.cooksd)
##  [1] 0.0668861786 0.2284452990 0.0193150251 0.0461479881 0.0005440446
##  [6] 0.0054108272 0.0970376280 0.0796340512 0.0006986207 0.0291503547
## [11] 0.0422197749 0.1367777095 0.3610134495
# Valoramos si hay alguna observación con distancia mayor que 1
abs(fitinf$.cooksd) > 1
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE

No hay ninguna observación que se pueda considerar influyente.

Ejemplo 2

# Ajuste del modelo
fit <- lm(cajas ~ vasos,data = ejemplo02)
fitinf <- fortify(fit)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","vasos",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method=mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

No se aprecia ninguna tendencia (recta horizontal) entre residuos y valores ajustado o variables predictoras. Además, los residuos presentan un comportamiento aleatorio. Ambos gráficos parecen verificar las hipótesis de varianza constante y linealidad.

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

A pesar de que hay pocos puntos para realizar un estudio más preciso, si parece apreciarse que los puntos de distribuyen a lo largo de la línea que marca la hipótesis de normalidad. Podemos asumir por tanto que se verifica la hipótesis de normalidad.

Pasamos ahora a la verificación de las hipótesis a través del los diferentes tests estadísticos asociados con las hipótesis del modelo.

# Varianza constante
bptest(fit)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 0.28704, df = 1, p-value = 0.5921
# Normalidad
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.94046, p-value = 0.4242
# Independencia
dwtest(fit)
## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 2.6522, p-value = 0.9099
## alternative hypothesis: true autocorrelation is greater than 0

Todos los tests resultan no significativos (pvalores mayores que 0.05) y por tanto tenemos evidencias para concluir que el modelo propuesto verifica todas las hipótesis necesarias.

Por último realizamos un análisis de influencia para conocer si hay alguna observación que se pueda considerar anómala:

# Calculamos los valores absolutos
abs(fitinf$.cooksd)
##  [1] 1.949301e-02 4.972583e-02 1.159003e-06 7.616037e-02 2.193083e-01
##  [6] 8.776774e-02 4.626570e-02 8.776774e-02 8.992655e-05 1.458694e-02
## [11] 1.184520e-01 4.386325e-01 1.949301e-02 1.100024e-03
# Valoramos si hay alguna observación con distancia mayor que 1
abs(fitinf$.cooksd) > 1
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE FALSE

No hay ninguna observación que se pueda considerar influyente.

Ejemplo 3

# Ajuste del modelo
fit <- lm(resistencia ~ madera,data = ejemplo03)
fitinf <- fortify(fit)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","madera",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method=mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Se aprecia una clara tendencia en los residuos. Se observa una parábola en los puntos o que indica que el modelo es inadecuado. En este tipo de situaciones se ha de incluir un modelo polinómico para corregir la falta de linealidad. EL modelo tendría la forma: \[ resistencia \sim madera + madera^2\]

El estudio de este tipo de modelos lo veremos más adelante. Por el momento la conclusión es que el modelo no es adecuado. No hace falta seguir con el proceso de diagnóstico ya que hemos detectado un problema que es necesario corregir.

Ejemplo 4

# Ajuste del modelo
fit <- lm(velocidad ~ distancia - 1,data = ejemplo04)
fitinf <- fortify(fit)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","distancia",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method=mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Se aprecia cierta tendencia en los datos que puede indicar falta de linealidad. Los residuos tienen un comportamiento aleatorio indicando que se puede verificar la hipótesis de varianza constante.

# Gráfico de normalidad de los residuos
ggplot(fitinf, aes(sample=.stdresid)) + 
  stat_qq() + 
  geom_abline() +
  theme_bw()

Parece apreciarse que los puntos de distribuyen a lo largo de la línea que marca la hipótesis de normalidad. Podemos asumir por tanto que se verifica la hipótesis de normalidad.

Pasamos ahora a la verificación de las hipótesis a través del los diferentes tests estadísticos asociados con las hipótesis del modelo.

# Varianza constante
# bptest(fit) No se puede utilizar este test cuando unicamente tenemos un coeficiente en el modelo 
# Normalidad
shapiro.test(fitinf$.stdresid)
## 
##  Shapiro-Wilk normality test
## 
## data:  fitinf$.stdresid
## W = 0.9773, p-value = 0.8413
# Independencia
dwtest(fit)
## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 2.0564, p-value = 0.5563
## alternative hypothesis: true autocorrelation is greater than 0

Las dos hipótesis quedan verificadas y por tanto podemos considerar como válido el modelo planteado.

Ejemplo 5

Probamos únicamente con el modelo I.

# Ajuste del modelo
fit <- lm(ynew ~ xnew,data = ejemplo05) # xnew = log(dia) 
fitinf <- fortify(fit)

Comenzamos con el análisis gráfico de los residuos

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","xnew",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method=mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

Se aprecia una clara curvatura en lo residuos indicando falta de linealidad. Dado que el modelo ajustado proviene de un modelo teórico no resulta posible ningún tipo de transformación para mejorar el modelo.

Realizamos un análisis de influencia para detectar posibles observaciones anómalas.

# Calculamos los valores absolutos
abs(fitinf$.cooksd)
##   [1] 5.413540e+00 9.264112e-01 3.224064e-02 1.086848e-01 4.296897e-02
##   [6] 1.059154e-03 2.177221e-03 2.605096e-02 4.884486e-03 5.949515e-04
##  [11] 8.207157e-03 3.062897e-03 4.408731e-03 3.657142e-03 1.474976e-03
##  [16] 1.751766e-02 7.808957e-03 2.256592e-06 9.509047e-03 3.830573e-03
##  [21] 2.050037e-03 5.148824e-03 2.483038e-03 3.453717e-03 3.811308e-03
##  [26] 1.233553e-03 8.938824e-03 6.141796e-03 2.423207e-03 9.591843e-03
##  [31] 6.437645e-03 4.216323e-03 7.381116e-03 7.046626e-03 6.613820e-03
##  [36] 6.004130e-03 5.436095e-03 3.867414e-03 3.247022e-03 2.986585e-03
##  [41] 2.357787e-03 3.035960e-03 1.419649e-03 7.449256e-03 5.140106e-03
##  [46] 1.391039e-03 2.189145e-03 2.013778e-03 3.677113e-03 2.310660e-03
##  [51] 2.647186e-03 2.183887e-03 2.073982e-03 3.846127e-03 1.703871e-03
##  [56] 3.778326e-03 5.103572e-04 2.337163e-03 1.754405e-03 1.770307e-03
##  [61] 1.857411e-03 1.010157e-03 1.087976e-03 5.761336e-04 5.114195e-04
##  [66] 1.435616e-03 1.497213e-03 1.117421e-03 1.788156e-03 2.861475e-03
##  [71] 6.225509e-04 8.081367e-06 1.596960e-03 4.114593e-04 3.924296e-04
##  [76] 1.432979e-03 5.242304e-04 1.355437e-04 7.061124e-04 5.000964e-04
##  [81] 5.398188e-04 7.089927e-04 3.134448e-04 7.825175e-04 3.224121e-04
##  [86] 5.398359e-04 9.378389e-04 5.265076e-04 1.966864e-04 8.456630e-04
##  [91] 7.165465e-04 4.572316e-04 3.021669e-04 5.230409e-05 7.757246e-04
##  [96] 3.777669e-05 1.222365e-03 7.615307e-04 6.521862e-05 2.048208e-06
## [101] 1.997789e-06 1.290655e-04 2.658470e-05 1.112869e-05 1.722802e-05
## [106] 2.775837e-05 1.369389e-05 2.126443e-04 3.399887e-07 1.679334e-04
## [111] 1.780475e-05 6.726814e-05 3.142151e-04 2.935191e-07 2.619769e-05
## [116] 3.983527e-06 8.124700e-07 1.020317e-04 1.758922e-04 2.575842e-04
## [121] 1.348879e-05 2.779081e-04 1.465376e-04 1.030050e-04 1.567684e-04
## [126] 2.642523e-04 4.557954e-05 8.644476e-05 4.942610e-04 1.568691e-04
## [131] 1.544151e-05 6.393376e-05 1.819769e-04 3.883620e-04 1.632558e-03
## [136] 1.027300e-05 1.201350e-03 8.399386e-05 1.130574e-05 1.361883e-03
## [141] 1.408270e-04 5.782935e-04 1.679950e-03 1.688715e-03 1.851244e-03
## [146] 9.301354e-04 1.842676e-03 3.020179e-04 2.980209e-06 1.867815e-03
## [151] 3.426100e-04 3.805474e-04 6.529610e-04 1.850860e-03 1.068752e-03
## [156] 1.260584e-03 6.714215e-04 1.466580e-03 5.045655e-04 1.597678e-03
## [161] 5.341060e-04 2.529161e-03 2.905025e-03 1.569756e-07 1.994630e-03
## [166] 1.325853e-03 1.088258e-03 2.307678e-03 1.303475e-03 1.504194e-03
## [171] 2.217072e-03 1.962798e-03 2.154374e-03 4.067611e-04 3.291106e-03
## [176] 3.965659e-03 2.384138e-03 2.127486e-03 2.054354e-03 3.352098e-03
## [181] 4.476106e-03 1.383748e-03 3.483050e-03 4.467336e-03 3.506282e-03
## [186] 3.545697e-03 1.867404e-03 3.304309e-03 2.448348e-03 4.499046e-03
## [191] 3.398315e-03 4.412632e-03 3.795086e-03 5.782256e-03 6.850507e-03
## [196] 1.527175e-03 3.386348e-03 7.160171e-03 5.849983e-03 7.306365e-03
# Valoramos si hay alguna observación con distancia mayor que 1
abs(fitinf$.cooksd) > 1
##   [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [166] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [188] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [199] FALSE FALSE

Hay por lo menos una observación de carácter influyente. La forma de proceder es eliminar dicha observación y ajustar de nuevo el modelo.

# Ajuste del modelo
datossel <- ejemplo05[-1,]
fit <- lm(ynew ~ xnew,data = datossel) # xnew = log(dia) 
fitinf <- fortify(fit)

Realizamos el diagnóstico gráfico:

# Gráfico de residuos versus valores ajustados, y versus variables predictoras
diagnosis <- fitinf[,c(".fitted","xnew",".stdresid")] 
datacomp = melt(diagnosis, id.vars='.stdresid')
ggplot(datacomp) +
  geom_jitter(aes(value,.stdresid, colour=variable),) + 
  geom_smooth(aes(value,.stdresid, colour=variable), method=mgcv::gam, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Residuos") +
  theme_bw()

A pesar de la eliminación de la observación influyente se sigue observando una fuerte tendencia en los residuos. El problema es sistemático y no por las observaciones recogidas. Es un indicativo claro de que el modelo teórico no es adecuado para estos datos. Ya habíamos visto anteriormente que el modelo II era mejor que este desde el punto de vista de la bondad de ajuste. Una alternativa puede ser, como en el ejemplo anterior, considerar un modelo polinómico: \[log(Y) = \beta_0 + \beta_1 log(X) + \beta_2 (log(X))^2\]

Realiza el diagnóstico para el modelo II.

Predicción

Una vez obtenido un modelo definitivo, la última fase de la modelización consiste en la predicción de la respuesta a partir de un nuevo conjunto de valores de la predictora o predictoras. Básicamente se trata de utilizar valores dentro del rango de la variable predictora para conocer el valor de la respuesta sin necesidad de realizar el diseño experimental. Una vez ajustado el modelo si consideramos una observación \(x_0\) dentro del rango de valores de \(X\) la predicción viene dada por: \[y_0 = \widehat{\beta_0} + \widehat{\beta_1} x_0.\]

Esto nos proporciona una estimación puntual del valor predicho, pero sin embargo es necesario proporcionar un intervalo de confianza para dicho valor para tener en cuenta la variabilidad observada en el modelo propuesto. Existen dos posibilidades de predicción en este sentido:

  • Predicción del valor medio de la respuesta. Se trata de predecir el valor medio de la respuesta para un valor especifico de la variable predictora (\(x_0\)). Esta es la herramienta de predicción habitual ya que tiene una menor variabilidad. La idea es que para un mismo valor de \(x_0\) obtendremos diferentes valores predichos de la respuesta, y por tanto, más que interesarnos la predicción de la respuesta, nos centramos en predecir la media de todos esos posibles valores de la respuesta.

  • Predicción del valor de la respuesta. Se trata de predecir el valor de la respuesta para un valor especifico de la variable predictora (\(x_0\)). Dado que estamos intentando predecir un único valor y no la media de un conjunto de valores el intervalo de confianza de predicción es mayor que en el caso anterior. Tenemos más variabilidad cuando queremos predecir un valor que cuando queremos predecir la media de un conjunto de valores.

Utilizaremos la función predict() para construir la predicción para un modelo dado. Veremos su aplicación en los diferentes ejemplos.

Predicción para los ejemplos

Pasamos a realizar el proceso de predicción de cada uno de los ejemplos. Veremos como obtener una predicción específica en las dos situaciones estudiadas, y obtendremos la representación gráfica para la predicción de todo el modelo obtenido.

Ejemplo 1

En primer lugar queremos predecir la pérdida de peso para diferentes valores de contenido en hierro inicial. Para ello construimos una secuencia de valores dentro del rango de valores del contenido de hierro.

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(peso ~ hierro, data = ejemplo01)
# Secuencia de valores de predicción
newdata <- data.frame(hierro = seq(min(ejemplo01$hierro),max(ejemplo01$hierro), length = 10))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = hierro, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo01, aes(x = hierro, y = peso)) +
  labs(x = "Contenido en Hierro", y = "Predicción del valor medio de la pérdida de peso") +
  theme_bw()

Para un valor de contenido de hierro de 1.31 se predice, con una confianza del 95%, una media de pérdida de peso entre 96.08 y 100.56. Además, en el gráfico de predicción, podemos ver que hay valores de la respuesta que no quedan dentro de la región de confianza de predicción.

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(peso ~ hierro, data = ejemplo01)
# Secuencia de valores de predicción
newdata <- data.frame(hierro = seq(min(ejemplo01$hierro),max(ejemplo01$hierro), length = 10))
# Predicción para una observación
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="prediction")) # Opción interval = "prediction"
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = hierro, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo01, aes(x = hierro, y = peso)) +
  labs(x = "Contenido en Hierro", y = "Predicción de la pérdida de peso") +
  theme_bw()

La diferencia entre ambas predicciones es que en el segundo caso la amplitud de los intervalos de confianza es mayor. En el primer gráfico las bandas de confianza tienden a encojerse cuando nos acercamos al punto medio definido por la variable predictora y la respuesta.

Ejemplo 2

En primer lugar queremos predecir el porcentaje de cajas defectuosas para diferentes valores del porcentaje de vasos defectuosos. Para ello construimos una secuencia de valores dentro del rango de valores.

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(cajas ~ vasos,data = ejemplo02)
# Secuencia de valores de predicción
newdata <- data.frame(vasos = seq(min(ejemplo02$vasos),max(ejemplo02$vasos), length = 10))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = vasos, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo02, aes(x = vasos, y = cajas)) +
  labs(x = "% vasos", y = "Predicción del valor medio del % cajas") +
  theme_bw()

# Ajuste del modelo y obtención de los residuos
ajuste <- lm(cajas ~ vasos,data = ejemplo02)
# Secuencia de valores de predicción
newdata <- data.frame(vasos = seq(min(ejemplo02$vasos),max(ejemplo02$vasos), length = 10))
# Predicción para un valor de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="prediction")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = vasos, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo02, aes(x = vasos, y = cajas)) +
  labs(x = "% vasos", y = "Predicción del valor medio del % cajas") +
  theme_bw()

Interpreta los resultados de los dos tipos de predicción obtenidos.

Ejemplo 4

En primer lugar queremos predecir la velocidad de retroceso de las galaxias para diferentes valores de la distancia. Para ello construimos una secuencia de valores dentro del rango de valores,

# Ajuste del modelo y obtención de los residuos
ajuste <-fit <- lm(velocidad ~ distancia - 1,data = ejemplo04)
# Secuencia de valores de predicción
newdata <- data.frame(distancia = seq(min(ejemplo04$distancia),max(ejemplo04$distancia), length = 10))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="confidence")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = distancia, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo04, aes(x = distancia, y = velocidad)) +
  labs(x = "Distancia", y = "Velocidad") +
  theme_bw()

Dado que el modelo no tiene término de interceptación, la predicción para valores bajos tiene rangos pequeños. El problema es que para distancias pequeñas es donde tenemos más datos observados. A pesar de que el modelo es válido (cumple con las hipótesis) su capacidad predictiva es muy reducida para valores pequeños de la distancia.

# Ajuste del modelo y obtención de los residuos
ajuste <-fit <- lm(velocidad ~ distancia - 1,data = ejemplo04)
# Secuencia de valores de predicción
newdata <- data.frame(distancia = seq(min(ejemplo04$distancia),max(ejemplo04$distancia), length = 10))
# Predicción para la media de la respuesta
newdata <- data.frame(newdata, predict(ajuste, newdata, interval="prediction")) # Opción interval = "confidence" 
newdata
# Gráfico resultante (predicción y bandas de predicción)
ggplot(newdata, aes(x = distancia, y = fit)) +
  geom_line() + 
  geom_ribbon(aes(ymax=upr, ymin=lwr), alpha=1/5) +
  geom_point(data = ejemplo04, aes(x = distancia, y = velocidad)) +
  labs(x = "Distancia", y = "Velocidad") +
  theme_bw()

Se puede ver la gran diferencia de este tipo de predicción con la anterior. Cuando deseamos predecir la velocidad para una distancia los intervalos de confianza se amplían de forma excesiva. Por ejemplo, para un valor de distancia de 1.0 la velocidad de mueve aproximadamente entre 0 y 800. En la predicción de la media dicho intervalo va entre 300 y 500 aproximadamente.

Bibliografía

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