Buscando una distribución a priori “no informativa” en el contexto de selección de modelos

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}\).

La distribución a priori g.

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:

Calculando el factor Bayes

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}}}\).

Selección del valor de \(g\)

Hay dos métodos típicos:

Ejemplo: prediciendo el consumo de gasolina de coches

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

Mirando la relación del consumo con las demás variables

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')

Inferencia y selección de modelos con métodos frecuentistas

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

Inferencia y selección con métodos bayesianos

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.