En el contexto de regresión lineal simple, se dice que un dato u observación está asociado a un valor extremo (outlier), si el valor \(y_{i}\) es inusual o muy extremo para el valor \(x_{i}\) asociado. Estos valores pueden identificarse fácilmente en el diagrama de dispersión inicial.

Análisis de algunos supuestos

Comencemos por graficar el scatterplot inicial

# Cargamos los datos
Datos = read.csv("ejemplo5.csv")

library(latex2exp)
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$"))

de los cuales puede observarse, para los valores de \(x=6\) y \(x=12\), outliers. Después observamos las gráficas que R trae por defecto para el análisis de los supuestos

# Ajuste del modelo
fit=lm(y~x, data=Datos)

#Gráficas de R para revisar algunos de los supuestos
par(mfrow=c(1,3)) 
par(mar=c(3, 4, 2, 1))
plot(fit, 1)
plot(fit, 3)
plot(fit, 2)

Donde en la gráfica para analizar el supuesto de la linealidad podemos ver cierto tipo de patrón y no la “nube” de puntos de la misma longitud que buscamos. Además pueden observarse dos valores en especifico que salen de la zona del resto; de igual manera, en la gráfica para analizar el supuesto de la homocedasticidad vemos que se marca cierto patrón. Finalmente, vemos que los puntos no caen en la diagonal de la última gráfica (que nos ayuda a analizar el supuesto de la normalidad). Así, concluimos que los datos como están son problemáticos.

Análisis de los Outliers

Recordemos que para analízar el potencial impacto en el ajuste utilizamos los hat values:

\[ h_{i}=\frac{1}{n}+\frac{(x_{i}-\bar{X})^{2}}{SSx} \]

los cuales indican la influencia (o leverage) de cada observación en las estimaciones considerando el valor de \(x_{i}\) asociado. Así, un outlier será potencialemnte influyente si tiene un valor alto \(h_{i}\) asociado, digamos \(h_{i}>2\bar{h}\) o \(h_{i}>3\bar{h}\) (con \(\bar{h}\) el número de parámetros en el modelo dividido entre el número de observaciones).

Estadística de Cook (Cook’s distance)

Esta estadística trata de identificar el impacto de una observación en la estimación de los parámetros.

Para esto, de forma similar a los errores estudentizados, para obtener la estadística \(D_{i}\) asociada a la observación \(i\), se ajusta un modelo sin considerar esa observación y los valores de los parámetros de ese ajuste se comparan con los correspondientes a la regresión que incluye todas las observaciones. Una regla “ad-hoc” (locución latina que significa literalmente “para esto”. Generalmente se refiere a una solución específicamente elaborada para un problema o fin preciso y, por tanto, no generalizable ni utilizable para otros propósitos) que comúnmente se usa para analizar estos valores es la siguiente:

se considera que una observación es influyente si

\[ D_{i}>\frac{4}{n-2} \]

donde el denominador corresponde a los grados de libertad asociados a la estimación de \(\sigma^{2}\).

En R viene implementada por defecto una gráfica asociada a la estadística de Cook

par(mfrow=c(1,1)) 
par(mar=c(4, 5, 2, 1))
plot(fit, 4)

del cual, mientras un valor tenga asociada una distancia de Cook mayor será más problable que dicho valor este cambiando el valor de las estimaciones para los betas. De acuerdo con la gráfica anterior tenemos que el valor en la observación 21 es la que nos esta causando problemas (la que nos pone alerta).

Luego, consideramos los valores

#Valores de referencia que pueden servir de apoyo
#Para la Cook's distance  4/(n-2)
auxiliar <- 4/(dim(Datos)[1]-2)

# Notamos que
(dim(Datos))
## [1] 21  2
# nos da un vector con el número de filas (observaciones) y el número columnas

#Para Leverage  2* hbar   o 3*hbar  con hbar=2/n
dos_h_bar <- 4/(dim(Datos)[1])
tres_h_bar <- 6/(dim(Datos)[1])

Después, sabemos que en la gráfica vista en clase

se incluye una curva para cierto valor de interés de \(D\), donde en los extremos superior e inferior derecho de la gráfica se encontrarán los potenciales valores influyentes. En R tenemos una gráfica por defecto que nos muestra el mismo gráfico que en la imagen anterior

plot(fit,5)

con valores por defecto para la distancia de Cook de 0.5 y 1. Luego, debemos incluir el valor de \(\frac{4}{n-2}\) (almacenado en la variable auxiliar) como un valor a considerar en la gráfica como distancia de Cook. Para ello

# Gráfica con valores para la distancia de Cook de 0.5, 1 y de 4/(n-2)
plot(fit, 5, cook.levels = c(auxiliar, 0.5, 1.0))

# agregamos algunas líneas:
# verticales
abline(v = dos_h_bar, col="blue", lty = 2)
abline(v = tres_h_bar, col="blue", lty = 2)

# horizontales
abline(h = 2, col="blue", lty = 2)
abline(h = -2, col="blue", lty = 2)

Por ende, el valor de la observación 21 queda en la zona de alerta y por tanto éste es un valor influyente.

Cálculo “manual” de la distancia de Cook 

library(broom)

# Cálculo de los errores
Datosfit=augment(fit)

# Veamos el dataframe resultante
head(Datosfit)
## # A tibble: 6 x 8
##       y     x .fitted  .resid .std.resid   .hat .sigma   .cooksd
##   <dbl> <int>   <dbl>   <dbl>      <dbl>  <dbl>  <dbl>     <dbl>
## 1  1.05     1    2.71 -1.66      -1.05   0.160    1.73 0.104    
## 2  1.82     2    2.95 -1.13      -0.696  0.118    1.76 0.0325   
## 3  3.75     3    3.19  0.561      0.339  0.0861   1.78 0.00540  
## 4  3.06     4    3.42 -0.366     -0.218  0.0636   1.78 0.00161  
## 5  3.60     5    3.66 -0.0551    -0.0327 0.0508   1.78 0.0000285
## 6  8.43     6    3.89  4.54       2.68   0.0478   1.40 0.180

y notamos en la última columna (.cooksd) el valor calculado de la distancia de Cook. Luego, recordemos que al observación 21 es en la que nos estamos enfonfocando, no obstante, puede ocurrir que los datos sean ordenados u el orden en que tenemos los datos cambie por alguna razón. De tal manera, así como tenemos los datos agregaremos una columna con un índice

Datosfit$index=1:length(Datosfit$x)
head(Datosfit)
## # A tibble: 6 x 9
##       y     x .fitted  .resid .std.resid   .hat .sigma   .cooksd index
##   <dbl> <int>   <dbl>   <dbl>      <dbl>  <dbl>  <dbl>     <dbl> <int>
## 1  1.05     1    2.71 -1.66      -1.05   0.160    1.73 0.104         1
## 2  1.82     2    2.95 -1.13      -0.696  0.118    1.76 0.0325        2
## 3  3.75     3    3.19  0.561      0.339  0.0861   1.78 0.00540       3
## 4  3.06     4    3.42 -0.366     -0.218  0.0636   1.78 0.00161       4
## 5  3.60     5    3.66 -0.0551    -0.0327 0.0508   1.78 0.0000285     5
## 6  8.43     6    3.89  4.54       2.68   0.0478   1.40 0.180         6

para después ordenar de mayor a menor los datos de acuerdo a la distancia de Cook asociada

DatosfitDi=Datosfit[order(-Datosfit$.cooksd),]
head(DatosfitDi)
## # A tibble: 6 x 9
##       y     x .fitted .resid .std.resid   .hat .sigma .cooksd index
##   <dbl> <int>   <dbl>  <dbl>      <dbl>  <dbl>  <dbl>   <dbl> <int>
## 1 0.525    12    5.31  -4.79     -3.16  0.234    1.23  1.52      21
## 2 8.43      6    3.89   4.54      2.68  0.0478   1.40  0.180      6
## 3 1.05      1    2.71  -1.66     -1.05  0.160    1.73  0.104      1
## 4 6.06      9    4.60   1.46      0.884 0.0972   1.74  0.0421    19
## 5 1.82      2    2.95  -1.13     -0.696 0.118    1.76  0.0325     2
## 6 5.43      6    3.89   1.53      0.907 0.0478   1.74  0.0207    16

Luego, volveremos a realizar el ajuste del modelo sin considerar a la observación cuyo índice es el 21. Para ello

# Graficamos los datos iniciales
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$") )

# Mostramos la línea correspondiente al ajuste del modelo
abline(fit, col= "red")

# Agregamos un índice en los datos
Datos$index=1:length(Datos$x)

# Realizamos el ajuste del modelo sin la observación con índice 21
fit2=lm(y~x, data=Datos[Datos$index!=21,])

# Mostramos la línea correspondiente al ajuste del modelo sin la observación 21
abline(fit2, col= "blue")

legend('topleft', legend = c("Ajuste con la obs 21", "Ajuste sin la obs 21"),
       lwd = 2, col=c("red", "blue"))

observándose así un cambio pronunciado en las pendientes de las rectas de los modelos. Luego, notemos que en

#Parece que el modelo no tiene "sentido"
summary(fit) 
## 
## Call:
## lm(formula = y ~ x, data = Datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7883 -0.6033 -0.2336  0.7669  4.5351 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   2.4765     0.7984   3.102  0.00587 **
## x             0.2364     0.1210   1.954  0.06562 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.734 on 19 degrees of freedom
## Multiple R-squared:  0.1673, Adjusted R-squared:  0.1235 
## F-statistic: 3.817 on 1 and 19 DF,  p-value: 0.06562

donde se está rechazando \(H_{0}\) en la prueba F asociada de la tabla anova, es decir, parece ser que el modelo no tiene sentido. Después haciendo el mismo análisis en el ajuste del modelo sin la observación 21 tenemos que

#Aquí el modelo sí tiene "sentido"
summary(fit2) 
## 
## Call:
## lm(formula = y ~ x, data = Datos[Datos$index != 21, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6521 -0.5932 -0.2631  0.2773  4.2014 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.67842    0.59336   2.829 0.011131 *  
## x            0.42505    0.09563   4.445 0.000313 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.228 on 18 degrees of freedom
## Multiple R-squared:  0.5233, Adjusted R-squared:  0.4968 
## F-statistic: 19.76 on 1 and 18 DF,  p-value: 0.0003129

Después, recordemos del scatterplot incial que la observación 6 parecía ser un outlier también. Así, consideramos ahora el modelo sin la observación 6 y sin la observación 21

# Graficamos los datos iniciales
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$") )

# Mostramos la línea correspondiente al ajuste del modelo
abline(fit, col= "red")

# Realizamos el ajuste del modelo sin la observación con índice 6 y 21
fit3=lm(y~x, data=Datos[Datos$index!=21 & Datos$index!=6,])

# Mostramos la línea correspondiente al ajuste del modelo sin la observación 21
abline(fit2, col= "blue")

# Mostramos la línea correspondiente al ajuste del modelo sin la observación 6 y 21
abline(fit3, col= "green")

legend('topleft', legend = c("Ajuste con todas las obs", "Ajuste sin la obs 21", "Ajuste sin la obs 6 y sin la obs 21"),
       lwd = 2, col=c("red", "blue", "green"))

Siendo así que no hay mucha diferencia entre el modelo sin la observación 21 y el modelos sin las observaciones 6 y 21, por lo que, en realidd, la observación 6 no es un valor influyente.Además

summary(fit3)
## 
## Call:
## lm(formula = y ~ x, data = Datos[Datos$index != 21 & Datos$index != 
##     6, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.39709 -0.40924 -0.02528  0.47134  1.42899 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.53077    0.34343   4.457 0.000346 ***
## x            0.41163    0.05525   7.450  9.5e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7092 on 17 degrees of freedom
## Multiple R-squared:  0.7655, Adjusted R-squared:  0.7517 
## F-statistic:  55.5 on 1 and 17 DF,  p-value: 9.497e-07

en esta ocasión el modelo sí tiene “sentido”.

Visualización de varios modelos

Dado que en la parte anterior ha sido necesario el trabajar con 3 modelos, tenemos una opción en R que nos permite visualizar y contrastar los resultados obtenidos en cada modelo. Para ello

# Librerías necesarias
library(huxtable)
library(jtools)

# Declaramos los modelos de los cuales deseamos contrastar sus resultados 
export_summs(fit, fit2, fit3, scale = FALSE)
Model 1Model 2Model 3
(Intercept)2.48 **1.68 *  1.53 ***
(0.80)  (0.59)   (0.34)   
x0.24   0.43 ***0.41 ***
(0.12)  (0.10)   (0.06)   
N21      20       19       
R20.17   0.52    0.77    
*** p < 0.001; ** p < 0.01; * p < 0.05.

Con lo anterior podemos notar que los valores de las betas para los modelos 2 (donde no se consideró a la observación 21) y 3 (donde no se consideraron a las observaciones 6 y 21) han cambiado muy poco. De esta manera tenemos que, incluir o no a la observación 6, no trae cambios significativos a la estimación de las betas en el modelo por lo que no lo quitaremos del modelo.

Regresamos de nuevo a revisar algunos supuestos

#Gráficas de R para revisar algunos de los supuestos
par(mfrow=c(1,3)) 
par(mar=c(4, 5, 2, 1))
plot(fit2, 1)
plot(fit2, 3)
plot(fit2, 2)

en donde, en las primeras dos gráficas, los patrones ya no son tan marcados como en las primeras gráficas del análisis de los supuestos. Después, volvemos a analizar

#Gráficas de R especificas para revisar outliers y valores influyentes

par(mfrow=c(1,2)) 
par(mar=c(4, 5, 2, 1))
plot(fit2, 4)

plot(fit2, 5, cook.levels = c(4/(dim(Datos)[1]-1), 0.5, 1.0))
abline(v = 4/(dim(Datos)[1]-1), col="blue", lty = 2)
abline(v = 6/(dim(Datos)[1]-1), col="blue", lty = 2)

que, después de quitar la observación 21, el análisis sobre los outliers nos alerta de la observación 6, sin embargo, podemos ver en la segunda gráfica como ésta no se encuentra en la zona de alerta. No obstante, anteriormente ya vimos que quitarlo o no, no cambia significativamente el valor de estimación de los betas.

Importante: Antes de quitar los outliers del modelo debemos ir al contexto del problema. Además, hay que apoyarse en este caso de los expertos y de la teoría que puediera haber detrás para guiarnos.