Métodos estadísticos I

Tarea 10 - Transformación Box-Cox

En este ejercicio replicaremos la transformación Box-Cox, que se utiliza para transformar la variable de regresión cuando no se cumple el supuesto de homocedasticidad. La transformación es la siguiente:

\[ z_i=\left\{\begin{matrix} \frac{y_i^\lambda -1}{\lambda \omega^{\lambda-1}} & & \lambda\neq 0\\ \\ \omega\ln(y_i)& &\lambda=0 \end{matrix}\right. \]

Donde: \( \omega=\exp\left[\frac{1}{n}\sum_{i=1}^n \ln(y_i)\right] \)

Se tienen que elegir una serie de \( \lambda \)'s propuestas, y la idea es optimizar la regresión. Se diseño la siguiente función para este propósito

nueva_box_cox <- function(modelo, lambda = seq(-2, 2, 1/10)) {
    x.reg <- as.matrix(modelo$model[, 2:ncol(modelo$model)])
    w <- exp(sum(log(modelo$model[, 1]))/length(modelo$model[, 1]))
    ssr <- rep(NA, length(lambda))
    for (i in 1:length(lambda)) {
        if (lambda[i] == 0) {
            yl <- w * log(modelo$model[, 1])
            modl <- lm(yl ~ x.reg)
            ssr[i] <- sum(modl$residuals^2)
        } else {
            yl <- ((modelo$model[, 1]^lambda[i]) - 1)/(lambda[i] * w^(lambda[i] - 
                1))
            modl <- lm(yl ~ x.reg)
            ssr[i] <- sum(modl$residuals^2)
        }
    }
    lam_op <- lambda[ssr == min(ssr)]
    plot(lambda, ssr, type = "l", xlab = expression(lambda), ylab = "Suma de cuadrados del error")
    abline(v = lam_op, lty = 2, col = 2)
    legend("top", bty = "n", legend = substitute("lambda" == lam_op))
}

Notese que se puede ingresas valores de lambda deseados, pero la función ya trae algunos por default.

Para evaluar la función propuesta se utilizan los siguientes datos de prueba:

x1 <- c(35, 35, 35, 35, 35, 35, 35, 35, 55, 55, 55, 55, 55, 55, 55, 55, 25, 
    65, 45, 45, 45, 45, 45, 45, 45, 45)
x2 <- c(0.3, 0.3, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.3, 0.3, 0.3, 0.3, 0.7, 0.7, 
    0.7, 0.7, 0.5, 0.5, 0.1, 0.9, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)
y1 <- c(76.5, 76, 79.9, 83.5, 89.5, 84.2, 85.7, 99.5, 89.4, 97.5, 103.2, 108.7, 
    115.2, 111.5, 102.3, 108.1, 80.2, 89.1, 77.2, 85.1, 71.5, 84.5, 77.5, 79.2, 
    71, 90.2)

Compararemos la función diseñada con la función en el paquete “MASS”, la propuesta es la primer gráfica y la segunda la que viene en “R”:

par(mfrow = c(1, 2))
library(MASS)
nueva_box_cox(lm(y1 ~ x1 + x2))
boxcox(lm(y1 ~ x1 + x2), lambda = seq(-2, 2, 1/10), plotit = TRUE)

Notemos que las gráficas no son iguales, pero el resultado es el mismo, lo que sucede es que en la función que diseñamos utilizamos para encontrar el valor óptimo de \( \lambda \) el minimizar la suma de cuadrados del error de la regresión, mientras que R utiliza el valor de \( \lambda \) que maximiza la log-verosimilitud, sin embargo dado que los estimadores de mínimos cuadrados de la regresión lineal son estimadores de máxima verosimilitud, el principio es el mismo, solo que la gráfica sale volteada.