Por medio de este documento se expone la aplicación de modelos de regresión lineal con un enfoque Bayesiano, es decir asumiendo que los parámetros del modelo no corresponden a estimaciones puntuales sino a distribuciones de probabilidad como tal; el ejemplo guía es tomado del capítulo 5 del libro “Learning Bayesian Models with R”, el cual se denomina “Bayesian Regression Models”. Además, se hace alusión también a la aplicación en R respecto a la selección de modelos Bayesianos, cuyo ejemplo fue tomado del manual del package \(BayesVarSel\).

Regresión lineal Bayesiana

Instalación y carga de paquetes

Para implementar los modelos de regresion Bayesiana es necesario cargar el paquete arm.

library(readr)
library(ggplot2)
library(gridExtra)
# install.packages("arm")
library(arm)

Archivo de datos

El dataset utilizado para este ejemplo contiene mediciones de variables relacionadas con la eficiencia energética de un conjunto de edificaciones. En total, hay 768 registros para ocho atributos o características (denotadas por la letra \(X\)) y dos respuestas o resultados (denotados por la letra \(Y\)):

  • \(X_{1}\): compactación relativa.
  • \(X_{2}\): área de la superficie.
  • \(X_{3}\): área de la pared.
  • \(X_{4}\): área del techo.
  • \(X_{5}\): altura total.
  • \(X_{6}\): orientación.
  • \(X_{7}\): área de acristalamiento.
  • \(X_{8}\): distribución del área de acristalamiento.
  • \(Y_{1}\): carga de calentamiento.
  • \(Y_{2}\): carga de refrigeración.

La siguiente es una muestra de dicho dataset:

datos <- read_delim("~/Noveno semestre/Seminario/Exposicion - Ronda 3/datos.csv", ";", escape_double = FALSE, trim_ws = TRUE) # Carga del archivo de datos

datos <- data.frame(datos) # Identificación de la base de datos como un data frame

head(datos) # Muestra
##     X1    X2    X3     X4 X5 X6 X7 X8    Y1    Y2
## 1 0.98 514.5 294.0 110.25  7  2  0  0 15.55 21.33
## 2 0.98 514.5 294.0 110.25  7  3  0  0 15.55 21.33
## 3 0.98 514.5 294.0 110.25  7  4  0  0 15.55 21.33
## 4 0.98 514.5 294.0 110.25  7  5  0  0 15.55 21.33
## 5 0.90 563.5 318.5 122.50  7  2  0  0 20.84 28.28
## 6 0.90 563.5 318.5 122.50  7  3  0  0 21.46 25.38

Dicho conjunto de datos fue extraído del website http://archive.ics.uci.edu/ml/datasets/Energy+efficiency#.

Regresión de la eficiencia energética

El objetivo en este caso es intentar predecir la carga de calentamiento (\(Y_{1}\)) en función de las características disponibles para las edificaciones, utilizando tanto regresión linear ordinaria como Bayesiana.

Análisis descriptivo

Con el objetivo de observar de manera preliminar qué variables podrían contribuir al modelo, se realiza un análisis descriptivo. Para esta cuestión, se presentan algunos gráficos bivariados entre la variable de respuesta y cada una de las variables predictoras, haciendo uso de la librería ggplot2:

attach(datos) # Permite referenciar los nombres de las variables de forma directa

# Gráficos bivariados
g1 <- ggplot(data = datos, aes(x = X1, y = Y1)) + geom_point()
g2 <- ggplot(data = datos, aes(x = X2, y = Y1)) + geom_point()
g3 <- ggplot(data = datos, aes(x = X3, y = Y1)) + geom_point()
g4 <- ggplot(data = datos, aes(x = X4, y = Y1)) + geom_point()
g5 <- ggplot(data = datos, aes(x = X5, y = Y1)) + geom_point()
g6 <- ggplot(data = datos, aes(x = X6, y = Y1)) + geom_point()
g7 <- ggplot(data = datos, aes(x = X7, y = Y1)) + geom_point()
g8 <- ggplot(data = datos, aes(x = X8, y = Y1)) + geom_point()

# Muestra de los gráficos
grid.arrange(g1,g2,g3,g4,g5,g6,g7,g8, nrow = 2, ncol = 4) 

detach(datos) # Desactiva la función attach

De los gráficos anteriores puede observarse que las variables \(X_{6}\) y \(X_{8}\) no presentan mucha influencia sobre la variable de respuesta. No obstante resulta necesario corroborar este hecho, ante lo cual se recurre al coeficiente de correlación de Spearman:

cor.spearman <- cor(datos[,1:8], datos[,9], method = "spearman")

cor.spearman
##            [,1]
## X1  0.622134697
## X2 -0.622134697
## X3  0.471457650
## X4 -0.804027000
## X5  0.861282577
## X6 -0.004163071
## X7  0.322860320
## X8  0.068343464

Debido a que los coeficientes de correlación permiten visualizar que la influencia de las mencionadas variables predictoras sobre la variable de respuesta es despreciable, se toma la decisión de no incluirlas dentro del modelo.

Regresión lineal ordinaria

Como primera medida, se eliminan del dataset las variables \(X_{6}\) y \(X_{8}\), dada la baja correlación presentada respecto a la variable de respuesta:

datos <- datos[, c(1,2,3,4,5,7,9)]

En segundo lugar, se divide la cantidad total de datos, con el objetivo de entrenar el modelo con algunos de ellos y, posteriormente, predecir valores de \(Y_{1}\) considerando la data de prueba:

set.seed(500) # Establecimiento de la semilla
muestra <- sample.int(n = nrow(datos), size = nrow(datos)*0.2, replace = F) # Tamaño de la muestra

datos.test <- datos[muestra,] # Test data
datos.train <- datos[-muestra,] # Training data

x.test <- datos.test[, 1:6] # Variables regresoras en el test data
y.test <- datos.test[, 7] # Variable de respuesta en el test data

Ya con estos resultados, se procede a calcular el modelo de regresión lineal ordinaria:

attach(datos.train)

mod.ols <- lm(Y1 ~ X1 + X2 + X3 + X4 + X5 + X7, data = datos.train) # Modelo lineal

summary(mod.ols) # Resumen del modelo
## 
## Call:
## lm(formula = Y1 ~ X1 + X2 + X3 + X4 + X5 + X7, data = datos.train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.2722  -1.3555   0.0594   1.3286   7.6166 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  97.843380  22.445134   4.359 1.53e-05 ***
## X1          -71.733244  12.108449  -5.924 5.25e-09 ***
## X2           -0.100353   0.020236  -4.959 9.19e-07 ***
## X3            0.066739   0.007997   8.345 4.78e-16 ***
## X4                  NA         NA      NA       NA    
## X5            3.919940   0.398574   9.835  < 2e-16 ***
## X7           20.627695   0.912929  22.595  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.013 on 609 degrees of freedom
## Multiple R-squared:  0.9115, Adjusted R-squared:  0.9108 
## F-statistic:  1254 on 5 and 609 DF,  p-value: < 2.2e-16

Puede observarse entonces que el modelo de regresión explica en gran medida la variabilidad del modelo, aún bajo una penalización por la cantidad de variables regresoras del modelo (\(R^2_{adj} = 0.9108\)); sin embargo, es importante observar que la variable \(X_{4}\) no tiene una estimación del coeficiente asociado dentro del modelo, lo cual se debe a una correlación de -0.9725 con la variable \(X_{5}\), hecho que conlleva al problema de multicolinealidad severa. Como tal correlación es muy alta, el cambio en el \(R^2_{adj}\) no es de importancia si se elimina esta variable dentro del modelo:

mod.ols <- lm(Y1 ~ X1 + X2 + X3 + X5 + X7, data = datos.train) # Modelo lineal

summary(mod.ols) # Resumen del modelo
## 
## Call:
## lm(formula = Y1 ~ X1 + X2 + X3 + X5 + X7, data = datos.train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.2722  -1.3555   0.0594   1.3286   7.6166 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  97.843380  22.445134   4.359 1.53e-05 ***
## X1          -71.733244  12.108449  -5.924 5.25e-09 ***
## X2           -0.100353   0.020236  -4.959 9.19e-07 ***
## X3            0.066739   0.007997   8.345 4.78e-16 ***
## X5            3.919940   0.398574   9.835  < 2e-16 ***
## X7           20.627695   0.912929  22.595  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.013 on 609 degrees of freedom
## Multiple R-squared:  0.9115, Adjusted R-squared:  0.9108 
## F-statistic:  1254 on 5 and 609 DF,  p-value: < 2.2e-16

Ahora bien, con el modelo obtenido de la manera tradicional, se proceden a realizar las predicciones sobre el dataset de prueba:

mod.coef <- mod.ols$coefficients # Coeficientes del modelo

ypred.ols <- predict.lm(mod.ols,x.test,interval = "prediction",se.fit = T) # Valores de predicción considerando el modelo lineal

head(ypred.ols$fit) # Muestra de los intervalos de predicción
##           fit       lwr      upr
## 359 37.336249 31.362331 43.31017
## 411 11.146684  5.204282 17.08909
## 431 15.637038  9.684876 21.58920
## 169  8.052529  2.105563 13.99950
## 236 11.931746  5.992607 17.87089
## 47  10.480114  4.512206 16.44802
yout.ols <- as.data.frame(cbind(y.test,ypred.ols$fit)) # Data frame con los valores de respuesta, la estimación y el intervalo de predicción
ols.upr <- yout.ols$upr # Límite superior del intervalo de predicción
ols.lwr <- yout.ols$lwr # Límite inferior del intervalo de predicción

Regresión lineal Bayesiana

La función bayesglm(), que se encuentra dentro del package arm, permite ajustar un modelo de regresión lineal Bayesiana cuando se especifica que la familia es Gaussiana y la función enlace es la identidad; además, es necesario establecer una distribución prior sobre los parámetros, para lo cual se dejará la configuración por defecto, como se muestra a continuación:

mod.bayes <- bayesglm(Y1 ~ X1 + X2 + X3 + X5 + X7, family=gaussian(link=identity), data=datos.train, prior.df = Inf, prior.mean = 0, prior.scale = NULL, maxit = 10000) # Modelo de regresión lineal Bayesiana

summary(mod.bayes)
## 
## Call:
## bayesglm(formula = Y1 ~ X1 + X2 + X3 + X5 + X7, family = gaussian(link = identity), 
##     data = datos.train, prior.mean = 0, prior.scale = NULL, prior.df = Inf, 
##     maxit = 10000)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -10.2693   -1.3745    0.0597    1.3607    7.6217  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  90.711583  21.843823   4.153 3.75e-05 ***
## X1          -67.829230  11.774667  -5.761 1.33e-08 ***
## X2           -0.094106   0.019724  -4.771 2.30e-06 ***
## X3            0.065193   0.007919   8.233 1.12e-15 ***
## X5            4.004814   0.393807  10.169  < 2e-16 ***
## X7           20.625599   0.912851  22.595  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 9.08027)
## 
##     Null deviance: 62464.0  on 614  degrees of freedom
## Residual deviance:  5529.9  on 609  degrees of freedom
## AIC: 3110
## 
## Number of Fisher Scoring iterations: 5

Se realizan posteriormente las predicciones considerando los datos de prueba:

ypred.bayes <- predict.glm(mod.bayes, newdata = x.test, se.fit = T) # Función para la predicción
head(ypred.bayes$fit) # Muestra de valores predichos con el modelo Bayesiano
##      359      411      431      169      236       47 
## 37.25313 11.10668 15.70432  8.01284 11.96227 10.54792

Ya que se pretende realizar una comparación entre las regresiones ordinaria y Bayesiana, se almacenan también los intervalos de predicción con este último método:

yout.bayes <- as.data.frame(cbind(y.test, ypred.bayes$fit)) # Data frame con los valores de respuesta y los valores predichos
names(yout.bayes) <- c("ytest","fit") # Cambio de nombre de las columnas del data frame
critval <- 1.96 # Valor crítico
bayes.upr <- ypred.bayes$fit + critval * ypred.bayes$se.fit # Límite superior de los intervalos
bayes.lwr <- ypred.bayes$fit - critval * ypred.bayes$se.fit # Límite inferior de los intervalos

Para la comparación de las predicciones mediante ambos modelos fueron realizados los siguientes plots:

p.ols <- ggplot(data = yout.ols, aes(x = y.test, y = fit)) + geom_point() + ggtitle("Regresión ordinaria") + labs(x = "Y-Prueba", y = "Y-Predicción") # Valores de prueba vs. valores predichos - OLS

p1 <- p.ols + geom_errorbar(ymin = ols.lwr, ymax = ols.upr) # Representación de los intervalos - OLS

p.bayes <- ggplot(data = yout.bayes, aes(x = y.test, y = fit)) + geom_point() + ggtitle("Regresión Bayesiana") + labs(x = "Y-Prueba", y = "Y-Predicción") # Valores de prueba vs. valores predichos - Bayes

p2 <- p.bayes + geom_errorbar(ymin = bayes.lwr, ymax = bayes.upr) # Representación de los intervalos - OLS

grid.arrange(p1,p2,ncol = 2)

Los gráficos anteriores permiten concluir que, si bien las estimaciones puntuales obtenidas son muy similares, los intervalos de predicción obtenidos con un nivel de confianza del 95% son mucho menos amplios que los obtenidos por regresión ordinaria. En el modelo Bayesiano, las predicciones son realizadas a partir de un conjunto de valores muestreados de la distribución posteriori y promediados para obtener la predicción y el intervalo final.

Simulación de la distribución posteriori

El package arm también permite encontrar la distribución posteriori de los parámetros del modelo, esto a través de la función sim.

sim(mod.bayes) # Funcionamiento de la función sim
## An object of class "sim"
## Slot "coef":
##        (Intercept)        X1          X2         X3       X5       X7
##   [1,]   103.39755 -73.70889 -0.10559772 0.06910849 3.699331 20.44815
##   [2,]    87.36962 -65.64984 -0.09178433 0.06641236 3.969437 20.32285
##   [3,]   105.10555 -75.37983 -0.10684450 0.06853857 3.813313 20.02428
##   [4,]    63.17525 -56.36951 -0.06583197 0.05135384 4.822557 20.59745
##   [5,]   119.28189 -83.09018 -0.12249297 0.08068719 3.528230 19.32815
##   [6,]    86.18015 -64.90446 -0.08757923 0.05551079 4.247401 20.06327
##   [7,]    77.41976 -61.02695 -0.08134565 0.05824920 4.395229 19.50667
##   [8,]   100.63567 -72.70757 -0.10307129 0.06785757 3.815561 21.35564
##   [9,]   115.70185 -79.10314 -0.11974058 0.07870961 3.338312 21.31690
##  [10,]    72.58481 -59.72585 -0.07768565 0.06168708 4.363099 21.08656
##  [11,]   107.79226 -78.36656 -0.10604241 0.06078075 4.093314 19.80828
##  [12,]   117.09395 -84.52900 -0.11388914 0.06592160 3.953655 18.84920
##  [13,]    95.33459 -66.11730 -0.10346094 0.07463592 3.512167 19.65830
##  [14,]   123.73394 -84.71907 -0.12219789 0.06821277 3.525705 21.19060
##  [15,]   107.96242 -75.31409 -0.11417173 0.08017160 3.450894 20.91814
##  [16,]    80.56332 -62.23542 -0.08395088 0.06151519 4.065712 20.29283
##  [17,]    72.83510 -58.08249 -0.07734452 0.05993022 4.177150 20.35756
##  [18,]    78.54909 -62.65066 -0.08436884 0.06978333 4.069071 20.04102
##  [19,]    68.44835 -52.16283 -0.07872180 0.06593740 3.925280 20.82341
##  [20,]    91.60098 -67.06400 -0.09719045 0.06935183 3.877980 19.81928
##  [21,]   118.69308 -82.08063 -0.11884395 0.06664869 3.769781 21.21989
##  [22,]   110.47621 -75.12668 -0.11432090 0.07224982 3.473581 21.47444
##  [23,]    52.43276 -47.50838 -0.06066807 0.05955312 4.405989 21.06610
##  [24,]   104.66089 -71.43924 -0.11040055 0.07452626 3.426702 19.85387
##  [25,]   104.14072 -76.39527 -0.10387589 0.06423721 4.033807 20.66715
##  [26,]    57.84678 -49.32094 -0.06499885 0.05675924 4.383209 20.91231
##  [27,]   110.86808 -78.59680 -0.11208530 0.07174360 3.652846 19.92449
##  [28,]    54.75569 -50.07884 -0.06105539 0.05648707 4.582410 20.01008
##  [29,]   131.47689 -92.28007 -0.12819143 0.06787651 3.868705 23.15749
##  [30,]    40.02921 -39.86709 -0.05169020 0.06159091 4.410089 20.64996
##  [31,]    94.51042 -69.40421 -0.09707749 0.06349404 3.982673 20.19606
##  [32,]   103.81661 -78.39739 -0.10150721 0.05664963 4.458221 21.15414
##  [33,]    82.20825 -65.67207 -0.08439815 0.06032038 4.343280 20.78743
##  [34,]    69.69380 -56.70807 -0.07624805 0.06210698 4.330251 20.20732
##  [35,]    85.16926 -67.19560 -0.08735371 0.06277583 4.257495 20.27847
##  [36,]   141.87627 -97.14008 -0.13783401 0.07445694 3.566160 20.35153
##  [37,]    93.02624 -66.26817 -0.09719604 0.06603815 3.703983 20.76477
##  [38,]    58.12360 -53.99218 -0.06098741 0.05174345 4.792691 19.91470
##  [39,]    90.83148 -66.32069 -0.09537419 0.06566086 3.872931 19.83975
##  [40,]    81.81662 -61.24598 -0.08761597 0.06341888 3.997390 21.41153
##  [41,]   100.37549 -71.58876 -0.10550508 0.07544842 3.565846 20.49204
##  [42,]    71.69836 -54.59910 -0.08120273 0.06785554 3.881149 21.01065
##  [43,]   111.54524 -77.20075 -0.11387100 0.07238131 3.471308 20.68894
##  [44,]   124.57094 -83.66557 -0.12693936 0.07875019 3.312875 19.27527
##  [45,]    31.64264 -37.07222 -0.03808879 0.04477801 4.929331 19.16890
##  [46,]    80.34541 -63.17517 -0.08589890 0.06702331 4.143430 20.69802
##  [47,]    98.43291 -74.41202 -0.09946197 0.06551471 4.189944 20.18731
##  [48,]    65.07077 -53.88567 -0.07210559 0.06106099 4.286568 20.73121
##  [49,]   108.41396 -80.43753 -0.10401619 0.05599640 4.334662 20.13302
##  [50,]   120.13454 -82.62189 -0.12252622 0.07599300 3.542614 20.70046
##  [51,]    98.82640 -74.19033 -0.09836526 0.06357960 4.048408 20.30950
##  [52,]    99.38544 -69.61544 -0.10512005 0.07251401 3.595367 20.88599
##  [53,]    96.98259 -70.24118 -0.10035498 0.06606245 3.868782 21.63091
##  [54,]    80.99536 -62.30385 -0.08510118 0.06196452 4.117937 20.28501
##  [55,]    92.17522 -64.09835 -0.10079115 0.07565769 3.361585 21.34414
##  [56,]    91.74237 -69.39976 -0.09544905 0.06815154 4.031597 20.15778
##  [57,]    64.94797 -50.89998 -0.07358319 0.06125589 4.046150 21.64161
##  [58,]    98.23457 -72.07294 -0.09997730 0.06672528 3.808502 21.85290
##  [59,]   116.77314 -83.59936 -0.11644362 0.07098140 3.771558 21.99024
##  [60,]    96.89332 -72.09307 -0.09782697 0.06069195 4.084443 22.89137
##  [61,]    37.94541 -37.17895 -0.05102165 0.05878559 4.502904 20.47600
##  [62,]    68.47065 -53.95056 -0.07718325 0.06583041 4.138351 18.27015
##  [63,]   143.85772 -95.39510 -0.14132589 0.07748839 3.205661 19.13487
##  [64,]   111.58667 -80.95348 -0.10861057 0.06208697 3.982040 21.16321
##  [65,]    64.98980 -55.96345 -0.06696448 0.05071172 4.632210 19.58580
##  [66,]    83.20272 -60.26517 -0.08952698 0.06111147 3.931935 21.55622
##  [67,]    88.63194 -67.15987 -0.09204807 0.06499191 4.073927 20.93287
##  [68,]    67.94466 -54.81373 -0.07528739 0.06550778 4.058482 19.28869
##  [69,]   117.86105 -80.89067 -0.12099472 0.07692585 3.402348 21.12234
##  [70,]    74.44077 -55.57573 -0.08467342 0.07396795 3.581348 21.25239
##  [71,]    70.13677 -56.30739 -0.07690707 0.06266176 4.204149 20.02449
##  [72,]    94.35499 -70.08041 -0.09877894 0.07113372 3.890819 20.13628
##  [73,]    60.30747 -53.37596 -0.06514193 0.05658430 4.486748 21.19865
##  [74,]    92.49116 -70.05256 -0.09142682 0.05426205 4.330332 20.29619
##  [75,]    54.50038 -45.73776 -0.06226916 0.05293534 4.335648 21.40574
##  [76,]   120.86646 -84.74433 -0.12147689 0.07593732 3.560156 20.43313
##  [77,]    97.74044 -73.43965 -0.09743820 0.06002337 4.202096 20.84473
##  [78,]    43.09671 -43.07523 -0.05056698 0.05082529 4.765482 20.01598
##  [79,]    58.11314 -50.22523 -0.06696895 0.06121705 4.379422 21.84887
##  [80,]    60.63978 -51.06433 -0.06685337 0.05605800 4.361983 19.70372
##  [81,]    70.80871 -57.15174 -0.07700156 0.06126981 4.255911 20.26774
##  [82,]    89.29438 -66.76669 -0.09302659 0.06677700 3.913590 20.21708
##  [83,]    81.05936 -62.94251 -0.08497644 0.06001351 4.281701 21.12668
##  [84,]    90.84181 -67.89544 -0.09373946 0.06508472 3.993434 19.39334
##  [85,]    87.66277 -61.95980 -0.09780054 0.07836488 3.388584 20.41318
##  [86,]    75.09826 -58.36332 -0.08003629 0.05973220 4.153280 19.90295
##  [87,]    66.44539 -56.97091 -0.06634106 0.04412930 4.769325 21.51168
##  [88,]   129.42294 -85.90055 -0.13147623 0.08022423 3.207304 19.12495
##  [89,]    99.12088 -74.76028 -0.10101533 0.06986586 4.006849 20.15396
##  [90,]    78.77698 -64.35359 -0.08260646 0.06355442 4.366933 21.37591
##  [91,]    92.37938 -65.03143 -0.10008457 0.07550425 3.421464 20.97958
##  [92,]    97.50105 -70.19962 -0.10064546 0.06797807 3.746567 21.25163
##  [93,]   102.37146 -72.75591 -0.10675967 0.07357297 3.616365 19.81418
##  [94,]    92.49471 -67.30911 -0.09896572 0.07295317 3.786947 18.97502
##  [95,]    77.51311 -61.35049 -0.08100444 0.05617383 4.399097 20.94769
##  [96,]    82.40414 -64.50502 -0.08408112 0.05827015 4.305937 19.22194
##  [97,]    94.23918 -69.40487 -0.09906945 0.07017039 3.868206 21.27457
##  [98,]   120.98056 -80.02374 -0.12719071 0.08349048 3.091233 21.50396
##  [99,]   103.29979 -75.16853 -0.10635787 0.07094130 3.904238 20.59513
## [100,]    50.10850 -47.05081 -0.05575877 0.04940676 4.733245 21.15752
## 
## Slot "sigma":
##   [1] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##   [9] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [17] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [25] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [33] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [41] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [49] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [57] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [65] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [73] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [81] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [89] 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213 3.013213
##  [97] 3.013213 3.013213 3.013213 3.013213
posterior.bayes <- as.data.frame(coef(sim(mod.bayes))) # Almacenamiento de los resultados a posteriori

head(posterior.bayes) # Muestra del data frame anterior
##   (Intercept)        X1          X2         X3       X5       X7
## 1    93.97765 -66.90550 -0.10020745 0.07243300 3.639024 18.57139
## 2    63.00396 -51.46422 -0.07013515 0.05770830 4.261702 20.86155
## 3    74.79306 -58.88003 -0.07988420 0.06090446 4.281314 18.93225
## 4    57.39910 -49.59851 -0.06374356 0.05708519 4.390134 19.23038
## 5    95.04840 -69.12629 -0.09849915 0.06563467 3.892951 20.61818
## 6    84.99291 -64.77459 -0.09097583 0.06990978 3.984047 21.21667

Finalmente, se muestran las distribuciones posteriores estimadas que resultan para los parámetros asociados a cada una de las covariables del modelo:

attach(posterior.bayes)

h0 <- ggplot(data = posterior.bayes, aes(x = posterior.bayes$`(Intercept)`)) + geom_histogram() + labs(x = "Intercepto")
h1 <- ggplot(data = posterior.bayes, aes(x = X1)) + geom_histogram()
h2 <- ggplot(data = posterior.bayes, aes(x = X2)) + geom_histogram()
h3 <- ggplot(data = posterior.bayes, aes(x = X3)) + geom_histogram()
h5 <- ggplot(data = posterior.bayes, aes(x = X5)) + geom_histogram()
h7 <- ggplot(data = posterior.bayes, aes(x = X7)) + geom_histogram()

grid.arrange(h0,h1,h2,h3,h5,h7, nrow = 2, ncol = 3)

detach(posterior.bayes)

De esta manera, se concluye que la regresión lineal Bayesiana permite mayor exactitud respecto a los intervalos de predicción obtenidos y, además, brinda información agregada en cuanto se obtienen distribuciones de probabilidad para los parámetros del modelo y no únicamente estimaciones puntuales.

Selección de modelos Bayesianos

Base de datos de análisis del crimen

Los criminólogos se interesan frecuentemente por el efecto que tienen los regímenes de castigo sobre los índices de delincuencia. Se trabajó entonces con la base de datos UScrime que contiene 47 estados de los EE.UU en los cuales se midieron 16 variables:

  • M: porcentaje de varones de 14 a 24 años.
  • So: variable indicadora para un estado del sur.
  • Ed: media de años de escolaridad.
  • Po1: los gastos de la policía en 1960.
  • Po2: los gastos de la policía en 1959.
  • LF: tasa de participación en la fuerza de trabajo.
  • M.F: número de hombres por cada 1000 mujeres.
  • Pop: población del estado.
  • NW: número de no blancos por cada 1000 personas.
  • U1: tasa de desempleo de los varones urbanos de 14 a 24 años.
  • U2: tasa de desempleo de los varones urbanos de 35 a 39 años.
  • GDP: producto interno bruto per cápita.
  • Ineq: desigualdad de ingresos.
  • Prob: probabilidad de encarcelamiento.
  • Time: tiempo medio de permanencia en las prisiones estatales.
  • y: tasa de delitos en una categoría particular por cabeza de la población.

Se tiene entonces como finalidad el predecir la tasa de delitos a partir de las demás covariables presentes en la base de datos.

En este sentido, se hará uso del paquete \(BayesVarSel\) para realizar la selección del mejor modelo a partir de un conjunto de modelos candidatos, se presentaron 2 escenarios, cuando las probabilidades a priori eran igual y cuando eran diferentes:

# install.packages("BayesVarSel")
library(BayesVarSel)

#Cargar los datos
data(UScrime)

attach(UScrime)

#Modelos candidatos para seleccionar el mejor

model1<- y ~ 1 + Prob
model2<- y ~ 1 + Prob + Time
model3<- y ~ 1 + Prob + + Time + Po1 + Po2
model4<- y ~ 1 + Prob + So
model5<- y ~ .


#Probabilidades a priori iguales
   
crime.BF<- Btest(models=list(basemodel=model1,
                             ProbTimemodel=model2,
                             ProbPolmodel=model3,
                             ProbSomodel=model4,
                             fullmodel=model5), data=UScrime)

crime.BF$BFi0
##     basemodel.to.basemodel ProbTimemodel.to.basemodel 
##                  1.0000000                  0.1306429 
##  ProbPolmodel.to.basemodel   ProbSomodel.to.basemodel 
##                305.5100305                  0.2518167 
##     fullmodel.to.basemodel 
##               8403.4002298
crime.BF$PostProbi
##     basemodel ProbTimemodel  ProbPolmodel   ProbSomodel     fullmodel 
##  1.148067e-04  1.499868e-05  3.507460e-02  2.891025e-05  9.647667e-01
#Probabilidades a priori diferentes

crime.BF2<- Btest(models=list(basemodel=model1, ProbTimemodel=model2,
ProbPolmodel=model3, ProbSomodel=model4, fullmodel=model5),
data=UScrime, prior.models = "User", priorprobs=list(basemodel=1/8,
ProbTimemodel=1/8, ProbPolmodel=1/2, ProbSomodel=1/8, fullmodel=1/8))

crime.BF2$BFi0
##     basemodel.to.basemodel ProbTimemodel.to.basemodel 
##                  1.0000000                  0.1306429 
##  ProbPolmodel.to.basemodel   ProbSomodel.to.basemodel 
##                305.5100305                  0.2518167 
##     fullmodel.to.basemodel 
##               8403.4002298
crime.BF2$PostProbi
##     basemodel ProbTimemodel  ProbPolmodel   ProbSomodel     fullmodel 
##  1.038764e-04  1.357072e-05  1.269412e-01  2.615782e-05  8.729152e-01

A partir de las probabilidades posteriores fue posible observar en los dos casos que el modelo completo es quien sería elegido, ya que presenta la mayor probabilidad entre todos los que se compararon.