Homocedasticidad

Continuemos con la práctica pasada y realicemos la verificación de este supuesto. Primero cargamos los paquetes que utilizaremos

después cargamos los datos

Datos=TolucaCompany
head(Datos)
##    x   y
## 1 80 399
## 2 30 121
## 3 50 221
## 4 90 376
## 5 70 361
## 6 60 224

Ahora bien, comenzamos por analizar un diagrama de dispersión \(x\) vs \(y\), en el cual se debe observar una nube de puntos sobre el eje \(y\) de la misma longitud sobre todos los valores de \(x\):

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

o alternativamente

ggplot(Datos, aes(x, y)) + geom_point() + geom_smooth()

por lo que podemos ver que la varianza es constante y lineal.

Después aplicaremos un método más general consistente en analizar los errores residuales. Para ello calcularemos los errores estandarizados

\[ e_{is}=\frac{y_{i}-\hat{y_{i}}}{\sqrt{\hat{\sigma^{2}}(1-h_{i})}} \]

Calculos “a mano” 

# Cálculos

xbar=mean(Datos$x)
SSx=sum((Datos$x-xbar)^2)
ybar=mean(Datos$y)
SSy=sum((Datos$y-ybar)^2)
SSxy=sum((Datos$y-ybar)*(Datos$x-xbar))
beta1=SSxy/SSx
beta0=ybar-beta1*xbar
n=length(Datos$x)
Datos$yhat=beta0+beta1*Datos$x
Datos$error=Datos$y-Datos$yhat
MSE=sum((Datos$error)^2)/(n-2)
# Cálculo de hi

Datos$hi=(1/n+ (Datos$x-xbar)^2/SSx)

# Cálculo de los errores estandarizados

Datos$errorSt=Datos$error/sqrt(MSE*(1-Datos$hi))

head(Datos)
##    x   y     yhat     error         hi    errorSt
## 1 80 399 347.9820  51.01798 0.04505051  1.0693154
## 2 30 121 169.4719 -48.47192 0.12080808 -1.0588176
## 3 50 221 240.8760 -19.87596 0.06020202 -0.4199365
## 4 90 376 383.6840  -7.68404 0.06020202 -0.1623473
## 5 70 361 312.2800  48.72000 0.04000000  1.0184611
## 6 60 224 276.5780 -52.57798 0.04505051 -1.1020124

en lo que sigue es visualizar el diagrama de dispersión de \(x\) vs \(e_{s}\) o bien \(\hat{y}\) vs \(e_{s}\).

# Gráficas

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

plot(Datos$yhat, Datos$errorSt, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$") )
plot(Datos$x, Datos$errorSt , xlab = TeX("$x$"), ylab=TeX("$e_{s}$"))
plot(Datos$x, sqrt(abs(Datos$errorSt)), xlab = TeX("$x$"), ylab=TeX("$\\sqrt{|e_{s}|}$") )

en este caso se debe presentar una nube de puntos sobre el eje \(e_{s}\) centrada en cero y de la misma longitud sobre los valores de \(x\). Además no podemos observar patrones o concluir algo de las primeras dos gráficas. Notemos que en la regresión lineal simple las primeras dos coinciden. En la tercer gráfica representa un caso extremo, exagerando si existe algún patrón. Esta última gráfica está implementada en R, en esta última gráfica se nos pueden presentar dudas sobre si el supuesto de la homocestacidad se está dando, además de que las primeras gráficas no me aportan mucho.

Por ende debemos recurrir a una prueba de hipótesis para complementar estas pruebas gráficas. Para ello constrastaremos

\[ H_{0}: \sigma^{2}=\gamma_{0}\ \ vs \ \ H_{a}: \sigma^{2}=\gamma_{0}+\gamma_{1}x \ \ \ \ \ \textrm{Test Breusch-Pagan} \]

Para ello realizamos el modelo lineal de regresión simple, donde buscamos no rechazar \(H_{0}\), es decir, que sea plausible asumir que la varianza no depende de forma lineal de \(x\), lo cual se consigue si el p-value es mayor a la significancia.

# Modelo lineal

fit=lm(y~x, data=Datos)

Opción 1: Uso de residuales studentilizados

# Librería necesaria
library(lmtest)


lmtest::bptest(fit)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 1.1326, df = 1, p-value = 0.2872

en el cual utilizamos el Test Breusch-Pagan con los valores studentilizados:

\[ e_{is}=\frac{y_{i}-\hat{y_{i}}}{\sqrt{\hat{\sigma^{2}_{(-i)}}(1-h_{i})}} \]

donde \(\hat{\sigma^{2}_{(-i)}}\) es el estimador de \(\sigma^{2}\) en el modelo de regresión obtenido al ajustar considerando los datos sin la \(i\)-ésima unidad. De la información obtenida tenemos que \(p-value>\alpha\) por lo que no se estaría rechazando \(H_{o}\), es decir, no se encontró evidencia contra la homocesticidad.

Opción 2: Uso de los residuales estandarizados

library(car)
car::ncvTest(fit)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.8209192, Df = 1, p = 0.36491

del cual podemos concluir igualmente que no se rechaza \(H_{0}\).

#R tiene una gráfica propia para verificar homocedasticidad (plot(fit, 3))
par(mfrow=c(1,2)) 
par(mar=c(4, 5, 1, 1))
plot(fit, 3)
plot(Datos$yhat, sqrt(abs(Datos$errorSt)), xlab = TeX("$\\widehat{y}$"), ylab=TeX("$\\sqrt{|e_{s}|}$") )

teniendo a la derecha el gráfico que nos da R y a la derecha la que habíamos planteado al inicio. Notemos que el gráfico de la derecha nos enseña una curva en rojo, la cual nos gustaría que fuera recta. Lo anterior nos podría implicar evidencia contra la homocestacidad, sin embargo por la prueba de hipótesis que realizamos si quitamos esa línea roja tampoco es evidente que se rompa la homcestacidad.

Finalmente, R tiene una función para crear errores de forma automatizada

library(broom)

Datosfit=augment(fit)
head(Datosfit)
## # A tibble: 6 x 8
##       y     x .fitted .resid .std.resid   .hat .sigma  .cooksd
##   <int> <int>   <dbl>  <dbl>      <dbl>  <dbl>  <dbl>    <dbl>
## 1   399    80    348.  51.0       1.07  0.0451   48.7 0.0270  
## 2   121    30    169. -48.5      -1.06  0.121    48.7 0.0770  
## 3   221    50    241. -19.9      -0.420 0.0602   49.7 0.00565 
## 4   376    90    384.  -7.68     -0.162 0.0602   49.9 0.000844
## 5   361    70    312.  48.7       1.02  0.04     48.8 0.0216  
## 6   224    60    277. -52.6      -1.10  0.0451   48.6 0.0286

y los gráficos

par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2)) 
plot(Datosfit$.fitted, Datosfit$.std.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$")   )
plot(Datos$yhat, Datos$errorSt, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$") )

en la cual contrastamos la gráfica obtenida de los errores que calculó en automático R contra la que nosotros calculamos.

Nota: Para efectuar un análisis podríamos visualizar la gráfica de los errores que R calcula junto con una prueba de hipótesis.