Regresión lineal múltiple y transformación Box-Cox

Introducción

En este trabajo se ajustará un modelo de regresión lineal múltiple para estudiar el comportamiento de los costos de seguros médicos a partir de diferentes variables explicativas.

Después de construir el modelo se verificarán tres supuestos fundamentales asociados a los errores: normalidad, homocedasticidad e independencia. Para cada supuesto se realizará una revisión gráfica y una prueba estadística formal.

Posteriormente se aplicará la transformación Box-Cox y se volverán a evaluar los mismos supuestos. De esta manera será posible comparar el comportamiento del modelo original con el modelo después de la transformación.

Librerías

Para comenzar, cargamos las librerías necesarias para importar la base de datos, construir el modelo de regresión y realizar las diferentes pruebas estadísticas.

library(readxl)
library(nortest)
library(lmtest)
library(car)
library(MASS)

Importación de los datos

A continuación importamos la base de datos desde Excel. La información se encuentra en la hoja denominada Datos.

datos <- read_excel(
  "C:/Users/carlo/Downloads/Base_seguros_Kaggle_simple.xlsx",
  sheet = "Datos"
)

Revisión de la base de datos

Antes de construir el modelo realizamos una revisión inicial de la base. Observamos el número de filas y columnas, los nombres de las variables y la estructura de los datos.

View(datos)
dim(datos)

## [1] 1338    8

names(datos)

## [1] "i"        "age"      "sex"      "bmi"      "children" "smoker"   "region"  
## [8] "charges"

str(datos)

## tibble [1,338 × 8] (S3: tbl_df/tbl/data.frame)
##  $ i       : num [1:1338] 1 2 3 4 5 6 7 8 9 10 ...
##  $ age     : num [1:1338] 19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr [1:1338] "female" "male" "male" "male" ...
##  $ bmi     : num [1:1338] 27.9 33.8 33 22.7 28.9 ...
##  $ children: num [1:1338] 0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr [1:1338] "yes" "no" "no" "no" ...
##  $ region  : chr [1:1338] "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num [1:1338] 16885 1726 4449 21984 3867 ...

Variables

Ahora definimos las variables que se utilizarán durante el análisis. La variable respuesta será costos, mientras que edad, sexo, índice de masa corporal, número de hijos, condición de fumador y región serán utilizadas como variables explicativas.

costos <- datos$charges
edad <- datos$age
sexo <- factor(datos$sex)
imc <- datos$bmi
hijos <- datos$children
fumador <- factor(datos$smoker)
region <- factor(datos$region)

La variable i corresponde únicamente a un identificador y no se utiliza como variable regresora.

Modelo de regresión lineal múltiple

Una vez definidas las variables, ajustamos el modelo de regresión lineal múltiple.

De manera general, un modelo de regresión lineal múltiple puede expresarse como:

\[Y_i = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki} + \varepsilon_i\]

En este caso buscamos explicar los costos médicos utilizando las variables edad, sexo, índice de masa corporal, número de hijos, condición de fumador y región.

modelo <- lm(
  costos ~ edad + sexo + imc + hijos + fumador + region
)

summary(modelo)

## 
## Call:
## lm(formula = costos ~ edad + sexo + imc + hijos + fumador + region)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11304.9  -2848.1   -982.1   1393.9  29992.8 
## 
## Coefficients:
##                 Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)     -11938.5      987.8 -12.086 < 0.0000000000000002 ***
## edad               256.9       11.9  21.587 < 0.0000000000000002 ***
## sexomale          -131.3      332.9  -0.394             0.693348    
## imc                339.2       28.6  11.860 < 0.0000000000000002 ***
## hijos              475.5      137.8   3.451             0.000577 ***
## fumadoryes       23848.5      413.1  57.723 < 0.0000000000000002 ***
## regionnorthwest   -353.0      476.3  -0.741             0.458769    
## regionsoutheast  -1035.0      478.7  -2.162             0.030782 *  
## regionsouthwest   -960.0      477.9  -2.009             0.044765 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6062 on 1329 degrees of freedom
## Multiple R-squared:  0.7509, Adjusted R-squared:  0.7494 
## F-statistic: 500.8 on 8 and 1329 DF,  p-value: < 0.00000000000000022

El resumen del modelo permite observar los coeficientes estimados y los principales resultados obtenidos para cada variable incluida en la regresión.

Residuales y valores ajustados

Para verificar los supuestos del modelo necesitamos trabajar con los residuales y con los valores ajustados.

Los residuales representan la diferencia entre los valores observados y los valores estimados por el modelo. Estos serán utilizados en los análisis de normalidad, homocedasticidad e independencia.

residuales <- modelo$residuals
ajustados <- modelo$fitted.values
n <- length(residuales)

Normalidad

El primer supuesto que vamos a evaluar es el supuesto de normalidad de los errores.

Para analizarlo realizaremos primero una verificación gráfica mediante un gráfico Q-Q y posteriormente una prueba formal de Lilliefors.

Verificación gráfica

El gráfico Q-Q permite comparar los cuantiles de los residuales con los cuantiles de una distribución normal teórica. Cuando los puntos siguen aproximadamente la línea de referencia, el comportamiento de los residuales se aproxima a la normalidad.

qqnorm(
  residuales,
  main = "Verificación gráfica del supuesto de normalidad",
  xlab = "Cuantiles teóricos",
  ylab = "Cuantiles de los residuales",
  xlim = c(-4, 4),
  ylim = c(-40000, 40000),
  pch = 19
)

qqline(
  residuales,
  col = "red",
  lwd = 2
)

Prueba de Lilliefors

Después de la revisión gráfica realizamos la prueba de Lilliefors. Esta prueba permite evaluar formalmente si los residuales presentan un comportamiento compatible con una distribución normal.

prueba_normalidad <- lillie.test(residuales)
prueba_normalidad

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuales
## D = 0.16081, p-value < 0.00000000000000022

KS <- 0.895/(sqrt(n) - 0.01 + (0.85/sqrt(n)))
KS

## [1] 0.02445895

A continuación se toma la decisión utilizando un nivel de significancia de 0.05.

if(prueba_normalidad$p.value < 0.05){
  print("No se cumple el supuesto de normalidad: p-value < 0.05")
}else{
  print("No se rechaza el supuesto de normalidad: p-value >= 0.05")
}

## [1] "No se cumple el supuesto de normalidad: p-value < 0.05"

Homocedasticidad

El segundo supuesto que vamos a evaluar es la homocedasticidad.

Este supuesto busca determinar si la varianza de los errores se mantiene aproximadamente constante a lo largo de los diferentes valores ajustados del modelo.

Primero realizaremos una revisión gráfica y posteriormente utilizaremos la prueba de Breusch-Pagan.

Verificación gráfica

En este gráfico representamos los valores ajustados en el eje horizontal y los residuales en el eje vertical. La línea horizontal en cero sirve como punto de referencia para observar la dispersión de los errores.

plot(
  x = ajustados,
  y = residuales,
  main = "Verificación del supuesto de homocedasticidad",
  xlab = "Valores ajustados",
  ylab = "Residuales",
  xlim = c(0, 60000),
  ylim = c(-40000, 40000),
  col = "blue",
  pch = 19
)

abline(
  h = 0,
  lty = 3,
  col = "red",
  lwd = 2
)

Prueba de Breusch-Pagan

Después de la revisión gráfica aplicamos la prueba de Breusch-Pagan para evaluar formalmente el supuesto de varianza constante.

prueba_homocedasticidad <- bptest(modelo)
prueba_homocedasticidad

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 121.74, df = 8, p-value < 0.00000000000000022

grados_libertad_bp <- as.numeric(
  prueba_homocedasticidad$parameter
)

chi_teorico <- qchisq(
  0.05,
  grados_libertad_bp,
  lower.tail = FALSE
)

chi_teorico

## [1] 15.50731

Finalmente se toma la decisión a partir del valor p obtenido.

if(prueba_homocedasticidad$p.value < 0.05){
  print("No se cumple el supuesto de varianza constante: p-value < 0.05")
}else{
  print("No se rechaza el supuesto de varianza constante: p-value >= 0.05")
}

## [1] "No se cumple el supuesto de varianza constante: p-value < 0.05"

Independencia

El tercer supuesto que vamos a analizar es el supuesto de independencia de los errores.

Para evaluar este supuesto observaremos primero los residuales según el orden de las observaciones y posteriormente aplicaremos la prueba de Durbin-Watson.

Verificación gráfica

El siguiente gráfico permite observar si existe algún patrón evidente en la secuencia de los residuales.

orden_temporal <- 1:n

plot(
  x = orden_temporal,
  y = residuales,
  xlab = "Orden de las observaciones",
  ylab = "Residuales",
  main = "Verificación del supuesto de independencia",
  xlim = c(0, 1400),
  ylim = c(-40000, 40000),
  type = "o",
  pch = 16,
  cex = 0.40,
  col = "blue"
)

abline(
  h = 0,
  lty = 3,
  col = "red"
)

Prueba de Durbin-Watson

Después de la inspección gráfica aplicamos la prueba de Durbin-Watson para evaluar formalmente la posible presencia de autocorrelación en los errores.

prueba_independencia <- durbinWatsonTest(
  modelo,
  alternative = "two.sided"
)

prueba_independencia

##  lag Autocorrelation D-W Statistic p-value
##    1     -0.04558149      2.088423   0.104
##  Alternative hypothesis: rho != 0

La decisión se toma nuevamente utilizando un nivel de significancia de 0.05.

if(prueba_independencia$p < 0.05){
  print("No se cumple el supuesto de independencia: p-value < 0.05")
}else{
  print("No existe evidencia de autocorrelación: p-value >= 0.05")
}

## [1] "No existe evidencia de autocorrelación: p-value >= 0.05"

Transformación Box-Cox

Después de evaluar los supuestos del modelo original, aplicamos la transformación Box-Cox.

El objetivo de este procedimiento es buscar un valor de \(\lambda\) que permita transformar la variable respuesta y posteriormente volver a evaluar el comportamiento de los supuestos.

La transformación se expresa como:

\[Y(\lambda)= \begin{cases} \dfrac{Y^\lambda-1}{\lambda}, & \lambda \neq 0 \\[6pt] \ln(Y), & \lambda=0 \end{cases}\]

Antes de realizar el procedimiento se verifica que los valores de la variable costos sean positivos.

Luego se prueban diferentes valores de \(\lambda\) entre -2 y 2 y se identifica el valor asociado con el máximo obtenido por Box-Cox.

stopifnot(all(costos > 0))

resultado_boxcox <- boxcox(
  modelo,
  lambda = seq(-2, 2, by = 0.01),
  plotit = TRUE
)

lambda_optimo <- resultado_boxcox$x[
  which.max(resultado_boxcox$y)
]

print(
  paste(
    "El valor óptimo de lambda es:",
    round(lambda_optimo, 4)
  )
)

## [1] "El valor óptimo de lambda es: 0.15"

Modelo después de Box-Cox

Una vez obtenido el valor óptimo de \(\lambda\), ajustamos nuevamente el modelo utilizando la transformación Box-Cox.

Se mantienen las mismas variables explicativas utilizadas en el modelo original, lo cual permite realizar posteriormente una comparación entre ambos modelos.

modelo_boxcox <- lm(
  ((costos^lambda_optimo) - 1)/lambda_optimo ~
    edad + sexo + imc + hijos + fumador + region
)

summary(modelo_boxcox)

## 
## Call:
## lm(formula = ((costos^lambda_optimo) - 1)/lambda_optimo ~ edad + 
##     sexo + imc + hijos + fumador + region)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7134 -0.8311 -0.2622  0.1921  8.4566 
## 
## Coefficients:
##                  Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)     11.590492   0.279354  41.490 < 0.0000000000000002 ***
## edad             0.127609   0.003365  37.923 < 0.0000000000000002 ***
## sexomale        -0.247118   0.094157  -2.625              0.00878 ** 
## imc              0.060078   0.008088   7.428    0.000000000000196 ***
## hijos            0.351572   0.038971   9.021 < 0.0000000000000002 ***
## fumadoryes       6.396904   0.116839  54.750 < 0.0000000000000002 ***
## regionnorthwest -0.231656   0.134690  -1.720              0.08568 .  
## regionsoutheast -0.564758   0.135373  -4.172    0.000032173982917 ***
## regionsouthwest -0.473899   0.135159  -3.506              0.00047 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.714 on 1329 degrees of freedom
## Multiple R-squared:  0.7766, Adjusted R-squared:  0.7752 
## F-statistic: 577.3 on 8 and 1329 DF,  p-value: < 0.00000000000000022

A continuación calculamos los residuales y valores ajustados correspondientes al modelo después de Box-Cox.

residuales_boxcox <- modelo_boxcox$residuals
ajustados_boxcox <- modelo_boxcox$fitted.values
n_boxcox <- length(residuales_boxcox)

Normalidad después de Box-Cox

Después de aplicar Box-Cox volvemos a revisar el supuesto de normalidad.

El objetivo es observar si la transformación produjo cambios en el comportamiento de los residuales con respecto al modelo original.

Verificación gráfica

Nuevamente utilizamos un gráfico Q-Q para comparar los residuales del modelo transformado con los cuantiles teóricos de una distribución normal.

qqnorm(
  residuales_boxcox,
  main = "Normalidad después de la transformación Box-Cox",
  xlab = "Cuantiles teóricos",
  ylab = "Cuantiles de los residuales",
  xlim = c(-4, 4),
  ylim = c(-10, 10),
  pch = 19
)

qqline(
  residuales_boxcox,
  col = "red",
  lwd = 2
)

Pruebas de normalidad después de Box-Cox

Aplicamos nuevamente la prueba de Lilliefors y calculamos el estadístico correspondiente. También se presentan las pruebas complementarias de Anderson-Darling y Shapiro-Wilk.

prueba_normalidad_boxcox <- lillie.test(
  residuales_boxcox
)

prueba_normalidad_boxcox

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuales_boxcox
## D = 0.22135, p-value < 0.00000000000000022

KS_boxcox <- 0.895/(sqrt(n_boxcox) - 0.01 + (0.85/sqrt(n_boxcox)))
KS_boxcox

## [1] 0.02445895

ad.test(residuales_boxcox)

## 
##  Anderson-Darling normality test
## 
## data:  residuales_boxcox
## A = 85.715, p-value < 0.00000000000000022

shapiro.test(residuales_boxcox)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_boxcox
## W = 0.81403, p-value < 0.00000000000000022

Posteriormente se toma la decisión utilizando el valor p de la prueba de Lilliefors.

if(prueba_normalidad_boxcox$p.value < 0.05){
  print("Después de Box-Cox no se cumple normalidad: p-value < 0.05")
}else{
  print("Después de Box-Cox no se rechaza normalidad: p-value >= 0.05")
}

## [1] "Después de Box-Cox no se cumple normalidad: p-value < 0.05"

Homocedasticidad después de Box-Cox

Ahora volvemos a evaluar el supuesto de homocedasticidad utilizando los residuales y valores ajustados del modelo transformado.

Esto permite comparar si la dispersión de los errores cambia después de aplicar Box-Cox.

Verificación gráfica

plot(
  x = ajustados_boxcox,
  y = residuales_boxcox,
  main = "Homocedasticidad después de Box-Cox",
  xlab = "Valores ajustados",
  ylab = "Residuales",
  xlim = c(0, 40),
  ylim = c(-10, 10),
  col = "blue",
  pch = 19
)

abline(
  h = 0,
  lty = 3,
  col = "red",
  lwd = 2
)

Prueba de Breusch-Pagan después de Box-Cox

Aplicamos nuevamente la prueba de Breusch-Pagan, esta vez sobre el modelo transformado.

prueba_homocedasticidad_boxcox <- bptest(
  modelo_boxcox
)

prueba_homocedasticidad_boxcox

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_boxcox
## BP = 53.949, df = 8, p-value = 0.000000007064

grados_libertad_bp_boxcox <- as.numeric(
  prueba_homocedasticidad_boxcox$parameter
)

chi_teorico_boxcox <- qchisq(
  0.05,
  grados_libertad_bp_boxcox,
  lower.tail = FALSE
)

chi_teorico_boxcox

## [1] 15.50731

Con el valor p obtenido se determina si después de Box-Cox se mantiene o no el supuesto de varianza constante.

if(prueba_homocedasticidad_boxcox$p.value < 0.05){
  print("Después de Box-Cox no se cumple varianza constante: p-value < 0.05")
}else{
  print("Después de Box-Cox no se rechaza varianza constante: p-value >= 0.05")
}

## [1] "Después de Box-Cox no se cumple varianza constante: p-value < 0.05"

Independencia después de Box-Cox

Finalmente volvemos a analizar el supuesto de independencia de los errores para el modelo transformado.

Primero se realiza una inspección gráfica y posteriormente se vuelve a aplicar la prueba de Durbin-Watson.

Verificación gráfica

orden_temporal_boxcox <- 1:n_boxcox

plot(
  x = orden_temporal_boxcox,
  y = residuales_boxcox,
  xlab = "Orden de las observaciones",
  ylab = "Residuales",
  main = "Independencia después de Box-Cox",
  xlim = c(0, 1400),
  ylim = c(-10, 10),
  type = "o",
  pch = 16,
  cex = 0.40,
  col = "blue"
)

abline(
  h = 0,
  lty = 3,
  col = "red"
)

Prueba de Durbin-Watson después de Box-Cox

prueba_independencia_boxcox <- durbinWatsonTest(
  modelo_boxcox,
  alternative = "two.sided"
)

prueba_independencia_boxcox

##  lag Autocorrelation D-W Statistic p-value
##    1     -0.03229471      2.062544   0.224
##  Alternative hypothesis: rho != 0

A continuación se toma la decisión sobre el supuesto de independencia después de Box-Cox.

if(prueba_independencia_boxcox$p < 0.05){
  print("Después de Box-Cox no se cumple independencia: p-value < 0.05")
}else{
  print("Después de Box-Cox no existe evidencia de autocorrelación: p-value >= 0.05")
}

## [1] "Después de Box-Cox no existe evidencia de autocorrelación: p-value >= 0.05"

Comparación final

Para finalizar, construimos una tabla que resume los valores p obtenidos en cada supuesto antes y después de aplicar la transformación Box-Cox.

Esta tabla permite comparar directamente los resultados de normalidad, homocedasticidad e independencia para ambos modelos.

comparacion <- data.frame(
  Supuesto = c(
    "Normalidad",
    "Homocedasticidad",
    "Independencia"
  ),
  p_valor_original = c(
    prueba_normalidad$p.value,
    prueba_homocedasticidad$p.value,
    prueba_independencia$p
  ),
  p_valor_boxcox = c(
    prueba_normalidad_boxcox$p.value,
    prueba_homocedasticidad_boxcox$p.value,
    prueba_independencia_boxcox$p
  )
)

knitr::kable(
  comparacion,
  digits = 8,
  caption = "Comparación de los supuestos antes y después de Box-Cox"
)

Comparación de los supuestos antes y después de Box-Cox

Supuesto	p_valor_original	p_valor_boxcox
Normalidad	0.000	0.00000000
Homocedasticidad	0.000	0.00000001
Independencia	0.104	0.22400000

Cierre del análisis

Después de aplicar la transformación Box-Cox, los supuestos que mostraron mejoría fueron la homocedasticidad y la independencia. Sin embargo, la homocedasticidad aún no se cumple completamente, mientras que la independencia se mantiene adecuada. El supuesto de normalidad no presentó mejoría con la transformación.