Partiendo del ejemplo visto en clases.
Carga de datos.
library(mlbench)
data(BostonHousing)
modelo_estimado<-lm(formula = medv~., data=BostonHousing)Previamente a realizar la simulación deben definirse las funciones (descomposición de Theil).
El paquete DescTools se utiliza para la descripción de datos.
library(DescTools)
# Bias Proportion, proporción de sesgo UM
Um<-function(pronosticado,observado){
library(DescTools)
((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado)
}
# Variance Proportion, proporción de avrianza US
Us<-function(pronosticado,observado){
library(DescTools)
((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
# Covariance Proportion, proporción de covarianza UC
Uc<-function(pronosticado,observado){
library(DescTools)
(2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}
#Coeficiente U de Theil
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) # Activa 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
# Formula a usar lm(formula = medv~.
df <- BostonHousing
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<- df$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 <- df[muestras[[j]], ]
Datos_Prueba <- df[-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)
)
} Para los datos de entrenamiento.
bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo BostonHousing",
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 |
Interpretación de performance.
Los regresores del modelo son capaces de explicar el 74.3% de la varianza de la variable endógena. Por otra parte, el ajuste mínimo es de 71.3% y el máximo poder explicativo es de 79.4%.
El error porcentual promedio (MAPE) del modelo es de 16.387%. El error mínimo porcentual promedio es de 14.813% y el máximo es de 17.691%. Los mejores modelos se encuentran en el 25% inferior de la distribución pues el error promedio va de 14.813% a 16.085%. El procentaje de modelos que tienen un error máximo promedio de 16.718 es de 75%.
Los mejores modelos se encuentran en la cuota inferior de la distribución dado que presentan los niveles de error más bajo.
En cuanto a las diferencias absolutas entre la predicción y la observación real, error absoluto medio (MAE), se tiene un error mínimo de $2.905 y un error máximo de $3.512, son valores bajos que indican una buena capacidad predictiva. De la misma manera, los mejores modelos se encuentran en la cuota inferior hasta el 25% de la distribución pues el error no sobrepasa los $3.204.
En términos de error cuadrático medio (RMSE) el modelo presenta pocos errores extremos pues en el mejor de los casos se tiene un error mínimo de $4.177 y un error máximo de $4.948.
El coeficiente de Theil expresa un buen desempeño pues en promedio es de 0.096, no existe mucha diferencia entre los datos reales y los pronosticados. En el extremo mínimo se tiene un valor de 0.097 y un valor máximo de 0.102.
El coeficiente de sesgo (Um) muestra que no hay diferencia entre el valor real y el pronosticado.
En términos de proporción de varianza (Us), los valores reales y estimados de la varianza muestran un nivel de semejanza alto 0.074, es decir la simulación es cercana a la realidad.
En términos de proporción de covarianza (Uc), el valor real y el pronosticado muestran una correlación de 0.928, se trata de una correlación alta.
Resultados para datos de simulación.
library(dplyr)
bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "html",
digits = 3,
table.placement = "H")| 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 |
Interpretación de performance.
El conjunto de regresores del modelo explica en promedio el 72.3% de la varianza de la variable endógena. En el mejor de los casos tiene un poder explicativo máximo de 84.0%, un mejor ajuste que para los datos de entrenamiento.
En términos de error cuadrático medio (RMSE) el modelo presenta pocos errores extremos pues en el mejor de los casos se tiene un error mínimo de $3.465 y un error máximo de $6.961. Sobrepasa a los valores de la data de entrenamiento.
En cuanto a las diferencias absolutas entre la predicción y la observación real, error absoluto medio (MAE), se tiene un error mínimo de $2.633 y un error máximo de $4.492, son valores mayores a los de la data de entrenamiento. De la misma manera, los mejores modelos se encuentran en la cuota inferior hasta el 25% de la distribución pues el error no sobrepasa los $3.216.
El error promedio del modelo es de 17.197%. En términos del error porcentual absoluto medio (MAPE), el poder predictivo del modelo ronda el 12.875% hasta el 23.137% de margen de error.
El coeficiente de sesgo (Um) muestra que la diferencia entre el valor real y el pronosticado es de 0.011, es bastante bajo.
En términos de proporción de varianza (Us), los valores reales y estimados de la varianza muestran un nivel de semejanza bastante alto 0.081, es decir la simulación es cercana a la realidad. Y es mayor al valor obtenido de la data de entrenamiento.
En términos de proporción de covarianza (Uc), el valor real y el pronosticado muestran una correlación de 0.918, se trata de una correlación alta.