#Modelo estimado
options(spicen= 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) 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#Pronostico
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)
#Intervalos de confianza del 95% y del 99%
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")| Ym | Li | Ls | |
| 95 | 26.116 | 15.558 | 36.673 |
| 99 | 26.116 | 12.221 | 40.010 |
#Predicción usando librería “forecast”
library(forecast)
library(kableExtra)
X2_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= X2_m,ts=FALSE)
kable(pronosticos,
caption = "Pronóstico e intervalos de confianza:",
digits = 2,format = "html")| Point Forecast | Lo 95 | Hi 95 | Lo 99 | Hi 99 |
|---|---|---|---|---|
| 26.12 | 15.56 | 36.67 | 12.22 | 40.01 |
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)))}#Simulación
options(scipen = 999999)
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(50)
numero_de_muestras<-1000
# 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),
Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$medv)
)
}
bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "html",
digits = 3)| 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 Simulación",
type = "html",
digits = 3)| 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 |