Estimación del Modelo

options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic) 
library(mlbench) 
data(BostonHousing)
modelo_boston<-lm(formula = medv~.,data=BostonHousing)
extract_eq(modelo_boston,wrap = TRUE)

\[ \begin{aligned} \operatorname{medv} &= \alpha + \beta_{1}(\operatorname{crim}) + \beta_{2}(\operatorname{zn}) + \beta_{3}(\operatorname{indus})\ + \\ &\quad \beta_{4}(\operatorname{chas}_{\operatorname{1}}) + \beta_{5}(\operatorname{nox}) + \beta_{6}(\operatorname{rm}) + \beta_{7}(\operatorname{age})\ + \\ &\quad \beta_{8}(\operatorname{dis}) + \beta_{9}(\operatorname{rad}) + \beta_{10}(\operatorname{tax}) + \beta_{11}(\operatorname{ptratio})\ + \\ &\quad \beta_{12}(\operatorname{b}) + \beta_{13}(\operatorname{lstat}) + \epsilon \end{aligned} \]

coeftest(modelo_boston)
## 
## t test of coefficients:
## 
##                 Estimate   Std. Error  t value              Pr(>|t|)    
## (Intercept)  36.45948839   5.10345881   7.1441    0.0000000000032834 ***
## crim         -0.10801136   0.03286499  -3.2865             0.0010868 ** 
## zn            0.04642046   0.01372746   3.3816             0.0007781 ***
## indus         0.02055863   0.06149569   0.3343             0.7382881    
## chas1         2.68673382   0.86157976   3.1184             0.0019250 ** 
## nox         -17.76661123   3.81974371  -4.6513    0.0000042456438076 ***
## rm            3.80986521   0.41792525   9.1161 < 0.00000000000000022 ***
## age           0.00069222   0.01320978   0.0524             0.9582293    
## dis          -1.47556685   0.19945473  -7.3980    0.0000000000006013 ***
## rad           0.30604948   0.06634644   4.6129    0.0000050705290227 ***
## tax          -0.01233459   0.00376054  -3.2800             0.0011116 ** 
## ptratio      -0.95274723   0.13082676  -7.2825    0.0000000000013088 ***
## b             0.00931168   0.00268596   3.4668             0.0005729 ***
## lstat        -0.52475838   0.05071528 -10.3471 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#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)))
}

SCRIPT DE SIMULACIÓN

options(scipen = 999999) 
library(dplyr) 
library(caret) 
library(DescTools) 
library(stargazer)  
set.seed(50)
numero_de_muestras<-1000 
muestras<- BostonHousing$medv %>%
  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<- BostonHousing[muestras[[j]], ]
Datos_Prueba<- BostonHousing[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = medv~.,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$medv),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        Datos_Entrenamiento$medv),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      Datos_Entrenamiento$medv),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       Datos_Entrenamiento$medv)*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv,type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$medv)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$medv,
                         type = 1), 
            Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$medv)
            )
}

Resultados de la Simulación

Medidas de Performance Datos del Modelo

bind_rows(Resultados_Performance_data_entrenamiento) %>% 
  stargazer(title = "Medidas de Performance Datos del Modelo",
            type = "html",
            digits = 3)
Medidas de Performance Datos del Modelo
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 1,000 0.743 0.013 0.713 0.734 0.751 0.794
RMSE 1,000 4.653 0.141 4.177 4.565 4.759 4.948
MAE 1,000 3.265 0.095 2.905 3.204 3.332 3.512
MAPE 1,000 16.387 0.464 14.813 16.085 16.718 17.691
THEIL 1,000 0.096 0.003 0.087 0.095 0.099 0.102
Um 1,000 0.000 0.000 0 0 0 0
Us 1,000 0.074 0.004 0.058 0.072 0.077 0.085
Uc 1,000 0.928 0.004 0.918 0.925 0.931 0.945

Medidas de Performance Simulación

bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación",
            type = "html",
            digits = 3)
Medidas de Performance Simulación
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 1,000 0.723 0.056 0.452 0.690 0.764 0.840
RMSE 1,000 4.862 0.575 3.465 4.450 5.226 6.961
MAE 1,000 3.411 0.281 2.633 3.216 3.598 4.492
MAPE 1,000 17.197 1.618 12.875 16.066 18.262 23.137
THEIL 1,000 0.101 0.012 0.073 0.092 0.108 0.148
Um 1,000 0.011 0.016 0.000 0.001 0.015 0.205
Us 1,000 0.081 0.066 0.00000 0.027 0.122 0.333
Uc 1,000 0.918 0.066 0.667 0.875 0.971 1.010