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.
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)
Consideraremos la prueba
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).
library(lmtest)
lmtest::bptest(fit)
##
## studentized Breusch-Pagan test
##
## data: fit
## BP = 9.8244, df = 1, p-value = 0.001722
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.
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