options(scipen = 999999)

library(mlbench)
library(lmtest)
library(stargazer)
library(equatiomatic)

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} \]

Ejemplo de simulación

#Proporción de sesgo
Um<-function(pronosticado,observado){
  library(DescTools)
  ((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado)
  }

#Proporción de varianza
Us<-function(pronosticado,observado){
  library(DescTools)
  ((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
  }

#Proporción de Covarianza
Uc<-function(pronosticado,observado){
  library(DescTools)
  (2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)
  }

#Coeficiente U de Theil
library(DescTools)
THEIL_U <-function(pronosticado,observado) {
   RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
  }

Scrip de simulación

Modificaciones

#DataFrame
DataFrame1 <- BostonHousing
#Formula del modelo
formula1 <- as.formula(medv~.) 
#Semilla
set.seed(50)
#Data de entrenamiento
DEntrenam <- 0.8
#Muestras
numero_de_muestras<-1000 
#Cambio de endógena más adelante
options(scipen = 999999)
library(dplyr) 
library(caret) 
library(DescTools)
library(stargazer) 
muestras<- DataFrame1$medv %>% createDataPartition(p = DEntrenam,
                      times = numero_de_muestras,
                      list = TRUE)
# Listas vacias
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)
for(j in 1:numero_de_muestras){
#Data verde
  Datos_Entrenamiento<- DataFrame1[muestras[[j]], ]
#Data roja
  Datos_Prueba<- DataFrame1[-muestras[[j]], ]

#Cambio de endógena
DatEntrenamiento <- Datos_Entrenamiento$medv
DatPrueba <- Datos_Prueba$medv

#Continua el scrip
Modelos_Entrenamiento[[j]]<-lm(formula = formula1 ,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,
                    DatEntrenamiento),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        DatEntrenamiento),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      DatEntrenamiento),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       DatEntrenamiento)*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         DatEntrenamiento,type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         DatEntrenamiento),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         DatEntrenamiento),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         DatEntrenamiento)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], DatPrueba),
            RMSE = RMSE(Pronostico_Prueba[[j]], DatPrueba),
            MAE = MAE(Pronostico_Prueba[[j]], DatPrueba),
            MAPE= MAPE(Pronostico_Prueba[[j]], DatPrueba)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], DatPrueba,
                         type = 1), 
            Um=Um(Pronostico_Prueba[[j]], DatPrueba),
            Us=Us(Pronostico_Prueba[[j]], DatPrueba),
            Uc=Uc(Pronostico_Prueba[[j]], DatPrueba)
            )
}

Resultados del modelo

Medidas de Performance Datos del Modelo

#(data verde)
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

#(data roja)
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