Covarianza cero

Recordemos que para analizar este supuesto debemos revisar el proceso de generación de los datos que se están utilizando en el análisis. Siendo lo deseable poder argumentar que la recolección de los datos proviene de elementos seleccionados de forma aleatoria, es decir, no se seleccionó el conjunto de datos en algún grupo particular que puediera ocasionar grupos de errores.

Linealidad, homocedasticidad y normalidad

Ahora bien, en este ejemplo continuaremos trabajando con los datos de la compañía Toluca

# Datos
library(ALSM)
Datos=TolucaCompany

# Scatterplot inicial 
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$"))

Pasamos después al análisis de los supuestos

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

# R tiene una función para obtener errores de forma automatizada
library(broom)
Datosfit=augment(fit)

par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2)) 
plot(Datosfit$.fitted, Datosfit$.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e$"))
plot(Datos$x, Datosfit$.resid, xlab = TeX("$X$"), ylab=TeX("$e$") )

y vemos las gráfica que R trae por defecto

#gráficas para linealidad-homocedasticidad-normalidad
par(mar=c(4, 5, 2, 1))
par(mfrow=c(1,3))
plot(fit, 1)
plot(fit, 3)
plot(fit, 2)

Supuesto de la covarianza cero

Para analizar este supuesto contamos con diversas herramientas. Comenzamos por observar el scatterplot de los residuales estandarizados (studentizados) versus el índice de las observaciones

#gráfica sobre el índice de los datos

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

donde, como tenemos 25 observaciones en nuestro dataframe tendremos hasta el índice de observación 25. De igual manera buscamos no ver algún patrón y/o tener una “nube” de puntos con la misma longitud. Debemos de tener cuidado al trabajar con esta gráfica pues, si ordenamos los índice de acuerdo a la \(x\), estamos analizando una gráfica que ya trabajamos para el supuesto de la homocedasticidad. Además, puede que al ver un patrón, donde el índice fue ordenado con base a la \(x\) (por ejemplo, cuando nos dan la base de datos previamente ordenada), este dándonos información sobre otras cosas (por ejemplo sobre la homocedasticidad) y no necesariamente sobre la aleatoriedad de los datos.

Después graficamos el autocorrelograma de los errores estandarizados utilizando acf() la cual viene incluida en R por defecto

#autocorrelograma de los errores

acf(Datosfit$.std.resid)

donde R nos ayuda con las líneas horizonatales punteadas y donde nos gustaría tener los valores de las correleaciones (de la correlación de Pearson) entre estas rectas (lo cual si pasa), lo que nos sirve para ver si podemos detectar algún cierto tipo de asociación monótona entre los errores que contraindicaría el supuesto.

Dada la gráfica anterior podemos intuir que el supuesto de la covarianza cero se está cumpliendo. Para complementar lo anterior pasamos a efectuar las pruebas de rachas:

#Prueba de rachas
library(lawstat)
lawstat::runs.test(Datosfit$.std.resid, plot.it = FALSE)
## 
##  Runs Test - Two sided
## 
## data:  Datosfit$.std.resid
## Standardized Runs Statistic = -1.015, p-value = 0.3101
library(randtests)
randtests::runs.test(Datosfit$.std.resid)
## 
##  Runs Test
## 
## data:  Datosfit$.std.resid
## statistic = -0.83485, runs = 11, n1 = 12, n2 = 12, n = 24, p-value =
## 0.4038
## alternative hypothesis: nonrandomness

en las cuales no se está rechazando \(H_{0}\), es decir, es plausible asumir \(H_{0}\) (las observaciones se puede considerar aleatorias). Notemos que en la prueba lawstat colocamos plot.it=FALSE, si cambiamos a plot.it=TRUE nos mostrará una gráfica ilustrativa

lawstat::runs.test(Datosfit$.std.resid, plot.it = TRUE)

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

en la cual nos indica, tomando como referencia la recta horizontal en cero, si los 25 eerores caen por encima o por debajo de dicha recta horizontal mostrándose así las rachas. Así, hasta ahora no hemos encontrado evidencia en contra de la covarianza cero.

Pasamos a realizar la prueba Durbin-Watson para contrastar la existencia de una posible correlación de orden 1

#Prueba para autocorrelación de orden 1

library(lmtest)
lmtest::dwtest(fit)
## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 1.4318, p-value = 0.08079
## alternative hypothesis: true autocorrelation is greater than 0

en la cual no se está rechazando \(H_{0}\) (i.e. es plausible asumir que los errores no están autocorrelacionados). Por otro lado, la librería car tiene una función para revisar más órdenes

library(car)  

# Probamos hasta el orden 5
durbinWatsonTest(fit, max.lag=5)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.25931930      1.431790   0.186
##    2      0.08274460      1.741432   0.696
##    3      0.01410007      1.844020   0.960
##    4      0.22115877      1.299156   0.232
##    5     -0.07531258      1.653310   0.910
##  Alternative hypothesis: rho[lag] != 0

en las cuales, cada una, no se está rechazando \(H_{0}\).

Nota importante: Se recomienda hacer la revisión del contexto al inicio del análisis. El uso de las herramientas de diagnóstico sería al final del modelado, pues usan los errores y estos podrían mostrar patrones por la falta del cumplimiento de los otros supuestos.

En caso de encontrar algún indicio del no cumplimiento de este supuesto se pueden usar otros modelos, por ejemplo: modelos de series de tiempo, modelos lineales mixtos, etc.