R: Bayesian model selection and averaging under block-diagonal...
postModeOrtho
R Documentation
Bayesian model selection and averaging under block-diagonal X'X for linear models.
Description
postModeOrtho is for diagonal X'X,
postModeBlockDiag for the more general block-diagonal X'X,
where X is the matrix with predictors.
Both functions return the model of highest posterior probability of any given
size using an efficient search algorithm. This sequence of models includes
the highest posterior probability model (HPM).
Posterior model probabilities, marginal variable inclusion probabilities
and Bayesian model averaging estimates are also provided.
The unknown residual variance is integrated out using an exact deterministic
algorithm of low computational cost (see details in reference).
Matrix with predictors. If an intercept is desired x should include
a column of 1's.
blocks
Factor or integer vector of length ncol(x) indicating the block
that each column in x belongs to.
priorCoef
Prior distribution for the coefficients. Object created
with momprior, imomprior, emomprior or zellnerprior.
priorDelta
Prior on model space. Use modelbbprior() for Beta-Binomial prior,
modelbinomprior(p) for Binomial prior with prior inclusion probability p,
or modelunifprior() for Uniform prior
priorVar
Inverse gamma prior on residual variance, created with igprior()
bma
Set to TRUE to obtain marginal inclusion probabilities and
Bayesian model averaging parameter estimates for each column of x.
maxvars
The search for the HPM is restricted to models with up to maxvars variables
(note: posterior model probabilities and BMA are valid regardless of maxvars)
momcoef
optional argument containing pre-computed coefficients needed to obtain
the marginal likelihood under the pMOM prior.
A first call to postModeBlockDiag returns these coefficients,
thus this argument is useful to speed up successive calls.
Details
The first step is to list a sequence of models with 0,...,maxvars variables which,
under fairly general conditions listed in Papaspiliopoulos & Rossell (2016),
is guaranteed to include the HPM.
Then posterior model probabilities are computed for all these models to determine
the HPM, evaluate the marginal posterior of the residual variance on a grid,
and subsequently compute the marginal density p(y) via adaptive quadrature.
Finally this adaptive grid is used to compute marginal inclusion probabilities
and Bayesian model averaging estimates.
For more details see Papaspiliopoulos & Rossell (2016).
Value
List with elemants
models
data.frame indicating the variables included in the sequence of models found
during the search of the HPM, and their posterior probabilities. The model with highest
posterior probability in this list is guaranteed to be the HPM.
phi
data.frame containing an adaptive grid of phi (residual variance) values and
their marginal posterior density p(phi|y).
logpy
log-marginal density p(y), i.e. normalization constant of p(phi|y).
bma
Marginal posterior inclusion probabilities and Bayesian model averaging estimates
for each column in x.
postmean.model
Coefficient estimates conditional on each of the models in models
momcoef
If a MOM prior was specified in priorCoef, momcoef stores some coefficients needed
to compute its marginal likelihood
Author(s)
David Rossell
References
Papaspiliopoulos O., Rossell D. Scalable Bayesian variable selection
and model averaging under block-orthogonal design. 2016
Examples
#Simulate data
set.seed(1)
p <- 500; n <- 510
x <- scale(matrix(rnorm(n*p),nrow=n,ncol=p),center=TRUE,scale=TRUE)
S <- cov(x)
e <- eigen(cov(x))
x <- t(t(x %*% e$vectors)/sqrt(e$values))
th <- c(rep(0,p-3),c(.5,.75,1)); phi <- 1
y <- x %*% matrix(th,ncol=1) + rnorm(n,sd=sqrt(phi))
#Fit
priorCoef=zellnerprior(tau=n); priorDelta=modelbinomprior(p=1/p); priorVar=igprior(0.01,0.01)
pm.zell <- postModeOrtho(y,x=x,priorCoef=priorCoef,priorDelta=priorDelta,priorVar=priorVar,
bma=TRUE)
#Best models
head(pm.zell$models)
#Posterior probabilities for sequence of models
nvars <- sapply(strsplit(as.character(pm.zell$models$modelid),split=','),length)
plot(nvars,pm.zell$models$pp,ylab='post prob',xlab='number of vars',ylim=0:1,xlim=c(0,50))
#Marginal posterior of phi
plot(pm.zell$phi,type='l',xlab='phi',ylab='p(phi|y)')
#Marginal inclusion prob & BMA estimates
plot(pm.zell$bma$margpp,ylab='Marginal inclusion prob')
plot(pm.zell$bma$coef,ylab='BMA estimate')