Verificación de supuestos

Análisis descriptivo

# Cargamos los datos
Datos <- read.csv("EjercicioVerSup.csv", header=TRUE)

Vemos una breve información sobre nuestras columnas

str(Datos)
## 'data.frame':    400 obs. of  4 variables:
##  $ y : num  18.9 14.3 13 14.3 12.8 ...
##  $ X3: num  0.202 0.66 0.648 0.661 0.875 ...
##  $ X4: num  2.1 2.14 2.04 2.14 2.06 ...
##  $ X1: chr  "A3" "A3" "A3" "A3" ...
# Hacemos que la variable X1 sea de tipo factor
Datos$X1<-as.factor(Datos$X1)

Algunas esdadísticas descriptivas

summary(Datos)
##        y              X3             X4        X1     
##  Min.   : 7.3   Min.   :0.10   Min.   :1.32   A1:100  
##  1st Qu.:12.6   1st Qu.:0.34   1st Qu.:1.86   A2:100  
##  Median :14.2   Median :0.54   Median :2.01   A3:100  
##  Mean   :14.9   Mean   :0.55   Mean   :2.00   A4:100  
##  3rd Qu.:16.6   3rd Qu.:0.77   3rd Qu.:2.13           
##  Max.   :28.5   Max.   :1.00   Max.   :2.64

Observemos el comportamiento de algunas estadísticas de los niveles de la variable categórica con respecto a la variable y:

Datos %>% group_by(X1) %>% 
  summarise(Observaciones = n(),
            Media = round(mean(y), 2),
            Mediana = round(median(y), 2),
            Varianza = round(var(y), 2))
## # A tibble: 4 x 5
##   X1    Observaciones Media Mediana Varianza
##   <fct>         <int> <dbl>   <dbl>    <dbl>
## 1 A1              100  14.8    14.6    10.9 
## 2 A2              100  15.1    14.4    10.7 
## 3 A3              100  15.2    14.3    12.6 
## 4 A4              100  14.5    13.7     9.83

Vemos algunos gráficos

ggpairs(Datos)

marginalmente podemos ver que el supuesto de la linealidad parece fallar, al menos en y contra X3 y de y contra X4. Veamos más a detalle los diagramas de dispersión

ggplot(Datos, aes(x=X3, y=y, color=X1))+
  geom_point()

ggplot(Datos, aes(x=X4, y=y, color=X1))+
  geom_point()

en lo cual vemos con más fuerza que el supuesto de la linealidad no se está cumpliendo; de igual manera parece que tenemos problemas sobre la homocedasticidad de los datos.

Linealidad

Nos auxiliaremos del gráfico que trae R base para el análisis de la linealidad

# Realizamos el ajuste
Ajuste <- lm(formula = y ~ X1 + X3 + X4, data = Datos)
summary(Ajuste)  
## 
## Call:
## lm(formula = y ~ X1 + X3 + X4, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.685 -1.098 -0.448  0.641  6.915 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.0973     0.7831   -0.12   0.9012    
## X1A2          0.5692     0.2174    2.62   0.0092 ** 
## X1A3         -0.4442     0.2184   -2.03   0.0427 *  
## X1A4          0.1321     0.2175    0.61   0.5440    
## X3           -8.8220     0.3021  -29.20   <2e-16 ***
## X4            9.8799     0.3802   25.98   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.5 on 394 degrees of freedom
## Multiple R-squared:  0.788,  Adjusted R-squared:  0.785 
## F-statistic:  292 on 5 and 394 DF,  p-value: <2e-16
# Valores ajustados (y-gorro) vs Residuales 
plot(Ajuste, 1)

del cual nos gustaría que la línea roja fuera lo más parecida a la línea horizontal en 0, y aunque no buscamos una coincidencia exacta, sí vemos que la línea está muy distorsionada.

Nos auxiliamos de lo siguiente:

# Analisis de los residuales por cada variable incluida 
# en el modelo

# Esperamos no rechazar (p-value > 0.05)
library(car)

residualPlots(Ajuste)

##            Test stat Pr(>|Test stat|)    
## X1                                       
## X3             30.71          < 2e-16 ***
## X4              1.45             0.15    
## Tukey test      7.51          5.8e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

donde buscamos no rechazar en cada prueba para asumir plausible el supuesto de la linealidad, no obstante, sólo para la variable X4 se consigue no rechazar.

Encontramos fuerte evidencia en contra del supuesto de linealidad. Al parecer nuestro modelo no esta funcionando

Varianza constante

Ahora nos auxiliamos de la gráfica 3 para analizar este supuesto

# Valores ajustados (y-gorro) vs (|residuales std|)^{1/2}
plot(Ajuste, 3)

vemos una curva en rojo, la cual nos gustaría que fuera recta. Lo anterior nos podría implicar evidencia contra la homocestacidad. Para complementar lo anterior realizamos una prueba de hipótesis (Breusch-Pagan)

library(lmtest)

# Esperamos no rechazar (p-value > 0.05)
lmtest::bptest(Ajuste)
## 
##  studentized Breusch-Pagan test
## 
## data:  Ajuste
## BP = 55, df = 5, p-value = 2e-10

se rechaza la hipótesis nula, es decir, tenemos fuerte evidencia en contra del supuesto de la homocedasticidad.

Encontramos fuerte evidencia en contra del supuesto de Homocedasticidad

Normalidad

Mostramos la Q-Q plot asociada

# Q-Q Plot
plot(Ajuste, 2)

vemos que gran cantidad de puntos en las colas se alejan de la recta. Lo anterior nos da indicios de que el supuesto de la normalidad tampoco se está cumpliendo. Corroboramos realizando pruebas de hipótesis:

# Utilizamos la librería broom para calcular los errores de forma automática
library(broom)

# Guardamos los errores en columnas de nuestro dataframe de los datos 
Datosfit <- augment(Ajuste)
shapiro.test(Datosfit$.std.resid)
## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfit$.std.resid
## W = 0.8, p-value <2e-16

se rechaza la hipótesis nula y por ende encontramos evidencia en contra del supuesto de al normalidad.

nortest::lillie.test(Datosfit$.std.resid) 
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Datosfit$.std.resid
## D = 0.1, p-value <2e-16
normtest::jb.norm.test(Datosfit$.std.resid)
## 
##  Jarque-Bera test for normality
## 
## data:  Datosfit$.std.resid
## JB = 654, p-value <2e-16

Encontramos fuerte evidencia en contra del supuesto de Normalidad

Aleatoriedad

# Grafica sobre el indice de los datos
par(mar=c(4, 5, 3, 1))
library(latex2exp)
plot(1:length(Datosfit$.std.resid), 
     Datosfit$.std.resid, 
     xlab=TeX("$i$"), 
     ylab=TeX("$e_s$")   )

tal parece que nuestros datos sí son aleatorios pues, a primera vista, no podemos detectar tendencias. Continuando

# autocorrelograma de los errores
acf(Datosfit$.std.resid)

la mayoría de los lags se quedan dentro del intervalo de confianza. Corroboramos realizando pruebas de rachas

# Prueba de rachas
library(lawstat)
lawstat::runs.test(Datosfit$.std.resid, plot.it = TRUE)

## 
##  Runs Test - Two sided
## 
## data:  Datosfit$.std.resid
## Standardized Runs Statistic = 0.7, p-value = 0.5

donde no rechazamos la hipótesis nula por lo cual es plausible asumir que nuestros datos son aleatorios.

No encontramos fuerte evidencia en contra del supuesto de Aleatoriedad, entonces podemos considerar plausible que se cumple con este supuesto.

Transformación

# Variables independientes
# No incluiremos a la variable X1
boxTidwell(y ~ X3 + X4, ~X1, data=Datos)
##    MLE of lambda Score Statistic (z) Pr(>|z|)    
## X3            -1                  46   <2e-16 ***
## X4             2                   5    2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## iterations =  3

vemos las lambdas recomendables para realizar las transformaciones.

Ajuste2 <- lm(y ~ X1 + I(X3^-1) + I(X4^2), data=Datos)
summary(Ajuste2)
## 
## Call:
## lm(formula = y ~ X1 + I(X3^-1) + I(X4^2), data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28867 -0.05779 -0.00364  0.06473  0.27427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.98917    0.02726   36.29   <2e-16 ***
## X1A2         0.49833    0.01372   36.33   <2e-16 ***
## X1A3        -0.51184    0.01378  -37.14   <2e-16 ***
## X1A4        -0.00847    0.01373   -0.62     0.54    
## I(X3^-1)     1.50261    0.00270  557.45   <2e-16 ***
## I(X4^2)      2.50155    0.00602  415.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.097 on 394 degrees of freedom
## Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 
## F-statistic: 9.31e+04 on 5 and 394 DF,  p-value: <2e-16

Linealidad

plot(Ajuste2, 1)

residualPlots(Ajuste2)

##            Test stat Pr(>|Test stat|)
## X1                                   
## I(X3^-1)        0.80             0.43
## I(X4^2)        -0.17             0.86
## Tukey test      0.47             0.64

vemos mejores resultados en los gráficos y en las pruebas de hipótesis.

Varianza constante

plot(Ajuste2, 3)

lmtest::bptest(Ajuste2)
## 
##  studentized Breusch-Pagan test
## 
## data:  Ajuste2
## BP = 2, df = 5, p-value = 0.8

notamos también una mejora entorno al supuesto de la homocedasticidad.

Normalidad

plot(Ajuste2, 2)

DatosAjuste2 <- augment(Ajuste2)
shapiro.test(DatosAjuste2$.std.resid)
## 
##  Shapiro-Wilk normality test
## 
## data:  DatosAjuste2$.std.resid
## W = 1, p-value = 0.6
nortest::lillie.test(DatosAjuste2$.std.resid)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  DatosAjuste2$.std.resid
## D = 0.03, p-value = 0.6

es plausible asumir que el supuesto de la normalidad se está cumpliendo.

Aleatoriedad

# Aleatoriedad

par(mar=c(4, 5, 3, 1))
library(latex2exp)
plot(1:length(DatosAjuste2$.std.resid), 
     DatosAjuste2$.std.resid,
     xlab=TeX("$i$"),
     ylab=TeX("$e_s$")   )

# autocorrelograma de los errores
acf(DatosAjuste2$.std.resid)

es plausible asumir que el supuesto de la aleatoriedad se está cumpliendo.

Interpretación

Observamos que el modelo

\[ E(y; X1, X3, X4) = \beta_0 + \beta_1\cdot X_1A_2 + \beta_2\cdot X_1A_3 + \beta_3\cdot X_1A_4 + \beta_4\cdot X_3^{-1} + \beta_5\cdot X_4^{2} \]

cumple con los supuestos de linealidad, homocedasticidad, normalidad y aleatoriedad.

summary(Ajuste2)
## 
## Call:
## lm(formula = y ~ X1 + I(X3^-1) + I(X4^2), data = Datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28867 -0.05779 -0.00364  0.06473  0.27427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.98917    0.02726   36.29   <2e-16 ***
## X1A2         0.49833    0.01372   36.33   <2e-16 ***
## X1A3        -0.51184    0.01378  -37.14   <2e-16 ***
## X1A4        -0.00847    0.01373   -0.62     0.54    
## I(X3^-1)     1.50261    0.00270  557.45   <2e-16 ***
## I(X4^2)      2.50155    0.00602  415.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.097 on 394 degrees of freedom
## Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 
## F-statistic: 9.31e+04 on 5 and 394 DF,  p-value: <2e-16

Interpretación de \(R^2\) y la prueba de la tabla ANOVA

  • \(R^2\) = 0.9992: 99.92% de la variabilidad en \(y\) es explicada por el modelo de RLM que incluye las variables X1, X3 y X4.

  • \(p-value\) = 2.2e-16 < 0.05 asociado a la prueba:

\[ H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = \beta_5 = 0 \ \ vs \ \ H_a: \beta_i \neq 0 \ \ \textrm{para alguna} \ \ i = 1,2,3,4,5 \]

Se rechaza \(H_0\), por lo tanto hay evidencia estadística para decir que el modelo tiene sentido, es decir, alguna de las variables \(X_1, X_3^{-1}\) o \(X_4^2\) ayuda a modelar \(E(Y)\).

Observación: Dado el summary(Ajuste2) es plausible considerar un modelo reducido donde se colapsen los niveles A1 y A4 de la variable X1.