options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic)
remotes::install_github("datalorax/equatiomatic")##
## * checking for file 'C:\Users\SANCHEZ\AppData\Local\Temp\RtmpkHJk8z\remotes25e436323009\datalorax-equatiomatic-0ac8f2b/DESCRIPTION' ... OK
## * preparing 'equatiomatic':
## * checking DESCRIPTION meta-information ... OK
## * checking for LF line-endings in source and make files and shell scripts
## * checking for empty or unneeded directories
## * building 'equatiomatic_0.3.3.tar.gz'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} \]
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)
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 = "text") 95 26.116 15.558 36.673 99 26.116 12.221 40.010 ———————–
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(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") | 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 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)))
}options(scipen = 999999) #No mostrar notación cientifica.
library(dplyr) # Para manejo de datos y activar el operador "pipe" %>%##
## Adjuntando el paquete: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(caret) # Permite Realizar muestreo sobre los data frame## Cargando paquete requerido: ggplot2## Cargando paquete requerido: latticelibrary(DescTools) # Contiene las funciones para calcular las medidas de performance##
## Adjuntando el paquete: 'DescTools'## The following objects are masked from 'package:caret':
##
## MAE, RMSE## The following object is masked from 'package:forecast':
##
## BoxCoxlibrary(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)
)
} #No olvidar este corchete ;)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 Max
## ---------------------------------------------
## R2 1,000 0.743 0.013 0.713 0.794
## RMSE 1,000 4.653 0.141 4.177 4.948
## MAE 1,000 3.265 0.095 2.905 3.512
## MAPE 1,000 16.387 0.464 14.813 17.691
## THEIL 1,000 0.096 0.003 0.087 0.102
## Um 1,000 0.000 0.000 0 0
## Us 1,000 0.074 0.004 0.058 0.085
## Uc 1,000 0.928 0.004 0.918 0.945
## ---------------------------------------------bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "text",
digits = 3)##
## Medidas de Performance Simulación
## ==============================================
## Statistic N Mean St. Dev. Min Max
## ----------------------------------------------
## R2 1,000 0.723 0.056 0.452 0.840
## RMSE 1,000 4.862 0.575 3.465 6.961
## MAE 1,000 3.411 0.281 2.633 4.492
## MAPE 1,000 17.197 1.618 12.875 23.137
## THEIL 1,000 0.101 0.012 0.073 0.148
## Um 1,000 0.011 0.016 0.000 0.205
## Us 1,000 0.081 0.066 0.00000 0.333
## Uc 1,000 0.918 0.066 0.667 1.010
## ----------------------------------------------load("C:/Users/SANCHEZ/Downloads/Gabriel Arturo Sánchez Henríquez - smoke.RData")
library(stargazer)
stargazer(desc,
title = "Descripción de las variables",
summary = FALSE,
type = "text",
rownames = FALSE)educ years of schooling
cigpric state cig. price, cents/pack white =1 if white
age in years
income annual income,
cigs cigs. smoked per day
restaurn =1 if rest. smk. restrictions lincome log(income)
agesq age2
lcigpric log(cigprice)
————————————–
stargazer(data[1:7,],
title = "Data frame",
summary=FALSE,
type="text",
rownames=FALSE)##
## Data frame
## ====================================================================
## educ cigpric white age income cigs restaurn lincome agesq lcigpric
## --------------------------------------------------------------------
## 16 60.506 1 46 20,000 0 0 9.903 2,116 4.103
## 16 57.883 1 40 30,000 0 0 10.309 1,600 4.058
## 12 57.664 1 58 30,000 3 0 10.309 3,364 4.055
## 13.500 57.883 1 30 20,000 0 0 9.903 900 4.058
## 10 58.320 1 17 20,000 0 0 9.903 289 4.066
## 6 59.340 1 86 6,500 0 0 8.780 7,396 4.083
## 12 57.883 1 35 20,000 0 0 9.903 1,225 4.058
## --------------------------------------------------------------------modelo<- lm(cigs ~ cigpric + lcigpric + income + lincome + age + agesq + educ + white + restaurn, data = data)
stargazer(modelo, title = "Modelo", type = "text", digits = 4) Dependent variable:
---------------------------
cigs | cigpric 2.0023 (1.4928) |
| lcigpric -115.2735 (85.4243) |
| income -0.00005 (0.0001) |
| lincome 1.4041 (1.7082) |
| age 0.7784*** (0.1606) |
| agesq -0.0092*** (0.0017) |
| educ -0.4948*** (0.1682) |
| white -0.5311 (1.4607) |
| restaurn -2.6442** (1.1300) |
| Constant 340.8044 (260.0156) |
Observations 807
R2 0.0552
Adjusted R2 0.0445
Residual Std. Error 13.4128 (df = 797)
F Statistic 5.1693*** (df = 9; 797)
=============================================== Note: p<0.1;
p<0.05; p<0.01
library(forecast)
library(kableExtra)
X_m1<-data.frame(cigpric=57.664,lcigpric=4.058424,income=20000,lincome=9.903487,age=31,agesq=7396,educ=10,white=1,restaurn=2)
confidense<-c(0.95,0.99)
pronosticos<-forecast(object = modelo,
level = confidense,
newdata = X_m1,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 |
|---|---|---|---|---|
| -52.9 | -87.29 | -18.5 | -98.14 | -7.65 |
options(scipen = 999999) #No mostrar notación cientifica.
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<-500 # 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<- data$cigs %>%
createDataPartition(p = 0.75,
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<- data[muestras[[j]], ]
Datos_Prueba<- data[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = cigs~.,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$cigs),
RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs),
MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs),
MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs)*100,
THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs,type = 1),
Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs),
Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs),
Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$cigs)
)
Resultados_Performance[[j]]<-data.frame(
R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$cigs)*100,
THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$cigs,
type = 1), # También se puede usar la función que creamos: THEIL_U
Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$cigs)
)
} #No olvidar este corchete ;)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 Max
## --------------------------------------------
## R2 500 0.058 0.007 0.040 0.083
## RMSE 500 13.304 0.198 12.629 13.828
## MAE 500 10.562 0.133 10.068 10.969
## MAPE 500 Inf.000 Inf Inf
## THEIL 500 0.522 0.006 0.505 0.541
## Um 500 0.000 0.000 0 0
## Us 500 0.613 0.020 0.554 0.669
## Uc 500 0.389 0.020 0.332 0.448
## --------------------------------------------bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "text",
digits = 3)##
## Medidas de Performance Simulación
## =============================================
## Statistic N Mean St. Dev. Min Max
## ---------------------------------------------
## R2 500 0.039 0.018 0.002 0.099
## RMSE 500 13.479 0.595 11.734 15.361
## MAE 500 10.728 0.290 9.834 11.626
## MAPE 500 Inf.000 Inf Inf
## THEIL 500 0.528 0.013 0.491 0.564
## Um 500 0.002 0.003 0.00000 0.021
## Us 500 0.597 0.051 0.476 0.751
## Uc 500 0.405 0.051 0.253 0.529
## ---------------------------------------------