Hemos visto en clase que en el modelo lineal \(y_i = \sum_{j=1}^k \beta_j x_{ij} + \epsilon_i\), donde \(\epsilon_i|\phi \sim \textrm{Normal}(0,1/\phi)\), podemos emplear la distribución a priori conjugada para \(\mathbf{\beta}\), es decir \(\mathbf{\beta}|\phi \sim \textrm{Normal}\left({\bf m},\frac{1}{\phi} {\bf V} \right)\) con una distribución a priori gamma para \(\phi\).
En el contexto de comparación de modelos, es razonable emplear una distribución impropia para \(\phi\), es decir: \(f(\alpha,\phi) \propto \frac{1}{\phi}\). Observamos que en este contexto, no invalida el factor Bayes, ya que el parámetro y la a priori son iguales en todos los modelos de regresión que empleamos.
Para fijar la a priori para \(\mathbf{\beta}\), es también natural suponer que \({\bf m} = {\bf 0}\) que significa que, en medio, pensamos que los regresores no van a tener efecto. El problema es ¿cómo elegimos la matriz \({\bf V}\).
Para la a priori g, fijamos \({\bf V} = g ({\bf X}^T{\bf X})^{-1}\).
Luego, la distribución a priori de \(\mathbf{\beta}|\phi,\textrm{datos}\) es \(\textrm{Normal}\left(\frac{g}{1+g}\hat{\mathbf{\beta}}, \frac{g}{(g+1) \phi}({\bf X}^T{\bf X})^{-1}\right)\).
Observamos que:
Para comparar un modelo \(\mathcal{M}\) con \(k\) parámetros con el modelo nulo sin regresores, \(\mathcal{M}_0\), se puede demostrar que si se emplea la a priori g, el factor Bayes es: \(B^k_0 = \frac{(1+g)^{\frac{n-k-1}{2}}}{(1+g (1 - R^2_k))^{\frac{n-1}{2}}}\).
Hay dos métodos típicos:
contiene datos de las características de varios coches y su consumo medio de gasolina.
rm(list=ls())
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
data(mtcars)
help(mtcars)
## starting httpd help server ...
## done
head(mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
La mayoria de las variables están linealmente relacionadas con mpg. Observamos que y son factores.
plot(mtcars$cyl,mtcars$mpg,xlab='cyl',ylab='mpg')
plot(mtcars$disp,mtcars$mpg,xlab='disp',ylab='mpg')
plot(mtcars$hp,mtcars$mpg,xlab='hp',ylab='mpg')
plot(mtcars$drat,mtcars$mpg,xlab='drat',ylab='mpg')
plot(mtcars$wt,mtcars$mpg,xlab='wt',ylab='mpg')
plot(mtcars$qsec,mtcars$mpg,xlab='qsec',ylab='mpg')
plot(mtcars$vs,mtcars$mpg,xlab='vs',ylab='mpg')
plot(mtcars$am,mtcars$mpg,xlab='am',ylab='mpg')
plot(mtcars$gear,mtcars$mpg,xlab='gear',ylab='mpg')
Si ajustamos el modelo completo, vemos que muchas variables no son significativas.
summary(lm(mpg ~ cyl+disp+hp+drat+wt+qsec+factor(vs)+factor(am)+gear+carb, data = mtcars))$coef
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.30337416 18.71788443 0.6573058 0.51812440
## cyl -0.11144048 1.04502336 -0.1066392 0.91608738
## disp 0.01333524 0.01785750 0.7467585 0.46348865
## hp -0.02148212 0.02176858 -0.9868407 0.33495531
## drat 0.78711097 1.63537307 0.4813036 0.63527790
## wt -3.71530393 1.89441430 -1.9611887 0.06325215
## qsec 0.82104075 0.73084480 1.1234133 0.27394127
## factor(vs)1 0.31776281 2.10450861 0.1509915 0.88142347
## factor(am)1 2.52022689 2.05665055 1.2254035 0.23398971
## gear 0.65541302 1.49325996 0.4389142 0.66520643
## carb -0.19941925 0.82875250 -0.2406258 0.81217871
Hay muchos posibles modelos. Podemos seleccionar variables con AIC o BIC y un procedimiento “stepwise”. Primero usamos AIC.
library(MASS)
## Warning: package 'MASS' was built under R version 3.5.2
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
fit <- lm(mpg ~ cyl+disp+hp+drat+wt+qsec+factor(vs)+factor(am)+gear+carb, data = mtcars)
stepA <- stepAIC(fit, direction="both", trace=FALSE)
summary(stepA)$coeff
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.617781 6.9595930 1.381946 1.779152e-01
## wt -3.916504 0.7112016 -5.506882 6.952711e-06
## qsec 1.225886 0.2886696 4.246676 2.161737e-04
## factor(am)1 2.935837 1.4109045 2.080819 4.671551e-02
El modelo seleccionado con AIC incluye , y .
Ahora usamos BIC, que implementamos incluyendo \(k=log(n)\) en el código anterior.
stepB <- stepAIC(fit, direction="both", trace=FALSE,k=dim(mtcars)[1])
summary(stepB)$coeff
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.285126 1.877627 19.857575 8.241799e-19
## wt -5.344472 0.559101 -9.559044 1.293959e-10
El modelo seleccionado con BIC incluye sólo
Podemos hacer lo mismo con el paquete . Primero lo hacemos con \(g = n\).
library(BAS)
## Warning: package 'BAS' was built under R version 3.5.3
n <- dim(mtcars)[1]
bfit <- bas.lm(mpg ~ cyl+disp+hp+drat+wt+qsec+factor(vs)+factor(am)+gear+carb, data = mtcars,prior="g-prior",alpha=n)
summary(bfit)
## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.0000000 1.00000 1.0000000 1.0000000 1.0000000
## cyl 0.3691409 1.00000 0.0000000 0.0000000 0.0000000
## disp 0.1529221 0.00000 0.0000000 0.0000000 0.0000000
## hp 0.3488517 0.00000 1.0000000 0.0000000 0.0000000
## drat 0.1403017 0.00000 0.0000000 0.0000000 0.0000000
## wt 0.9231081 1.00000 1.0000000 1.0000000 1.0000000
## qsec 0.3524140 0.00000 0.0000000 1.0000000 1.0000000
## factor(vs)1 0.1314339 0.00000 0.0000000 0.0000000 0.0000000
## factor(am)1 0.2414594 0.00000 0.0000000 0.0000000 1.0000000
## gear 0.1375086 0.00000 0.0000000 0.0000000 0.0000000
## carb 0.2065967 0.00000 0.0000000 0.0000000 0.0000000
## BF NA 1.00000 0.7686289 0.7474264 0.8418717
## PostProbs NA 0.12820 0.0985000 0.0958000 0.0405000
## R2 NA 0.83020 0.8268000 0.8264000 0.8497000
## dim NA 3.00000 3.0000000 3.0000000 4.0000000
## logmarg NA 21.84769 21.5845436 21.5565712 21.6755631
## model 5
## Intercept 1.000000
## cyl 1.000000
## disp 0.000000
## hp 1.000000
## drat 0.000000
## wt 1.000000
## qsec 0.000000
## factor(vs)1 0.000000
## factor(am)1 0.000000
## gear 0.000000
## carb 0.000000
## BF 0.487542
## PostProbs 0.023400
## R2 0.843100
## dim 4.000000
## logmarg 21.129312
Ahora lo hacemos con una distribución a priori bayesiana jerárquica para \(g\).
bfit2 <- bas.lm(mpg ~ cyl+disp+hp+drat+wt+qsec+factor(vs)+factor(am)+gear+carb, data = mtcars,prior="JZS")
summary(bfit2)
## P(B != 0 | Y) model 1 model 2 model 3 model 4
## Intercept 1.0000000 1.00000 1.0000000 1.0000000 1.0000000
## cyl 0.3665729 1.00000 0.0000000 0.0000000 0.0000000
## disp 0.1259293 0.00000 0.0000000 0.0000000 0.0000000
## hp 0.3300561 0.00000 1.0000000 0.0000000 0.0000000
## drat 0.1133176 0.00000 0.0000000 0.0000000 0.0000000
## wt 0.9288863 1.00000 1.0000000 1.0000000 1.0000000
## qsec 0.3320197 0.00000 0.0000000 1.0000000 1.0000000
## factor(vs)1 0.1096095 0.00000 0.0000000 0.0000000 0.0000000
## factor(am)1 0.1992327 0.00000 0.0000000 0.0000000 1.0000000
## gear 0.1105780 0.00000 0.0000000 0.0000000 0.0000000
## carb 0.1695350 0.00000 0.0000000 0.0000000 0.0000000
## BF NA 1.00000 0.7544591 0.7322344 0.6464464
## PostProbs NA 0.16130 0.1217000 0.1181000 0.0391000
## R2 NA 0.83020 0.8268000 0.8264000 0.8497000
## dim NA 3.00000 3.0000000 3.0000000 4.0000000
## logmarg NA 21.50092 21.2191645 21.1892641 21.0646538
## model 5
## Intercept 1.0000000
## cyl 0.0000000
## disp 0.0000000
## hp 0.0000000
## drat 0.0000000
## wt 1.0000000
## qsec 0.0000000
## factor(vs)1 0.0000000
## factor(am)1 0.0000000
## gear 0.0000000
## carb 0.0000000
## BF 0.0442609
## PostProbs 0.0321000
## R2 0.7528000
## dim 2.0000000
## logmarg 18.3832651
Observamos que en ambos casos, el mejor modelo es distinto al modelo seleccionado con el AIC, que queda como el cuarto modelo entre los bayesianos. El modelo seleccionado con BIC es el quinto modelo con la a priori bayesiana jerárquica.