options(scipen = 999999)
library(equatiomatic) # optativo remotes::install_github("datalorax/equatiomatic")
library(mlbench) #esta librería tiene el data frame BostonHousing
data(BostonHousing)options(scipen = 999999)
library(lmtest)
library(stargazer)
#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) #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} \]
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") #Poner results='asis' en opciones del chunk| Ym | Li | Ls | |
| 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)
#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") #Poner results='asis' en opciones del chunk| 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" %>%##
## Attaching package: '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## Loading required package: ggplot2## Loading required package: latticelibrary(DescTools) # Contiene las funciones para calcular las medidas de performance##
## Attaching package: '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 = "html",
digits = 3)| 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 = "html",
digits = 3)| 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/abyme/Downloads/Marta Abigail Meza Robles - smoke.RData")options(scipen = 9999999)
library(lmtest)
library(stargazer)
library(equatiomatic)
library(mlbench)
modelo<-lm(formula=cigs~.,data = data)
extract_eq(modelo,wrap = TRUE)\[ \begin{aligned} \operatorname{cigs} &= \alpha + \beta_{1}(\operatorname{educ}) + \beta_{2}(\operatorname{cigpric}) + \beta_{3}(\operatorname{white})\ + \\ &\quad \beta_{4}(\operatorname{age}) + \beta_{5}(\operatorname{income}) + \beta_{6}(\operatorname{restaurn}) + \beta_{7}(\operatorname{lincome})\ + \\ &\quad \beta_{8}(\operatorname{agesq}) + \beta_{9}(\operatorname{lcigpric}) + \epsilon \end{aligned} \]
coeftest(modelo)t test of coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 340.804374604 260.015587269 1.3107 0.190334
educ -0.494780616 0.168180198 -2.9420 0.003356 ** cigpric 2.002267667
1.492831189 1.3413 0.180220
white -0.531051635 1.460721806 -0.3636 0.716287
age 0.778359013 0.160555612 4.8479 0.0000015001 income
-0.000046194 0.000133491 -0.3460 0.729402
restaurn -2.644241351 1.129998690 -2.3400 0.019528
lincome 1.404061178 1.708165841 0.8220 0.411340
agesq -0.009150353 0.001749292 -5.2309 0.0000002158 * lcigpric
-115.273464445 85.424315195 -1.3494 0.177585
— Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05
‘.’ 0.1 ’ ’ 1
#Bias Proportion
Um<-function(pronosticado_1,observado_1){
library(DescTools)
((mean(pronosticado_1)-mean(observado_1))^2)/MSE(pronosticado_1,observado_1)
}
#Variance Proportion
Us<-function(pronosticado_1,observado_1){
library(DescTools)
((sd(pronosticado_1)-sd(observado_1))^2)/MSE(pronosticado_1,observado_1)
}
#Covariance Proportion
Uc<-function(pronosticado_1,observado_1){
library(DescTools)
(2*(1-cor(pronosticado_1,observado_1))*sd(pronosticado_1)*sd(observado_1))/MSE(pronosticado_1,observado_1)}
#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U<-function(pronosticado_1,observado_1){
library(DescTools)
RMSE(pronosticado_1,observado_1)/(sqrt(mean(pronosticado_1^2))+sqrt(mean(observado_1^2)))
}options(scipen = 999999)
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(50)
numero_muestras<-500
muestras_1<- data$cigs %>%
createDataPartition(p = 0.75,
times = numero_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_1<-vector(mode = "list",
length = numero_muestras)
Pronostico_Prueba_1<-vector(mode = "list",
length = numero_muestras)
Resultados_Performance_data_entrenamiento_1<-vector(mode = "list",
length = numero_muestras)
Resultados_Performance_1<-vector(mode = "list",
length = numero_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_muestras){
Datos_Entrenamiento_1<- data[muestras_1[[j]], ]
Datos_Prueba_1<- data[-muestras_1[[j]], ]
Modelos_Entrenamiento_1[[j]]<-lm(formula=cigs~cigpric+lcigpric+income+lincome+age+agesq+educ+white+restaurn,data=Datos_Entrenamiento_1)
Pronostico_Prueba_1[[j]]<-Modelos_Entrenamiento_1[[j]] %>% predict(Datos_Prueba_1)
Resultados_Performance_data_entrenamiento_1[[j]]<-data.frame(
R2_1 = R2(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs),
RMSE = RMSE(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs),
MAE = MAE(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs),
MAPE = MAPE(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs)*100,
THEIL=TheilU(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs,type = 1),
Um=Um(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs),
Us=Us(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs),
Uc=Uc(Modelos_Entrenamiento_1[[j]]$fitted.values,
Datos_Entrenamiento_1$cigs)
)
Resultados_Performance_1[[j]]<-data.frame(
R2_3 = R2(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs),
RMSE = RMSE(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs),
MAE = MAE(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs),
MAPE= MAPE(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs)*100,
THEIL=TheilU(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs,
type = 1), # También se puede usar la función que creamos: THEIL_U
Um=Um(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs),
Us=Us(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs),
Uc=Uc(Pronostico_Prueba_1[[j]], Datos_Prueba_1$cigs)
)
} #No olvidar este corchete ;)Resultados_Performance_data_entrenamiento_1<-bind_rows(Resultados_Performance_data_entrenamiento_1)
Resultados_Performance_data_entrenamiento_1<-
Resultados_Performance_data_entrenamiento_1[complete.cases(Resultados_Performance_data_entrenamiento_1),]
Resultados_Performance_data_entrenamiento_1%>%
stargazer(data,title = "Medidas de Performance Datos del Modelo",
type = "html",
digits = 3, summary.stat=c("n","mean","sd","min","p25","p75","max"))| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| R2_1 | 500 | 0.058 | 0.007 | 0.040 | 0.054 | 0.063 | 0.083 |
| RMSE | 500 | 13.304 | 0.198 | 12.629 | 13.175 | 13.439 | 13.828 |
| MAE | 500 | 10.562 | 0.133 | 10.068 | 10.474 | 10.654 | 10.969 |
| MAPE | 500 | Inf.000 | Inf | Inf | Inf | Inf | |
| THEIL | 500 | 0.522 | 0.006 | 0.505 | 0.517 | 0.526 | 0.541 |
| Um | 500 | 0.000 | 0.000 | 0 | 0 | 0 | 0 |
| Us | 500 | 0.613 | 0.020 | 0.554 | 0.601 | 0.625 | 0.669 |
| Uc | 500 | 0.389 | 0.020 | 0.332 | 0.376 | 0.401 | 0.448 |
| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| educ | 807 | 12.471 | 3.057 | 6.000 | 10.000 | 13.500 | 18.000 |
| cigpric | 807 | 60.300 | 4.738 | 44.004 | 58.141 | 63.179 | 70.129 |
| white | 807 | 0.879 | 0.327 | 0 | 1 | 1 | 1 |
| age | 807 | 41.238 | 17.027 | 17 | 28 | 54 | 88 |
| income | 807 | 19,304.830 | 9,142.958 | 500 | 12,500 | 30,000 | 30,000 |
| cigs | 807 | 8.686 | 13.722 | 0 | 0 | 20 | 80 |
| restaurn | 807 | 0.247 | 0.431 | 0 | 0 | 0 | 1 |
| lincome | 807 | 9.687 | 0.713 | 6.215 | 9.433 | 10.309 | 10.309 |
| agesq | 807 | 1,990.135 | 1,577.166 | 289 | 784 | 2,916 | 7,744 |
| lcigpric | 807 | 4.096 | 0.083 | 3.784 | 4.063 | 4.146 | 4.250 |
bind_rows(Resultados_Performance_1)%>%
stargazer(data,title = "Medidas de Performance de Simulacion",
type = "html",
digits = 3, summary.stat=c("n","mean","sd","min","p25","p75","max"))| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| R2_3 | 500 | 0.039 | 0.018 | 0.002 | 0.027 | 0.051 | 0.099 |
| RMSE | 500 | 13.479 | 0.595 | 11.734 | 13.064 | 13.853 | 15.361 |
| MAE | 500 | 10.728 | 0.290 | 9.834 | 10.528 | 10.904 | 11.626 |
| MAPE | 500 | Inf.000 | Inf | Inf | Inf | Inf | |
| THEIL | 500 | 0.528 | 0.013 | 0.491 | 0.519 | 0.536 | 0.564 |
| Um | 500 | 0.002 | 0.003 | 0.00000 | 0.0003 | 0.003 | 0.021 |
| Us | 500 | 0.597 | 0.051 | 0.476 | 0.563 | 0.636 | 0.751 |
| Uc | 500 | 0.405 | 0.051 | 0.253 | 0.367 | 0.441 | 0.529 |
| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| educ | 807 | 12.471 | 3.057 | 6.000 | 10.000 | 13.500 | 18.000 |
| cigpric | 807 | 60.300 | 4.738 | 44.004 | 58.141 | 63.179 | 70.129 |
| white | 807 | 0.879 | 0.327 | 0 | 1 | 1 | 1 |
| age | 807 | 41.238 | 17.027 | 17 | 28 | 54 | 88 |
| income | 807 | 19,304.830 | 9,142.958 | 500 | 12,500 | 30,000 | 30,000 |
| cigs | 807 | 8.686 | 13.722 | 0 | 0 | 20 | 80 |
| restaurn | 807 | 0.247 | 0.431 | 0 | 0 | 0 | 1 |
| lincome | 807 | 9.687 | 0.713 | 6.215 | 9.433 | 10.309 | 10.309 |
| agesq | 807 | 1,990.135 | 1,577.166 | 289 | 784 | 2,916 | 7,744 |
| lcigpric | 807 | 4.096 | 0.083 | 3.784 | 4.063 | 4.146 | 4.250 |