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 ' ' 1Um<-function(Pronosticado,observado){
library(DescTools)
((mean(Pronosticado)-mean(observado))˄2)/MSE(Pronosticado,observado)
}Us<-function(pronosticado,observado){
library(DescTools)
((sd(pronosticado)-sd(observado))˄2)/MSE(pronosticado,observado)
}Uc<-function(pronosticado,observado){
library(DescTools)
(2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/
MSE(pronosticado,observado)}THEIL_U<-function(pronosticado,observado){
library(DescTools)
RMSE(pronosticado,observado)/(sqrt(mean(pronosticado˄2))+sqrt(mean(observado˄2)))
}options(scipen = 999999)
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(50)
numero_de_muestras<-1000
## Lista con las 1000 muestras (indica la posicion de la fila en cada dataframe)
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 prueba
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)
## Estimacion de los modelos lineales para cada muestra, los pronosticos y calculos de los estadisticos 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 Simulacion
bind_rows(Resultados_Performance_data_entrenamiento)%>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "text",
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
## ---------------------------------------------------------------bind_rows(Resultados_Performance)%>%
stargazer(title = "Medidas de Performance de Simulacion",
type = "text",
digits = 3)##
## Medidas de Performance de Simulacion
## ================================================================
## 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
## ----------------------------------------------------------------