Replica de Ejercicio de Simulacion

options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic) 
library(mlbench) 
data(BostonHousing)
bostonmodel<-lm(formula = medv~.,data=BostonHousing)
extract_eq(bostonmodel,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(bostonmodel)

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

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

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)

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