ParteMiguel
Miguel Ángel López Ayala
2026-06-05
library(readr)
insurance <- read_csv("C:/Users/MINEDUCYT/Downloads/insurance.csv")
View(insurance, "insurance")Modelo Estimado
library(stargazer)
modelo_seguro <- lm(charges ~ age + bmi + children + smoker + region, data = insurance)
stargazer(modelo_seguro, type = "html", title = "Modelo de Regresión Lineal - Seguro de Salud")| Dependent variable: | |
| charges | |
| age | 256.974*** |
| (11.891) | |
| bmi | 338.665*** |
| (28.559) | |
| children | 474.566*** |
| (137.740) | |
| smokeryes | 23,836.300*** |
| (411.856) | |
| regionnorthwest | -352.182 |
| (476.120) | |
| regionsoutheast | -1,034.360** |
| (478.537) | |
| regionsouthwest | -959.375** |
| (477.778) | |
| Constant | -11,990.270*** |
| (978.762) | |
| Observations | 1,338 |
| R2 | 0.751 |
| Adjusted R2 | 0.750 |
| Residual Std. Error | 6,060.178 (df = 1330) |
| F Statistic | 572.697*** (df = 7; 1330) |
| Note: | p<0.1; p<0.05; p<0.01 |
library(fitdistrplus)
fit_normal<-fitdist(data = modelo_seguro$residuals,distr = "norm")
plot(fit_normal)
PRUEBAS DE NORMALIDAD
# Extración de Residuos
residuos <- residuals(modelo_seguro)
options(scipen = 99) # Evitar notación científica en los resultados
# A) Prueba de Shapiro-Wilk
# Calculo de W
salida_SW <- shapiro.test(modelo_seguro$residuals)
print(salida_SW)##
## Shapiro-Wilk normality test
##
## data: modelo_seguro$residuals
## W = 0.89909, p-value < 0.00000000000000022# Resultado final de Wn
Wn_salida<-qnorm(salida_SW$p.value,lower.tail = FALSE)
print(Wn_salida)## [1] 11.07031# B) Prueba de Jarque-Bera
library(tseries)
options(scipen = 99) # Evitar notación científica en los resultados
jb_test <- jarque.bera.test(residuos)
print(jb_test)##
## Jarque Bera Test
##
## data: residuos
## X-squared = 720.52, df = 2, p-value < 0.00000000000000022options(scipen = 99) # Evitar notación científica en los resultados
# C) Prueba de Kolmogorov-Smirnov con residuos estandarizados
library(nortest)
prueba_KS<-lillie.test(modelo_seguro$residuals)
prueba_KS##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo_seguro$residuals
## D = 0.16124, p-value < 0.00000000000000022REPRESENTACIÓN GRÁFICA
library(fastGraph)
# Parámetros del nivel de significancia y estadísticos
alpha_sig <- 0.05
JB_stat <- jb_test$statistic
gl <- jb_test$parameter
# Cálculo del Valor Crítico Teórico
VC_jb <- qchisq(1 - alpha_sig, gl, lower.tail = TRUE)
# Grafico
shadeDist(VC_jb,
ddist = "dchisq",
parm1 = gl,
lower.tail = FALSE,
xmin = 0,
xmax = 10,
sub = paste("VC:", round(VC_jb, 2), " ", "JB:", round(JB_stat, 2)))
Conclusion: El comportamiento de los residuos es idéntico: todas las pruebas estadísticas arrojan un p-valor menor a 0.05, confirmando el rechazo de la hipótesis nula de normalidad.
Prueba de Heterocedasticidad
Prueba de White
library(lmtest)
# Ejecutar la prueba de White para el nuevo modelo
PWhite <- bptest(
modelo_seguro,
~ I(age^2) + I(bmi^2) + I(children^2) +
age * bmi + age * children + age * smoker + age * region +
bmi * children + bmi * smoker + bmi * region +
children * smoker + children * region +
smoker * region,
data = insurance
)
# Mostrar los resultados
print(PWhite)##
## studentized Breusch-Pagan test
##
## data: modelo_seguro
## BP = 139.81, df = 28, p-value < 0.00000000000000022Como p-value < nivel de signficancia(0.05). Se rechaza la hipótesis nula; por lo tanto, existe evidencia estadística de la presencia de heterocedasticidad en la varianza de los residuos
library(skedastic)
options(scipen = 99) # Evitar notación científica en los resultados
skedastic::white(modelo_seguro, interactions = FALSE) ## # A tibble: 1 × 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 128. 1.92e-20 14 White's Test greaterConclusión. Como p-value < nivel de signficancia(0.05). Se rechaza la hipótesis nula; por lo tanto, existe evidencia estadística de la presencia de heterocedasticidad en la varianza de los residuos.
Prueba del multiplicador de Lagrange (Breusch-Godfrey)
Prueba de Breusch-Godfrey de segundo orden
library(lmtest)
bgtest(modelo_seguro,order = 2)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_seguro
## LM test = 5.1285, df = 2, p-value = 0.07698Conclusión. Como p value(0.07)> nivel de significancia(0.05) No se rechaza la hipotesis nula, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden “2”.
Prueba de Breusch-Godfrey de primer orden
library(lmtest)
bgtest(modelo_seguro,order = 1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelo_seguro
## LM test = 2.8251, df = 1, p-value = 0.0928Conclusión. Como p value(0.07)> nivel de significancia(0.05) No se rechaza la hipotesis nula, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden “1”.
Mas pruebas de autocorrelación.
library(lmtest)
dwtest(modelo_seguro,alternative = "two.sided",iterations = 1000)##
## Durbin-Watson test
##
## data: modelo_seguro
## DW = 2.089, p-value = 0.1035
## alternative hypothesis: true autocorrelation is not 0library(car)
durbinWatsonTest(modelo_seguro,simulate = TRUE,reps = 1000)## lag Autocorrelation D-W Statistic p-value
## 1 -0.04582739 2.088964 0.102
## Alternative hypothesis: rho != 0Estimadores HAC
options(scipen = 99999)## Warning in options(scipen = 99999): invalid 'scipen' 99999, used 9999library(lmtest)
#Modelo Sin corregir:
coeftest(modelo_seguro)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11990.270 978.762 -12.2505 < 0.00000000000000022 ***
## age 256.974 11.891 21.6101 < 0.00000000000000022 ***
## bmi 338.665 28.559 11.8584 < 0.00000000000000022 ***
## children 474.566 137.740 3.4454 0.000588 ***
## smokeryes 23836.301 411.856 57.8753 < 0.00000000000000022 ***
## regionnorthwest -352.182 476.120 -0.7397 0.459618
## regionsoutheast -1034.360 478.537 -2.1615 0.030834 *
## regionsouthwest -959.375 477.778 -2.0080 0.044846 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Modelo Corregido
options(scipen = 99999)## Warning in options(scipen = 99999): invalid 'scipen' 99999, used 9999library(lmtest)
library(sandwich)## Warning: package 'sandwich' was built under R version 4.5.3#Corregido
#HC0 Corrige Sólo Heterocedasticidad, use HC1 para corregir también Autocorrelación de Primer Orden
estimacion_omega<-vcovHC(modelo_seguro,type = "HC0")
coeftest(modelo_seguro,vcov. = estimacion_omega)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11990.27 1034.02 -11.5958 < 0.00000000000000022 ***
## age 256.97 11.86 21.6664 < 0.00000000000000022 ***
## bmi 338.67 31.53 10.7410 < 0.00000000000000022 ***
## children 474.57 129.80 3.6563 0.0002659 ***
## smokeryes 23836.30 571.14 41.7344 < 0.00000000000000022 ***
## regionnorthwest -352.18 482.84 -0.7294 0.4658903
## regionsoutheast -1034.36 499.59 -2.0704 0.0386042 *
## regionsouthwest -959.38 459.57 -2.0875 0.0370287 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1options(scipen = 99)
library(robustbase)
library(stargazer)
modelo_seguro_robust<-lmrob(charges ~ age + bmi + children + smoker + region, data = insurance)
# print(summary(modelo_seguro_robust))
stargazer(modelo_seguro,modelo_seguro_robust,type = "html",title = "comparativa")| Dependent variable: | ||
| charges | ||
| OLS | MM-type | |
| linear | ||
| (1) | (2) | |
| age | 256.974*** | 267.595*** |
| (11.891) | (1.819) | |
| bmi | 338.665*** | 10.938*** |
| (28.559) | (3.814) | |
| children | 474.566*** | 442.761*** |
| (137.740) | (17.081) | |
| smokeryes | 23,836.300*** | 33,153.790*** |
| (411.856) | (207.147) | |
| regionnorthwest | -352.182 | -182.454*** |
| (476.120) | (60.985) | |
| regionsoutheast | -1,034.360** | -544.654*** |
| (478.537) | (66.774) | |
| regionsouthwest | -959.375** | -564.381*** |
| (477.778) | (61.272) | |
| Constant | -11,990.270*** | -3,988.614*** |
| (978.762) | (132.853) | |
| Observations | 1,338 | 1,338 |
| R2 | 0.751 | 0.995 |
| Adjusted R2 | 0.750 | 0.995 |
| Residual Std. Error (df = 1330) | 6,060.178 | 893.972 |
| F Statistic | 572.697*** (df = 7; 1330) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
Pronosticos
## Pronosticos
options(scipen = 99)
library(lmtest)
library(stargazer)
library(equatiomatic) # optativo remotes::install_github("datalorax/equatiomatic")## Warning: package 'equatiomatic' was built under R version 4.5.3##
## Adjuntando el paquete: 'equatiomatic'## The following object is masked from 'package:datasets':
##
## penguins#Modelo estimado medv~. indica "medv" en función del resto de variables del dataframe
#Adaptado a tu modelo_seguro justificado con tus variables
modelo_seguro <- lm(formula = charges ~ age + bmi + children + smoker + region, data = insurance)
extract_eq(modelo_seguro,wrap = TRUE) #optativo\[ \begin{aligned} \operatorname{charges} &= \alpha + \beta_{1}(\operatorname{age}) + \beta_{2}(\operatorname{bmi}) + \beta_{3}(\operatorname{children})\ + \\ &\quad \beta_{4}(\operatorname{smoker}_{\operatorname{yes}}) + \beta_{5}(\operatorname{region}_{\operatorname{northwest}}) + \beta_{6}(\operatorname{region}_{\operatorname{southeast}}) + \beta_{7}(\operatorname{region}_{\operatorname{southwest}})\ + \\ &\quad \epsilon \end{aligned} \]
Usando Predict
library(stargazer)
options(scipen = 2)
#Data para la predicción X'm basado en un perfil de cliente para insurance
X_m<-data.frame(age=35, bmi=28.5, children=2, smoker="no", region="southwest")
# Intervalos de Confianza del 95% y del 99%
confidense<-c(0.95,0.99)
#Predicción usando predict (calculando ambos niveles para la matriz)
pred_95 <- predict(object = modelo_seguro, newdata = X_m, interval = "prediction", level = 0.95)
pred_99 <- predict(object = modelo_seguro, newdata = X_m, interval = "prediction", level = 0.99)
predicciones <- list()
predicciones$fit <- rbind(pred_95, pred_99)
rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym","Li","Ls")
stargazer(predicciones$fit,
title = "Pronósticos e intervalos de confianza",
type = "html") #Poner results='asis' en opciones del chunk| Ym | Li | Ls | |
| 95 | 6,645.506 | -5,265.556 | 18,556.570 |
| 99 | 6,645.506 | -9,016.512 | 22,307.520 |
Predicción usnado libreria forescast
library(forecast)
library(kableExtra)
options(scipen = 2)
#Data para la predicción X'm
X_m<-data.frame(age=35, bmi=28.5, children=2, smoker="no", region="southwest")
#Nivel de confianza para el intervalo de confianza
confidense<-c(0.95,0.99)
#Realizando el pronóstico con forecast
pronosticos<-forecast(object = modelo_seguro,
level = confidense,
newdata = X_m,ts = FALSE)
kable(pronosticos,
caption = "Pronóstico e intervalos de confianza:",
digits = 1,format = "html") #Poner results='asis' en opciones del chunk| Point Forecast | Lo 95 | Hi 95 | Lo 99 | Hi 99 |
|---|---|---|---|---|
| 6645.5 | -5265.6 | 18556.6 | -9016.5 | 22307.5 |
#Simulación.
#Bias Proportion
Um<-function(pronosticado,observado){
library(DescTools)
((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado)
}
#Variance Proportion
Us<-function(pronosticado,observado){
library(DescTools)
((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
#Covariance Proportion
Uc<-function(pronosticado,observado){
library(DescTools)
(2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}
#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U<-function(pronosticado,observado){
library(DescTools)
RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}
options(scipen = 99) #No mostrar notación cientifica.
library(dplyr) # Para manejo de datos y activar el operador "pipe" %>%
library(caret) # Permite Realizar muestreo sobre los data frame
library(DescTools) # Contiene las funciones para calcular las medidas de performance
library(stargazer) # Para dar formato, y obtener resumen estadistico de las simulaciones
set.seed(50) # Permite fijar la semilla aleatoria, para reproducir los resultados obtenidos en esta clase
# SOLICITADO: Simulación con 5,000 réplicas / iteraciones
numero_de_muestras<-5000
# Se crea la lista con las 5000 muestras sobre los cargos de seguros (80% entrenamiento)
muestras<- insurance$charges %>%
createDataPartition(p = 0.8,
times = numero_de_muestras,
list = TRUE)
# Listas vacias, que contendran los datos de entrenamiento, los pronosticos para los datos de prueba, y para las estadisticas de cada muestra
Modelos_Entrenamiento<-vector(mode = "list",
length = numero_de_muestras)
Pronostico_Prueba<-vector(mode = "list",
length = numero_de_muestras)
Resultados_Performance_data_entrenamiento<-vector(mode = "list",
length = numero_de_muestras)
Resultados_Performance<-vector(mode = "list",
length = numero_de_muestras)
#Estimación de los modelos lineales para cada muestra, los pronósticos y cálculo de las estadisticas de performance.
for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- insurance[muestras[[j]], ]
Datos_Prueba<- insurance[-muestras[[j]], ]
# Modelo adaptado a tu estructura justificada
Modelos_Entrenamiento[[j]]<-lm(formula = charges ~ age + bmi + children + smoker + region, data = Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Resultados_Performance_data_entrenamiento[[j]]<-data.frame(
R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges),
RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges),
MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges),
MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges)*100,
THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges,type = 1),
Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges),
Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges),
Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$charges)
)
Resultados_Performance[[j]]<-data.frame(
R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$charges),
RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$charges),
MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$charges),
MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$charges)*100,
THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$charges,
type = 1), # También se puede usar la función que creamos: THEIL_U
Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$charges),
Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$charges),
Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$charges)
)
} #No olvidar este corchete ;)bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo (Entrenamiento 80%)",
type = "html",
digits = 3)| Statistic | N | Mean | St. Dev. | Min | Max |
| R2 | 5,000 | 0.751 | 0.008 | 0.725 | 0.778 |
| RMSE | 5,000 | 6,033.374 | 73.045 | 5,735.630 | 6,256.792 |
| MAE | 5,000 | 4,166.144 | 51.019 | 3,977.420 | 4,327.530 |
| MAPE | 5,000 | 42.069 | 0.886 | 38.878 | 44.981 |
| THEIL | 5,000 | 0.173 | 0.003 | 0.163 | 0.182 |
| Um | 5,000 | 0.000 | 0.000 | 0 | 0 |
| Us | 5,000 | 0.072 | 0.003 | 0.063 | 0.080 |
| Uc | 5,000 | 0.929 | 0.003 | 0.921 | 0.938 |
bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación (Predictivo - 5,000 Réplicas)",
type = "html",
digits = 3)| Statistic | N | Mean | St. Dev. | Min | Max |
| R2 | 5,000 | 0.749 | 0.031 | 0.625 | 0.839 |
| RMSE | 5,000 | 6,092.745 | 292.925 | 5,109.142 | 7,164.056 |
| MAE | 5,000 | 4,211.507 | 156.204 | 3,699.289 | 4,826.858 |
| MAPE | 5,000 | 42.452 | 2.845 | 33.856 | 54.501 |
| THEIL | 5,000 | 0.174 | 0.010 | 0.141 | 0.215 |
| Um | 5,000 | 0.004 | 0.006 | 0.000 | 0.059 |
| Us | 5,000 | 0.078 | 0.043 | 0.00000 | 0.282 |
| Uc | 5,000 | 0.922 | 0.043 | 0.706 | 1.002 |