sampler: sampler

In certifiedwaif/blma: Bayesian Linear Model Averaging

sampler

Description

sampler

Usage

sampler(iterations, vy, mX, prior = "BIC", modelprior = "uniform",
modelpriorvec_in=NULL, cores=1)

Arguments

iterations	The number of iterations to run the MCMC sampler for
vy_in	Vector of responses
mX_in	The matrix of covariates which may or may not be included in each model
prior	– the choice of mixture $g$-prior used to perform Bayesian model averaging. The choices available include: "BIC"– the Bayesian information criterion obtained by using the cake prior of Ormerod et al. (2017). "ZE"– special case of the prior structure described by Maruyama and George (2011). "liang_g1"– the mixture g-prior of Liang et al. (2008) with prior hyperparameter a=3 evaluated directly using Equation (10) of Greenaway and Ormerod (2018) where the Gaussian hypergeometric function is evaluated using the gsl library. Note: this option can lead to numerical problems and is only meant to be used for comparative purposes. "liang_g2"– the mixture g-prior of Liang et al. (2008) with prior hyperparameter a=3 evaluated directly using Equation (11) of Greenaway and Ormerod (2018). "liang_g_n_appell"– the mixture g/n-prior of Liang et al. (2008) with prior hyperparameter a=3 evaluated using the appell R package. "liang_g_approx"– the mixture g/n-prior of Liang et al. (2008) with prior hyperparameter eqna=3 using the approximation Equation (15) of Greenaway and Ormerod (2018) for model with more than two covariates and numerical quadrature (see below) for models with one or two covariates. "liang_g_n_quad"– the mixture g/n-prior of Liang et al. (2008) with prior hyperparameter eqna=3 evaluated using a composite trapezoid rule. "robust_bayarri1"– the robust prior of Bayarri et al. (2012) using default prior hyper choices evaluated directly using Equation (18) of Greenaway and Ormerod (2018) with the gsl library. "robust_bayarri2"– the robust prior of Bayarri et al. (2012) using default prior hyper choices evaluated directly using Equation (19) of Greenaway and Ormerod (2018).
modelprior	The model prior to use. The choices of model prior are "uniform", "beta-binomial" or "bernoulli". The choice of model prior dictates the meaning of the parameter modelpriorvec.
cores	The number of cores to use. Defaults to 1.
modelpriorvec	If modelprior is "uniform", then the modelpriorvec is ignored and can be null. If modelprior is "beta-binomial" then modelpriorvec should be length 2 with the first element containing alpha hyperparameter for the beta prior and the second element containing the beta hyperparameter for beta prior. If modelprior is "bernoulli", then modelpriorvec must be of the same length as the number columns in mX. Each element i of modelpriorvec contains the prior probability of the the ith covariate being included in the model.

Value

The object returned is a list containing:

• "mGamma"– An iterations by p binary matrix containing the sampled models.

• "vinclusion_prob"– The vector of inclusion probabilities.

• "vlogBF"– The vector of logs of the Bayes Factors of the models in mGamma.

References

Bayarri, M. J., Berger, J. O., Forte, A., Garcia-Donato, G., 2012. Criteria for Bayesian model choice with application to variable selection. Annals of Statistics 40 (3), 1550-1577.

Greenaway, M. J., J. T. Ormerod (2018) Numerical aspects of Bayesian linear models averaging using mixture g-priors.

Liang, F., Paulo, R., Molina, G., Clyde, M. a., Berger, J. O., 2008. Mixtures of g priors for Bayesian variable selection. Journal of the American Statistical Association 103 (481), 410-423.

Ormerod, J. T., Stewart, M., Yu, W., Romanes, S. E., 2017. Bayesian hypothesis tests with diffuse priors: Can we have our cake and eat it too?

Examples

mD <- MASS::UScrime
notlog <- c(2,ncol(MASS::UScrime))
mD[,-notlog] <- log(mD[,-notlog])

for (j in 1:ncol(mD)) {
  mD[,j] <- (mD[,j] - mean(mD[,j]))/sd(mD[,j])
}

varnames <- c(
  "log(AGE)",
  "S",
  "log(ED)",
  "log(Ex0)",
  "log(Ex1)",
  "log(LF)",
  "log(M)",
  "log(N)",
  "log(NW)",
  "log(U1)",
  "log(U2)",
  "log(W)",
  "log(X)",
  "log(prison)",
  "log(time)")

vy <- mD$y
mX <- data.matrix(cbind(mD[1:15]))
colnames(mX) <- varnames
K <- 100
p <- ncol(mX)
sampler_result <- sampler(10000, vy, mX, prior = "BIC",
				  modelprior = "uniform",
                                  modelpriorvec_in=NULL)

certifiedwaif/blma documentation built on July 8, 2023, 10:29 p.m.