options(scipen=999999)
library(lmtest)
library(stargazer)
library(equatiomatic)# optativo remotes::install_github("datalorax/equatiomatic")
library(mlbench)# esta libreria tiene el dataframe BostonHousing
data(BostonHousing)
#Modelo estimado medv~. indica "medv" en funcion del resto de variables del dataframe
modelo_boston<-lm(formula=medv~.,data= BostonHousing)
extract_eq(modelo_boston,wrap=TRUE)#optativo\[ \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} \]
library(stargazer)
#Data para la prediccion X'm
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)
#Prediccion usando predict
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=" Pronosticos e intervalos de confianza",
type="text")# Poner results="asis" en opciones del chunk##
## Pronosticos e intervalos de confianza
## =======================
## Ym Li Ls
## -----------------------
## 95 26.116 15.558 36.673
## 99 26.116 12.221 40.010
## -----------------------library(forecast)
library(kableExtra)
#Data para la prediccion Xm
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)
#Nivel de confianza para el intervalo de confianza
confidense<-c(0.95, 0.99)
#Realizando el pronostico con forecast
pronosticos<-forecast(object=modelo_boston,
level= confidense,
newdata= X_m, ts=FALSE)
kable(pronosticos,
caption="Pronosticos e intervalos de confianza:",
digits=2, format="html")#Poner results="asis" en opciones del chuck| Point Forecast | Lo 95 | Hi 95 | Lo 99 | Hi 99 |
|---|---|---|---|---|
| 26.12 | 15.56 | 36.67 | 12.22 | 40.01 |
# 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 proporcion
Uc<-function(pronosticado,observado){
library(DescTools)
(2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/
MSE(pronosticado,observado)}# coeficiente U de Theil(Tambien aparece en la libreria Desctools)
THEIL_U<-function(pronosticado,observado){library(Desctools)
RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))}options(scipen = 999999)# No mostrar notacion cientifica
library(dplyr)# Para manejo de datos y activar el operador "pipe" %>%
library(caret)# Permite realizar muestreo sobre Las 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
# Numerode muestrasque se optendran del dataframe
# Se crea la listacon las 1000 muestras (indica la posicion de la fila en cada data frame)
muestras<-BostonHousing$medv %>%
createDataPartition(p=0.8,
times = numero_de_muestras,
list=TRUE)
#listas vacias, que contendranlos datos de entrenamiento,los pronosticos para los datos de pruueba 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)
# Estimacion de los modelos lineales de cada muestra, los pronosticos y calculos 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),
# tambien se puede usar la funcion 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) )}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
## ----------------------------------------------------------------