Transformaciones de tipo Box-Cox

Trabajaremos con los datos de una universidad (empresa) y los costos de la colegiatura. Para ello, los datos los importaremos de la siguiente librería

# Librería con los datos
library(Ecdat)

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

# almacenamos los datos
Datos <- University

donde la variable \(y\) está denominada por nassets (posible costo económico que le daríamos a la universidad como empresa) y la variable \(x\) por stfees (costo de la colegiatura), de donde veremos si la univeridad como empresa y cuán atractiva es determina el costo de las colegiaturas.

Análisis de los supuestos

comencemos por ver el scatterplot de los datos

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

de donde el supuesto de la linealidad claramente no se cumple. Podemos utilizar el scatterplot por defecto que trae R para hacer más notorio lo anterior

# Realizamos el ajuste
fit=lm(nassets~stfees, data=Datos)

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

y también la gráfica para revisar la homocedasticidad

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

Con lo cual podemos intuir que ambos supuestos no se están cumpliendo, siendo necesaria transformar a nuestra variables \(y\). Recordemos que la gráfica que nos brinda R es una especie de exageración en cuanto al patrón que puede presentarse. Podemos hallar una gráfica para ver si el supuesto de la homcedasticidad se cumple usando la toría ya trabajada, donde utilizaremos la librería broom para el cálculo de los errores residuales estandarizados

# Librería para obtener los errores estandarizados
library(broom)

# Nos agrega una columna con los errores (resid) y los errores estandarizados (std.resid)
Datosfit=augment(fit)

# Graficamos nuevamente con estos errores obtenidos
par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,1)) 
plot(Datosfit$.fitted, Datosfit$.std.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$")   )

del cual, nuevamente, se puede ver que el supuesto de la homocedasticidad no se está cumpliendo.

Comparemos ambas gráficas obtenidas, una con los errores estandarizados y otra con la gráfica que R tiene implementada

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

# Gráfica con los errores estandarizados 
plot(Datosfit$.fitted, Datosfit$.std.resid, xlab = TeX("$\\widehat{y}$"), ylab=TeX("$e_{s}$"))

# Gráfica que R trae por defecto
plot(fit, 3)

Pruebas de hipótesis

Consideraremos la prueba

  • H0: varianza no depende de forma lineal en x vs Ha: varianza depende de forma lineal en x.

de donde se busca no rechazar, es decir, que sea plausible asumir que la varianza no depende de forma lineal de \(x\) (i.e. p-value mayor a significancia).

Prueba de hipotesis con los errores studentizados

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

Prueba de hipotesis con los residuales estandarizados

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

Siendo así que en ambas pruebas rechazamos \(H_{0}\) por lo que es más evidente el problema de la homocedasticidad.

Transformación de los datos

Primero se trabajará cuando la variable “y” es estrictamente positiva. Recordemos que se prefiere un valor de \(\lambda\) conocido para no complicar mucho la interpretación de los datos. Luego

summary(powerTransform(fit))
## bcPower Transformation to Normality 
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1    0.1894        0.33       0.0483       0.3305
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                            LRT df    pval
## LR test, lambda = (0) 6.473425  1 0.01095
## 
## Likelihood ratio test that no transformation is needed
##                            LRT df       pval
## LR test, lambda = (1) 107.1042  1 < 2.22e-16

El intervalo sugiere que \(\lambda=.2\) podría ser una opción (notar que se rechaza \(\lambda=0\) y \(\lambda=1\)). Después transformamos los datos

Datos$nassets_lambda0=bcPower(Datos$nassets, .2)
Datos$nassetsBC=(Datos$nassets^(.2)-1)/.2

y graficamos

# Realizamos el nuevo ajuste con los datos transformados
fit2=lm(nassetsBC~stfees, data=Datos)

# Realizamos la gráfica
par(mar=c(4, 5, 1, 1))
par(mfrow=c(1,1)) 
plot(Datos$stfees, Datos$nassetsBC, xlab = TeX("$stfees$"), ylab=TeX("$(y^{.2}-1)/.2$"))

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

Realizamos nuevamente la prueba de homocedasticidad con la gráfica de R

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

y realizamos la prueba de hipótesis

lmtest::bptest(fit2)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit2
## BP = 0.16697, df = 1, p-value = 0.6828
car::ncvTest(fit2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.2670924, Df = 1, p = 0.60529

con lo cual conseguimos no rechazar \(H_{0}\). Así, en conjunto de la gráfica y la prueba de hipótesis podemos suponer que el supuesto de la homocedasticidad ya se cumple para estos datos transformados.

Veamos el supuestos de la linealidad

plot(fit2, 1)

el cual parece tener un mejor ajuste a la recta horizontal en 0, siendo así que este supuesto parece cumplirse en este modelo ajustado de los datos transformados.

Para datos negativos o con ceros se puede usar una constante \(\gamma\) positiva para trabajar con valores positivos.

#powerTransform(,  family="bcnPower")
summary(powerTransform(fit,  family="bcnPower"))
## bcnPower transformation to Normality 
## 
## Estimated power, lambda
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1   -0.0558           0      -0.9142       0.8026
## 
## Estimated location, gamma
##    Est gamma Std Err. Wald Lower Bound Wald Upper Bound
## Y1  16677.21 26052.79                0         67740.68
## 
## Likelihood ratio tests about transformation parameters
##                               LRT df      pval
## LR test, lambda = (0)   0.5455977  1 0.4601221
## LR test, lambda = (1) 122.5897026  1 0.0000000

Lo cual inidca que podríamos tomar \(\lambda=0\) y \(\gamma=16677.21\). De ahí que

# Transformamos los datos
Datos$nassetsBCalt=bcnPower(Datos$nassets, lambda=0, gamma=16677.21)

# Ajustamos el modelo a los datos transformados
fit3=lm(nassetsBCalt~stfees, data=Datos)

# Graficamos
plot(Datos$stfees, Datos$nassetsBCalt, xlab = TeX("$stfees$"), ylab=TeX("$ln(z)$") )
abline(fit3$coefficients, col="red")

y también

plot(fit3, 3)

finalmente la prueba de hipotesis

lmtest::bptest(fit3)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit3
## BP = 0.31854, df = 1, p-value = 0.5725
car::ncvTest(fit3)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.4741867, Df = 1, p = 0.49107

Una opción es transformar los datos sumándole la constante positiva antes de usar la transformación BoxCox simple