options(scipen = 999999)## Warning in options(scipen = 999999): invalid 'scipen' 999999, used 9999library(lmtest)## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(equatiomatic) # optativo remotes::install_github("datalorax/equatiomatic")##
## Adjuntando el paquete: 'equatiomatic'## The following object is masked from 'package:datasets':
##
## penguinslibrary(mlbench) #esta librería tiene el data frame BostonHousing
data(BostonHousing)
#Modelo estimado medv~. indica "medv" en función del resto de variables del dataframe
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) 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 ' ' 1library(stargazer)
#Data para la predicción 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)
#Predicción 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 = "Pronósticos e intervalos de confianza",
type = "html")##
## <table style="text-align:center"><caption><strong>Pronósticos e intervalos de confianza</strong></caption>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>Ym</td><td>Li</td><td>Ls</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">95</td><td>26.116</td><td>15.558</td><td>36.673</td></tr>
## <tr><td style="text-align:left">99</td><td>26.116</td><td>12.221</td><td>40.010</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr></table>library(forecast)
library(kableExtra)
#Data para la predicción 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)
#Nivel de confianza para el intervalo de confianza
confidense<-c(0.95,0.99)
#Realizando el pronóstico con forecast
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 |
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.## Warning in options(scipen = 999999): invalid 'scipen' 999999, used 9999library(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)
)
}bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "html",
digits = 3)##
## <table style="text-align:center"><caption><strong>Medidas de Performance Datos del Modelo</strong></caption>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">R2</td><td>1,000</td><td>0.743</td><td>0.013</td><td>0.713</td><td>0.794</td></tr>
## <tr><td style="text-align:left">RMSE</td><td>1,000</td><td>4.653</td><td>0.141</td><td>4.177</td><td>4.948</td></tr>
## <tr><td style="text-align:left">MAE</td><td>1,000</td><td>3.265</td><td>0.095</td><td>2.905</td><td>3.512</td></tr>
## <tr><td style="text-align:left">MAPE</td><td>1,000</td><td>16.387</td><td>0.464</td><td>14.813</td><td>17.691</td></tr>
## <tr><td style="text-align:left">THEIL</td><td>1,000</td><td>0.096</td><td>0.003</td><td>0.087</td><td>0.102</td></tr>
## <tr><td style="text-align:left">Um</td><td>1,000</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">Us</td><td>1,000</td><td>0.074</td><td>0.004</td><td>0.058</td><td>0.085</td></tr>
## <tr><td style="text-align:left">Uc</td><td>1,000</td><td>0.928</td><td>0.004</td><td>0.918</td><td>0.945</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "html",
digits = 3)##
## <table style="text-align:center"><caption><strong>Medidas de Performance Simulación</strong></caption>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">R2</td><td>1,000</td><td>0.723</td><td>0.056</td><td>0.452</td><td>0.840</td></tr>
## <tr><td style="text-align:left">RMSE</td><td>1,000</td><td>4.862</td><td>0.575</td><td>3.465</td><td>6.961</td></tr>
## <tr><td style="text-align:left">MAE</td><td>1,000</td><td>3.411</td><td>0.281</td><td>2.633</td><td>4.492</td></tr>
## <tr><td style="text-align:left">MAPE</td><td>1,000</td><td>17.197</td><td>1.618</td><td>12.875</td><td>23.137</td></tr>
## <tr><td style="text-align:left">THEIL</td><td>1,000</td><td>0.101</td><td>0.012</td><td>0.073</td><td>0.148</td></tr>
## <tr><td style="text-align:left">Um</td><td>1,000</td><td>0.011</td><td>0.016</td><td>0.000</td><td>0.205</td></tr>
## <tr><td style="text-align:left">Us</td><td>1,000</td><td>0.081</td><td>0.066</td><td>0.00000</td><td>0.333</td></tr>
## <tr><td style="text-align:left">Uc</td><td>1,000</td><td>0.918</td><td>0.066</td><td>0.667</td><td>1.010</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>