# Cargar librerías necesarias
library(MASS)
library(ggplot2)
library(mosaic)
library(knitr)
library(kableExtra)
library(modeest) # Necesario para la función mfv()
library(e1071)
library(nortest)
library(lmtest)
library(alr4)4. Transformaciones de BoxCox
1 Introducción
En el análisis de regresión lineal simple, el modelo se asume de la forma: \[Y = \beta_0 + \beta_1 X + \epsilon\] Para que las inferencias (pruebas de hipótesis, intervalos de confianza y predicciones) sean válidas, los errores (\(\epsilon_i\)) deben cumplir con los supuestos fundamentales: normalidad, homocedasticidad (varianza constante) e independencia.
Cuando estos supuestos no se cumplen, frecuentemente evidenciados por curvaturas en los gráficos de dispersión o embudos en los residuales, la transformación de Box-Cox ofrece una herramienta formal de máxima verosimilitud para encontrar la potencia óptima \(\lambda\) en la variable respuesta \(Y\) que restaure los supuestos del modelo.
2 Transformación de BoxCox
Sea \(y_i\) la variable de respuesta positiva. La transformación de Box-Cox se define como:
\[\begin{equation} y_{i}^{(\lambda)} = \begin{cases} \frac{y_{i}^{\lambda} - 1}{\lambda}, & \text{si } \lambda \neq 0 \\ \ln(y_i), & \text{si } \lambda = 0 \end{cases} \end{equation}\]
La transformación Box-Cox proporciona un método formal para encontrar el valor de \(\lambda\) que mejor normaliza los datos (solo funciona con datos positivos \(y_i>0\)).
El valor de \(\lambda\) se determina mediante estimación de máxima verosimilitud del modelo transformado para encontrar el valor que aproxime lo máximo posible la distribución normal de los datos.
\[L^{(\lambda)}=-\frac{n}{2}ln(\hat \sigma^2_{(\lambda)}) + (\lambda-1)\displaystyle \sum_{i=1}^nln(y_i)\]
El valor de \(\lambda\) encontrado dice qué transformación se está aplicando. A menudo se redondea al valor “común” más cercano para facilitar la interpretación.
\[\lambda \in [-5,5]\]
2.1 Transformaciones comúnes
| \(\lambda\) | Transformación |
|---|---|
| -3 | \[\frac{1}{Y^3}\] |
| -2 | \[\frac{1}{Y^2}\] |
| -1 | \[\frac{1}{Y}\] |
| -0.5 | \[\frac{1}{\sqrt Y}\] |
| 0 | \[ln(Y)\] |
| 0.5 | |
| 1 | \[\sqrt Y\] |
| 2 | \[Y^2\] |
| 3 | \[Y^3\] |
Box-Cox no es robusto ante valores atípicos.
Se aplica principalmente a la transformación de la variable respuesta \(Y\).
Todos los valores de la variable a transformar deberán ser positivos.
Después de aplicar el método es importante verificar de nuevo los supuestos del modelo.
No usar el método si se cumplen los supuestos del modelo.
Es posible encontrar un IC aproximado para el parámetro de transformación \(\lambda\). Este IC puede ser útil para seleccionar el valor final de \(\lambda\).
Un intervalo de confianza aproximado del \(100(1−α)\%\) para \(\lambda\) se define como el conjunto de valores que cumplen:
\[L(\hat \lambda)-L(\lambda)\leq \frac{1}{2} \chi^2_{\alpha,1}\]
\(\lambda\) El parámetro de la transformación de potencia de Box-Cox.
L(𝜆) Log-Verosimilitud.
\(\hat\lambda\) Estimación puntual (EMV).
\(L_{max}(\lambda)-\frac{1}{2} \chi^2_{\alpha,1}\) Es el umbral para construir un intervalo de confianza.
\(\lambda_1\) Es el límite inferior para el Intervalo de Confianza.
\(\lambda_2\) Es el límite superior para el Intervalo de Confianza.
A pesar del valor que minimiza la suma de cuadrados de los residuos es \(\hat \lambda\), los valores que se encuentran entre \(𝜆_1\) y \(𝜆_2\) podrían utilizarse de forma eficiente para realizar una transformación.
Un valor de \(\lambda\) de 1 es equivalente a utilizar los datos originales. Por lo tanto, si el intervalo de confianza de λ óptimo incluye 1, entonces no es necesaria una transformación.
Si incluye 0, el logaritmo natural puede considerarse una buena transformación.
Si está “alejado” de 0 y 1, los datos necesitan una transformación de potencia específica.
3 Ejemplo 1
Estaremos siguiendo los siguientes pasos para obtener un modelo optimo en los dos ejemplos:
1. Descripción de los datos.
2. Determinar modelo inicial.
3. Revisión de supuestos.
4. Determinar la lambda estimada óptima para la transformación Box-Cox.
5. Determinar modelo con la transformación.
6. Revisión de supuestos.
7. Comparación entre modelos. Nota: Para todas las pruebas de hipótesis se usará un nivel de significancia del 0.05
Para este laboratorio utilizaremos R y el paquete MASS, que incluye la función boxcox().
x_kwh <- c(679, 292, 1012, 493, 582, 1156, 997, 2189, 1097, 2078, 1818, 1700, 747, 2030, 1643, 414, 354, 1276, 745, 435, 540, 874, 1543, 1029, 710, 1434,
837, 1748, 1381, 1428, 1255, 1777, 370, 2316, 1130, 463, 770, 724, 808, 790, 783, 406, 1242, 658, 1746, 468, 1114, 413, 1787, 3560, 1495, 2221, 1526)
y_kw <- c(0.79, 0.44, 0.56, 0.79, 2.70, 3.64, 4.73, 9.50, 5.34, 6.85, 5.84, 5.21, 3.25, 4.43, 3.16, 0.50, 0.17, 1.88, 0.77, 1.39, 0.56, 1.56, 5.28, 0.64, 4.00, 0.31,
4.20, 4.88, 3.48, 7.58, 2.63, 4.99, 0.59, 8.19, 4.79, 0.51, 1.74, 4.10, 3.94, 0.96, 3.29, 0.44, 3.24, 2.14, 5.71, 0.64, 1.90, 0.51, 8.33, 14.94, 5.11, 3.85, 3.93)
# Crear el data frame
datos_electricos <- data.frame(X = x_kwh, Y = y_kw)
# Primeras observaciones
head(datos_electricos) X Y
1 679 0.79
2 292 0.44
3 1012 0.56
4 493 0.79
5 582 2.70
6 1156 3.64modelo_original <- lm(Y ~ X, data = datos_electricos)
summary(modelo_original)
Call:
lm(formula = Y ~ X, data = datos_electricos)
Residuals:
Min 1Q Median 3Q Max
-4.1399 -0.8275 -0.1934 1.2376 3.1522
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.8313037 0.4416121 -1.882 0.0655 .
X 0.0036828 0.0003339 11.030 4.11e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.577 on 51 degrees of freedom
Multiple R-squared: 0.7046, Adjusted R-squared: 0.6988
F-statistic: 121.7 on 1 and 51 DF, p-value: 4.106e-15## Diagnóstico inicial
plot(y_kw ~ x_kwh, data = datos_electricos)
plot(modelo_original)
Observa la dispersión de los residuales y busca curvatura, patrones sistemáticos y posibles observaciones influyentes.
3.1 Supuesto de Normalidad de los residuales
3.1.1 Planteamiento de las Hipótesis
Para cualquier prueba de normalidad sobre los residuales del modelo (\(\varepsilon_i\)):
\(\text{H}_0\) (Hipótesis Nula): Los residuales del modelo siguen una distribución normal (\(\varepsilon_i \sim N(0, \sigma^2)\)).
\(\text{H}_1\) (Hipótesis Alternativa): Los residuales del modelo no siguen una distribución normal.
shapiro.test(modelo_original$residuals)
Shapiro-Wilk normality test
data: modelo_original$residuals
W = 0.97435, p-value = 0.3087Como \(p\text{-valor}=0.3087 > 0.05\) no se rechaza \(\text{H}_0\). No existen evidencias de que los residuales no se distribuyen de forma normal.
3.2 Supuesto de Homocedasticidad
3.2.1 Planteamiento de las Hipótesis
\(\text{H}_0\) (Hipótesis Nula): Los residuales son homocedásticos (la varianza de los errores es constante, \(\sigma_i^2 = \sigma^2\)).
\(\text{H}_1\) (Hipótesis Alternativa): Los residuales son heterocedásticos (la varianza de los errores no es constante).
# 2. Prueba formal de Breusch-Pagan
prueba_bp <- bptest(modelo_original)
print(prueba_bp)
studentized Breusch-Pagan test
data: modelo_original
BP = 6.8012, df = 1, p-value = 0.00911Como \(p\text{-valor}=0.0091 < 0.05\) se rechaza \(\text{H}_0\). Los residuales son heterocedásticos (la varianza de los errores no es constante).
3.3 Supuesto de Independencia
3.3.1 Planteamiento de las Hipótesis
Para la prueba de Durbin-Watson:
\(\text{H}_0\) (Hipótesis Nula): \(\rho = 0\) (Los errores son independientes / No existe autocorrelación en los residuales).
\(\text{H}_1\) (Hipótesis Alternativa): \(\rho \neq 0\) (Existe autocorrelación en los residuales).
# 2. Prueba de Durbin-Watson
prueba_dw <- dwtest(modelo_original)
print(prueba_dw)
Durbin-Watson test
data: modelo_original
DW = 2.0946, p-value = 0.6262
alternative hypothesis: true autocorrelation is greater than 0# 3. Gráfico de autocorrelación (ACF) y residuales en función del orden
datos_indep <- data.frame(
Orden = 1:nrow(datos_electricos),
Residuales = residuals(modelo_original)
)
# Gráfico de residuales seguidos por orden de observación
ggplot(datos_indep, aes(x = Orden, y = Residuales)) +
geom_point(color = "#2b5c8f", size = 2.5) +
geom_line(color = "#2b5c8f", linetype = "dashed", alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "#d9534f", linewidth = 1) +
labs(
title = "Residuales por Orden de Observación",
subtitle = "Evaluación de Independencia de Errores",
x = "Orden de la Observación",
y = "Residuales"
) +
theme_minimal()
Como Si \(p\text{-valor} > 0.05\): No se rechaza \(\text{H}_0\). Se cumple el supuesto de independencia.
3.4 Transformación de CoxBox
# Calcular Box-Cox
bc <- boxcox(modelo_original, lambda = seq(-2, 2, by = 0.1))
# Encontrar el lambda óptimo exacto
lambda_optimo <- bc$x[which.max(bc$y)]
cat("El valor óptimo estimado de lambda es:", lambda_optimo, "\n")El valor óptimo estimado de lambda es: 0.5454545 Aplicaremos la transformación sugerida (raíz cuadrada):
# Aplicamos transformación raíz cuadrada (lambda = 0.5) o logarítmica según el resultado de Box-Cox
# Usaremos transformación logarítmica si lambda está cerca de 0, o raíz cuadrada si está cerca de 0.5
# Para este caso de estudio eléctrico, probemos con log(Y) o raíz de Y:
datos_electricos$Y_trans <- (datos_electricos$Y)^0.5
modelo_trans <- lm(Y_trans ~ X, data = datos_electricos)
summary(modelo_trans)
Call:
lm(formula = Y_trans ~ X, data = datos_electricos)
Residuals:
Min 1Q Median 3Q Max
-1.39185 -0.30576 -0.03875 0.25378 0.81027
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.822e-01 1.299e-01 4.481 4.22e-05 ***
X 9.529e-04 9.824e-05 9.699 3.61e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.464 on 51 degrees of freedom
Multiple R-squared: 0.6485, Adjusted R-squared: 0.6416
F-statistic: 94.08 on 1 and 51 DF, p-value: 3.614e-13shapiro.test(modelo_trans$residuals)
Shapiro-Wilk normality test
data: modelo_trans$residuals
W = 0.96602, p-value = 0.1353prueba_bp <- bptest(modelo_trans)
print(prueba_bp)
studentized Breusch-Pagan test
data: modelo_trans
BP = 0.012597, df = 1, p-value = 0.9106prueba_dw <- dwtest(modelo_trans)
print(prueba_dw)
Durbin-Watson test
data: modelo_trans
DW = 2.3042, p-value = 0.8642
alternative hypothesis: true autocorrelation is greater than 0## Diagnóstico inicial
plot(y_kw ~ x_kwh, data = datos_electricos)
plot(Y_trans ~ x_kwh, data = datos_electricos)
plot(modelo_trans)
3.5 Detección de Outliers
n <- nrow(datos_electricos)
p <- 2 # Intercepto y pendiente
# 1. Leverage
hat_vals <- hatvalues(modelo_trans)
umbral_hat <- 3 * (p / n)
# 2. Distancia de Cook
cook_vals <- cooks.distance(modelo_trans)
umbral_cook <- 4 / (n - p)
# 3. RStudent
r_student <- rstudent(modelo_trans)
# Identificar filas problemáticas
atipicos_rstudent <- which(abs(r_student) > 3)
atipicos_cook <- which(cook_vals > umbral_cook)
cat("Observaciones con RStudent > 3:", atipicos_rstudent, "\n")Observaciones con RStudent > 3: 26 plot(modelo_trans, which = 4, col = "blue", pch = 16)
abline(h = umbral_cook, col = "red", lty = 2)
# Identificamos y filtramo (la de mayorCook)
electric_filtrado <- datos_electricos[-atipicos_rstudent, ]
electric_filtrado X Y Y_trans
1 679 0.79 0.8888194
2 292 0.44 0.6633250
3 1012 0.56 0.7483315
4 493 0.79 0.8888194
5 582 2.70 1.6431677
6 1156 3.64 1.9078784
7 997 4.73 2.1748563
8 2189 9.50 3.0822070
9 1097 5.34 2.3108440
10 2078 6.85 2.6172505
11 1818 5.84 2.4166092
12 1700 5.21 2.2825424
13 747 3.25 1.8027756
14 2030 4.43 2.1047565
15 1643 3.16 1.7776389
16 414 0.50 0.7071068
17 354 0.17 0.4123106
18 1276 1.88 1.3711309
19 745 0.77 0.8774964
20 435 1.39 1.1789826
21 540 0.56 0.7483315
22 874 1.56 1.2489996
23 1543 5.28 2.2978251
24 1029 0.64 0.8000000
25 710 4.00 2.0000000
27 837 4.20 2.0493902
28 1748 4.88 2.2090722
29 1381 3.48 1.8654758
30 1428 7.58 2.7531800
31 1255 2.63 1.6217275
32 1777 4.99 2.2338308
33 370 0.59 0.7681146
34 2316 8.19 2.8618176
35 1130 4.79 2.1886069
36 463 0.51 0.7141428
37 770 1.74 1.3190906
38 724 4.10 2.0248457
39 808 3.94 1.9849433
40 790 0.96 0.9797959
41 783 3.29 1.8138357
42 406 0.44 0.6633250
43 1242 3.24 1.8000000
44 658 2.14 1.4628739
45 1746 5.71 2.3895606
46 468 0.64 0.8000000
47 1114 1.90 1.3784049
48 413 0.51 0.7141428
49 1787 8.33 2.8861739
50 3560 14.94 3.8652296
51 1495 5.11 2.2605309
52 2221 3.85 1.9621417
53 1526 3.93 1.9824228# Reajustar modelo sin el outlier
modelo_filtrado <- lm(Y_trans ~ X, data = electric_filtrado)
summary(modelo_filtrado)
Call:
lm(formula = Y_trans ~ X, data = electric_filtrado)
Residuals:
Min 1Q Median 3Q Max
-0.82253 -0.29841 -0.06996 0.26894 0.77846
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.884e-01 1.188e-01 4.952 8.76e-06 ***
X 9.708e-04 8.998e-05 10.789 1.17e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4243 on 50 degrees of freedom
Multiple R-squared: 0.6995, Adjusted R-squared: 0.6935
F-statistic: 116.4 on 1 and 50 DF, p-value: 1.174e-14shapiro.test(modelo_filtrado$residuals)
Shapiro-Wilk normality test
data: modelo_filtrado$residuals
W = 0.96105, p-value = 0.08691prueba_bp <- bptest(modelo_filtrado)
print(prueba_bp)
studentized Breusch-Pagan test
data: modelo_filtrado
BP = 0.43153, df = 1, p-value = 0.5112prueba_dw <- dwtest(modelo_filtrado)
print(prueba_dw)
Durbin-Watson test
data: modelo_filtrado
DW = 1.854, p-value = 0.2869
alternative hypothesis: true autocorrelation is greater than 0