Homocedasticidad: regresión ponderada

Trabajaremos con la relación entre el salario y los años de experiencia de los trabajadores en la compañía Initech. Para ello descargamos los datos

Datos = read.csv("initech.csv")

donde consideramos a \(y\) como salario (salary) y \(x\) a los años de experiencia (years).

Después, comenzamos por el scatterplot de los datos

# Librería para poder escribir las etiquetas en Látex
library(latex2exp)

# Ajuste del modelo
fit <- lm(salary ~ years, data=Datos)

# Scatterplot
par(mar=c(4, 5, 1, 1))
plot(Datos$years, Datos$salary, xlab = TeX("$X=years$"), ylab=TeX("$salary$"))

abline(fit$coefficients, col="red")

donde podemos ver que el supuesto de la linealidad no se está cumpliendo. Recordemos que al transformar \(x\) sólo arreglaremos el problema de la linealidad pero nunca el problema de la homocedasticidad. Utilicemos la siguiente transformación (sacada bajo la manga) en \(x\) para arreglar el problema de la linealidad. Después

# Transformamos los datos
Datos$Xprima=(Datos$years+10)^2

# Ajustamos el modelo a los datos transformados
fit2=lm(salary ~ Xprima, data=Datos)

# Graficamos
plot(Datos$Xprima, Datos$salary, xlab = TeX("$X'=(years+10)^2$"), ylab=TeX("$salary$"))

abline(fit2$coefficients, col="red")

y utilizamos la gráfica que trae por defecto R para notar el problema de la linealidad

plot(fit2, 1)

con lo cual parece indicar que el problema de la linealidad se ha arreglado. Después veamos que el supuesto de la homocedasticidad no se cuemple

# Gráfica que trae R para ver el problema de la homocedasticidad
plot(fit2, 3)

en conjunto con las pruebas de hipótesis

lmtest::bptest(fit2)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit2
## BP = 27.895, df = 1, p-value = 1.281e-07
car::ncvTest(fit2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 30.97662, Df = 1, p = 2.6116e-08

que nos conducen a rechazar a \(H_{0}\), es decir, rechazar que la varianza no depende de forma lineal en \(x\), no obstante recordemos que buscamos en este caso no rechazar.

Después, recordmos que la gráfica que trae R incluida exagera el comportamiento que puedan tener nuestros datos para resaltar posibles patrones. Sin embargo, podemos visualizar los siguientes diagramas de dispersión que no exageran tanto los posibles patrones. Estos son los diagramas de dispersión contrastando los residuales estandarizados versus \(\hat{y}\) y versus \(x'\) (Xprima), que recordemos que para el caso de la regresión lineal simple éstos coinciden. Para ello

# Cargamos la librería broom para el cálculo de los errores
library(broom)

# Se agregan estos errores al dataframe de los datos
Datosfit=augment(fit2)

# Graficamos
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(Datosfit$Xprima, Datosfit$.std.resid, xlab = TeX("$X'$"), ylab=TeX("$e_{s}$"))

Notemos que nos gustaría ver “una nube” de puntos casi de la misma longitud sobre el eje de las \(x\).

Después, usamos los errores al cuadrado para identificar la gráfica patrón

# Usamos los residuales (errores observados) al cuadrado
Datosfit$res2=Datosfit$.resid^2

par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2))
plot(Datosfit$Xprima, Datosfit$res2, xlab = TeX("$X'$"), ylab=TeX("$e^{2}$") )
plot(Datosfit$Xprima^2, Datosfit$res2, xlab = TeX("$X'^2$"), ylab=TeX("$e^{2}$") )

intentanto ver si podemos identificar algún tipo de función patrón en ellas, lo cual es muy difícil de observar pues podemos notar mucho “ruido” en nuestros datos.

En general, para observar el comportamiento es mejor agrupar sobre intervalos de \(x\) y sacar promedios. En este caso, hay varios valores de una misma \(x\), así que sólo promediaremos. Para conseguir lo anterior utilizaremos DatosfitAgrup= Datosfit %>% group_by(Xprima) %>% summarise(res2media=mean(res2)) del cual, primero se agruparán los valores de la Xprima en los que sean iguales

head(Datosfit)
## # A tibble: 6 x 9
##   salary Xprima .fitted .resid .std.resid   .hat .sigma   .cooksd       res2
##    <int>  <dbl>   <dbl>  <dbl>      <dbl>  <dbl>  <dbl>     <dbl>      <dbl>
## 1  41504    121  30888. 10616.     0.439  0.0296 24679. 0.00294   112708658.
## 2  32619    121  30888.  1731.     0.0715 0.0296 24703. 0.0000781   2997862.
## 3  44322    121  30888. 13434.     0.555  0.0296 24664. 0.00470   180484000.
## 4  40038    144  35246.  4792.     0.198  0.0277 24698. 0.000557   22959930.
## 5  46147    144  35246. 10901.     0.450  0.0277 24678. 0.00288   118824217.
## 6  38447    144  35246.  3201.     0.132  0.0277 24701. 0.000248   10244174.

por ejemplo para aquellos en que Xprima es igual a 144 y promediaremos sus respectivos valores de los residuales al cuadrado. Así

# Cargamos la librería necesaria
library(tidyverse)

# Agrupamos los datos:
# group_by: agrupamos de acuerdo a los datos en Xprima
# summarise: obtendremos el promedio de los residuales al cuadrado
DatosfitAgrup= Datosfit %>% group_by(Xprima) %>% summarise(res2media=mean(res2))

# Graficamos
par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2)) 
plot(DatosfitAgrup$Xprima, DatosfitAgrup$res2media, xlab = TeX("$X'$"), ylab=TeX("$e^{2}$") )
plot(DatosfitAgrup$Xprima^2, DatosfitAgrup$res2media, xlab = TeX("$X'^2$"), ylab=TeX("$e^{2}$") )

Del cual, a primera vista, en el gráfico de la izquierda podemos ver que el “ruido” en los datos ha desaparecido. Luego, analizando ambas gráficas podemos pensar que el patrón puede ser una recta o una curva cuadrática. Probaremos pues el caso de la recta idéntidad y el de una curva cuádratica.

Procederemos después ha efectuar el ajuste de los datos con la regresión ponderada para el caso en que tomamos la recta identidad, para ello

fitpond <- lm(salary~Xprima,weights = 1/(Xprima), data=Datos)
summary(fitpond)
## 
## Call:
## lm(formula = salary ~ Xprima, data = Datos, weights = 1/(Xprima))
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -1991.49  -662.08   -22.32   647.60  2211.16 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12385.625   3107.348   3.986 0.000129 ***
## Xprima        181.855      6.627  27.444  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 933.1 on 98 degrees of freedom
## Multiple R-squared:  0.8849, Adjusted R-squared:  0.8837 
## F-statistic: 753.1 on 1 and 98 DF,  p-value: < 2.2e-16

Donde una vez que identificamos a la función \(g()\), que en este caso es la idéntidad, los pesos (weights) se definen como \[w_{i}=1/g(x_{i})\] y como \(g\) es la idéntidad tenemos que los pesos están dados por 1/Xprima.

Procedemos a hacer las pruebas respectivas a este nuevo modelo:

# Gráfica
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(fitpond, 3)

de donde aún puede notarse cierto tipo de patrón y además no podemos ver la “nube de puntos” con la misma longitud que nos gustaría tener. Procedemos a la prueba de hipótesis

### Nota: bptest no considera los errores ponderados, no usar lmtest::bptest(fitpond) cuando la regresión es ponderada

car::ncvTest(fitpond)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 7.047207, Df = 1, p = 0.0079389

y nuevamente se rechaza \(H_{o}\), es decir, el modelo continúa sin funcionar. Recurriendo nuevamente a la gráfica de los errores

# Cargamos la librería
library(broom)

# Calculamos los errores y los agregamos al dataframe
Datosfitpond=augment(fitpond)

# Graficamos
par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2)) 
plot(Datosfitpond$.fitted, Datosfitpond$.std.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$")   )
plot(Datosfitpond$Xprima, Datosfitpond$.std.resid, xlab = TeX("$X'$"), ylab=TeX("$e_{s}$")   )

de los cuales aún podemos ver cierto tipo de patrón en los datos. Procemos entonces a tomar a \(g\) como una funció cuadrática.

# Ajuste de los datos con g una función cuadrática
fitpond2 <- lm(salary~Xprima,weights = 1/(Xprima^2), data=Datos)

# Graficamos
par(mfrow=c(1,1)) 
par(mar=c(4, 5, 1, 1))
plot(fitpond2, 3)

y la prueba de hipótesis

car::ncvTest(fitpond2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.115539, Df = 1, p = 0.29088

no rechazándose \(H_{0}\), por consiguiente, y gracias a la gráfica, tenemos en este caso que el problema de la homocedasticidad se ha resuelto. Finalmente las gráficas de los errores

# Agregamos los errores del nuevo modelo al dataframe
Datosfitpond2=augment(fitpond2)

par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2)) 
plot(Datosfitpond2$.fitted, Datosfitpond2$.std.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$")   )
plot(Datosfitpond2$Xprima, Datosfitpond2$.std.resid, xlab = TeX("$X'$"), ylab=TeX("$e_{s}$")   )

teniendo que los puntos van por “cierta” banda, que es justo lo que buscamos.