Librerías

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

Importación de los datos

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

Revisión de la base de 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 ...
summary(datos)
##        i               age               sex            bmi       
##  Min.   :   1.0   Min.   :18.00   Length   :1338   Min.   :15.96  
##  1st Qu.: 335.2   1st Qu.:27.00   N.unique :   2   1st Qu.:26.30  
##  Median : 669.5   Median :39.00   N.blank  :   0   Median :30.40  
##  Mean   : 669.5   Mean   :39.21   Min.nchar:   4   Mean   :30.66  
##  3rd Qu.:1003.8   3rd Qu.:51.00   Max.nchar:   6   3rd Qu.:34.69  
##  Max.   :1338.0   Max.   :64.00                    Max.   :53.13  
##     children           smoker           region        charges     
##  Min.   :0.000   Length   :1338   Length   :1338   Min.   : 1122  
##  1st Qu.:0.000   N.unique :   2   N.unique :   4   1st Qu.: 4740  
##  Median :1.000   N.blank  :   0   N.blank  :   0   Median : 9382  
##  Mean   :1.095   Min.nchar:   2   Min.nchar:   9   Mean   :13270  
##  3rd Qu.:2.000   Max.nchar:   3   Max.nchar:   9   3rd Qu.:16640  
##  Max.   :5.000                                     Max.   :63770
colSums(is.na(datos))
##        i      age      sex      bmi children   smoker   region  charges 
##        0        0        0        0        0        0        0        0

Variables

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

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

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

Parámetros estimados

coef(modelo)
##     (Intercept)            edad        sexomale             imc           hijos 
##     -11938.5386        256.8564       -131.3144        339.1935        475.5005 
##      fumadoryes regionnorthwest regionsoutheast regionsouthwest 
##      23848.5345       -352.9639      -1035.0220       -960.0510
coef(summary(modelo))
##                    Estimate Std. Error     t value
## (Intercept)     -11938.5386  987.81918 -12.0857530
## edad               256.8564   11.89885  21.5866552
## sexomale          -131.3144  332.94544  -0.3944020
## imc                339.1935   28.59947  11.8601306
## hijos              475.5005  137.80409   3.4505546
## fumadoryes       23848.5345  413.15335  57.7232020
## regionnorthwest   -352.9639  476.27579  -0.7410914
## regionsoutheast  -1035.0220  478.69221  -2.1621870
## regionsouthwest   -960.0510  477.93302  -2.0087563
##                                                                                                          Pr(>|t|)
## (Intercept)     0.00000000000000000000000000000005579044345927094398877865710772994134458713233470916748046875000
## edad            0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000007783217
## sexomale        0.69334751915996695181831910304026678204536437988281250000000000000000000000000000000000000000000
## imc             0.00000000000000000000000000000064981939262579103801578672694461147330002859234809875488281250000
## hijos           0.00057696824232809972108487750475092070701066404581069946289062500000000000000000000000000000000
## fumadoryes      0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## regionnorthwest 0.45876893258606588066328413333394564688205718994140625000000000000000000000000000000000000000000
## regionsoutheast 0.03078173928092441807846668666570622008293867111206054687500000000000000000000000000000000000000
## regionsouthwest 0.04476492951783491575090678793458209838718175888061523437500000000000000000000000000000000000000

Residuales y valores ajustados

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

Normalidad

Verificación gráfica

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

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
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

Verificación gráfica

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

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
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

Verificación gráfica

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

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

prueba_independencia
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.04558149      2.088423   0.106
##  Alternative hypothesis: rho != 0
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

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

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"

Transformación de la variable respuesta

if(abs(lambda_optimo) < 0.0001){
  costos_boxcox <- log(costos)
}else{
  costos_boxcox <- ((costos^lambda_optimo) - 1)/lambda_optimo
}

Modelo transformado

modelo_boxcox <- lm(
  costos_boxcox ~ edad + sexo + imc + hijos + fumador + region
)

summary(modelo_boxcox)
## 
## Call:
## lm(formula = costos_boxcox ~ 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
coef(modelo_boxcox)
##     (Intercept)            edad        sexomale             imc           hijos 
##     11.59049155      0.12760860     -0.24711829      0.06007766      0.35157189 
##      fumadoryes regionnorthwest regionsoutheast regionsouthwest 
##      6.39690367     -0.23165577     -0.56475795     -0.47389917
coef(summary(modelo_boxcox))
##                    Estimate  Std. Error   t value
## (Intercept)     11.59049155 0.279353886 41.490354
## edad             0.12760860 0.003364978 37.922566
## sexomale        -0.24711829 0.094156506 -2.624548
## imc              0.06007766 0.008087890  7.428100
## hijos            0.35157189 0.038970805  9.021417
## fumadoryes       6.39690367 0.116839193 54.749639
## regionnorthwest -0.23165577 0.134690129 -1.719916
## regionsoutheast -0.56475795 0.135373490 -4.171850
## regionsouthwest -0.47389917 0.135158793 -3.506240
##                                                                                                                                                                                                                                                                   Pr(>|t|)
## (Intercept)     0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004860288
## edad            0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000068635781507748265290076716826916936
## sexomale        0.00877586816222560364697535817413154290989041328430175781250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## imc             0.00000000000019627470680956672851452671224592450016643851995468139648437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## hijos           0.00000000000000000063325551131876240207665973436235162807861343026161193847656250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## fumadoryes      0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## regionnorthwest 0.08568043839034020225930987635365454480051994323730468750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## regionsoutheast 0.00003217398291721979610444451247452946063276613131165504455566406250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## regionsouthwest 0.00046961534929182733666755411583437762601533904671669006347656250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
residuales_boxcox <- modelo_boxcox$residuals
ajustados_boxcox <- modelo_boxcox$fitted.values
n_boxcox <- length(residuales_boxcox)

Normalidad después de Box-Cox

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
)

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
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

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_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
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

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_independencia_boxcox <- durbinWatsonTest(
  modelo_boxcox,
  alternative = "two.sided"
)

prueba_independencia_boxcox
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.03229471      2.062544   0.244
##  Alternative hypothesis: rho != 0
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

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 = 6,
  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.000
Homocedasticidad 0.000 0.000
Independencia 0.106 0.244