Ejercicio Simulación vista en clase
Estimación Modelo
library(lmtest)
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) 3.6459e+01 5.1035e+00 7.1441 3.283e-12 ***
## crim -1.0801e-01 3.2865e-02 -3.2865 0.0010868 **
## zn 4.6420e-02 1.3727e-02 3.3816 0.0007781 ***
## indus 2.0559e-02 6.1496e-02 0.3343 0.7382881
## chas1 2.6867e+00 8.6158e-01 3.1184 0.0019250 **
## nox -1.7767e+01 3.8197e+00 -4.6513 4.246e-06 ***
## rm 3.8099e+00 4.1793e-01 9.1161 < 2.2e-16 ***
## age 6.9222e-04 1.3210e-02 0.0524 0.9582293
## dis -1.4756e+00 1.9945e-01 -7.3980 6.013e-13 ***
## rad 3.0605e-01 6.6346e-02 4.6129 5.071e-06 ***
## tax -1.2335e-02 3.7605e-03 -3.2800 0.0011116 **
## ptratio -9.5275e-01 1.3083e-01 -7.2825 1.309e-12 ***
## b 9.3117e-03 2.6860e-03 3.4668 0.0005729 ***
## lstat -5.2476e-01 5.0715e-02 -10.3471 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Predicciones Stargazer
library(stargazer)
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)
predicciones<-predict(object = modelo_boston, newdata = X_m,
interval = "prediction", level = confidense,
se.fit =TRUE)
rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym", "Li", "Ls")
stargazer(predicciones$fit, type = "html", digits = 4,
title = "Pronóstico - Intervalo de Confianza")
Pronóstico - Intervalo de Confianza
|
|
|
|
Ym
|
Li
|
Ls
|
|
|
|
95
|
26.1157
|
15.5582
|
36.6732
|
|
99
|
26.1157
|
12.2211
|
40.0103
|
|
|
Librería
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)
#Pronóstico (Forecast)
pronosticos<-forecast(object = modelo_boston, level = confidense,
newdata = X_m, ts = FALSE)
kable(pronosticos, caption = "**Pronóstico - Intervalo de Confianza:",
digits = 2, format = "html")
**Pronóstico - Intervalo de Confianza:
|
Point Forecast
|
Lo 95
|
Hi 95
|
Lo 99
|
Hi 99
|
|
26.12
|
15.56
|
36.67
|
12.22
|
40.01
|
Simulación
#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
THEIL_U <- function(pronosticado,observado){
library(DescTools)
RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}
Simulación
library(dplyr) # Para manejo de datos y activar el operador "pipe" %>%
library(caret) # Permite Realizar muestreo sobre los data frame
library(DescTools) # Contiene las funciones para calcular las medidas de performance
library(stargazer) # Para dar formato, y obtener resumen estadistico de las simulaciones
set.seed(50) # Permite fijar la semilla aleatoria, para reproducir los resultados obtenidos en esta clase
numero_de_muestras<-1000 # Numero de muestras que se optendran del data frame
# Se crea la lista con las 1000 muestras (indica la posición de la fila en cada data frame)
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), # También se puede usar la función que creamos: THEIL_U
Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$medv)
)
}
Simulacion vista de resultados
# Resultados 1
library(stargazer)
bind_rows(Resultados_Performance) %>%
stargazer(title = "Performance Simulación", type = "html", digits = 2)
Performance Simulación
|
|
|
Statistic
|
N
|
Mean
|
St. Dev.
|
Min
|
Pctl(25)
|
Pctl(75)
|
Max
|
|
|
|
R2
|
1,000
|
0.72
|
0.06
|
0.45
|
0.69
|
0.76
|
0.84
|
|
RMSE
|
1,000
|
4.86
|
0.57
|
3.46
|
4.45
|
5.23
|
6.96
|
|
MAE
|
1,000
|
3.41
|
0.28
|
2.63
|
3.22
|
3.60
|
4.49
|
|
MAPE
|
1,000
|
17.20
|
1.62
|
12.88
|
16.07
|
18.26
|
23.14
|
|
THEIL
|
1,000
|
0.10
|
0.01
|
0.07
|
0.09
|
0.11
|
0.15
|
|
Um
|
1,000
|
0.01
|
0.02
|
0.00
|
0.001
|
0.01
|
0.20
|
|
Us
|
1,000
|
0.08
|
0.07
|
0.0000
|
0.03
|
0.12
|
0.33
|
|
Uc
|
1,000
|
0.92
|
0.07
|
0.67
|
0.88
|
0.97
|
1.01
|
|
|
#Resultado 2
library(stargazer)
bind_rows(Resultados_Performance) %>%
stargazer(title = "Performance Simulación", type = "html", digits = 2)
Performance Simulación
|
|
|
Statistic
|
N
|
Mean
|
St. Dev.
|
Min
|
Pctl(25)
|
Pctl(75)
|
Max
|
|
|
|
R2
|
1,000
|
0.72
|
0.06
|
0.45
|
0.69
|
0.76
|
0.84
|
|
RMSE
|
1,000
|
4.86
|
0.57
|
3.46
|
4.45
|
5.23
|
6.96
|
|
MAE
|
1,000
|
3.41
|
0.28
|
2.63
|
3.22
|
3.60
|
4.49
|
|
MAPE
|
1,000
|
17.20
|
1.62
|
12.88
|
16.07
|
18.26
|
23.14
|
|
THEIL
|
1,000
|
0.10
|
0.01
|
0.07
|
0.09
|
0.11
|
0.15
|
|
Um
|
1,000
|
0.01
|
0.02
|
0.00
|
0.001
|
0.01
|
0.20
|
|
Us
|
1,000
|
0.08
|
0.07
|
0.0000
|
0.03
|
0.12
|
0.33
|
|
Uc
|
1,000
|
0.92
|
0.07
|
0.67
|
0.88
|
0.97
|
1.01
|
|
|