Transformación de tipo Box-Tidwell

No transformación de la x

Continuaremos con los datos de la compañía Toluca. Para ello cargamos los datos y graficaremos el scatterplot para ver el comportamiento de los datos. Además realizaremos el ajuste del modelo y graficaremos la línea de dicho ajuste

# Cargamos las librerías necesarias:
# Datos
library(ALSM)

# Etiquetas en Látex
library(latex2exp)

# Guardamos los datos
Datos=TolucaCompany

# Realizamos el ajuste del modelo
Fit=lm(y~x, data=Datos)

# Graficamos
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$coefficients, col="red")

Después procedemos ha analizar las gráficas de los residuales verses la \(\hat{y}\) y versus \(x\). Recordemos que para el caso de loa homocedasticidad se analizan las gráficas de los residuales estandarizados contra la \(\hat{y}\) y la \(x\).

# Librería necesaria
library(broom)

# Cálculo de los errores, lo cuales se agregan a nuestro dataframe
DatosFit=augment(Fit)

# Gráficas
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$") )

Notemos que las gráficas son las mismas pues estamos en el caso de la regresión lineal simple. De las gráfica anteriores queremos no encontrar algún patrón que resalte, lo cual es indicio de que no hay evidencia en contra de la linealidad . Luego, veamos la gráfica que R trae por defecto para analizar el supuesto de la linealidad

plot(Fit, 1)

del cual nos gustaría que la línea roja fuera lo más parecida a la línea horizontal en 0, sin embargo no buscamos una coincidencia exacta.

Después, utilizamos una herramienta para argumentar la no necesidad de transformación de la \(x\) a partir de la familia de transformaciones Box-Tidwell

library(car)
boxTidwell(y~x, data=Datos)
##  MLE of lambda Score Statistic (z) Pr(>|z|)
##         1.2396              0.4932   0.6219
## 
## iterations =  2

donde el p-value de 0.6219 corresponde a la prueba de hipótesis

\[ H_{0}: \lambda=1\ \ vs \ \ H_{a}: \lambda\neq 1 \]

  • Si no se rechaza, entonces es plausible considerar el modelo con \(x\) directamente
  • Si se rechaza, entonces conviene usar el estimador de \(\lambda\) para transformar a \(x\) u otra transformación.

En nuestro ejemplo estamos no rechazando \(H_{0}\) por lo que es plausible no realizar ninguna transformación.

Transformación de la x

Comenzamos por cargar los datos, realizar el ajuste y graficar el scatterplot inicial

# Datos
datos = read.csv("ejemplo4.csv")

# Ajuste
fit <- lm(y~x, data=datos)

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

del cual claramente podemos ver que la linealidad no se está cumpliendo. De tal manera:

  1. Podemos analizar las gráficas de los errores contra \(\hat{y}\) y contra \(x\)

    # Cálculo de los errores, los cuales son agregados al dataframe
datosfit=augment(fit)

# Gráfica
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$"))

    la cual evidencia de forma contundente el patrón.

  2. Podemos analizar la gráfica que trae R por defecto para exagerar los posibles patrones

    plot(fit, 1)

De ahí que debamos buscar algún tipo de transformación para la \(x\). Para ello

boxTidwell(y~x, data=datos)
##  MLE of lambda Score Statistic (z)  Pr(>|z|)    
##         1.9776              31.325 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## iterations =  3

donde se está rechazando \(H_{0}\), es decir, se rechaza que \(\lambda=1\). Un valor de \(\lambda\) propuesto es \(\lambda=1.97\) que podemos redondear a 2. De tal manera

# Transformamos los datos
datos$Xprima=datos$x^2

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

Y checamos las gráficas de contraste con los errores

datosfit2=augment(fit2)

# Graficas
par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,2)) 
plot(datosfit2$.fitted, datosfit2$.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e$")   )
plot(datos$Xprima, datosfit2$.resid, xlab = TeX("$X^2$"), ylab=TeX("$e$"))

donde ya no podemos observar un patrón evidente. También en

plot(fit2, 1)

tampoco se ve un patrón evidente. Podemos completar utilizando la prueba de hipótesis

boxTidwell(y~Xprima, data=datos)
##  MLE of lambda Score Statistic (z) Pr(>|z|)
##        0.98879             -0.7474   0.4548
## 
## iterations =  1

con la cual no se rechaza \(H_{0}\), es decir, es plausible asumir que \(\lambda=1\)

Finalmente vemos si el supuesto de la homocedasticidad se está cumpliendo

#Homocedasticidad

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

del cual podemos ver la “nube” de puntos de igual longitud. Para estar seguros hacemos prueba de hipótesis

library(lmtest)
lmtest::bptest(fit2)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit2
## BP = 0.05566, df = 1, p-value = 0.8135

del cual no se rechaza \(H_{0}\), es decir que se puede considerar a la varianza como constante cuando la comparamos con un tipo de relación lineal con respecto a la variable Xprima. Además

summary(powerTransform(fit2))
## bcPower Transformation to Normality 
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1    1.0084           1        0.978       1.0387
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                            LRT df       pval
## LR test, lambda = (0) 403.2822  1 < 2.22e-16
## 
## Likelihood ratio test that no transformation is needed
##                             LRT df   pval
## LR test, lambda = (1) 0.2894144  1 0.5906

nos indica que no hay necesidad de transformar a la \(y\). En conclusión, los supuesto de linealidad y homocedasticidad se están cumpliendo.

Para acabar, tomemos el valor de \(\hat{\beta_{0}}\) y \(\hat{\beta_{1}}\) obtenidos del ajuste del modelo para la Xprima

fit2$coefficients
## (Intercept)      Xprima 
##   1.0567288   0.5008747

con lo cual podemos graficar la curva asociada a la recta ajustada

# Graficamos el scatterplot de los datos no transformados

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

# Línea de la recta del ajuste
abline(fit, col="red")

# ------------------------------------------------------------------------

# Función para gráficar la curva
f <- function(x){
  fit2$coefficients[1] + fit2$coefficients[2] * x ** 2
  }

# Graficamos la curva 
curve(f, from=min(datos$x), to = max(datos$x), col="purple", add=T)

legend("topleft", legend = c("curva asociada a la recta", "recta ajustada"),
       lwd = 2, col=c("purple", "red"))