Estimación del modelo
library(lmtest)
library(stargazer)
library(equatiomatic)
library(mlbench)
data(BostonHousing)
modelo_boston <- lm(medv ~ ., data = BostonHousing)
extract_eq(model = 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
Preducción usando “Predict”
X_m<-data.frame(crim=0.05,zn=15,indus=2,chas="0",nox=0.004,
rm=5,age=85,dis=5.56,rad=2,tax=300,ptratio=17,b=0.00005,lstat=5)
confidense<-c(0.95,0.99)
predict(object = modelo_boston,
newdata = X_m,
interval = "prediction",
level = confidense,
se.fit =TRUE)->predicciones
rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym","Li","Ls")
stargazer(predicciones$fit,
title = "Pronósticos e intervalos de confianza",
type = "html", digits = 2)
Pronósticos e intervalos de confianza
|
|
|
|
Ym
|
Li
|
Ls
|
|
|
|
95
|
26.12
|
15.56
|
36.67
|
|
99
|
26.12
|
12.22
|
40.01
|
|
|
Predicción Usando forecast
library(forecast)
library(kableExtra)
X_m<-data.frame(crim=0.05,zn=15,indus=2,chas="0",nox=0.004,
rm=5,age=85,dis=5.56,rad=2,tax=300,ptratio=17,b=0.00005,lstat=5)
confidense<-c(0.95,0.99)
pronosticos<-forecast(object = modelo_boston,
level = confidense,
newdata = X_m,ts = FALSE)
kable(pronosticos,
caption = "Pronóstico e intervalos de confianza:",
digits = 2,format = "html")
Pronóstico e intervalos de confianza:
|
Point Forecast
|
Lo 95
|
Hi 95
|
Lo 99
|
Hi 99
|
|
26.12
|
15.56
|
36.67
|
12.22
|
40.01
|
Script de Simulación
Um<-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)))
}
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(50)
numero_de_muestras<-500
muestras<- BostonHousing$medv %>%
createDataPartition(p = 0.75,
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)
)
}
Resultados de simulación
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
|
Max
|
|
|
|
R2
|
500
|
0.743
|
0.016
|
0.704
|
0.797
|
|
RMSE
|
500
|
4.646
|
0.172
|
4.084
|
4.957
|
|
MAE
|
500
|
3.263
|
0.116
|
2.944
|
3.503
|
|
MAPE
|
500
|
16.414
|
0.546
|
14.578
|
17.849
|
|
THEIL
|
500
|
0.096
|
0.003
|
0.085
|
0.103
|
|
Um
|
500
|
0.000
|
0.000
|
0
|
0
|
|
Us
|
500
|
0.074
|
0.005
|
0.057
|
0.088
|
|
Uc
|
500
|
0.928
|
0.005
|
0.915
|
0.946
|
|
|
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
|
Max
|
|
|
|
R2
|
500
|
0.725
|
0.051
|
0.524
|
0.832
|
|
RMSE
|
500
|
4.868
|
0.518
|
3.788
|
6.462
|
|
MAE
|
500
|
3.411
|
0.257
|
2.703
|
4.208
|
|
MAPE
|
500
|
17.060
|
1.528
|
12.186
|
21.680
|
|
THEIL
|
500
|
0.101
|
0.011
|
0.078
|
0.134
|
|
Um
|
500
|
0.009
|
0.012
|
0.000
|
0.081
|
|
Us
|
500
|
0.087
|
0.064
|
0.00001
|
0.295
|
|
Uc
|
500
|
0.912
|
0.066
|
0.692
|
1.008
|
|
|