
\documentclass{article}
\usepackage[utf8]{inputenc}

\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amscd}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage[margin = 3cm]{geometry}

\usepackage{enumerate}


% environments
\usepackage{framed}
\usepackage[ruled,vlined,norelsize]{algorithm2e}

% graphics
\usepackage{float}
\usepackage{graphicx}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}

% captions
\usepackage{caption}
\usepackage{subcaption}
\captionsetup[figure]{font=small,labelfont={small, bf}}
\captionsetup[subfigure]{font=footnotesize,labelfont={footnotesize, bf}}

\usepackage{xcolor}

\usepackage{hyperref}

\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt}

% vectors
\newcommand{\va}{\mathbf{a}}
\newcommand{\vb}{\mathbf{b}}
\newcommand{\vc}{\mathbf{c}}
\newcommand{\vd}{\mathbf{d}}
\newcommand{\ve}{\mathbf{e}}
\newcommand{\vf}{\mathbf{f}}
\newcommand{\vg}{\mathbf{g}}
\newcommand{\vh}{\mathbf{h}}
\newcommand{\vi}{\mathbf{i}}
\newcommand{\vj}{\mathbf{j}}
\newcommand{\vk}{\mathbf{k}}
\newcommand{\vl}{\mathbf{l}}
\newcommand{\vm}{\mathbf{m}}
\newcommand{\vn}{\mathbf{n}}
\newcommand{\vo}{\mathbf{o}}
\newcommand{\vp}{\mathbf{p}}
\newcommand{\vq}{\mathbf{q}}
\newcommand{\vr}{\mathbf{r}}
\newcommand{\vs}{\mathbf{s}}
\newcommand{\vt}{\mathbf{t}}
\newcommand{\vu}{\mathbf{u}}
\newcommand{\vv}{\mathbf{v}}
\newcommand{\vw}{\mathbf{w}}
\newcommand{\vx}{\mathbf{x}}
\newcommand{\vy}{\mathbf{y}}
\newcommand{\vz}{\mathbf{z}}

\newcommand{\valpha}{\boldsymbol{\alpha}}
\newcommand{\vbeta}{\boldsymbol{\beta}}
\newcommand{\vgamma}{\boldsymbol{\gamma}}
\newcommand{\vmu}{\boldsymbol{\mu}}
\newcommand{\vtheta}{\boldsymbol{\theta}}
\newcommand{\veta}{\boldsymbol{\eta}}
\newcommand{\vepsilon}{\boldsymbol{\varepsilon}}
\newcommand{\vsigma}{\boldsymbol{\sigma}}

\newcommand{\vzero}{\mathbf{0}}
\newcommand{\vone}{\mathbf{1}}

% matrices

\newcommand{\mA}{\mathbf{A}}
\newcommand{\mB}{\mathbf{B}}
\newcommand{\mC}{\mathbf{C}}
\newcommand{\mD}{\mathbf{D}}
\newcommand{\mE}{\mathbf{E}}
\newcommand{\mF}{\mathbf{F}}
\newcommand{\mG}{\mathbf{G}}
\newcommand{\mH}{\mathbf{H}}
\newcommand{\mI}{\mathbf{I}}
\newcommand{\mJ}{\mathbf{J}}
\newcommand{\mK}{\mathbf{K}}
\newcommand{\mL}{\mathbf{L}}
\newcommand{\mM}{\mathbf{M}}
\newcommand{\mN}{\mathbf{N}}
\newcommand{\mO}{\mathbf{O}}
\newcommand{\mP}{\mathbf{P}}
\newcommand{\mQ}{\mathbf{Q}}
\newcommand{\mR}{\mathbf{R}}
\newcommand{\mS}{\mathbf{S}}
\newcommand{\mT}{\mathbf{T}}
\newcommand{\mU}{\mathbf{U}}
\newcommand{\mV}{\mathbf{V}}
\newcommand{\mW}{\mathbf{W}}
\newcommand{\mX}{\mathbf{X}}
\newcommand{\mY}{\mathbf{Y}}
\newcommand{\mZ}{\mathbf{Z}}

\newcommand{\msigma}{\boldsymbol{\Sigma}}
\newcommand{\mgamma}{\boldsymbol{\Gamma}}

\newcommand{\defeq}{\mathrel{\mathop:}=}
\newcommand{\eqdef}{\mathrel{\mathop=}:}

\DeclareMathOperator{\logit}{logit}
\DeclareMathOperator{\expit}{expit}

\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}

% expected value
\newcommand{\E}{\mathbb{E}}

% natural numbers, real numbers
\newcommand{\N}{\mathbb{N}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}



\begin{document}

<<setup, include=FALSE, cache=FALSE>>=
library(knitr)
# set global chunk options
opts_chunk$set(fig.align='center')
options(formatR.arrow=TRUE,width=90)
@

%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{EMVS Vignette}
%\VignetteEncoding{UTF-8}

\title{EMVS Vignette}
\author{Veronika Rockova and Gemma Moran}
\maketitle
\section{Introduction}

This vignette describes the algorithm ``EMVS", the EM approach to Bayesian Variable Selection (Rockova and George, 2014) and its usage in the R package \texttt{EMVS}. This R package implementation of EMVS has two options for prior specification:

1. A ``conjugate" or ``scale-invariant" prior on the regression coefficients, as detailed in Rockova and George (2014);

2. An ``independent" prior where the regression coefficients and error variance are treated as independent \emph{a priori}, which is recommended by Moran, Rockova and George (2018).

This vignette details the EMVS algorithm where the second, independent prior formulation is used, and provides an example of its usage in the package.

\section{Model}

Consider the classical linear regression model
\begin{equation}
\mathbf{Y}= \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim N_n(0, \sigma^2\mathbf{I}_n) \label{regression}
\end{equation}
where $\mathbf{Y}\in\mathbb{R}^n$ is a vector of responses, $\mathbf{X} = [\mathbf{X}_1,\dots, \mathbf{X}_p] \in \mathbb{R}^{n\times p}$ is a fixed regression matrix of $p$ potential predictors, $\boldsymbol{\beta} = (\beta_1,\dots, \beta_p)^T \in \mathbb{R}^p$ is a vector of unknown regression coefficients and $\boldsymbol{\epsilon} \in \mathbb{R}^n$ is the noise vector of independent normal random variables with $\sigma^2$ as their unknown common variance. 

The goal of EMVS is to find which predictors $\mathbf{x}_i$ should be included in the model. In the Bayesian paradigm, this is facilitated by the introduction of binary latent variables $\boldsymbol{\gamma} = (\gamma_1,\dots, \gamma_p)^T$, $\gamma_i \in \{0,1\}$, where $\gamma_i = 1$ indicates that $\mathbf{x}_i$ is to be included in the model. 

The hierarchical model is given by:
\begin{align}
\pi(\boldsymbol{\beta}|\boldsymbol{\gamma}, v_0, v_1) = N_p(\mathbf{0}, \mathbf{D}_{\mathbf{\gamma}}) \label{beta_prior}
\end{align}
where $\mathbf{D} = \text{diag}(d_1, \dots, d_p)$ with $d_i = (1 - \gamma_i)v_0 + \gamma_iv_1$ for $0 \leq v_0 < v_1$. Following George and McCulloch (1997), Rockova and George (2014) recommend setting the hyperparameters $v_0$ and $v_1$ to be small and large fixed values respectively; this yields the canonical ``spike" and ``slab" of Bayesian variable selection.  Note that the above prior does not depend on the error variance $\sigma^2$, unlike the original formulation in Rockova and George (2014).

The prior on the error variance is an inverse gamma:
\begin{align}
\pi(\sigma^2) = \text{IG}(\nu/2, \nu\lambda/2).\label{sigma_prior}
\end{align}

To incorporate uncertainty regarding which variables should be included in the model, EMVS additionally specifies a prior for the latent indicator variables $\boldsymbol{\gamma}$. This prior is iid Bernoulli:
\begin{align}
\pi(\boldsymbol{\gamma}|\theta) = \theta^{|\boldsymbol{\gamma}|}(1-\theta)^{p - |\boldsymbol{\gamma}|}
\end{align}
where $\theta \in [0, 1]$ and $|\boldsymbol{\gamma}| = \sum_{i=1}^p \gamma_i$.  The proportion of non-zero regression coefficients, $\theta$, is unknown: this parameter is assigned a beta prior:
\begin{align}
\pi(\theta) \propto \theta^{a-1}(1-\theta)^{b-1}, \quad a, b > 0.
\end{align}

\section{Algorithm}

The EMVS algorithm treats the latent indicators $\boldsymbol{\gamma}$ as ``missing data" and indirectly maximizes the posterior $\pi(\boldsymbol{\beta}, \theta, \sigma^2|\mathbf{Y})$ by iteratively maximizing the following objective function:
\begin{align}
\mathcal{Q}(\boldsymbol{\beta}, \theta, \sigma|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}) = E_{\boldsymbol{\gamma}|\cdot}[\log\pi(\boldsymbol{\beta}, \theta, \sigma^2, \boldsymbol{\gamma}|\mathbf{Y})|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}, \mathbf{Y}] \label{Q}
\end{align}
where $E_{\boldsymbol{\gamma}|\cdot}(\cdot)$ denotes the conditional expectation $E_{\boldsymbol{\gamma}|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}, \mathbf{Y}}(\cdot)$.  At the $k$th iteration, given $(\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)})$, an E-step is first applied, which computes the expectation of the right side of \eqref{Q} to obtain $\mathcal{Q}$. This is followed by an M-step, which maximizes $\mathcal{Q}$ over $(\boldsymbol{\beta}, \theta, \sigma^2)$ to yield the values of $(\boldsymbol{\gamma}^{(k+1)}, \theta^{(k+1)}, \sigma^{2(k+1)})$.

The objective function is of the form:
\begin{align}
\mathcal{Q}(\boldsymbol{\beta}, \theta, \sigma^2|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}) = C + \mathcal{Q}_1(\boldsymbol{\beta}, \sigma^2|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}) + \mathcal{Q}_2(\theta|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)})
\end{align}
where
\begin{align}
\mathcal{Q}_1(\boldsymbol{\beta}, \sigma^2|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}) &= - \frac{1}{2\sigma^2}\lVert \mathbf{Y} - \mathbf{X}\boldsymbol{\beta}\rVert^2 -\frac{n+\nu + 2}{2}\log(\sigma^2) - \frac{\nu\lambda}{2\sigma^2}\\
&\quad - \frac{1}{2}\sum_{i=1}^p\beta_i^2 E_{\boldsymbol{\gamma}|\cdot}\left[ \frac{1}{v_0(1-\gamma_i) + v_1\gamma_i}\right],
\end{align}
and
\begin{align}
\mathcal{Q}_2(\theta|\boldsymbol{\beta}^{(k)}, \theta^{(k)}, \sigma^{2(k)}) = \sum_{i=1}^p \log\left(\frac{\theta}{1-\theta}\right) E_{\boldsymbol{\gamma}|\cdot}\gamma_i + (a-1)\log(\theta) + (p+b-1)\log(1-\theta).
\end{align}

\subsection{E-Step}
We have:
\begin{align}
E_{\boldsymbol{\gamma}|\cdot}\gamma_i = P(\gamma_i = 1|\boldsymbol{\beta}^{(k)}, \theta^{(k)}) = p_i^*
\end{align}
where
\begin{align}
p_i^* = \frac{\theta^{(k)}\phi_{v_1}(\beta_i^{(k)})}{\theta^{(k)}\phi_{v_1}(\beta_i^{(k)}) + (1-\theta^{(k)})\phi_{v_0}(\beta_i^{(k)})}
\end{align}
where $\phi_v(x) = \frac{1}{\sqrt{2\pi v}}\exp(-x^2/2v)$ is the normal density with zero mean and variance, $v$. 

To complete the E-step, we have
\begin{align}
E_{\boldsymbol{\gamma}|\cdot}\left[ \frac{1}{v_0(1-\gamma_i) + v_1\gamma_i}\right] &= \frac{E_{\boldsymbol{\gamma}|\cdot}(1-\gamma_i)}{v_0} + \frac{E_{\boldsymbol{\gamma}|\cdot}\gamma_i}{v_1} \\
&= \frac{1-p_i^*}{v_0} + \frac{p_i^*}{v_1}\\
& \equiv d_i^*.
\end{align}
We denote the matrix $\mathbf{D}^* = \text{diag}(d_1^*,\dots, d_p^*)$.

\subsection{M-Step}
The objective function \eqref{Q} yields closed form updates for each of ($\boldsymbol{\beta}^{(k+1)}, \theta^{(k+1)}, \sigma^{2(k+1)})$. These are:
\begin{align}
\boldsymbol{\beta}^{(k+1)} &= [\mathbf{X}^T\mathbf{X} + \sigma^{2(k)}\mathbf{D}^*]^{-1}\mathbf{X}^T\mathbf{Y}  \label{beta_m_step}, \\
\sigma^{2(k+1)} &= \frac{\lVert \mathbf{Y} - \mathbf{X}\boldsymbol{\beta}^{(k+1)}\rVert^2 + \nu\lambda}{n + \nu + 2}, \\
\theta^{(k+1)} &= \frac{\sum_{i=1}^pp_i^* + a -1 }{a + b- p - 2}.
\end{align}

In problems where $p \gg n$, the calculation cost of \eqref{beta_m_step} is substantially reduced by using the Sherman-Morrison-Woodbury formula to obtain:
\begin{align}
\boldsymbol{\beta}^{(k+1)} = \sigma^2[\mathbf{D}^{*-1}-\mathbf{D}^{*-1}\mathbf{X}^T(\sigma^2\mathbf{I}_{n\times n} + \mathbf{X}\mathbf{D}^{*-1}\mathbf{X}^T)^{-1}\mathbf{X}\mathbf{D}^{*-1}]\mathbf{X}^T\mathbf{Y}.
\end{align}

The EMVS algorithm iterates between the above steps until convergence. The default convergence criterion is: $\lVert \boldsymbol{\beta}^{(k+1)} - \boldsymbol{\beta}^{(k)}\rVert^2 \leq 10^{-5}$.

\subsection{Thresholding the EM Output for Variable Selection}

Once we have obtained our MAP estimates $(\widehat{\boldsymbol{\beta}}, \widehat{\theta},\widehat{\sigma}^2)$, we can find the most probable $\boldsymbol{\gamma}$ given those values. This is obtained by setting
\begin{align}
\widehat{\gamma}_i = 
\begin{cases}
1 &\text{if } P(\gamma_i = 1|\widehat{\boldsymbol{\beta}}, \widehat{\theta}, \widehat{\sigma}^2) \geq 0.5 \\
0 &\text{otherwise.}
\end{cases}
\end{align}
This selection of $\widehat{\boldsymbol{\gamma}}$ is equivalent to thresholding the $\widehat{\beta}_i$ (Rockova and George, 2014).  This thresholding occurs at the intersection points $\pm \beta_i^*$ of the $P(\gamma_i = 1|\widehat{\boldsymbol{\beta}}, \widehat{\theta})$ weighted mixture of the spike-and-slab priors, namely,
\begin{align}
\pm \beta_i^*(v_0, v_1, \widehat{\theta}) = \pm\sqrt{2v_0 \log(\omega_ic)c^2/(c^2 - 1)}
\end{align}
where $c^2 = v_1/v_0$ and $\omega_i = (1-\widehat{\theta})/\widehat{\theta}$. Then, the thresholding rule is
\begin{align}
\widehat{\gamma}_i = 
\begin{cases}
1 &\text{if } |\widehat{\beta}_i| \geq \beta_i^*(v_0, v_1, \widehat{\theta}). \\
0 &\text{otherwise.}
\end{cases}
\end{align}


\section{Dynamic Posterior Exploration}

The speed of the EM algorithm allows for posterior modes to be found for multiple values of the spike parameter $v_0$. This facilitates a \emph{dynamic posterior exploration} strategy where the slab parameter $v_1$ is held fixed and $v_0$ is gradually decreased to approximate the ideal point mass spike prior.  This is akin to an annealing strategy; when $v_0 = v_1$, the posterior is convex (in this case, it is equivalent to ridge regression) and as $v_0$ is decreased, it becomes multimodal. By starting with large $v_0$ and using the resultant solution as a ``warm start" for the next (smaller) value of $v_0$, the procedure can more easily find modes of the posterior.  Decreasing $v_0$ also serves to shrink smaller coefficients to zero and find a sparse solution to the regression problem. The original paper (Rockova and George, 2014) illustrates a ``forward'' strategy for EMVS by starting with small $v_0$ and gradually increasing it. However, our recommendation is to use a ``backward" strategy wherein we start with large $v_0$ and gradually \emph{decrease} it.  In this ``backward" strategy, the algorithm also stabilizes at a solution for small $v_0$, in many cases eliminating the need to choose a particular ``best" spike parameter.

The \texttt{EMVS} function in the R package has three options for dynamic posterior exploration: \texttt{direction = c("backward", "foward", "null")}. These are described below.

\begin{itemize}
\item \texttt{direction = "backward"} \emph{(default, recommended)}: this option starts with the largest value of \texttt{v0} and finds a solution to the EMVS algorithm. This solution is used as the inital value for the algorithm with the second largest value of \texttt{v0}. This process is repeated until the smallest value of \texttt{v0}.
\item \texttt{direction = "forward"}: this option is similar to the \texttt{backward} option, but starts with the smallest value of \texttt{v0} and then uses the resultant solution as the initial value for the next largest value of \texttt{v0}, repeating the process until the largest value of \texttt{v0}.
\item \texttt{direction = "null"}: this option uses $\mathtt{beta\_init}$ as the initial value for each value of \texttt{v0} and as such is not a dynamic posterior exploration strategy.
\end{itemize}

\section{Prior Specification}

As mentioned in the introduction, the \texttt{EMVS} function in the R package has two options for the prior specification: \texttt{independent = c(TRUE, FALSE)}. These are described below.

\begin{itemize}
\item \texttt{independent = TRUE} \emph{(default, recommended)}: this option implements EMVS with the ``independent'' prior  detailed in Section 2 of this vignette. 
\item \texttt{independent = FALSE}: this option implements EMVS with the ``conjugate'' or ``scale-invariant'' prior on the regression coefficients detailed in the original paper (Rockova and George, 2014).
\end{itemize}

The reason we recommend the ``independent'' prior for EMVS is that it yields better error variance estimates (Moran et al., 2018). Intuitively, the ``conjugate'' prior adds $p$ ``pseudo-observations'' of the error variance $\sigma^2$, which can result in severe underestimation of the error variance. For more details, see Moran et al. (2018).

\section{Example}

In this section, we demonstrate the basic usage of \texttt{EMVS} for both the independent and conjugate prior implementations. We conclude with a comparison of the error variance estimates from the two, highlighing the benefit of the independent prior formulation.

We begin by loading the package:
<<>>=
library(EMVS)
@

We create a toy dataset with $n = 100$ and $p = 1000$:
<<>>=
set.seed(12022018)
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)
@


\subsection{Independent Prior}

We set the parameters for EMVS: the ``ladder" of values of $v_0$, the slab parameter $v_1$, the initial $\boldsymbol{\beta}$ and the hyperparameters of the beta distribution, $a$, $b$. 

<< results = "hide">>=
# 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
@

We then run \texttt{EMVS} using the independent prior formulation described in Section 2 of this vignette (\texttt{independent = TRUE}). The output is stored in $\mathtt{result\_ind}$.

<< results = "hide" >>=
# independent prior on regression coefficients and variance
result_ind = EMVS(Y = Y, X = X, v0 = v0, v1 = v1, type = "betabinomial", 
                  independent = TRUE, beta_init = beta_init, 
                  sigma_init = sigma_init,
                  a = a, b = b, log_v0 = TRUE)
@

The function \texttt{EMVSplot} plots the estimates of the regression coefficients, $\widehat{\vbeta}$, over all the values of \texttt{v0} (i.e. the regularization plot). Here the $v_0$ are plotted on the log scale - this is because we set the option $\mathtt{log\_v0 = TRUE}$ when running \texttt{EMVS}.
<< fig.width = 2.75, fig.height = 3 , dev.args = list(pointsize =8) >>=
EMVSplot(result_ind, "both", FALSE)
@

Note that the default option \texttt{direction = "backward"} results in the EMVS algorithm stabilizing for small $v_0$. This eliminates the need to choose a $v_0$; we take the coefficients at the smallest $v_0$ as our solution.

\subsection{Conjugate Prior}

We now run \texttt{EMVS}, using the conjugate prior formulation as outlined in the original paper (\texttt{independent = FALSE}). The output is stored in $\mathtt{result\_conj}$. We re-initialize both the slab parameter \texttt{v1} and the ladder of spike values \texttt{v0} as the scale of the variance is different for the independent and conjugate formulations.
<< results = "hide">>=
v0 = seq(0.1, 2, length.out = 20)
v1 = 1000
result_conj = EMVS(Y, X, v0 = v0, v1 = v1, type = "betabinomial", 
                   independent = FALSE, beta_init = beta_init, 
                   sigma_init = sigma_init, a = a, b = b)
@

The function \texttt{EMVSbest} shows the maximum value of the $\mathtt{log\_g}$ function over all the \texttt{v0} values (Rockova and George, 2014) and the non-zero indices of the model with the highest value of $\mathtt{log\_g}$ (the marginal posterior).
<<>>=
EMVSbest(result_conj)
@

We plot both the regularization path and the values of the $\mathtt{log\_g}$ function.
<< fig.width = 5.5, fig.height= 3, dev.args = list(pointsize =8) >>=
EMVSplot(result_conj, "both", FALSE)
@

Note that the $\mathtt{log\_g}$ function is unavailable for the \texttt{independent = TRUE} implementation as the independent prior formulation does not yield a closed form for the marginal posterior. An approximation will be implemented in a future update to the \texttt{EMVS} package.  As discussed earlier, however, the ``backward'' strategy stabilizes for small values $v_0$ and so we recommend taking the coefficients found at the smallest value of $v_0$ as the solution.  In many cases this eliminates the need for the $\mathtt{log\_g}$ function as a criterion. 

\subsection{Variance Estimation}

We now compare the error variance estimates from the independent and conjugate prior implementations of EMVS.

The below plot shows the estimates of the standard deviation of the error ($\widehat{\sigma}$) for each value of $v_0$ for both the independent and conjugate prior formulations, as well as the regularization plots for the coefficients, $\vbeta$.

<< echo = F, warning = F, fig.width = 5.5, fig.height = 6, dev.args = list(pointsize =8)  >>=
par(mfrow = c(2, 2))

# independent regularization plot
select <- apply(result_ind$prob_inclusion, 2, function(x){as.numeric(x > 0.5)})
matplot(log(result_ind$v0), result_ind$betas, xlab = expression(log(v[0])), ylab = expression(hat(beta)), lwd = 1, col = "grey", lty = 2, type = "l")
matpoints(log(result_ind$v0), result_ind$betas * select, xlab = expression(log(v[0])), ylab = expression(hat(beta)), lwd = 1, col = 4, lty = 2, pch = 19)
matpoints(log(result_ind$v0), result_ind$betas*(1-select), xlab = expression(log(v[0])), ylab = expression(hat(beta)), lwd = 1, col = 2, lty = 2, pch = 19)
par(xpd = T)
labels = paste("X", 1:ncol(result_ind$betas), sep = "")
labels[select[length(result_ind$v0), ] == 0] <-""
text(max(log(result_ind$v0)) + 0.75, result_ind$betas[length(result_ind$v0), ], labels = labels, col = 4)
title("Regularization Plot (Independent)")

# independent sigma estimates
plot(log(result_ind$v0), result_ind$sigmas, xlab = expression(log(v[0])), ylab = expression(hat(sigma)), lwd = 1, col = "grey", lty = 2, type = "l", ylim = c(0, 1.1))
points(log(result_ind$v0), result_ind$sigmas, xlab = expression(log(v[0])), ylab = expression(hat(sigma)), lwd = 1, lty = 2, pch = 19)
par(xpd = F)
abline(h = 1, col = "red")
title(expression(paste(bold(hat(sigma)), bold(" (Independent)"))))
ind_sig_max = round(max(result_ind$sigmas), 4)


# conjugate regularization plot
select <- apply(result_conj$prob_inclusion, 2, function(x){as.numeric(x > 0.5)})
matplot(result_conj$v0, result_conj$betas, xlab = expression(v[0]), ylab = expression(hat(beta)), lwd = 1, col = "grey", lty = 2, type = "l")
matpoints(result_conj$v0, result_conj$betas * select, xlab = expression(v[0]), ylab = expression(hat(beta)), lwd = 1, col = 4, lty = 2, pch = 19)
matpoints(result_conj$v0, result_conj$betas*(1-select), xlab = expression(v[0]), ylab = expression(hat(beta)), lwd = 1, col = 2, lty = 2, pch = 19)
par(xpd = T)
labels = paste("X", 1:ncol(result_conj$betas), sep = "")
labels[select[length(result_conj$v0), ] == 0] <-""
text(max(result_conj$v0) * (1.1), result_conj$betas[length(result_conj$v0), ], labels = labels, col = 4)
title("Regularization Plot (Conjugate)")

# conjugate sigma estimates
plot(result_conj$v0, result_conj$sigmas, xlab = expression(v[0]), ylab = expression(hat(sigma)), lwd = 1, col = "grey", lty = 2, type = "l")
points(result_conj$v0, result_conj$sigmas, xlab = expression(v[0]), ylab = expression(hat(sigma)), lwd = 1, lty = 2, pch = 19)
par(xpd = F)
abline(h = 1, col = "red")
title(expression(paste(bold(hat(sigma)), bold(" (Conjugate)"))))
conj_sig_max = round(max(result_conj$sigmas), 4)

@

We can see that the conjugate prior formulation severely underestimates the true error variance ($\sigma = 1$) with the estimate $\widehat{\sigma}_{conj} = \Sexpr{conj_sig_max}$. In contrast, the estimate of $\sigma$ at $\mathtt{v0 = exp(-10)}$ for the independent case is $\widehat{\sigma}_{ind} = \Sexpr{ind_sig_max}$, much closer to the true value. 

\section{Conclusion}

In this vignette, we described the EMVS algorithm of Rockova and George (2014) with an independent prior formulation. We demonstrated how this algorithm can be applied using the \texttt{EMVS} R package.

\section{References}

George, E. I. and McCulloch, R. E. (1997), ``Approaches for Bayesian Variable Selection", \emph{Statistica Sinica}, 7, 339-373

Moran, G. E., Rockova, V. and George, E. I., (2018) ``On variance estimation for Bayesian variable selection", [arXiv:1801.03019]

Rockova, V. and George, E. I. (2014), ``EMVS: The EM approach to Bayesian variable selection," \emph{Journal of the American Statistical Association}, 109, 828-846

\end{document}
