Evidencia Empirica

library(readr)
insurance <- read_csv("insurance.csv")
head(insurance, 10)
## # A tibble: 10 × 7
##      age sex      bmi children smoker region    charges
##    <dbl> <chr>  <dbl>    <dbl> <chr>  <chr>       <dbl>
##  1    19 female  27.9        0 yes    southwest  16885.
##  2    18 male    33.8        1 no     southeast   1726.
##  3    28 male    33          3 no     southeast   4449.
##  4    33 male    22.7        0 no     northwest  21984.
##  5    32 male    28.9        0 no     northwest   3867.
##  6    31 female  25.7        0 no     southeast   3757.
##  7    46 female  33.4        1 no     southeast   8241.
##  8    37 female  27.7        3 no     northwest   7282.
##  9    37 male    29.8        2 no     northeast   6406.
## 10    60 female  25.8        0 no     northwest  28923.

Modelo Estimado

library(lmtest)
## Cargando paquete requerido: zoo
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(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_seguro <- lm(log(charges) ~ age + children + smoker * bmi, data = insurance)
extract_eq(modelo_seguro, wrap = TRUE)

\[ \begin{aligned} \operatorname{log(charges)} &= \alpha + \beta_{1}(\operatorname{age}) + \beta_{2}(\operatorname{children}) + \beta_{3}(\operatorname{smoker}_{\operatorname{yes}})\ + \\ &\quad \beta_{4}(\operatorname{bmi}) + \beta_{5}(\operatorname{smoker}_{\operatorname{yes}} \times \operatorname{bmi}) + \epsilon \end{aligned} \]

library(stargazer)
modelo_seguro <- lm(log(charges) ~ age + children + smoker * bmi, data = insurance)
stargazer(modelo_seguro, type = "html", title = "Modelo de Regresión Lineal - Seguro de Salud")
Modelo de Regresión Lineal - Seguro de Salud
Dependent variable:
log(charges)
age 0.035***
(0.001)
children 0.102***
(0.010)
smokeryes 0.186
(0.148)
bmi 0.001
(0.002)
smokeryes:bmi 0.044***
(0.005)
Constant 7.272***
(0.074)
Observations 1,338
R2 0.777
Adjusted R2 0.776
Residual Std. Error 0.435 (df = 1332)
F Statistic 927.462*** (df = 5; 1332)
Note: p<0.1; p<0.05; p<0.01

Verificación de los supuestos del Modelo Clásico de Regresión Lineal Múltiple

Verifiación de los ajuste de los residuos

library(fitdistrplus)
fit_normal<-fitdist(data = modelo_seguro$residuals,distr = "norm")
plot(fit_normal)

Verificación de la Normalidad

\(H_0 : \varepsilon \sim N(0, \sigma^2)\) “Los Residuos siguen una distribución normal con media cero y varianza constante”

\(H_A : \varepsilon \not\sim N(0, \sigma^2)\) “Los Residuos no siguen una distribución normal con media cero y varianza constante”

Prueba de Shapiro-Wilk

# 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.81858, p-value < 0.00000000000000022
# Resultado final de Wn
Wn_salida<-qnorm(salida_SW$p.value,lower.tail = FALSE)
print(Wn_salida)
## [1] 12.54004

Prueba de Jaque-Bera

# 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 = 1973.8, df = 2, p-value < 0.00000000000000022

Prueba de Kolmogorov-Smirnov con residuos estandarizados

options(scipen = 99) # Evitar notación científica en los resultados
# C) Prueba de Kolmogorov-Smirnov
library(nortest)
prueba_KS<-lillie.test(modelo_seguro$residuals)
prueba_KS
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_seguro$residuals
## D = 0.24721, p-value < 0.00000000000000022

Resultado en Forma Grafica

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

## Verificación de Multicolinealidad ### Indice de Condición

\(\kappa(\mathbf{x}) = \sqrt{\frac{\lambda_{max}}{\lambda_{min}}}\)

Interpretación del Índice de Condición:

Si \(\kappa(\mathbf{x}) \leq 20\), la multicolinealidad es leve y no representa un problema significativo.

Cuando \(20 < \kappa(\mathbf{x}) < 30\), se considera que la multicolinealidad es de nivel moderado.

Si \(\kappa(\mathbf{x}) \geq 30\), la multicolinealidad es severa y puede comprometer los resultados del modelo.

# OPCION 1 
#  Cálculo del Indice de Condición usando librería “mctest”
library(mctest)
X_mat<-model.matrix(modelo_seguro)
mctest(mod = modelo_seguro)
## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0387         0
## Farrar Chi-Square:      4338.8326         1
## Red Indicator:             0.3125         0
## Sum of Lambda Inverse:    53.9039         1
## Theil's Method:           -0.9406         0
## Condition Number:         19.7272         0
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

Factores Inflacionarios de la Varianza (FIV)

VIF Interpretación

Igual a 1 - Sin multicolinealidad

Menor a 5 - Baja o moderada, generalmente aceptable

Menor a 10 - Posible problema de multicolinealidad

Mayor a 10. Problema serio de multicolinealidad

Nota: El criterio puede cambiar dependiendo del investigador.

# CÁLCULO DE LAS VIFs usando performance
library(performance)
VIFs<-multicollinearity(x = modelo_seguro,verbose = FALSE)
VIFs
## # Check for Multicollinearity
## 
## Low Correlation
## 
##      Term  VIF        VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##       age 1.02 [ 1.00,     1.56]     1.01      0.98     [0.64, 1.00]
##  children 1.00 [ 1.00, 2.72e+08]     1.00      1.00     [0.00, 1.00]
##       bmi 1.30 [ 1.22,     1.39]     1.14      0.77     [0.72, 0.82]
## 
## High Correlation
## 
##        Term   VIF        VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##      smoker 25.14 [22.66,    27.91]     5.01      0.04     [0.04, 0.04]
##  smoker:bmi 25.45 [22.93,    28.25]     5.04      0.04     [0.04, 0.04]

Resultados en Forma Gráfica

# CÁLCULO DE LAS VIFs usando mctest
library(mctest)
mc.plot(mod = modelo_seguro,vif = 2)

Verificación de la Homocedasticidad

\(\hat{\delta}^2_u = \delta_0 + \delta_1 x_1 + \delta_2 x_2 + \ldots + \delta_k x_k + \alpha_2 x_1^2 + \alpha_2 x_2^2 \ldots + \alpha_k x_k^2 + \theta_1 x_1 x_2 + \theta_1 x_1 x_3 + \theta_k C2 \alpha_i x_j\)

Hipótesis

\(H_0\): \(\delta_1 = \delta_2 = \ldots = \delta_k = \alpha_1 = \alpha_2 = \ldots = \alpha_k = \theta_1 = \theta_2 = \ldots = \theta_k C_2 = 0\)

Existe evidencia estadística de que los residuos son homocedásticos (la varianza del error es constante).

\(H_1\): \(\delta_1 = \delta_2 = \ldots = \delta_k = \alpha_1 = \alpha_2 = \ldots = \alpha_k = \theta_1 = \theta_2 = \ldots = \theta_k C_2 \neq 0\)

Existe evidencia estadística de que los residuos son heterocedásticos (la varianza del error no es constante).

Prueba de White

library(lmtest)
options(scipen = 99) # Evitar notación científica en los resultados
# 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  + bmi * children + bmi * smoker  +  children * smoker , data = insurance)

# Mostrar los resultados
print(PWhite)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_seguro
## BP = 128.98, df = 13, p-value < 0.00000000000000022

Usando Libreraria Skedastic (Otra alternativa)

library(skedastic)
## Warning: package 'skedastic' was built under R version 4.5.3
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      113. 1.59e-19        10 White's Test greater

Verificación de Autocorrelación

Prueba de Durbin-Watson

Se asume que los residuos presentan la siguiente estructura:

\(u_i = \rho u_{i-1} + v_i\)

\(H_0 : \rho = 0\)

No existe autocorrelación de primer orden en los residuos del modelo estimado

\(H_1 : \rho \neq 0\)

Se detecta presencia de autocorrelación de primer orden en los residuos del modelo estimado

library(lmtest)
dwtest(modelo_seguro,alternative = "two.sided",iterations = 1000)
## 
##  Durbin-Watson test
## 
## data:  modelo_seguro
## DW = 2.0473, p-value = 0.3876
## alternative hypothesis: true autocorrelation is not 0
library(car)
## Warning: package 'car' was built under R version 4.5.3
## Cargando paquete requerido: carData
durbinWatsonTest(modelo_seguro,simulate = TRUE,reps = 1000)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.02468678      2.047337   0.386
##  Alternative hypothesis: rho != 0

Multiplicadores de Lagrange (Breusch-Godfrey)

\(u_j = \rho_1 u_{i-1} + \rho_2 u_{i-2} + \ldots + \rho_m u_{i-m} + u_j\)

\(H_0: \rho_1 = \rho_2 = \rho_3 = \ldots = \rho_m = 0\)

Los residuos del modelo no presentan dependencia entre sí y carecen de autocorrelación hasta el orden \(m\).

\(H_1: \rho_1 = \rho_2 = \rho_3 = \ldots = \rho_m \neq 0\)

Los residuos del modelo presentan dependencia entre sí y exhiben autocorrelación hasta el orden \(m\).

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 = 0.83555, df = 2, p-value = 0.6585
library(lmtest)
bgtest(modelo_seguro,order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo_seguro
## LM test = 0.81855, df = 1, p-value = 0.3656

Forma Gráfica

# Extraemos los residuos de tu modelo
residuos <- residuals(modelo_seguro)

# Graficamos con ZOOM
acf(residuos, 
    lag.max = 5,                      # Zoom en X: Muestra solo hasta el lag 5
    ylim = c(-0.2, 0.2),              # Zoom en Y: Corta el gráfico entre -0.2 y 0.2
    main = "Correlograma de los Residuos",
    xlab = "Lag", 
    ylab = "ACF")

# Estimadores HAC

options(scipen = 9999)
library(lmtest)
#Modelo Sin corregir:
coeftest(modelo_seguro)
## 
## t test of coefficients:
## 
##                 Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   7.27183076 0.07427125 97.9091 <0.0000000000000002 ***
## age           0.03500144 0.00085352 41.0082 <0.0000000000000002 ***
## children      0.10228988 0.00988305 10.3500 <0.0000000000000002 ***
## smokeryes     0.18624931 0.14782454  1.2599              0.2079    
## bmi           0.00085930 0.00222231  0.3867              0.6991    
## smokeryes:bmi 0.04421299 0.00471952  9.3681 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo Corregido
options(scipen = 9999)
library(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)   7.27183076 0.07610536 95.5495 <0.0000000000000002 ***
## age           0.03500144 0.00098966 35.3671 <0.0000000000000002 ***
## children      0.10228988 0.00887285 11.5284 <0.0000000000000002 ***
## smokeryes     0.18624931 0.13388648  1.3911              0.1644    
## bmi           0.00085930 0.00214688  0.4003              0.6890    
## smokeryes:bmi 0.04421299 0.00436163 10.1368 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
options(scipen = 99)
library(robustbase)
library(stargazer)
modelo_seguro_robust<-lmrob(log(charges) ~ age + bmi + children + smoker  + smoker:bmi, data = insurance)
# print(summary(modelo_seguro_robust))
stargazer(modelo_seguro,modelo_seguro_robust,type = "html",title = "comparativa")
comparativa
Dependent variable:
log(charges)
OLS MM-type
linear
(1) (2)
age 0.035*** 0.042
(0.001)
children 0.102*** 0.111
(0.010)
smokeryes 0.186 0.414
(0.148)
bmi:smokeryes 0.032
bmi 0.001 -0.002
(0.002)
smokeryes:bmi 0.044***
(0.005)
Constant 7.272*** 6.964
(0.074)
Observations 1,338 1,338
R2 0.777 0.969
Adjusted R2 0.776 0.969
Residual Std. Error (df = 1332) 0.435 0.129
F Statistic 927.462*** (df = 5; 1332)
Note: p<0.1; p<0.05; p<0.01

Pronosticos

options(scipen = 99)
library(lmtest)
library(stargazer)
library(equatiomatic) # optativo remotes::install_github("datalorax/equatiomatic")

#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 = log(charges) ~ age + bmi + children + smoker + smoker:bmi, data = insurance)
extract_eq(modelo_seguro,wrap = TRUE) #optativo

\[ \begin{aligned} \operatorname{log(charges)} &= \alpha + \beta_{1}(\operatorname{age}) + \beta_{2}(\operatorname{bmi}) + \beta_{3}(\operatorname{children})\ + \\ &\quad \beta_{4}(\operatorname{smoker}_{\operatorname{yes}}) + \beta_{5}(\operatorname{bmi} \times \operatorname{smoker}_{\operatorname{yes}}) + \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")
# 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)

# Reconvertir de log a dolares aplicando exp()
pred_95_exp <- exp(pred_95)
pred_99_exp <- exp(pred_99)

predicciones <- list()
predicciones$fit <- rbind(pred_95_exp, pred_99_exp)

rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym","Li","Ls")
stargazer(predicciones$fit,
          title = "Pronosticos e intervalos de confianza (en dolares)",
          type = "html") #Poner results='asis' en opciones del chunk
## 
## <table style="text-align:center"><caption><strong>Pronosticos e intervalos de confianza (en dolares)</strong></caption>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>Ym</td><td>Li</td><td>Ls</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">95</td><td>6,160.733</td><td>2,621.688</td><td>14,477.170</td></tr>
## <tr><td style="text-align:left">99</td><td>6,160.733</td><td>2,003.236</td><td>18,946.660</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr></table>

Predicción usnado libreria forescast

library(forecast)
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.3
options(scipen = 2)
#Data para la predicción X'm
X_m<-data.frame(age=35, bmi=28.5, children=2, smoker="no")
#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)

# Reconvertir de log a dolares aplicando exp()
pronosticos_df <- as.data.frame(pronosticos)
pronosticos_exp <- exp(pronosticos_df)

kable(pronosticos_exp,
      caption = "Pronostico e intervalos de confianza (en dolares):",
      digits = 1,format = "html") #Poner results='asis' en opciones del chunk
Pronostico e intervalos de confianza (en dolares):
Point Forecast Lo 95 Hi 95 Lo 99 Hi 99
6160.7 2621.7 14477.2 2003.2 18946.7

#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 = log(charges) ~ age + bmi + children + smoker  + smoker:bmi, 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,
                    log(Datos_Entrenamiento$charges)),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        log(Datos_Entrenamiento$charges)),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      log(Datos_Entrenamiento$charges)),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       log(Datos_Entrenamiento$charges))*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         log(Datos_Entrenamiento$charges),type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                  log(Datos_Entrenamiento$charges)),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                  log(Datos_Entrenamiento$charges)),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                  log(Datos_Entrenamiento$charges))
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            RMSE = RMSE(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            MAE = MAE(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            MAPE= MAPE(Pronostico_Prueba[[j]], log(Datos_Prueba$charges))*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], log(Datos_Prueba$charges),
                         type = 1), # También se puede usar la función que creamos: THEIL_U
            Um=Um(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            Us=Us(Pronostico_Prueba[[j]], log(Datos_Prueba$charges)),
            Uc=Uc(Pronostico_Prueba[[j]], log(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)
Medidas de Performance Datos del Modelo (Entrenamiento 80%)
Statistic N Mean St. Dev. Min Max
R2 5,000 0.777 0.008 0.753 0.803
RMSE 5,000 0.434 0.007 0.408 0.455
MAE 5,000 0.263 0.005 0.244 0.282
MAPE 5,000 2.886 0.061 2.660 3.100
THEIL 5,000 0.024 0.0004 0.022 0.025
Um 5,000 0.000 0.000 0 0
Us 5,000 0.063 0.002 0.055 0.071
Uc 5,000 0.938 0.002 0.930 0.946 
bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación (Predictivo - 5,000 Réplicas)",
            type = "html",
            digits = 3)
Medidas de Performance Simulación (Predictivo - 5,000 Réplicas)
Statistic N Mean St. Dev. Min Max
R2 5,000 0.777 0.031 0.673 0.872
RMSE 5,000 0.436 0.028 0.343 0.527
MAE 5,000 0.266 0.014 0.216 0.321
MAPE 5,000 2.911 0.155 2.364 3.495
THEIL 5,000 0.024 0.002 0.019 0.029
Um 5,000 0.004 0.005 0.000 0.047
Us 5,000 0.068 0.034 0.0004 0.223
Uc 5,000 0.932 0.034 0.773 1.000