Add text here.

Analysis settings

These variables specify the number of samples to simulate (“n”), the correlation between the two simulated variables, X1 and X2 (“xcor”), and the regression coefficients used to simulate the outcomes (“beta”).

n    <- 40
xcor <- 0.99
beta <- c(3,0)

These two variables specify how the variational parameters are initialized in the three varbvs runs.

alpha0 <- rbind(c(0,1),
                c(1,0),
                c(1,1))
mu0    <- rbind(c(0,3),
                c(2.5,0),
                c(1,1.5))

For creating the surface plot of the variational objective, the objective is evaluated at all grid points specified by “mu_grid”.

mu_grid <- seq(-1,3.5,0.01)

Set up environment

Load a few packages as well as some functions used in the analysis below.

library(MASS)
library(varbvs)
library(ggplot2)
library(cowplot)
source("visualize_varbvs_surface_functions.R")

Initialize the sequence of pseudorandom numbers.

set.seed(1)

Generate data set

Simulate the \(n \times 2\) matrix of predictors, in which the two predictors are strongly correlated.

S <- rbind(c(1,xcor),c(xcor,1))
X <- mvrnorm(n,c(0,0),S)
X <- scale(X,center = TRUE,scale = FALSE)

Simulate the continuously-valued outcomes.

y <- drop(X %*% beta + rnorm(n))
y <- y - mean(y)

Fit model

Run the co-ordinate ascent updates with three different initial estimates of the variational parameters.

fit1 <- varbvs(X,NULL,y,sigma = 1,sa = 1,logodds = 0,tol = 0,
               alpha = alpha0[1,],mu = mu0[1,],verbose = FALSE)
fit2 <- varbvs(X,NULL,y,sigma = 1,sa = 1,logodds = 0,tol = 0,
               alpha = alpha0[2,],mu = mu0[2,],verbose = FALSE)
fit3 <- varbvs(X,NULL,y,sigma = 1,sa = 1,logodds = 0,tol = 0,
               alpha = alpha0[3,],mu = mu0[3,],verbose = FALSE)
pts <- as.data.frame(rbind(fit1$beta,fit2$beta,fit3$beta))
print(pts)
#        X1      X2
# 1 0.02278 2.54849
# 2 2.61971 0.02711
# 3 1.25072 1.40955

The co-ordinate ascent updates should converge to three different stationary points of the K-L divergence (equivalently, the variational lower bound).

Plot variational objective surface

Compute the value of the variational objective (the K-L divergence) at each grid point \((\mu_1, \mu_2)\), where \(\mu_j\) gives the (approximate) posterior mean of the jth regression coefficient given that the coefficient is nonzero.

pdat <- expand.grid(mu1 = mu_grid,mu2 = mu_grid)
pdat <- compute_kl_grid(X,y,pdat)

Create a contour plot to visualize the 2-d surface of the Kullback-Divergence divergence with respect to the variational estimates of the posterior mean coefficients.

p1 <- create_contour_plot(pdat,pts)
print(p1)

This contour plot confirms that there are indeed three stationary points (local minima) of the variational objective: at two of these stationary points, one of the regression coefficients is zero, whereas both coefficients are nonzero, and nearly equal, at the third stationary point.

A second look

To better appreciate this result, let’s plot the objective along a “slice” that cuts through the three modes:

pdat2 <- data.frame(mu1 = mu_grid)
pdat2 <- transform(pdat2,mu2 = -0.9275*mu1 + 2.5696)
pdat2 <- compute_kl_grid(X,y,pdat2)
pdat2 <- transform(pdat2,beta1 = alpha1*mu1)
p2 <- ggplot(as.data.frame(pdat2),aes(x = beta1,y = KL)) +
      geom_line(color = "dodgerblue",size = 1) +
      labs(x = "posterior mean of beta1",y = "K-L divergence")
print(p2)

From this plot we clearly see the three stationary points of the variational objective. The global mean is reached at \(\mu_1 = 1.25, \mu_2 = 1.41\), but the local mode at \(\mu_1 = 0.02, \mu_2 = 2.55\) achieves nearly the same value.

Session info

This is the version of R and the packages that were used to generate these results.

sessionInfo()
# R version 3.4.3 (2017-11-30)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.6
# 
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] cowplot_0.9.4 ggplot2_3.1.0 varbvs_2.6-5  MASS_7.3-48  
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.0          compiler_3.4.3      pillar_1.3.1       
#  [4] RColorBrewer_1.1-2  plyr_1.8.4          tools_3.4.3        
#  [7] digest_0.6.17       evaluate_0.11       tibble_2.1.1       
# [10] gtable_0.2.0        lattice_0.20-35     pkgconfig_2.0.2    
# [13] rlang_0.3.1         Matrix_1.2-12       yaml_2.2.0         
# [16] withr_2.1.2         stringr_1.3.1       dplyr_0.8.0.1      
# [19] knitr_1.20          rprojroot_1.3-2     grid_3.4.3         
# [22] tidyselect_0.2.5    glue_1.3.0          R6_2.2.2           
# [25] rmarkdown_1.10      latticeExtra_0.6-28 purrr_0.2.5        
# [28] magrittr_1.5        backports_1.1.2     scales_0.5.0       
# [31] htmltools_0.3.6     assertthat_0.2.0    colorspace_1.4-0   
# [34] labeling_0.3        nor1mix_1.2-3       stringi_1.2.4      
# [37] lazyeval_0.2.1      munsell_0.4.3       crayon_1.3.4