EMVS: Bayesian Variable Selection using EM Algorithm
In EMVS: The Expectation-Maximization Approach to Bayesian Variable Selection

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/EMVS.R

EMVS is a fast deterministic approach to identifying sparse high posterior models for Bayesian variable selection under spike-and-slab priors in linear regression. EMVS performs dynamic posterior exploration, which outputs a solution path computed at a grid of values for the spike variance parameter v0.

1
2
3

EMVS(Y, X, v0, v1, type = c("betabinomial", "fixed"), independent = TRUE,  
      beta_init, sigma_init, epsilon = 10^(-5), temperature, theta, a, b, v1_g,
    direction=c("backward", "forward", "null"), standardize = TRUE, log_v0 = FALSE)

`Y`	Vector of continuous responses (n x 1). The responses are expected to be centered.
`X`	Matrix of regressors (n x p). Continous predictors are expected to be standardized to have mean zero and standard deviation one.
`v0`	Spike variance parameters. Either a numeric value for a single run or a sequence of increasing values for dynamic posterior exploration.
`v1`	Slab variance parameter. Needs to be greater than v0.
`type`	Type of the prior distribution over the model space: `type="betabinomial"` for the betabinomial prior with shape parameters `a` and `b`, `type="fixed"` for the Bernoulli prior with a fixed inclusion probability `theta`.
`independent`	If TRUE, the regression coefficients and the error variance are taken to be independent a priori (default). If FALSE, a conjugate prior is used as in Rockova and George (2014).
`beta_init`	Vector (p x 1) of initial values for the regression parameters `beta`. If missing, a default vector of starting values obtained as a limiting case of deterministic annealing used beta^0=[X'X+0.5(1/v1+1/v0)I_p]^{-1}X'Y.
`sigma_init`	Initial value for the residual variance parameter.
`epsilon`	Convergence margin parameter. The computation at each `v0` is terminated when \|\|beta^{k+1}-beta^{k}\|\|_2<epsilon.
`temperature`	Temperature parameter for deterministic annealing. If missing, a default value `temperature=1` used.
`theta`	Prior inclusion probability for `type="fixed"`.
`a,b`	Scale parameters of the beta distribution for `type="betabinomial"`.
`v1_g`	Slab variance parameter value for the g-function. If missing, a default value `v1` used.
`direction`	Direction of the sequential reinitialization in dynamic posterior exploration. The default is `direction="backward"` - this initializes the first computation at `beta_init` using the largest value of `v0` and uses the resulting output as a warm start for the next largest value `v0` in a backward direction (i.e. from the largest to the smallest value of `v0`). The option `direction="forward"` proceeds from the smallest value of `v0` to the largest value of `v0`, using the output from the previous solution as a warm start for the next. `direction = "null"` re-initializes at `beta_init` for each `v0`.
`standardize`	If `TRUE` (default), the design matrix X is standardized (mean zero and variance `n`).
`log_v0`	If TRUE, the v0s are displayed on the log scale in `EMVSplot`.

An EM algorithm is applied to find posterior modes of the regression parameters in linear models under spike and slab priors. Variable selection is performed by threshodling the posterior modes to obtain models gamma with high posterior probability P(gamma|Y) . The spike variance v0 can be altered to obtain models with various degrees of sparsity. The slab variance is set to a fixed value v1>v0. The thresholding is based on the conditional posterior probabilities of inclusion, which are outputed of the procedure. Variables are included as long as their inclusion probability is above 0.5. Dynamic exploration is achieved by considering a sequence of increasing spike variance parameters v0. For each v0, a candidate model is obtained. For the conjugate prior case, the best model is then picked according to a criterion ("log g-function"), which equals to the log of the posterior model probability up to a constant

log g(gamma)=log P(gamma| Y) + C.

Independent and sequential initializations are implemented. Sequential initialization uses previously found modes are warm starts in both forward and backward direction of the given sequence of v0 values.

A list object, for which EMVSplot and EMVSbest functions exist.

`betas`	Matrix of estimated regression coefficients (posterior modal estimates) of dimension (L x p), where L is the length of v0.
`log_g_function`	Vector (L x 1) of log posterior model probabilities (up to a constant) of subsets found for each v0. (Only available for `independent = FALSE`).
`intersects`	Vector (L x 1) of posterior weighted intersection points between spike and slab components.
`sigmas`	Vector (L x 1) of estimated residual variances.
`v1`	Slab variance parameter values used.
`v0`	Spike variance parameter values used.
`niters`	Vector (L x 1) of numbers of interations until convergence for each v0
`prob_inclusion`	A matrix (L x p) of conditional inclusion probabilities. Each row corresponds to a single v0 value.
`type`	Type of the model prior used.
`type`	Type of initialization used, `type="null"` stands for the default cold start.
`theta`	Vector (L x 1) of estimated inclusion probabilities for `type="betabinomial"`.

Veronika Rockova

Rockova, V. and George, E. I. (2014) EMVS: The EM Approach to Bayesian Variable Selection. Journal of the American Statistical Association

EMVSplot, EMVSsummary, EMVSbest

# Linear regression with p>n variables
library(EMVS)

n = 100
p = 1000
X = matrix(rnorm(n * p), n, p)
beta = c(1.5, 2, 2.5, rep(0, p-3))
Y = X[,1] * beta[1] + X[,2] * beta[2] + X[,3] * beta[3] + rnorm(n)

# conjugate prior on regression coefficients and variance
v0 = seq(0.1, 2, length.out = 20)
v1 = 1000
beta_init = rep(1, p)
sigma_init = 1
a = b = 1
epsilon = 10^{-5}

result = EMVS(Y, X, v0 = v0, v1 = v1, type = "betabinomial", 
independent = FALSE, beta_init = beta_init, sigma_init = sigma_init,
epsilon = epsilon, a = a, b = b)

EMVSplot(result, "both", FALSE)

EMVSbest(result)

# independent prior on regression coefficients and variance
v0 = exp(seq(-10, -1, length.out = 20))
v1 = 1
beta_init = rep(1,p)
sigma_init = 1
a = b = 1
epsilon = 10^{-5}

result = EMVS(Y, X, v0 = v0, v1 = v1, type = "betabinomial", 
independent = TRUE, beta_init = beta_init, sigma_init = sigma_init,
epsilon = epsilon, a = a, b = b, log_v0 = TRUE)

EMVSplot(result, "both", FALSE)

EMVSbest(result)

Loading Data 
EMVS Begins
v0 = 2
Epsilon 1000.9
Epsilon 7.41599
Epsilon 3.25743
Epsilon 2.86853e-09
v0 = 1.9
Epsilon 11.648
Epsilon 7.43358
Epsilon 2.34771
Epsilon 1.40644e-09
v0 = 1.8
Epsilon 11.6311
Epsilon 7.58138
Epsilon 2.01765
Epsilon 1.16163e-09
v0 = 1.7
Epsilon 11.606
Epsilon 7.77628
Epsilon 1.67192
Epsilon 9.25894e-10
v0 = 1.6
Epsilon 11.5701
Epsilon 8.02749
Epsilon 1.31945
Epsilon 7.0738e-10
v0 = 1.5
Epsilon 11.5193
Epsilon 8.34255
Epsilon 0.973738
Epsilon 5.12886e-10
v0 = 1.4
Epsilon 11.4478
Epsilon 8.72343
Epsilon 0.653857
Epsilon 3.47554e-10
v0 = 1.3
Epsilon 11.3462
Epsilon 9.15808
Epsilon 0.383471
Epsilon 2.15069e-10
v0 = 1.2
Epsilon 11.1998
Epsilon 9.60577
Epsilon 0.184809
Epsilon 1.17704e-10
v0 = 1.1
Epsilon 10.984
Epsilon 9.98079
Epsilon 0.0668746
Epsilon 5.52867e-11
v0 = 1
Epsilon 10.6572
Epsilon 10.1512
Epsilon 0.0159243
Epsilon 2.29775e-11
v0 = 0.9
Epsilon 10.1474
Epsilon 9.96809
Epsilon 0.00211387
Epsilon 1.05124e-11
v0 = 0.8
Epsilon 9.33921
Epsilon 9.30834
Epsilon 0.00015402
Epsilon 6.70881e-12
v0 = 0.7
Epsilon 8.09927
Epsilon 8.10766
Epsilon 6.89928e-06
v0 = 0.6
Epsilon 6.51005
Epsilon 6.52244
Epsilon 2.94502e-07
v0 = 0.5
Epsilon 5.2486
Epsilon 5.25781
Epsilon 4.25061e-07
v0 = 0.4
Epsilon 4.46095
Epsilon 4.4666
Epsilon 1.51772e-06
v0 = 0.3
Epsilon 2.6948
Epsilon 2.69871
Epsilon 8.72089e-07
v0 = 0.2
Epsilon 0.297059
Epsilon 0.298025
Epsilon 4.05292e-08
v0 = 0.1
Epsilon 1.51324e-05
Epsilon 6.80023e-07
Done! 
[1] "Best Model Found"
$log_g_function
[1] -277.4459

$indices
[1] 1 2 3

Loading Data 
EMVS Begins
v0 = 0.367879
Epsilon 1000.9
Epsilon 3.34086e-05
Epsilon 6.00443e-07
v0 = 0.22908
Epsilon 6.67257e-05
Epsilon 3.76853e-05
Epsilon 1.02919e-05
Epsilon 3.48915e-06
v0 = 0.142649
Epsilon 0.000397206
Epsilon 0.000146879
Epsilon 4.10709e-05
Epsilon 7.11398e-06
v0 = 0.0888281
Epsilon 0.00123288
Epsilon 0.000502512
Epsilon 0.000120913
Epsilon 1.05823e-05
Epsilon 8.39699e-07
v0 = 0.0553136
Epsilon 0.00363265
Epsilon 0.00143932
Epsilon 0.000383756
Epsilon 1.61988e-05
Epsilon 6.53279e-07
v0 = 0.034444
Epsilon 0.0111381
Epsilon 0.00288489
Epsilon 0.00251479
Epsilon 5.6918e-05
Epsilon 6.35677e-07
v0 = 0.0214484
Epsilon 0.0440164
Epsilon 0.216169
Epsilon 2.67671
Epsilon 5.30661e-06
v0 = 0.013356
Epsilon 0.0211353
Epsilon 0.00875779
Epsilon 0.000742659
Epsilon 4.38609e-06
v0 = 0.00831683
Epsilon 0.0410251
Epsilon 0.0908389
Epsilon 1.81721
Epsilon 0.000169428
Epsilon 1.11577e-07
v0 = 0.00517892
Epsilon 0.0393829
Epsilon 1.42075
Epsilon 0.0401635
Epsilon 3.55991e-06
v0 = 0.00322494
Epsilon 0.00845322
Epsilon 0.00460163
Epsilon 0.00032534
Epsilon 1.34582e-05
Epsilon 2.90744e-07
v0 = 0.00200818
Epsilon 0.0141653
Epsilon 0.00472028
Epsilon 0.000732522
Epsilon 0.000120404
Epsilon 1.57296e-05
Epsilon 1.66377e-06
v0 = 0.0012505
Epsilon 0.021615
Epsilon 0.00357136
Epsilon 0.000756838
Epsilon 0.000219508
Epsilon 7.07498e-05
Epsilon 2.35991e-05
Epsilon 7.92357e-06
v0 = 0.000778692
Epsilon 0.0239493
Epsilon 0.00190486
Epsilon 0.000347045
Epsilon 8.69339e-05
Epsilon 2.47781e-05
Epsilon 7.54036e-06
v0 = 0.000484895
Epsilon 0.0128673
Epsilon 0.000742589
Epsilon 8.14548e-05
Epsilon 1.14052e-05
Epsilon 1.71157e-06
v0 = 0.000301946
Epsilon 0.00481121
Epsilon 0.00023011
Epsilon 1.21327e-05
Epsilon 7.79699e-07
v0 = 0.000188023
Epsilon 0.00168616
Epsilon 6.34437e-05
Epsilon 1.40928e-06
v0 = 0.000117083
Epsilon 0.000591897
Epsilon 1.74732e-05
Epsilon 1.51869e-07
v0 = 7.29077e-05
Epsilon 0.00021581
Epsilon 5.11008e-06
v0 = 4.53999e-05
Epsilon 7.8134e-05
Epsilon 1.53696e-06
Done! 
[1] "Log posterior not available for independent prior case"