mc3_pep {PEPBVS}

Bayesian variable selection with MC3 algorithm

Description

Given a response vector and an input data matrix, performs Bayesian variable selection using the MC3 algorithm. Normal linear models are assumed for the data with the prior distribution on the model parameters (beta coefficients and error variance) being the PEP or the intrinsic. The prior distribution on the model space can be the uniform on the model space or the uniform on the model dimension (special case of the beta-binomial prior).

Usage

mc3_pep(
  x,
  y,
  intrinsic = FALSE,
  reference.prior = TRUE,
  beta.binom = TRUE,
  ml_constant.term = FALSE,
  burnin = 1000,
  itermc3 = 11000
)

Arguments

x

A matrix of numeric (of size nxp), input data matrix. This matrix contains the values of the p explanatory variables without an intercept column of 1's.

y

A vector of numeric (of length n), response vector.

intrinsic

Boolean, indicating whether the PEP (FALSE) or the intrinsic - which is a special case of it - (TRUE) should be used as prior on the regression parameters. Default value=FALSE.

reference.prior

Boolean, indicating whether the reference prior (TRUE) or the dependence Jeffreys prior (FALSE) is used as baseline. Default value=TRUE.

beta.binom

Boolean, indicating whether the beta-binomial distribution (TRUE) or the uniform distribution (FALSE) should be used as prior on the model space. Default value=TRUE.

ml_constant.term

Boolean, indicating whether the constant (marginal likelihood of the null/intercept-only model) should be included in computing the marginal likelihood of a model (TRUE) or not (FALSE). Default value=FALSE.

burnin

Non-negative integer, the burnin period for the MC3 algorithm. Default value=1000.

itermc3

Positive integer (larger than burnin), the (total) number of iterations for the MC3 algorithm. Default value=11000.

Details

The function works when p<=n-2 where p is the number of explanatory variables and n is the sample size.

It is suggested to use this function (i.e. MC3 algorithm) when p is larger than 20.

The reference model is the null model (i.e. intercept-only model).

The case of missing data (i.e. presence of NA's either in the input matrix or the response vector) is not currently supported.

The intercept term is included in all models.

If p>1, the input matrix needs to be of full rank.

The reference prior as baseline corresponds to hyperparameter values d0=0 and d1=0, while the dependence Jeffreys prior corresponds to model-dependent-based values for the hyperparameters d0 and d1, see Fouskakis and Ntzoufras (2022) for more details.

The MC3 algorithm was first introduced by Madigan and York (1995) while its current implementation is described in the Appendix of Fouskakis and Ntzoufras (2022).

When ml_constant.term=FALSE then the log marginal likelihood of a model in the output is shifted by -logC1 (logC1: marginal likelihood of the null model).

When the prior on the model space is beta-binomial (i.e. beta.binom=TRUE), the following special case is used: uniform prior on model size.

Value

mc3_pep returns an object of class pep, as this is described in detail in full_enumeration_pep. The difference is that here the number of rows of the first list element is not 2^p but the number of unique models ‘visited’ by the MC3 algorithm. Further, the posterior probability of a model corresponds to the estimated posterior probability as this is computed by the relative Monte Carlo frequency of the ‘visited’ models by the MC3 algorithm.

References

Fouskakis, D. and Ntzoufras, I. (2022) Power-Expected-Posterior Priors as Mixtures of g-Priors in Normal Linear Models. Bayesian Analysis, 17(4): 1073-1099. doi:10.1214/21-BA1288

Madigan, D. and York, J. (1995) Bayesian Graphical Models for Discrete Data. International Statistical Review, 63(2): 215–232. doi:10.2307/1403615

See Also

full_enumeration_pep

Examples

data(UScrime_data)
y <- UScrime_data[,"y"]
X <- UScrime_data[,-15]
set.seed(123)
res <- mc3_pep(X,y,itermc3=3000)
resu <- mc3_pep(X,y,beta.binom=FALSE,itermc3=3000)
resj <- mc3_pep(X,y,reference.prior=FALSE,burnin=500,itermc3=2200)