Covarianza cero

Comenzamos por cargar los datos

library(latex2exp)
Datos = read.csv("ejemplo4.csv")

y después realizamos el scatterplot inicial

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

# Gráfica
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$"))

abline(fit, col="red")

Linealidad, homocedasticidad y normalidad

Dado que este ejemplo ya fue trabajado, para recordar que los estos supuestos se cumplen cuando transformamos a la variable \(x\) bastará con ver las gráficas que R trae por defecto. Primero

#Se transforma X a X^2 para lograr linealidad
Datos$Xprima=Datos$x^2

# Se realiza el ajuste con los datos transformados
fit2=lm(y~Xprima, data=Datos)

# Gráfica
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(Datos$Xprima, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$"))

abline(fit2, col="red")

Y finalmente

#gráficas para linealidad-homocedasticidad-normalidad

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

Supuesto de la covarianza cero

  • Comenzamos por ver el gráfico sobre el índice de los datos

    #gráfica sobre el índice de los datos

# Librería para calcular los errores
library(broom)
Datosfit2=augment(fit2)

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

    del cual no puede observarse algún patrón evidente.

  • Seguimos con el autocorrelograma

    #autocorrelograma de los errores

acf(Datosfit2$.std.resid)

    donde tampoco vemos algún valor que salga de los valores de las rectas horizontales punteadas. Así, hasta el momento no hemos encontrado evidencia en contra del supuesto de la covarianza cero.

  • Pasamos a realizar las pruebas de rachas

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

    ## 
##  Runs Test - Two sided
## 
## data:  Datosfit2$.std.resid
## Standardized Runs Statistic = -0.40204, p-value = 0.6877
    library(randtests)
randtests::runs.test(Datosfit2$.std.resid)
    ## 
##  Runs Test
## 
## data:  Datosfit2$.std.resid
## statistic = -0.40204, runs = 49, n1 = 50, n2 = 50, n = 100, p-value =
## 0.6877
## alternative hypothesis: nonrandomness
    siendo así que no se está rechanzado \(H_{0}\), lo cual es justamente lo que buscamos.

  • Finalmente la prueba de Durbin-Watson

    #Prueba para autocorrelación de orden 1 (Durbin-Watson)
lmtest::dwtest(fit2)
    ## 
##  Durbin-Watson test
## 
## data:  fit2
## DW = 2.0289, p-value = 0.5368
## alternative hypothesis: true autocorrelation is greater than 0
    library(car)  
# La librería car tiene una función para revisar más órdenes
durbinWatsonTest(fit2, max.lag=5)
    ##  lag Autocorrelation D-W Statistic p-value
##    1     -0.02506064      2.028857   0.898
##    2     -0.10994695      2.196005   0.274
##    3      0.14725154      1.630055   0.112
##    4     -0.09643929      2.064546   0.488
##    5     -0.01099229      1.887869   0.956
##  Alternative hypothesis: rho[lag] != 0
    no rechazándose \(H_{0}\) en todas las pruebas. Concluimos pues que no hay evidencia suficiente en contra del supuesto de la covarianza cero.

Otro caso

Analicemos este supuesto en el modelo antes de transformar los datos

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

# Gráfica original
par(mar=c(4, 5, 3, 1))
par(mfrow=c(1,3))
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$") )

# Gráfica sobre el índice de los datos
plot(1:length(Datosfit$.std.resid), Datosfit$.std.resid, xlab = TeX("$i$"), ylab=TeX("$e_s$"))

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

## 
##  Runs Test - Two sided
## 
## data:  Datosfit$.std.resid
## Standardized Runs Statistic = -5.6285, p-value = 1.818e-08
durbinWatsonTest(fit, max.lag=5)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.6951288     0.5710322   0.000
##    2       0.1907604     1.5756213   0.050
##    3      -0.3047335     2.5662998   0.004
##    4      -0.6735497     3.2966808   0.000
##    5      -0.7799410     3.4816793   0.000
##  Alternative hypothesis: rho[lag] != 0

donde en cada prueba se está rechazando \(H_{0}\), es decir, es plausible asumir que los errores están autocorrelacionados (\(H_{a}\)).

Finalmente comparemos las gráficas:

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

# Datos transformados
lawstat::runs.test(Datosfit2$.std.resid, plot.it = TRUE)
## 
##  Runs Test - Two sided
## 
## data:  Datosfit2$.std.resid
## Standardized Runs Statistic = -0.40204, p-value = 0.6877
# Datos no transformados
lawstat::runs.test(Datosfit$.std.resid, plot.it = TRUE)

## 
##  Runs Test - Two sided
## 
## data:  Datosfit$.std.resid
## Standardized Runs Statistic = -5.6285, p-value = 1.818e-08

siendo la gráfica de la derecha en la cual se hace más visible las rachas.