getBmaSamples: Get posterior samples from the Bayesian Model Average (BMA)

In hypergsplines: Bayesian model selection with penalised splines and hyper-g prior

Description

Get posterior samples from the Bayesian Model Average (BMA)

Usage

  getBmaSamples(config, logPostProbs, nSamples, modelData,
    mcmc = getMcmc(), computation = getComputation())

Arguments

config	the data frame/matrix with model specifications, e.g. the result from exhaustive
logPostProbs	vector of log posterior probabilites (will be exponentiated and normalized within the function) for the weighting of the models in config
nSamples	number of samples to simulate
modelData	the data necessary for model estimation, which is the result from modelData or glmModelData
mcmc	MCMC options produced by getMcmc, only matters for generalised response models. Then, the burn-in and thinning parameters will be applied for each sampled model.
computation	computation options produced by getComputation, only matters for generalised response models.

Value

A list with samples from the shrinkage hyperparameter, the regression variance, and the (linear and spline) coefficients, analogous to the return value from getSamples or glmGetSamples. The only difference is that “linearCoefs” and “splineCoefs” contain zeroes for samples where the model did not contain that covariate linearly or smoothly. This is necessary to ensure compatibility with getFunctionSamples and getFitSamples. Moreover, the model specifications matrix is appended with columns “postProb” and “sampleFreq”, containing the posterior probability and the sampling frequency, respectively.

Author(s)

Daniel Sabanes Bove [e-mail]

Examples

## get some data
attach(longley)

## get model data
md <- modelData(y=Employed,
                X=cbind(GNP, Armed.Forces),
                gPrior="hyper-g/n")

## get models table
tab <- exhaustive(modelData=md)$models
tab

## get posterior samples from the BMA assuming
## a flat model prior
res <- getBmaSamples(config=tab,
                     logPostProbs=tab$logMargLik,
                     nSamples=1000L,
                     modelData=md)
str(res)

summary(res$t)
hist(res$t, nclass=100)

## now for generalised response:

## get the model data
md <- glmModelData(y=as.numeric(Employed > 64),
                   X=cbind(GNP, Armed.Forces),
                   family=binomial)

## get models table
tab <- exhaustive(modelData=md,
                  computation=
                  getComputation(higherOrderCorrection=FALSE))$models

## get posterior samples from the BMA assuming
## a flat model prior
res <- getBmaSamples(config=tab,
                     logPostProbs=tab$logMargLik,
                     nSamples=1000L,
                     modelData=md,
                     mcmc=
                     getMcmc(burnin=10L,
                             step=1L),
                     computation=
                     getComputation(higherOrderCorrection=FALSE))
str(res)

hist(res$t, nclass=100)
res$config

hypergsplines documentation built on May 2, 2019, 6:14 p.m.