\section{Introduction}

Much current research in genetics focuses on combining heterogeneous data for the same samples. This is prompted by the increasing availability of diverse molecular data types from a single study and should lead to a more complete understanding of biological systems, as combined inference for data from different molecular layers might unravel regulatory interactions within and across layers \citep{civelek2014systems}. An example is \emph{protein quantitative trait locus} (pQTL) analyses to detect associations between hundreds of thousands of genetic variants and hundreds of proteomic expression levels. When associated with disease genetic variants, proteins are often regarded as intermediate phenotypes or molecular proxies for the disease of interest, as they may provide direct insights on biological processes underlying the clinical condition, and likewise for other molecules such as genes or metabolites involved in so-called eQTL or mQTL analyses.

An important goal of molecular QTL studies is to detect genetic variants with \emph{pleiotropic effects}, i.e., variants that regulate the expression levels of several molecules, as the genome location where they lie may initiate essential functional mechanisms \citep{breitling2008genetical}. As pleiotropy has been acknowledged as a central property of genetic variants causing phenotypic variation, efforts are being made to uncover its activity patterns and gain understanding of the shared biological processes it induces \citep{sivakumaran2011abundant, solovieff2013pleiotropy}.
It is also of interest to identify those cases where the same molecule is simultaneously influenced by several genetic variants. This dual task requires a model that allows flexible selection, ideally performed jointly on the genetic variants and the molecular outcomes.

Although in practice univariate analyses still dominate, several proposals for joint modelling of multiple outcomes in genetic association problems have recently been made. \citet{flutre2013statistical} and \citet{zhou2014efficient} model the outcomes as multivariate responses having a matrix-variate distribution. The former rely on a Bayes factor framework to uncover the associations between a given genetic variant and any subset of outcomes, whereas the latter propose a linear mixed model whose random effect accounts for relatedness among individuals. While these methods show improvements over the fully marginal regression approach, they are restricted to problems with a few outcomes, because their models involve unstructured covariance matrices. \citet{o2012multiphen}'s MultiPhen method reverses the classical regression setup and fits a succession of models where each genetic variant is regressed on several outcomes. This eliminates the need to model large covariance matrices, but also entails a marginal treatment of the genetic variants. Moreover, MultiPhen does not penalize model complexity, which may cause instabilities when many outcomes are modelled. Molecular QTL problems are particularly complex, because in addition to the so-called $p \gg n$ paradigm, whereby the number of covariates (genetic variants) $p$ greatly exceeds the number of samples $n$, there is also a ``large $d$'' characteristic, as the number of responses (expression levels) $d$ is usually large. Methodologies to accommodate this are needed, especially when joint response modelling is sought.

Two-stage procedures are natural approaches to association problems with both $p$ and $d$ large. 2HiGWAS \citep{jiang20152higwas} is essentially an implementation of the screening strategy of \citet{fan2008sure} in the context of longitudinal outcomes. It consists of a dimension reduction step, where each genetic variant is tested against each outcome, followed by a penalized regression recast into functional mapping in which only the genetic variants remaining after screening are involved. While the second stage of 2HiGWAS is an interesting approach to joint covariate and response modelling, the method is not tailored to typical molecular QTL analysis, as it is designed for outcomes measured over time. Also, the fully marginal first stage screening, if too stringent, may cause important predictors to drop out of the analysis and thus lead to false negatives. Instead of pruning the covariate set, \citet{wang2016block} summarize information at outcome level in the context of eQTL analyses. They precluster the expression levels using a block-mixture model and test for association between each genetic variant and the resulting clustering. This approach requires that the stability of the clustering and its functional relevance are carefully checked, as the group pattern chosen is critical to subsequent analysis. More generally, two-stage procedures often rely on ad-hoc thresholding decisions at the first stage which influence the conclusions of the second stage.

The approaches sketched above use very diverse strategies to model predictor and outcome variables in high-dimensional settings. Trade-offs between realistic modelling and computational efficiency are inevitable, but two components seem critical to practical and powerful analyses: involving all variables in a single multivariate model, and maintaining interpretable inference. Unified selection of genetic variants and associated outcomes is then possible, unlike for most existing methods, whose focus is on selecting either predictors or outcomes. The Bayesian framework seems particularly suitable, as it offers flexible modelling possibilities in which biological beliefs may be naturally incorporated.
The hierarchical regression approach of \citet{richardson2010bayesian}, HESS (hierarchical evolutionary stochastic search), is an appealing example of such approaches; it can identify associations between hundreds of covariates and up to a few thousand responses from a single model. \citet{jia2007mapping} and \citet{scott2012integrated} propose methods called BAYES and iBMQ (integrated Bayesian hierarchical model for eQTL mapping) based on models similar to that of HESS, but a major drawback of all three approaches is the lack of scalability of their Markov chain Monte Carlo inference procedures. Problems whose size corresponds to actual genome-wide association studies with molecular outcomes (with hundreds of thousands of genetic variants and hundreds to thousands of outcomes and individuals) are out of reach even for HESS, which is based on adaptive parallel tempering/evolutionary Monte Carlo techniques. We are unaware of any fully multivariate approach that can deal with such data within a reasonable time.

In this paper, we propose a Bayesian inference strategy that avoids sampling, via an algorithm that is fast and whose convergence is easy to monitor, while having performance comparable with Markov chain Monte Carlo approaches. We describe a variational inference procedure for a model similar to that of \citet{richardson2010bayesian}, comprising a series of parallel linear regressions, one for each response, combined in a hierarchical manner to leverage shared information. A spike-and-slab prior \citep{ishwaran2005spike} is used to induce sparsity of the regression coefficients, and the probability that a given covariate affects any response is modelled through a parameter that is shared across responses. Interpretable posterior quantities, such as the probability of association of each covariate-response pair, are produced. An efficient algorithm for this model is crucial, as its number of parameters can be very large. Our variational approach can update the parameters jointly in a tractable manner. \citet{carbonetto2012scalable} provide a good discussion of variational inference in the context of genetic association and propose a variational regression method called ``varbvs'', which can be seen as a single-response counterpart of our approach. In our multiple response setting, we show by simulation that our procedure is accurate and reliable, and is better than existing methods at selecting variables for very large problems. Hence, it offers clear added-value in practice: it enables complex and flexible Bayesian inference based on large batches of genetic data, without having to prune them beforehand.

The paper is organized as follows. Section \ref{SecMod} describes the model, discusses its relation to earlier proposals, and presents a procedure to allow for sparsity control at both covariate and response levels. Section \ref{SecVB} gives an overview of variational Bayes approaches and describes our inference strategy. Section \ref{SecEQ} compares variational and Markov chain Monte Carlo inferences on the same model, also using direct approximations of posterior quantities. Section \ref{SecNS} describes numerical experiments for larger problems, comparing our method with several predictor selection methods, including the varbvs approach of  \citet{carbonetto2012scalable}, and with methods performing combined covariate and response selection, namely HESS \citep{richardson2010bayesian} and iBMQ \citep{scott2012integrated}. The section also presents a permutation-based comparison of our method with varbvs on a real mQTL problem. Section \ref{SecConcl} summarizes the discussion and highlights further possible developments.

Although the applicability of our method is not restricted to any particular context, all the numerical experiments presented in this paper use settings tailored to genome-wide association studies. The data-generation schemes are designed to embody common biological assumptions and are described in Section \ref{SecPred}. The method is implemented in the publicly available R package \texttt{locus}.

\section{Model and earlier proposals}\label{SecMod}

Let $y=(y_1, \ldots, y_d)$ be a $n \times d$ matrix of $d$ centered responses and let $X$ be a $n \times  p$ matrix of $p$ covariates, for each of $n$ samples. The covariates are $p$ genetic variants, more precisely single nucleotide polymorphisms, SNPs, and the responses might represent $d$ gene, protein or metabolite expression levels for $n$ individuals, depending on whether an eQTL, pQTL or mQTL problem is considered. Our model is intended to accommodate all the constraints entailed by molecular QTL analyses; it is adapted from that of HESS \citep{richardson2010bayesian}. Suppose that
\begin{eqnarray*}
y_t &\mid& \beta_t, \tau_t \sim \mathcal{N}_n\left(X\beta_t, \tau_t^{-1} I_n\right), \qquad\qquad\qquad \tau_t \overset{\text{ ind}}{\sim} \textup{Gamma}(\eta_t, \kappa_t)\,, \qquad\qquad\qquad t = 1, \ldots, d\,,\\
\beta_{st} &\mid& \gamma_{st}, \tau_t, \sigma^2 \sim \gamma_{st}\,\mathcal{N}\left(0, \sigma^2\,\tau_{t}^{-1}\right) + (1-\gamma_{st})\,\delta_0 \,,\qquad\qquad\qquad \sigma^{-2} \sim \textup{Gamma}(\lambda, \nu)\,,\quad\\
\gamma_{st}&\mid& \omega_{s} \overset{\text{ iid}}{\sim} \text{Bernoulli}\left(\omega_{s}\right), \qquad\qquad\qquad\qquad \omega_s  \overset{\text{ ind}}{\sim} \textup{Beta}(a_s, b_s)\,,\qquad\qquad\qquad\quad   s=1, \ldots, p\,,
\end{eqnarray*}
where $\delta_0$ is the Dirac distribution. Each response, $y_t$, is related linearly to the covariates and has a specific precision, $\tau_t$. The responses are conditionally independent across the regressions, but dependence among responses associated with the same covariates is captured through the prior specification of parameters $\omega_s$ and $\sigma^{-2}$, which are shared across the responses. This formulation circumvents modelling the covariance between the responses, which is infeasible when $d$ is large. Each covariate-response pair has its own regression parameter, $\beta_{st}$, for which sparsity is induced using a spike-and-slab prior. The binary parameter $\gamma_{st}$ acts as a ``pair selection'' indicator; covariate $X_s$ is associated with response $y_t$ if and only if $\gamma_{st}=1$. The parameter $\sigma$ represents the typical size of nonzero effects and is modulated by the residual variance, $\tau_t^{-1}$, of the response concerned by the effect. The parameters $\gamma_{s1},\ldots, \gamma_{sd}$ specify the response(s) associated with $X_s$ and are identically distributed as Bernoulli with common parameter $\omega_s$. Thus, $\omega_s$ controls the proportion of responses associated with covariate $X_s$. The goal of inference is variable selection. Selection of predictors can be performed by ranking the posterior means of $\{\omega_s\}$ and selection of covariate-response pairs can be performed by ranking the posterior probabilities of inclusion, i.e., the posterior means of~$\{\gamma_{st}\}$.

Our model differs from that of \citet{richardson2010bayesian} in two respects. One concerns the treatment of the regression coefficient parameters $\beta_{st}$: we use independent priors, whereas \citeauthor{richardson2010bayesian} rely on g-priors \citep{zellner1986assessing}. The main motivation for our choice is that the effects of genetic variants on a given outcome can be understood as causal, since no retroactive process can affect the variants, and they can take place at locations of the genome that are far apart, so their correlation structure need not reflect the spatial correlation of the SNPs; see \cite{guan2011bayesian}.
\citet{jia2007mapping} also rely on independent priors for the regression coefficients of BAYES, but they model the latter with a mixture of two normal distributions rather than a spike-and-slab prior and impose a residual variance parameter that is common to all responses. This stringent assumption may represent a weakness of their proposal.

The second difference concerns the third level of the model. \citet{richardson2010bayesian} opt for a quite complex specification, in which
\begin{equation}\label{EqPriorHess}\omega_{st}=\rho_s\omega_t\,, \qquad\quad \omega_t \sim \textup{Beta}(a_t, b_t)\,, \qquad\quad  \rho_s \sim \textup{Gamma}(c_s, d_s)\,,\qquad\quad 0 \leq \omega_{st} \leq 1\,.\end{equation}
In their case, the inclusion probability of $X_s$ for response $y_t$ is modelled through $\omega_t$; it is specific to that response but can be regulated using the parameter $\rho_s$, common to all responses.
\citet{jia2007mapping} and \citet{scott2012integrated} propose other variants for this prior. The former choose a treatment similar to ours, with $\omega_{st} \equiv \omega_s \sim \text{Dirichlet}(1,1)$, and the latter consider an additional level of hierarchy,
\begin{equation}\label{EqPriorIbqm}
\omega_{st} \mid  a_s, b_s, \pi_s \sim \pi_s \delta_0 + (1-\pi_s)\text{Beta}(a_s, b_s)\,, \qquad\qquad\quad \pi_s \sim \text{Beta}(a_0, b_0)\,,
\end{equation}
with $a_s \sim \text{Exp}(\lambda_a)$ and $b_s \sim \text{Exp}(\lambda_b)$. Our choice $\omega_{st}\equiv\omega_s \sim \textup{Beta}(a_s, b_s)$ is partly driven by our wish to design a simpler model and partly by practical considerations, since it ensures a closed form for our variational algorithm, unlike with (\ref{EqPriorHess}). While such a formulation was mentioned by \citet{richardson2010bayesian} and by \citet{scott2012integrated}, they did not pursue it because of concerns regarding its ability to control for multiplicity. Indeed, our model inherently enforces sparse associations as the number of responses, $d$, increases, but no control is achieved when the number of covariates, $p$, grows. We address this below by providing a procedure to induce a correction through the prior of $\omega_s$.

Part of the flexibility of our model comes from the fact that the hyperparameters, $a$, $b \in \mathbb{R}^p_+$ (for $\omega$'s Beta prior), $\lambda$, $\nu \in \mathbb{R}_+$ (for $\sigma^{-2}$'s Gamma prior) and $\eta$, $\kappa \in \mathbb{R}^d_+$ (for $\tau$'s Gamma prior) are readily interpreted. One option is to set them based on external information regarding the likelihood of given associations, if available. For instance, to favour associations with covariate $X_s$, one can set $a_s$ and $b_s$ so that the prior proportion of responses affected by $X_s$, $\text{E}(\omega_s)$, is large. The use of such assumptions may be very efficient, but it may also skew the inference towards existing knowledge. In the simulations presented in this paper, we assume that the regression and variance parameters are exchangeable, i.e., that all covariates, and responses have the same prior propensity to be involved in associations, by selecting a single value for all components of $a$, $b$, $\eta$ and $\kappa$. Without favouring any covariate or response, however, we can control signal sparsity at the level of covariates by specifying (possibly through cross-validation) a prior average number of covariates, $p^*$, expected to be included in the model. Setting
 \begin{equation}\label{EqHyperAdjust}
 a_s \equiv 1\,, \qquad\quad b_s \equiv d(p-p^*)/p^*\,, \qquad\quad 0 < p^* < p\,,
 \end{equation}
 the prior probability that $X_s$ is associated with at least one response is
$$p\left( \cup_{t=1}^d\{\gamma_{st}=1\} \right) = 1 - \frac{\prod_{j=1}^d (b_s + d -j )}{\prod_{j=1}^d ( a_s + b_s + d -j )} = \frac{p^*}{p}\,,$$
and simpler models are favoured as $p$ increases. To see this, one can consider the prior odds ratio representing the support for a model
to have an additional response associated with~$X_s$, i.e.,
\begin{equation}\label{EqPOR}
\text{POR}(q_s-1:q_s)=\frac{p\left( \sum_{t=1}^d\gamma_{st}=q_s-1\right)}{p\left( \sum_{t=1}^d\gamma_{st}=q_s\right)}= \frac{b_s + d - q_s }{ a_s + q_s -1 }\,,\qquad\qquad q_s=1, \ldots, d\,.
\end{equation}
Clearly, penalties arise and increase with the total number of responses in the model, $d$. Figure~\ref{FigPOR} displays $(\ref{EqPOR})$ for $q_s=1, \ldots, 5$ as a function of $p$ and indicates that, when $a_s$ and $b_s$ are specified as in $(\ref{EqHyperAdjust})$, the penalties also increase with the total number of covariates, $p$, therefore naturally adjusting for multiplicity. Moreover, the penalties are not uniform when moving from one to two responses associated with $X_s$, or from four to five, for instance.

\begin{figure
\centering
  \noindent\includegraphics[scale=1]{./VB_paper_Fig1-eps-converted-to.pdf}
  \caption{\small Prior odds ratios, $\text{POR}(q_s-1:q_s)$, for $q_s=1,\ldots, 5$, $a_s$ and $b_s$ as in (\ref{EqHyperAdjust}), $d=100$, $p^*=2$, and for a total number of covariates ranging from $p=5$ to $5,000$; see \citet{scott2010bayes} for a similar visualization of prior odds ratios in a single response context.}\label{FigPOR}
\end{figure}

The experiment reported in Table \ref{TabAdj} confirms that adjustment takes place in practice. It considers problems with $p_0=20$ ``active'' covariates, i.e., those associated with at least one response, and an increasing number of ``noise'' covariates and it compares the regime with $a_s$ and $b_s$ set according to (\ref{EqHyperAdjust}) to an ``uncorrected'' regime with $a_s\equiv1$, $b_s\equiv 2d-1$, so that the prior mean number of responses associated with $X_s$ is $0.5$, i.e., $\text{E}(\omega_s) \equiv (2d)^{-1}$. The number of false positives, based on a posterior probability of inclusion greater than $0.5$, grows linearly with $p$ when the uncorrected model is used but remains roughly constant close to zero with correction (\ref{EqHyperAdjust}), giving a clear multiplicity adjustment. Other experiments confirmed strong sparsity control; the reported findings on real data should therefore be plausible when (\ref{EqHyperAdjust}) is used.

\begin{table
\begin{center}
\captionof{table}{\small Multiplicity adjustment at covariate level. The mean numbers of false positives (FP) and true positives (TP) obtained with the uncorrected and corrected regimes are compared for $p_0=20$ active covariates and an increasing number of noise covariates, $p-p_0$. The total number of responses is $d=25$. $64$ replicates were performed; standard errors are in parentheses.\label{TabAdj}}
\footnotesize
\begin{tabular}{rccccc}
  \hline
 $p$ & $50$ & $250$ & $500$ & $1,000$ & $2,500$ \\
  \hline
 Mean \# of FP \\Uncorrected & $0.61$ ($0.66$) & $5.23$ ($2.24$) & $10.58$ ($3.16$) & $22.38$ ($4.48$) & $52.61$ ($8.39$) \\
  Corrected & $0.77$ ($0.81$) & $0.70$ ($0.99$) & $0.61$ ($0.77$) & $0.39$ ($0.58$) & $0.44$ ($0.59$)\\
  Mean \# of TP \\
  Uncorrected & $19.95$ ($0.21$) & $20.00$ ($0.00$) & $19.98$ ($0.12$) & $20.00$ ($0.00$) & $19.94$ ($0.24$) \\
  Corrected & $19.97$ ($0.18$)& $19.91$ ($0.34$) & $19.81$ ($0.43$) & $19.77$ ($0.56$) & $19.38$ ($0.86$)\\
   \hline
\end{tabular}
\end{center}
\end{table}

\section{Variational inference}\label{SecVB}

Section \ref{SecMod} described some differences between our model and those of \citet[][BAYES]{jia2007mapping}, \citet[][HESS]{richardson2010bayesian} and \citet[][iBMQ]{scott2010bayes}, but a more fundamental distinction concerns the inference procedure. The three earlier methods rely on Markov chain Monte Carlo (MCMC) techniques and require massive computing resources when the dimensionality of the problem is large. We instead employ a variational inference procedure, which is deterministic and hence can be much cheaper \citep{ormerod2010explaining}.

Instead of sampling from the joint posterior probability $p(\theta \mid y)$ of the parameter vector of interest, $\theta$,
variational approaches proceed by replacing it by a tractable analytical approximation, $q(\theta)$. We focus on so-called \emph{mean-field} variational formulations \citep{xing2002generalized, attias2000variational} to construct such a class, i.e., we assume that $q(\theta)$ factorizes over some partition of $\theta$, $\{\theta_j\}_{j=1,\ldots, J}$,
$$q(\theta) = \prod_{j=1}^J q_j(\theta_j)\,;$$
no further assumption is made about the distribution, and in particular no constraint is imposed on the functional forms of the $q_j(\theta_j)$. Here, we consider the factorization
\begin{eqnarray}\label{EqFact}
q\left(\beta, \gamma, \tau, \sigma^{-2}, \omega\right) =  \left\{ \prod_{s=1}^p\prod_{t=1}^d  q(\beta_{st}, \gamma_{st}) \right\}\left\{\prod_{s=1}^p q(\omega_s)\right\}\left\{\prod_{t=1}^{d} q(\tau_t) \right\} q\left(\sigma^{-2}\right)\,,
\end{eqnarray}
and turn the inference into an optimization problem where $q(\theta)$ is obtained by minimizing its Kullback--Leibler divergence $\text{KL}\left(q \middle\| p\right)$ from the target distribution, $p(\theta \mid y)$.
Because the marginal log-likelihood may be written as
\begin{equation}\label{EqKL}
\log p(y) = \mathcal{L}(q) + \text{KL}\left(q \middle\| p\right)\,,
\end{equation}
where
$$\mathcal{L}(q) = \int q(\theta) \log\left\{\frac{p(y, \theta)}{q(\theta)}\right\}\text{d}\theta\,,\qquad \text{KL}\left(q \middle\| p\right) = - \int q(\theta) \log\left\{ \frac{p(\theta\mid y)}{q(\theta)}\right\} \text{d}\theta\,, $$
minimizing the Kullback--Leibler divergence amounts to maximizing $\mathcal{L}(q)$, which represents a lower bound for $\log p(y)$. To this end, we observe that
\begin{eqnarray}
\mathcal{L}(q) &=& \int \prod_{k=1}^J q_k(\theta_k) \left\{ \log p(y,\theta) - \sum_{k=1}^J\log q_k(\theta_k) \right\} \text{d}\theta_1\cdots\text{d}\theta_J\nonumber\\
&=& \int q_j(\theta_j) \left\{ \int \log p(y,\theta) \prod_{k\neq j} q_k(\theta_k) \text{d}\theta_k - \log q_j(\theta_j) \right\} \text{d}\theta_j + \text{cst} \nonumber\\
&=& \int q_j(\theta_j) \log\left\{\frac{p_{-j}(\theta_j; y)}{q_j(\theta_j)}\right\} \text{d}\theta_j + \text{cst}\,, \qquad\qquad\qquad\qquad\qquad\quad j =1, \ldots, J\,,\label{EqExpKL}
\end{eqnarray}
where $\text{cst}$ is constant with respect to $\theta_j$ and where we introduced the distribution
$$p_{-j}(\theta_j; y) = \text{cst} \times \exp\left[ \text{E}_{-j}\left\{ \log p(y, \theta) \right\}\right]\,,$$
with $\text{E}_{-j}\{\cdot\}$ denoting the expectation with respect to the distributions $q_k$ over all variables $\theta_k$, $k\neq j$. The right-hand side of (\ref{EqExpKL}) corresponds to the negative Kullback--Leibler divergence between $q_j(\theta_j)$ and $p_{-j}(\theta_j; y)$, plus a constant. Hence, assuming that the $q_k(\theta_k)$, $k\neq j$, are fixed, the distribution $q_j(\theta_j)$ which maximizes $\mathcal{L}(q)$ is $q_j(\theta_j) = p_{-j}(\theta_j; y)$, i.e., the maximum of $\mathcal{L}(q)$ occurs when
\begin{equation}\label{EqQj}
\log q_j(\theta_j) = \text{E}_{-j} \{\log p(y, \theta)\} + \text{cst}\,,\qquad\qquad j=1, \ldots, J\,.
\end{equation}
The relations $(\ref{EqQj})$ give rise to cyclic dependencies among the densities $q_j(\theta_j)$. This suggests an iterative algorithm whose convergence can easily be monitored by evaluating changes in the lower bound $\mathcal{L}(q)$. Our choice (\ref{EqFact}) ensures that the coordinate updates can be derived in closed form; in particular, the semi-conjugacy of our model implies that the prior densities of all parameters are preserved by the variational densities. For instance, a spike-and-slab distribution with modified parameters is recovered at posterior level, $q(\beta_{st}, \gamma_{st}) = q(\beta_{st} \mid \gamma_{st}) q(\gamma_{st})\,,$ with
$$\beta_{st} \mid \gamma_{st}=1, y\; \sim \mathcal{N}\left(\mu_{\beta, st}, \sigma^2_{\beta, st}\right)\,,\qquad \beta_{st} \mid \gamma_{st}=0, y\; \sim \delta_0\,,\qquad \gamma_{st} \mid y \sim \text{Bernoulli}\left(\gamma^{(1)}_{st}\right)\,,$$
where the \emph{variational parameters} $\mu_{\beta, st}$, $\sigma^2_{\beta, st}$, $\gamma^{(1)}_{st}$ are to be updated iteratively. Convergence is ensured by the convexity of $\mathcal{L}(q)$ in each of the $q_j(\theta_j)$ \citep[\S\S\,3.1.5, 3.2.4, 3.2.5]{boyd2004convex}. The algorithm and its derivation are given in Appendix \ref{AppB}.

\section{Empirical quality assessment of the variational approximation}\label{SecEQ}

\subsection{Tightness of the marginal log-likelihood lower bound}\label{SecLB}

In this section we evaluate the closeness of the variational density $q$ to the target posterior distribution by approximating the Kullback--Leibler divergence $\text{KL}\left(q \middle\| p\right)$. Because of relation $(\ref{EqKL})$, this amounts to assessing the tightness of the variational lower bound for the marginal log-likelihood, $\mathcal{L}(q)$.
For small problems, the likelihood $p(y)$ may be accurately approximated using simple Monte Carlo sums.
We have
$$
p(y) = \int\cdots\int \mathrm{d}\omega\,\mathrm{d}\sigma^{-2}  \, \left\{ \prod_{s=1}^p\,p( \omega_s)\right\} \,p\left(\sigma^{-2}\right)\prod_{t=1}^d \left\{\sum_{\gamma_t \in \{0,1\}^p}  p\left(y_t\mid\gamma_t, \sigma^{-2}\right)  \prod_{s=1}^p\,p\left(\gamma_{st} \mid\omega_s\right) \right\}\,,
$$
with
{
$$p\left(y_t\mid\gamma_t, \sigma^{-2}\right) = \begin{cases} \displaystyle (2\pi)^{-n/2}\, \Gamma\left(\frac{n}{2} + \eta_t\right) \frac{\kappa_t^{\eta_t}}{\Gamma(\eta_t)} \left(\kappa_t + \frac{\Vert y_t \Vert^2}{2}\right)^{-n/2-\eta_t}, \hspace{2cm} q_{\gamma_t}=0\,,
\\\text{}\\\displaystyle(2\pi)^{-n/2} \begin{vmatrix} V_{\gamma_t,\sigma^{-2}}\end{vmatrix}^{-1/2} \Gamma\left(\frac{n}{2} + \eta_t\right) \frac{\kappa_t^{\eta_t}}{\Gamma(\eta_t)} \left(\kappa_t +\frac{S^2_{\gamma_t}}{2} \right)^{-n/2-\eta_t} \left(\sigma^{-2}\right)^{q_{\gamma_t}/2}, \\ \hspace{10.1cm}\text{otherwise}\,,
\end{cases}$$
}
where $$q_{\gamma_t} = \sum_{s=1}^p \gamma_{st}\,,\qquad\quad V_{\gamma_t,\sigma^{-2}} = X_{\gamma_t}^TX_{\gamma_t} + \sigma^{-2}I_{q_{\gamma_t}}\,, \qquad\quad
S^2_{\gamma_t,\sigma^{-2}} = \Vert y_t\Vert^2 - y_t^T X_{\gamma_{t}} V_{\gamma_t,\sigma^{-2}}^{-1} X_{\gamma_{t}}^T y_t\,;$$
see Appendix \ref{AppC} for details. As no closed form is available for the remaining integrals, we use
$$
p(y) \approx \frac{1}{I} \sum_{i=1}^I \prod_{t=1}^d \left\{ \sum_{\gamma_t \in \{0,1\}^p}  p\left(y_t\mid\gamma_t, \left(\sigma^{-2}\right)^{(i)}\right)  \prod_{s=1}^p\,p\left(\gamma_{st} \mid\omega_s^{(i)}\right)\right\}\,,
$$
where we independently generate
\begin{equation}\label{EqSMC}\left(\sigma^{-2}\right)^{(i)} \sim \text{Gamma}(\lambda, \nu)\,, \qquad \omega_s^{(i)} \sim \text{Beta}(a_s, b_s)\,,\qquad s=1, \ldots, p\,,\qquad i=1, \ldots, I\,.\end{equation}

Figure \ref{FigLB} displays the relative difference $\left\{\log p(y) - \mathcal{L}(q)\right\}/\log p(y)$ for problems with $p=5$ covariates, $d=6$ responses and increasing sample sizes, $n$. In the left panel, the covariates are independent of each other, and so are the responses. In the right panel, the covariates are equicorrelated with correlation coefficient $\rho=0.75$, and so are the responses. In both cases, the mean relative difference is below $1\%$ with $n=50$ and seems to decrease as $n$ grows. Although we are not aware of any such study with which to benchmark our results, these values seem very small, suggesting that our variational distribution $q$ adequately reflects the target distribution $p$, at least for small problems. Likewise, the variational lower bound $\mathcal{L}(q)$ may be used as a proxy for the marginal log-likelihood when performing model selection; this use will be illustrated in Section \ref{SecReal}. The fact that the variational lower bound remains tight in the correlated data case is reassuring, as it suggests that the independence assumptions underlying the mean-field factorization of $q$ may only weakly impact the quality of the approximation.

\begin{figure
\centering
  \noindent\includegraphics[scale=1]{./VB_paper_Fig2_left-eps-converted-to.pdf}\hspace{0.4cm}
  \includegraphics[scale=1]{./VB_paper_Fig2_right-eps-converted-to.pdf}
    \caption{\small $\text{Log}_{10}$ relative difference between the marginal log-likelihood and the variational lower bound. Left: independent covariates and responses. Right: correlated covariates and responses, $\rho=0.75$. Problems with $p=5$ covariates, of which $p_0=3$ randomly selected as ``active'' (associated with at least one response), and $d=6$ responses, of which $d_0=3$ ``active'' (associated with at least one covariate). Each active covariate is associated with an additional active response with probability $0.25$ and explains on average $3.5\%$ of the variance of its corresponding response(s). The number of draws for the simple Monte Carlo approximations is $I=50,000$; the number of replicates for each sample size is $150$.}\label{FigLB}
\end{figure}

\subsection{Comparison with Markov Chain Monte Carlo}\label{SecComp}

\begin{table
\begin{center}
\captionof{table}{\small Variational Bayes (VB), MCMC and simple Monte Carlo estimates for $\beta$ and $\omega$ (components corresponding to noise averaged). Standard errors are in parentheses. \label{TabBetaOmega}}
\footnotesize
\begin{tabular}{rccccccc}
  \hline
  & & & Active & & & Inactive\\
 $10 \times$ & $\beta_{1,2}$ & $\beta_{2,1}$ & $\beta_{3,2}$ & $\beta_{4,1}$ & $\beta_{4,2}$ & $\beta_{\text{rest}}$ (avg) \\
  \hline
Truth & $-1.75$ & $2.87$ & $2.37$ & $3.73$ & $-4.76$ & $0.00$ \\
  VB & $-1.74$ ($0.01$) & $1.86$ ($0.01$) & $1.70$ ($0.02$) & $2.26$ ($0.01$) & $-3.48$ ($0.01$) & $0.02$ ($0.04$) \\
  MCMC & $-1.74$ ($0.33$) & $1.86$ ($0.32$) & $1.69$ ($0.34$) & $2.26$ ($0.32$) & $-3.48$ ($0.34$) & $0.02$ ($0.14$) \\
   \hline
      \hline
    & & &  Active & & & Inactive\\
 & $\omega_1$ & $\omega_2$ & $\omega_3$ & $\omega_4$ & &$\omega_{\text{rest}}$ (avg) \\
  \hline
  True prop. of active resp. & $0.2$ & $0.2$ & $0.2$ & $0.4$ & & $0$ \\
VB & $0.21$ ($0.14$) & $0.25$ ($0.15$) & $0.19$ ($0.14$) & $0.33$ ($0.17$) & & $0.03$ ($0.06$) \\
  MCMC & $0.23$ ($0.18$) & $0.26$ ($0.18$) & $0.21$ ($0.16$) & $0.35$ ($0.18$) & & $0.04$ ($0.08$) \\
  Simple Monte Carlo & $0.25$ & $0.26$ & $0.21$ & $0.35$ & & $0.05$ \\
   \hline
\end{tabular}
\end{center}
\end{table}

We complement our quality assessment by comparing several variational posterior quantities with those for MCMC inference on problems of moderate size. A fair comparison is not straightforward, as these two types of inference rely on stopping rules and convergence diagnostics of very different natures. While the convergence criterion for variational inference comes down to a tolerance to be prescribed, the ability of MCMC sampling to adequately explore the model space for a given chain length can be difficult to evaluate, and usually varies greatly with the problem size. To alleviate the risk of inaccurate MCMC inference, we run $10^5$ iterations and discard the first half. We also support our comparison with selected quantities approximated by simple Monte Carlo sums, namely, the posterior probability of inclusion of a covariate $X_s$ for a response~$y_t$,
\begin{eqnarray}
p( \gamma_{st}=1 \mid y ) &=& \frac{1}{p(y)} \;\frac{1}{I}\sum_{i=1}^I \left[ \prod_{t'\neq t} \left\{\sum_{\gamma_{t'} \in \{0,1\}^p}  p\left(y_{t'}\mid\gamma_{t'}, \left(\sigma^{-2}\right)^{(i)}\right)  \prod_{s'=1}^p\,p\left(\gamma_{s't'} \mid\omega_{s'}^{(i)}\right) \right\} \right.
\nonumber\\
&&\hspace{0.15cm}\times\left. \left\{\sum_{\gamma_{t} \in \{0,1\}^p:\; \gamma_{st}=1}  p\left(y_{t}\mid\gamma_{t}, \left(\sigma^{-2}\right)^{(i)}\right)  \prod_{s'=1}^p\,p\left(\gamma_{s't} \mid\omega_{s'}^{(i)}\right) \right\}\right]\,,\nonumber
\end{eqnarray}
and the posterior mean of $\omega_s$, controlling the proportion of responses associated with covariate~$X_s$,
\begin{eqnarray}
\text{E}( \omega_{s} \mid y ) &=&  \frac{1}{p(y)} \;\frac{1}{I}\sum_{i=1}^I \;\omega_{s}^{(i)}\;\prod_{t=1}^d \left\{\sum_{\gamma_t \in \{0,1\}^p}  p\left(y_t\mid\gamma_t, \left(\sigma^{-2}\right)^{(i)}\right)  \prod_{s'=1}^p\,p\left(\gamma_{s't} \mid\omega_{s'}^{(i)}\right) \right\}\,,\nonumber
\end{eqnarray}
 with the samples $\left(\sigma^{-2}\right)^{(i)}$ and $\{\omega_s^{(i)}\}$ generated as in $(\ref{EqSMC})$, with $I=2 \times 10^5$ draws.

Table $\ref{TabBetaOmega}$ reports the variational, MCMC and simple Monte Carlo estimates of $\beta$ and $\omega$ for a problem with $p=8$ covariates, $d=5$ responses for $n=250$ samples, and with each nonzero association explaining on average $13.5\%$ of response variance.
The estimates all agree closely. Those of the five active regression coefficients, $\beta_{1,2}$, $\beta_{2,1}$, $\beta_{3,2}$, $\beta_{4,1}$ and $\beta_{4,2}$, are significantly different from zero, unlike the average estimate of the inactive coefficients. A plot of the MCMC and variational posterior densities (the latter obtained in closed form), given in Figure \ref{FigPostDens} of Appendix \ref{AppSMCp}, shows that the posterior modes of the inactive coefficients are all zero. Moreover, in this case the variational distributions are usually solely made up of a clear spike at zero, whereas the MCMC histograms correspond roughly to a centered Gaussian distribution with average standard deviation $0.014$. Table $\ref{TabBetaOmega}$ also indicates a shrinkage effect for both variational and MCMC posterior means of the nonzero $\beta$ compared to the true values. This is a consequence of the spike-and-slab prior but does not seem to hamper the detection of the association signals, since the posterior probabilities of inclusion of the true nonzero associations are concentrated around $1$, while those corresponding to noise are usually much lower, whether obtained by MCMC, variational or simple Monte Carlo procedures; see Figure \ref{FigPPI} of Appendix \ref{AppSMCp}.
Finally, the estimates of $\{\omega_s\}$ in Table $\ref{TabBetaOmega}$ provide a fair approximation to the actual proportion of responses associated with a given covariate.

Two additional numerical experiments comparing variational and MCMC posterior quantities are provided in Appendix \ref{AppSMCp}. One compares the estimates of $\omega$ and $\tau$ with the true values when the data are generated from the model with $p=100$ covariates and $d=~10$ responses. It also provides receiver operating characteristic (ROC) curves assessing the pairwise variable selection performance for both inference types. The other simulation gathers the observed values, $y$, and the estimated posterior means of $X\beta$ obtained by variational and MCMC procedures and an oracle. Both experiments indicate equivalent performance for MCMC and variational inferences.

\section{Statistical performance}\label{SecNS}

\subsection{Predictor selection}\label{SecPred}

The problems considered in Section \ref{SecComp} were small enough to allow accurate and tractable MCMC inference. In this section, we assess the performance of our approach on larger problems by comparing it to popular variable selection methods; i.e., with joint modelling of outcomes and covariates (elastic net for multivariate Gaussian responses), with joint modelling of covariates only (Bayesian multiple regression based on MCMC inference, ``BAS'', or variational inference, ``varbvs''), or with fully marginal modelling (univariate ordinary least squares and ``lmBF'' Bayesian regressions). Complete descriptions and references are in Appendix \ref{AppMeths}. The methods are compared by measuring their ability to detect the active covariates, i.e., to determine which covariates are associated with at least one response. For our variational approach, this task is achieved by ranking the posterior means of the $\omega_s$, which control the proportion of responses associated with a given covariate.

Our data-generation design is based on generally accepted principles of population genetics. We simulate SNPs under Hardy--Weinberg equilibrium from a binomial distribution with probabilities corresponding to minor allele frequencies of common variants, chosen in the interval $(0.05, 0.5)$ uniformly at random, and we generate outcomes from Gaussian distributions with specific error variances. The dependence structure of these variables is either enforced block-wise with preselected auto- or equicorrelation coefficients or chosen to be that of real data. The labels of the active SNPs and outcomes are picked randomly, and each active SNP is associated to one (randomly selected) active outcome and to each of the remaining active outcomes with a prescribed probability; some outcomes are therefore under pleiotropic control. The proportions of outcome variance explained per SNP are simulated for all associations from a positively skewed Beta distribution to favour the generation of smaller effects, and they are then rescaled to match a given average proportion. To mimic the result of natural selection, the effect sizes are inversely related to the SNP minor allele frequencies. For more details, see Appendix~\ref{AppGen}.

 \begin{figure
\centering
  \noindent\includegraphics[scale=0.72]{./VB_paper_Fig3_left-eps-converted-to.pdf}\quad
  \includegraphics[scale=0.72]{./VB_paper_Fig3_middle-eps-converted-to.pdf}\quad
 \includegraphics[scale=0.72]{./VB_paper_Fig3_right-eps-converted-to.pdf}
    \caption{\small Truncated average receiver operating characteristic curves with $95\%$ confidence intervals for predictor selection obtained from $48$ replications. The competing methods are used in three studies of different sizes, based on their computational tractability. Left: $p=5,000$ covariates spatially autocorrelated with correlation coefficient $\rho_X = 0.75$, $d=50$ outcomes equicorrelated by blocks with four blocks of equal sizes and correlation coefficients $\rho_Y=0.8, 0.3, 0.2$ and $0.5$, $p_0=100$ active covariates, $d_0=40$ active outcomes, $n=250$ observations, probability of association with an additional outcome $p_{\text{add}} = 0.15$, average outcome variance percentage explained by the active covariates $p_{\text{ve}} = 30.0\%$. Middle: $p=500$ independent covariates, $d=20,000$ outcomes equicorrelated by blocks of size $10$ with $\rho_X \in \{0.5, \ldots, 0.8\}$, $p_0 = 300$, $d_0 = 12,500$, $n = 300$, $p_{\text{add}} = 0.01$, $p_{\text{ve}} = 55.8\%$. Right: $p=150,000$ covariates autocorrelated by blocks of size $100$ with $\rho_X \in \{0.5, \ldots, 0.9\}$, $d = 200$ outcomes with same correlation structure than real protein expression levels \citep[Diogenes study,][see Appendix \ref{AppPre}]{2010_Diogenes_Obesity_rev}, $p_0 = 500$, $d_0=150$, $n=200$, $p_{\text{add}} = 0.05$, $p_{\text{ve}} = 62.6\%$. The univariate ordinary least squares and varbvs curves overlap.
}\label{FigROCAll}
\end{figure}

We perform $48$ replications for each of three simulation configurations. The first configuration has moderate numbers of covariates ($p=5,000$) and outcomes ($d=50$), and allows time-consuming methods to run within hours. The second has many outcomes ($d=20,000$) and the third has many covariates ($p=150,000$); these numbers approach those encountered in molecular QTL studies. The remaining settings (numbers of active outcomes and covariates, of observations, effect sizes, etc) are detailed in the caption to Figure \ref{FigROCAll}.

The ROC curves in Figure \ref{FigROCAll} indicate that our approach outperforms the other methods. It is appreciably more powerful for low false positive rates, which are of particular interest for the highly sparse scenarios typically expected for genome-wide association studies. Despite the correlation among the covariates and the outcomes, our method does not seem to suffer from the independence assumptions implied by the mean-field approximation, as suggested by the results of Section \ref{SecLB}. The marginal ordinary least squares and marginal lmBF regressions appear to miss many associations because of their univariate modelling of covariates, but jointly accounting for the covariates may not suffice, as suggested by the rather poor performances of the Bayesian multiple regression approaches, BAS and varbvs, which apply separate multiple linear regressions for each outcome. It appears that the ability of our approach to exploit the similarity across outcomes yields more power to detect their shared associations. Finally, even though the multivariate elastic net models jointly the covariate and outcome variables, its inference suffers from the assumption that to each covariate corresponds a single regression coefficient, shared for all responses. As a consequence, regression estimates of covariates with weak or few associations with the responses may be shrunk to zero.

\subsection{Combined selection of predictors and outcomes}

\begin{figure
\centering
  \noindent \includegraphics[scale=0.97]{./VB_paper_Fig4_left-eps-converted-to.pdf}
      \includegraphics[scale=0.97]{./VB_paper_Fig4_right-eps-converted-to.pdf}
          \caption{
  \small Marginal posterior probabilities of inclusion (PPI) obtained by our approach, and those of HESS (left) and those of iBMQ (right), for a problem with $p=250$ covariates of which $p_0=50$ are active, with $d=100$ outcomes, of which $d_0= 50$ are active, and $n=250$ samples.}\label{FigVBHESSiBMQ}
\end{figure}

Unlike the classical variable selection methods used as comparators in Section \ref{SecPred}, our approach and those of \citet{richardson2010bayesian}, HESS, and \citet{scott2012integrated}, iBMQ, are tailored to molecular QTL problems: they quantify the associations between each covariate-response pair in a single model, and thus provide flexible and unified frameworks for detecting pairs of associated SNP-molecules, as well as pleiotropic SNPs associated with many molecular outcomes.
In this section, we compare the three approaches in terms of the posterior quantities used to perform such selection.
As both HESS and iBMQ rely on MCMC sampling, we consider smaller problems than in Section \ref{SecPred} in order to ensure convergence within a reasonable time. The simulated datasets have $p=250$ covariates, of which $p_0=50$ are active, and $d=100$ outcomes, of which $d_0=50$ are active, the probability of association being $0.05$, for $n=250$ samples. On average, the active covariates account for $22\%$ of the variance of an outcome with which they are associated. HESS was run with three MCMC chains, the number selected by the authors for their simulations but with $50,000$ iterations of which $25,000$ were discarded as burn-in. For iBMQ, $50,000$ iterations were saved after removal of $50,000$ burn-in samples, as suggested in the package documentation for a problem of comparable dimensions. Inference for one replication took on average $10$ seconds with our method, around $21$ minutes with iBMQ and $4$ hours with HESS (GPU computation option disabled, since no GPU was available to us) on an Intel Xeon CPU at 2.60 GHz with 64 GB RAM.

\begin{table
\begin{center}
\captionof{table}{\small Mean true positive rate (TPR) and true negative rate (TNR) for our approach, HESS and iBMQ based on median probability models. Settings: $p=250$, $p_0=50$, $d=100$, $d_0=50$, $n=250$, $48$ replicates. Standard errors are in parentheses.\label{TabTPR}}
\footnotesize
  \begin{tabular}{rcc}
  \hline
 $100 \times$ & TPR  & TNR \\
  \hline
  VB  & $58.9$ ($5.0$) & $99.9$ ($0.0$) \\
  HESS & $57.9$ ($5.3$) & $99.9$ ($0.0$)   \\
  iBMQ & $0.1$ ($0.2$) & $99.8$ ($0.0$)  \\
   \hline
\end{tabular}
\end{center}
\end{table}

Figure \ref{FigVBHESSiBMQ} compares the marginal posterior probabilities of inclusion obtained by our method with those of HESS and iBMQ. We observe a strong correlation between our approach and HESS, with a quite good ability to discriminate between active and inactive covariate-response pairs. There is a discrepancy at the zero ordinate, where HESS signals a series of false positives and few true positives. The comparison with iBMQ is more contrasted, as the values of its posterior probabilities of inclusion for many true associations are below $0.1$ and indistinguishable from noise. The same conclusions are reached when running the three methods on $47$ additional datasets, as suggested by Table \ref{TabTPR}, which gathers sensitivity and specificity measures based on median probability models \citep{barbieri2004optimal} (consisting of those covariate-response pairs whose posterior inclusion probability is higher than $0.5$).  As discussed in Section~\ref{SecMod}, control of signal sparsity can be induced through the prior for $\omega_s$, still, rather than median probability models, one may prefer to use a data-driven false discovery threshold in order to prescribe a desired level of false discoveries.

Figure \ref{FigVBHESS} compares the patterns recovered by HESS and by our method, again based on the marginal posterior probabilities of inclusion from the first replicate. Visual comparison of the true positive rates suggests that the abilities of the two approaches to detect the true associations are very similar. Our approach indicates the presence of associations in the region of active covariates only, whereas the HESS pattern is blurrier in regions of inactive covariates. The posterior means of $\{\omega_s\}$ from our approach discriminate quite well between active and inactive covariates, and so do, for HESS, the posterior probabilities $\text{pr}(\rho_s >1 \mid y)$ ($s=1, \ldots, p$), described by \citet{richardson2010bayesian} as capturing the propensity for a given covariate to influence several responses simultaneously.

\begin{figure
\centering
    \noindent\includegraphics[scale=0.96]{./VB_paper_Fig5_left-eps-converted-to.pdf}
      \includegraphics[scale=0.96]{./VB_paper_Fig5_right-eps-converted-to.pdf}
    \caption{\small Posterior quantities for detection of associations with HESS (left) and with our approach (right), for a simulated dataset with $p=250$ independent covariates (here only the first $100$ are displayed), of which $p_0=50$ are active and placed first, with $d=100$ responses ($50$ active) and $n=250$ individuals. Marginal posterior probabilities of inclusion (central panel), true positive rates for predictor and response selection based on posterior probability of inclusion being $> 0.5$ (bottom and right panels), posterior probability $\text{pr}(\rho_s > 1 \mid y)$ for HESS and posterior mean $\text{E}_q(\omega_s\mid y)$ for our approach (left panel). The simulated associations are shown by red crosses. \label{FigVBHESS}}
\end{figure}

\subsection{Application to a real mQTL dataset}\label{SecReal}

We end these numerical experiments by illustrating our approach on data from a large multicenter dietary intervention study called Diogenes \citep{2010_Diogenes_Obesity_rev}. The study contains a series of genomic data types collected at different stages of a dietary treatment provided to the cohort. Its goal is to uncover molecular mechanisms underlying the metabolic status of overweight individuals and improve understanding of the factors predisposing weight regain after a diet. Here, we perform a metabolite quantitative trait locus (mQTL) analysis; in this context, the metabolites may be viewed as proxies for the clinical condition of interest, weight maintenance. We also use this illustration on real data to further highlight the benefits of modelling the outcomes jointly via an extensive permutation-based comparison with the single-response variational method varbvs \citep{carbonetto2012scalable}.

After quality control, the data consist of $p=215,907$ tag SNPs and $d=125$ metabolite expression levels, adjusted for age, center and gender, for $n=317$ individuals. The SNPs were collected on Illumina HumanCore arrays and the metabolites were quantified in plasma using liquid chromatography-mass spectrometry (LC-MS). They span cholesterol esters (CholE), phosphatidylcholines (PC), phosphatidylethanolamines (PE), sphingomyelins (SM), di- (DG) and triglycerides (TG). Appendix \ref{AppPre} provides more details.

\begin{table
\begin{center}
\captionof{table}{\small Number of associations declared by our method and by varbvs, and number of signals in common at selected permutation-based false discovery rates. For each case, the number of associations also declared by univariate screening at Benjamini--Hochberg FDR of $25\%$ is in parentheses. \label{TabFDR}}
\footnotesize
\begin{tabular}{cccc}
  \hline
  &\# declared:\\
Permutation-based FDR (\%) & VB & varbvs & VB $\cap$ varbvs \\
  \hline
  5 & 21 & 19 & 8 \\
  10 & 26 & 19 & 8 \\
  15 & 47 & 21 & 10 \\
  20 & 76 & 31 & 12 \\
  25 & 89 (48 univ.) & 47 (19 univ.) & 14 (13 univ.) \\
   \hline
\end{tabular}
\end{center}
\end{table}

\begin{figure
\centering
  \noindent\includegraphics[scale=0.61]{./VB_paper_Fig6_left-eps-converted-to.pdf}\quad
    \includegraphics[scale=0.91]{./VB_paper_Fig6_right-eps-converted-to.pdf}
    \caption{\small SNPs and pairwise associations declared by our approach for the Diogenes study. Left: Manhattan plot for SNP association and evidence of pleiotropy. Right: posterior probabilities of inclusion for SNP-metabolite associations declared at estimated FDR of $25\%$ and overlap with the associations declared by the varbvs method at same FDR level (pink crosses) and declared by univariate screening at Benjamini--Hochberg FDR of $25 \%$ (orange underscores).}\label{FigReal}
\end{figure}

In order to adjust for multiplicity, we specify the hyperparameters for $\omega$ according to the discussion of Section~\ref{SecMod} and choose the prior average number of active SNPs, $p^*$, by grid search within a $3$-fold cross-validation procedure that maximizes the variational lower bound. After hyperparameter selection, the algorithm converged in $83$ iterations, taking about $10$ hours on an Intel Xeon CPU at 2.60 GHz with 512 GB RAM. The posterior means $\text{E}_q(\omega_s \mid y)$ suggest the presence of several active SNPs, spread across the chromosomes (Figure \ref{FigReal}), but we use the marginal posterior inclusion probabilities, $\text{E}_q(\gamma_{st} \mid y)$, to declare pairwise associations and active SNPs.

We compare varbvs and our method on real data based on the number of associations declared by each method at specific false discovery rates estimated by permutations. We apply Efron's Bayesian interpretation of the false discovery rate \citep{efron2008microarrays} to posterior probabilities of inclusion, and use an empirical null distribution based on $B=400$ permutations to compute the estimate
\begin{equation}\label{EqFDR}\widehat{\textup{FDR}}(\tau) = \frac{\textup{median}_{b = 1, \ldots, B} \#\{\textup{PPI}_{st}^{(b)} > \tau\}}{\#\{\textup{PPI}_{st} > \tau\}}\,,\qquad 0 <\tau < 1\,,\end{equation}
for a grid of thresholds $\tau$; we then fit a cubic spline to the resulting false discovery rates to find thresholds for specific rates. The analysis suggests that our method is more powerful, with $89$ associations declared at an estimated FDR of $25 \%$, against $47$ for varbvs; the superiority of our method is further highlighted by Table $\ref{TabFDR}$. Figure~$\ref{FigReal}$ displays the associations declared by our method and the overlap with those declared by varbvs at estimated FDR of $25\%$; the associations detected also largely agree with those obtained with marginal screening at Benjamini--Hochberg FDR of $25\%$.

Database searches on the functional relevance of the detected associations give hints of promising biological functions related to metabolic activities for $12$ of the $25$ SNPs declared as active by our procedure. For instance, the most outstanding SNP Figure \ref{FigReal}, $rs4316911$, shows many associations with triglyceride levels, and turns out to be located less than $150 kb$ from the protein coding gene ITGA6 known to be linked to diabetic kidney disease \citep{iyengar2015genome}. The second most prominent pleiotropic SNP, $rs174535$, is declared by our approach to be associated with phospholipids, more precisely with $14$ different phosphatidylcholine levels, of which four are ether-linked/plasmalogen (PC-O). Interestingly, this latter SNP has been recently reported to be related to metabolite levels; among others, it was found to be associated with trans fatty acid levels and plasma phospholipid levels \citep{mozaffarian2015genetic}, in line with our findings. Moreover, it was found to be an eQTL for the fatty acid desaturase genes FADS1 and FADS2. The SNP $rs3903703$ too has been identified as associated with very long-chain fatty acid levels \citep{lemaitre2015genetic}. This seems to agree with our findings, in which $rs3903703$ exhibits associations with sphingomyelin, a type of lipid containing fatty acids of different chain lengths. The complete subset of SNPs with metabolism-related links found by our procedure is given in Appendix \ref{AppE}.  Additional details on this real data study, as well as on its replication using simulated data, are also provided there.

\section{Conclusions}\label{SecConcl}

We have described a scalable and efficient approach to joint variable selection from large numbers of candidate predictor and outcome variables. As it exploits the similarity across outcomes through a flexible hierarchical structure, our procedure outperforms the most popular predictor selection approaches in high-dimensional set-ups. The variational approximation on which our approach relies provides accurate posterior quantities, with reduced computational effort relative to MCMC procedures; in particular, it yields inferences comparable to those of the MCMC procedure HESS \citep{richardson2010bayesian}. Convergence control is automatic, whereas convergence assessment for MCMC algorithms can be difficult, especially in high dimensions. Our simulations also show that variable selection remains powerful when the predictors and outcomes are correlated, notwithstanding the independence assumptions underlying the mean-field factorization.

The key added-value of our approach is its applicability to molecular QTL datasets without the need for prior dimension reduction. In an application, our approach recovered several previously reported SNP-metabolite associations, and declared more associations than the single-outcome method ``varbvs'' \citep{carbonetto2012scalable} at prescribed false discovery rates, thus highlighting the benefits of jointly modelling the outcomes. To the best of our knowledge, no competing Bayesian approach for joint inference on two high-dimensional sets of variables can deal with the problem sizes typically encountered in molecular QTL analyses.

Bioinformatics is moving towards whole-genome analyses, for which several million genetic variants need to be considered, so it seems worthwhile to consider further speed-up strategies for our approach. One option is to use new optimization procedures. So-called natural gradient methods, which rely on the Riemannian structure of variational approximate distributions, seem particularly attractive, as they can be orders of magnitude faster than conventional gradient algorithms \citep{honkela2007natural}. Another possibility is a ``split-and-merge'' strategy, i.e., first partitioning the variable space and then inferring a global variational distribution on the aggregated dataset. \citet{tran2016parallel} designed a variant of this approach for sample space partitioning. At the recombination step, they proposed to ``merge'' the variational distributions by exploiting the independence assumptions of the mean-field formulation. Both strategies could lead to significant computational gains.

\section{Software}

The algorithm and the data-generation functions used in this paper are implemented in the publicly available R package \texttt{locus}.

\section*{Acknowledgements}
The authors are grateful to the Associate Editor and the two referees for their very helpful comments that improved the paper. We thank Jérôme Carayol, Loris Michel and Armand Valsesia for their valuable comments. We also thank James Holzwarth, Bruce O'Neel and Jaroslaw Szymczak for giving us access to computing resources.

\section*{Funding}

The Diogenes trial was supported by the European Commission (FP6-2005-513946) and Nestlé Institute of Health Sciences.

\bibliographystyle{plainnat}


\section*{Abstract}
We propose a simple extension of the Standard Model (SM) which has a viable dark matter (DM) candidate, as well as can explain the generation of tiny neutrino masses. The DM is an electroweak (EW) singlet scalar $S$, odd under an imposed exact $Z_2$ symmetry, interacting to SM through `Higgs-portal' coupling, while all other particles are even under $Z_2$. The model also has an EW isospin $3/2$ scalar, $\Delta$ and a pair of EW isospin vector, $\Sigma$ and $\bar{\Sigma}$, responsible for generating tiny neutrino mass via the effective dimension seven operator. Thanks to the additional interactions with $\Delta$, the scalar singlet DM $S$ survives a large region of parameter space by relic density constraints from WMAP/PLANCK and direct search bounds from updated LUX data. Constraints on the model from Large Hadron Collider (LHC) has also been discussed.
\end{abstract}

\maketitle

\newpage
\section{Introduction}
\label{sec:intro}

Evidence of Physics beyond the Standard Model (SM) have essentially come from two discoveries: the existence of non-zero neutrino masses and cosmological evidence for the existence of the Dark matter (DM). Although it is not established whether the DM is of astrophysical origin, or of particle physics, many models have been proposed as an extension of SM to accommodate a stable weakly interacting massive particle (WIMP) which can satisfy DM constraints. Similarly, several models of neutrino mass generation have been proposed to satisfy the observed neutrino masses and mixing. However, to bring them under one umbrella seems harder if not impossible \cite{Davoudiasl:2004be, Bhattacharya:2016lts}. The main incentive for addressing such a cause is that both neutrinos and DMs necessarily have to be coupled weakly to the SM.

In this note we try to address both the issues of neutrino masses and DM together with a minimal possible extension of the SM. The simplest way of accommodating DM is to assume the existence of a singlet scalar ($S$) which is coupled to the SM through Higgs portal coupling. The stability of such a DM is ensured with imposing a $Z_2$ symmetry under which the DM is odd, while the SM is even. The phenomenology of such a case have been discussed in different contexts for the simplicity and predictability \cite{singlet,singlet-FIMP,singlet-Z3,singlet-other}. However, non-observation of DM in direct search experiments is pushing this model under tight constraints. Excepting for Higgs resonance ($m_S=m_H/2$) region, the singlet scalar DM is essentially ruled out from the direct detection constraints \cite{Feng:2014vea, Bhattacharya:2016ysw} to a very large DM mass.

However, we demonstrate here that scalar singlet DM with Higgs-portal interaction can still survive without incorporating semi-annihilation \cite{semi-annihilation}, or multi-component feature \cite{Bhattacharya:2016ysw} if we can think of additional interactions of such DM to annihilate to some non-SM particles and produce right amount of relic density as observed by WMAP \cite{Hinshaw:2012aka} or PLANCK \cite{Ade:2013zuv}. In such a situation the effective DM-SM portal coupling required to generate correct density is {\it reduced} and so is the direct search cross-section; keeping the DM alive after strong bounds of LUX \cite{Akerib:2013tjd} and XENON100 \cite{Aprile:2012nq}. While the presence of additional annihilation channels for the DM to those beyond the SM particles has already been discussed in literature, for example, in MSSM, we highlight the fact here that such a feature helps alleviating the pressure from non-observation of the DM candidate in direct search experiments. This is first possible attempt to bring such a phenomena correlated to the cause of neutrino mass generation mechanism in one model framework to the best of our knowledge.

While we know that neutrino masses can effectively be generated though the existence of additional EW quadruplet scalar ($\Delta$) \cite{Babu:2009aq}, the singlet scalar $S$ can couple to it and can annihilate to them whenever the mass of the DM is larger than $\Delta$ ($m_S>m_{\Delta}$). This effectively then reduces annihilation cross-sections for $SS \to SM$ and keeps the model alive in a much larger region of parameter space from direct search experiment. We also point out that constraints coming from neutrino masses do not affect too much to the dark sector.

After the sad demise of the 750 GeV diphoton excess at the Large Hadron Collider (LHC) \cite{ATLAS:2016eeo,CMS:2016crm}, we point out the modified bound on the EW Quadruplet and the DM in absence of any signal beyond SM coming from the LHC experiments.

The paper is organized as follows. In Section \ref{sec:model}, we discuss the extensions beyond SM to accommodate DM and neutrino masses. In Sec. \ref{sec:dm}, we analyze the constraints on the model from relic density and direct detection searches of the dark matter. Numerical simulations for the LHC signatures of our model is discussed  in Section \ref{sec:collider}. Finally in section \ref{sec:conclusions}, we give our conclusions.

\section{Model and Formalism}
\label{sec:model}

Our model is based on the SM symmetry group $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$ supplemented by an unbroken discrete $Z_2$ symmetry. In the fermion sector, in addition to the SM fermions, we add two vector-like $SU(2)$ triplet leptons, $\Sigma $ and $\bar{\Sigma}$. In the scalar sector, in addition to the usual SM Higgs doublet, $H$, we introduce an isospin $3/2$ scalar, $\Delta$, and an EW singlet, $S$. The Singlet scalar, $S$ is odd under $Z_2$, while all other particles in the model are even under $Z_2$.

Let us mention that the extra particles $\Delta$, $\Sigma$ and $\bar{\Sigma}$ are introduced to generate tiny neutrino masses via the dimension seven operators \cite{Babu:2009aq,Ghosh:2016lnu}. We will discuss more on this after we introduce the scalar potential. This model can explain neutrino masses for reasonable choice of parameters, and can accommodate both normal and inverted hierarchy for the neutrino masses. The singlet scalar, $S$ along with the $Z_2$ symmetry introduced in this framework, provides with a viable candidate for DM, and is one of the major motivations of this work. The particle contents along with their quantum numbers are shown in the Table~\ref{Table1}.

\begin{table}[htb]
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
&$SU(3)_{C} \times SU(2)_{L} \times U(1)_{Y}$\\\hline
\small Matter &${\begin{pmatrix} u \\ d \end{pmatrix}}_L\sim(3,2,\frac{1}{3}), u_R\sim (3,1,\frac{4}{3}), d_R \sim (3,1,-\frac{2}{3})$ \\
&$ {\begin{pmatrix} \nu_e \\ e \end{pmatrix}}_L\sim (1,2,-1), e_R\sim (1,1,-2), \nu_R\sim (1,1,-2)$ \\&$ {\begin{pmatrix} \Sigma^{++} \\ \Sigma^{+} \\ \Sigma^{0} \end{pmatrix}}\sim (1,3,2)$, $ {\begin{pmatrix} \bar{\Sigma}^{0} \\ \bar{\Sigma}^{-} \\ \bar{\Sigma}^{--} \end{pmatrix}}\sim (1,3,-2)$ \\
\hline
\small Gauge & $G^\mu_{a,a=1-8}, A^\mu_{i, i=1-3}, B^\mu$ \\
\hline
\small Higgs & ${\begin{pmatrix} H^{+} \\ H^{0} \end{pmatrix}}\sim(1,2,1)$, ${\begin{pmatrix} \Delta^{+++} \\ \Delta^{++} \\ \Delta^{+} \\\Delta^{0}\end{pmatrix}}\sim(1,4,3)$, $ S\sim(1,1,0)$\\
\hline
\end{tabular}
\\
\end{center}
\caption{Matter, gauge and Higgs contents of the model.}
\label{Table1}
\end{table}

The most general renormalizable  scalar potential  consistent with {scalar spectrum of this} model is given by,

\begin{equation}
\begin{split}
V ( H, \Delta)= - \mu_H^{2}H^{\dagger}H + \mu_\Delta^{2}\Delta^{\dagger}\Delta +   \frac{\lambda_{1}}{2}(H^{\dagger}H)^{2} + \frac{\lambda_{2}}{2}(\Delta^{\dagger}\Delta)^{2} \\+  \lambda_{3}(H^{\dagger}H)(\Delta^{\dagger}\Delta) + \lambda_{4}(H^{\dagger}\tau_{a}H)(\Delta^{\dagger}T_{a}\Delta)  + \lbrace\lambda_{5}H^{3}\Delta^{\star} + h.c. \rbrace ,
\end{split}
\end{equation}
\begin{equation}
\begin{split}
V ( H,\Delta,S) = V ( H, \Delta)+\mu_S^{2}S^{2}  + \frac{\lambda_{7}}{2}S^{4} + \lambda_{8}(H^{\dagger}H)S^{2}  \\+ \lambda_{9}(\Delta^{\dagger}\Delta)S^{2},
\end{split}
\end{equation}
where  $\tau_{a}$ and $T_{a}$  are the generators of $SU(2)$ in the doublet {and} four-plet representations, respectively.

As was shown in \cite{Babu:2009aq}, even with positive ${\mu_{\Delta}}^2$, due to the $\lambda_5$ term in the potential, and the fields $\Sigma$ and $\bar\Sigma$, the neutral component of $\Delta$ acquires an induced VEV  at the tree level, $v_{\Delta} = - \lambda_5 v_H^3 / M_{\Delta}^2$, where $\langle H \rangle =v_H/\sqrt 2$ is the usual EW VEV. This gives rise to effective dimension seven operator $LLHH(H^{\dagger} H)/ M^{3}$, and generate tiny neutrino masses \cite{Babu:2009aq}. The additional singlet S that we have introduced gets no VEV ($v_s = 0$) to keep the $Z_2$ symmetry intact. Hence, we impose the condition $\mu_S^2>0$.

The mass of the neutral member of the quadruplet is given by  \cite{Babu:2009aq,Ghosh:2016lnu}
\begin{equation}
 M_{\Delta}^{2} = \mu_\Delta^{2} + \lambda_{3}v_{H}^{2} + \frac{3}{4}\lambda_{4}v_{H}^{2},
 \end{equation}

The mass splittings between the members of $\Delta$ are given by
\begin{equation}
M_{i}^{2} = M_{\Delta}^{2} - q_{i}\frac{\lambda_{4}}{4}v_{H}^{2},
\label{spectrum}
\end{equation}
where $q_{i}$ is the (non-negative) electric charge of the respective field. The mass
splittings are equally spaced and there are two possible mass orderings. For $\lambda_{4}$ positive,
we have the ordering $M_{\Delta^{+++}} < M_{\Delta^{++}} < M_{\Delta^{+}} < M_{\Delta^{0}}$ and for $\lambda_{4}$ negative, we have the ordering $M_{\Delta^{+++}} > M_{\Delta^{++}} > M_{\Delta^{+}} > M_{\Delta^{0}}$.

Later, it will turn out that that the mass splitting $\Delta M$ plays an important role in the decays, specially that of $\Delta^{++}$. So let us make some comments on the allowed vales of $\Delta M$.  As can be seen from the above equation, this mass splitting is arbitrary depending on the value of $\lambda_4$. however, as shown in \cite{Babu:2009aq}, there is an upper limit of $38$ GeV  on $\Delta M$ coming from the constraint on the $\rho$ parameter. There is also a theoretical lower limit of $1.4$ GeV on $\Delta M $\cite{Babu:2009aq}. In our analysis in this paper, we satisfy both limits.

The gauge singlet scalar $S$ which is odd under a $Z_2$ symmetry provides with a simplest DM candidate which has portal interactions with $\Delta$ (through $\lambda_9$) in addition to the SM Higgs (through $\lambda_8$).

The SM scalar singlet $S$ acquires mass through EW symmetry breaking as
\begin{equation}
M_S^2=\mu_S^2+\lambda_8 v_H^2/2
\end{equation}

Note here that $\mu_S^2>0$ implies following inequality

\begin{equation}
M_S^2>\lambda_8 v_H^2/2, ~{\rm or}~ \lambda_8< \frac{2 M_S^2}{v_H^2}
\end{equation}

Being singlet, the scalar-$S$ does not couple to the SM gauge bosons at tree level. The Yukawa interactions involving $S$ and the SM fermions are also forbidden by the EW as well as $Z_2$ symmetry.

We now address an important point regarding the stability of the DM particle $S$ due to the added $Z_2$ symmetry. It is well known that a discrete symmetry is vulnerable to Planck scale physics due to anomalies unless it is of gauge origin, and satisfies the discrete anomaly-free conditions \cite{Krauss, Banks-etc, Ibanez}. For example, in minimal supersymmetric Standard Model (MSSM), we introduce the discrete symmetry, matter parity, $P_M$ which is of gauge origin and satisfy the discrete anomaly conditions. This prevents the existence of dimension four baryon and lepton number violating operators in the superpotential and guarantees the stability of stability of proton (the proton decay is still possible through dimension five operator in MSSM, however the decay rate within the limit of proton life time). Following the work of Ibanez and Ross \cite{Ibanez}, if the discrete symmetry is $Z_N$, and $q_i$ are the charges of the fermions of the theory under $Z_N$ should satisfy the following condition:

\begin{eqnarray}
\sum_i q_i^3=mN+\eta n\frac{N^3}{8}
\label{eq:anomaly-cond}
\end{eqnarray}

where $\eta=0,1 ~ \rm{for}~N=\rm{even},~\rm{odd}$ respectively and $m$ and $n$ are integers. In our model, all the SM particles are even under $Z_2$, while the singlet scalar $S$ is odd. Thus, the SM fermions present in the model $(Q,u_L^c,d_L^c, L,e_L^c)$ have $Z_2$ quantum numbers $(1,1,1,1,1)$. We now show that our model satisfy the discrete anomaly-free conditions. For the cubic $Z_2^3$ anomaly, we find
\begin{eqnarray}
\sum_i q_i^3=\sum_i (1)^3=15=2m+n
\end{eqnarray}
using Eq.\ref{eq:anomaly-cond} with $m$ and $n$ being integers. This is easily satisfied with for example, with $m=7$ and $n=1$. For mixed gravitational anomaly, we get
\begin{eqnarray}
\sum_i q_i=\sum_i (1)=15=2p+q
\end{eqnarray}
where $p,q$ are integers and once again can easily be satisfied with $p=7, ~q=1$.
For mixed anomaly, for example, $Z_2-SU(2)-SU(2)$, we get
\begin{eqnarray}
\sum_{\rm{doublet}} q-i= 4 (1)=4=2r
\end{eqnarray}
where $r$ is an integer and is also easily satisfied. Thus our model satisfies all the anomaly free conditions for the imposed $Z_2$ symmetry to be of gauge origin leading to the stability of the DM.

\section{Dark Matter Analysis}
\label{sec:dm}
Scalar singlet extension of SM to accommodate DM through Higgs portal interaction is under tension as the allowed region of relic density space has been ruled out to a very large DM mass excepting for the Higgs resonance by non-observation in direct search experiment, especially the LUX data \cite{Akerib:2013tjd, Bhattacharya:2016ysw}. Possibilities to evade direct search bound for a DM component is an important question and present day DM research has to answer to that query. Here we present one such phenomena that successfully demonstrates a case for scalar singlet DM which can evade the direct search bound allowing the DM valid through a large region of parameter space.

The scalar singlet $S$ introduced here interacts with the scalar quadruplet $\Delta$ and can annihilate through $SS \to \Delta^0 \Delta^0, \Delta^+ \Delta^-, \Delta^{++} \Delta^{--}, \Delta^{+++} \Delta^{---}$ on top of annihilations to SM particles through Higgs portal interactions. Relic density of the DM in the present universe is obtained by the annihilation cross-section of the DM as

\begin{equation}
\Omega h^2 =\frac{0.1{\rm pb}}{\langle \sigma v \rangle}
\label{eq:Omega-sigmav}.
\end{equation}

The thermally averaged annihilation cross-section $\langle \sigma v \rangle $ can be written in terms of the usual cross-section as

\begin{equation}
\langle \sigma v \rangle_{ab \to cd} =\frac{T g_a g_b}{2 (2\pi)^4 n_a^{eq} n_b^{eq}} \int_{s_0}^{\infty} ds \frac{\lambda(s,m_a^2,m_b^2)}{\sqrt{s}} K_1(\frac{\sqrt{s}}{T}) \sigma_{ab \to cd}.
\label{eq:sigmav}
\end{equation}

Here $g_{a,b}$ corresponds to the degrees of freedom of annihilating particles, $K_1(x)$ is the first Bessel function, $s$ corresponds to the center-of-mass energy available for the process with $s_0=(m_a+m_b)^2$, $T$ is the temperature and $\lambda(a,b,c)=(a^2+b^2+c^2-2ab-2bc-2ca)^{\frac{1}{2}}$. $n_{a,b}^{eq}$ represents the equilibrium distributions of the DMs annihilating, which we assume to be non-relativistic and given by
\begin{equation}
n_X =  \int \frac{g_X d^3 p}{(2 \pi)^3 2E}  \tilde{f}_X,  \hspace{.4 cm}
n_X^{eq} =  \int \frac{g_X d^3 p}{(2 \pi)^3 2E}  \tilde{f}_X^{EQ},    \hspace{.4 cm}   \tilde{f}_X^{eq}  = \frac{1}{e^{E/T}\pm1}
\label{eq:neq}
\end{equation}

In the following analysis, there are two major contributions to the annihilation cross-section of $S$ as has already been mentioned and can be written as

\begin{equation}
\langle \sigma v \rangle=\langle \sigma v \rangle_{SS \to SM}+\langle \sigma v \rangle_{SS \to \Delta\D}
\end{equation}
The first part of the cross-section is well known and the corresponding annihilation cross-sections to fermions, gauge bosons and the the SM Higgs boson, $h$, and the singlet Higgs boson, $S$ can be written as \cite{Bhattacharya:2016ysw}

\begin{eqnarray}
(\sigma v_{rel})_{s s \rightarrow f \overline f}&=&\frac{1}{4\pi s \sqrt s} \frac{N_c\lambda_8^2 m_f^2}{(s-m_h^2)^2+m_h^2 \Gamma_h^2}(s-4m_f^2)^\frac{3}{2}\nonumber\\
{(\sigma v_{rel})}_{s s \rightarrow W^+ W^-}&=&\frac{\lambda_8^2}{8\pi} \frac{s}{(s-m_h^2)^2+m_h^2 \Gamma_h^2}(1+\frac{12m_W^4}{s^2}-\frac{4m_W^2}{s})(1-\frac{4m_W^2}{s})^\frac{1}{2}\nonumber\\
{(\sigma v_{rel})}_{s s \rightarrow Z Z}&=&\frac{\lambda_8^2}{16\pi} \frac{s}{(s-m_h^2)^2+m_h^2 \Gamma_h^2}(1+\frac{12m_Z^4}{s^2}-\frac{4m_Z^2}{s})(1-\frac{4m_Z^2}{s})^\frac{1}{2}\nonumber\\
{(\sigma v_{rel})}_{s s \rightarrow h h}&=&\frac{\lambda_8^2}{16\pi s}[1+\frac{3m_h^2}{(s-m_h^2)}-\frac{4\lambda_8 v^2}{(s-2m_h^2)}]^2 (1-\frac{4m_h^2}{s})^\frac{1}{2}.
\label{eq:annihilationSM}
\end{eqnarray}
Where $N_c=3$ is the color factor for quark and $N_c=1$ for leptons, $m_h=125$ GeV is the Higgs mass and $\Gamma_h$ is Higgs decay width at resonance ($\Gamma_{h\to SM}=4.07$ MeV).
The cross-section to $ss \to \Delta\D$ can be written as

\begin{equation}
{(\sigma v_{rel})}_{s s \rightarrow \Delta \Delta}=\frac{4 \sqrt{s-4m_{\Delta}^2}}{8\pi s \sqrt s}\lambda_9^2
\label{eq:annihilationDM}
\end{equation}

Where we assumed that the mass of all the charged and neutral components of $\Delta$ have the same mass. The upper limit on this mass splitting is 38 GeV \cite{Babu:2009aq}. The factor 4 is essentially indicating four different annihilations $SS \to \Delta^0 \Delta^0, \Delta^+ \Delta^-, \Delta^{++} \Delta^{--}, \Delta^{+++} \Delta^{---}$ which contribute equally in absence of a mass difference between them as we have assumed here.

We have inserted the model in micrOMEGAs \cite{Belanger:2014vza} and scan over the DM parameter space. The relevant parameter space of this model is spanned by the two mass parameters: the DM mass,  $M_s$, the common quadruplet mass, $M_{\Delta}$, and the two couplings, the Higgs portal coupling $\lambda_8$  and $\lambda_9$, the one connecting DM and quadruplet and are given by the set

\begin{equation}
\{M_s,~M_{\Delta}, ~\lambda_8, ~\lambda_9\}
\label{eq:parameters}
\end{equation}

In the following we vary the parameters given in Eq.~\ref{eq:parameters} and find the allowed region of correct relic abundance for the DM, $S$  satisfying WMAP~\cite{Hinshaw:2012aka} constraint ~\footnote{The range we use corresponds to the
WMAP results; the PLANCK constraints $0.112 \leq \Omega_{\rm DM} h^2 \leq 0.128$~\cite{Ade:2013zuv}, though more stringent, do not lead to
significant changes in the allowed regions of parameter space.}
\begin{eqnarray}
0.094 \leq \Omega_{\rm DM} h^2 \leq 0.128 \,.
\label{eq:wmap.region}
\end{eqnarray}

In Fig.~\ref{fig:Omega-M1}, we show the variation in relic density of $S$ with respect to the DM mass for fixed value of $M_{\Delta}=$ 400 GeV on the left and 700 GeV on the right. The couplings $\lambda_8=\lambda_9$ is varied in a long range and indicated through different color codes as follows: $\{0.01-0.1\} ({\rm Blue}),~ \{0.1-0.5\} ({\rm Green}),~ \{0.5-1\} ({\rm Purple}),~\{1-2\} ({\rm Orange})$. The correct density as in Eq.~\ref{eq:wmap.region} is indicated by the red horizontal lines.

 \begin{figure}[htb!]
$$
 \includegraphics[height=5.5cm]{Omega-M1-Del400.png}
 \includegraphics[height=5.5cm]{Omega-M1-Del700.png}
 $$
 \caption{$\Omega h ^2$ versus $M_s$ for different choices of $\lambda_8=\lambda_9:~ \{0.01-0.1\} ({\rm Blue}),~ \{0.1-0.5\} ({\rm Green}),~ \{0.5-1\} ({\rm Purple}),~\{1-2\} ({\rm Orange})$. Left: $M_{\Delta}=400~{\rm GeV}$, Right: $M_{\Delta}=700~{\rm GeV}$ are chosen for illustration. The correct density is indicated through the red horizontal lines.}
\label{fig:Omega-M1}
\end{figure}

\begin{figure}[htb!]
$$
 \includegraphics[height=5.5cm]{M1-L1-Del400.png}
 \includegraphics[height=5.5cm]{M1-L1-Del700.png}
 $$
\caption{Correct relic density region in $M_s-\lambda_8$ plane. We have assumed $\lambda_8=\lambda_9$ for simplicity here. In the left we have fixed $M_{\Delta}=400~{\rm GeV}$, on the right: $M_{\Delta}=700~{\rm GeV}$ chosen for illustration.}
\label{fig:M1-l}
\end{figure}

Next we turn to relic density allowed parameter space of the model. In the simplest scan as shown in Fig.~\ref{fig:Omega-M1}, the allowed region of parameter space can be depicted in $M_s-\lambda_8$ plane with the assumption of $\lambda_8=\lambda_9$ for constant $M_{\Delta}$. This is shown in Fig.~\ref{fig:M1-l} for fixed values of quadruplet mass as $M_{\Delta}=400~{\rm GeV}$ on the left and  $M_{\Delta}=700~{\rm GeV}$ on the right. We have chosen a value so that the chances of conflicting with direct LHC search bound is less. Both in Fig.~\ref{fig:Omega-M1} and Fig.~\ref{fig:M1-l}, we see a clear bump and a subsequent drop in relic density for $M_s>M_\Delta$ which is 400 GeV on the left and 700 GeV on the right. This is simply due to the additional cross-section of $SS \to \Delta\D$.

 \begin{figure}[htb!]
$$
 \includegraphics[height=5.5cm]{Omega-M1-del500.png}
 $$
 $$
 \includegraphics[height=5.5cm]{Omega-M1-del500-l8.png}
  \includegraphics[height=5.5cm]{Omega-M1-del500-l8l9.png}
 $$
 \caption{Top: $\Omega h ^2$ versus $M_s$ for different choices of $\lambda_8:~ \{0.001-0.01\} ({\rm Blue}),~ \{0.01-0.05\} ({\rm Green}),~ \{0.05-0.1\} ({\rm Purple}),~\{0.1-0.2\} ({\rm Orange})$ while $\lambda_9 :\{0.01-2.0\}$ varies. Bottom Left: $\Omega h ^2$ versus $M_s$ for $\lambda_8:~ \{0.001-0.1\} ({\rm Blue})$; Bottom right: $\Omega h ^2$ versus $M_s$ with $\lambda_8:~ \{0.001-0.1\} $ and different choices of $\lambda_9:~ \{0.001-0.1\} ({\rm Purple}),~ \{0.1-0.5\} ({\rm Green}),~ \{0.5-1\} ({\rm Orange}),~\{1-2\} ({\rm Brown})$. $M_{\Delta}=500~{\rm GeV}$ is chosen for illustration. The correct density is indicated through the red horizontal lines.}
\label{fig:Omega-M2}
\end{figure}

The situation gets even more interesting when we relax the condition imposed on couplings as $\lambda_8=\lambda_9$ and vary them freely as independent parameters. We show such an example in Fig.~\ref{fig:Omega-M2} where we choose the $M_{\Delta}=500~{\rm GeV}$ for illustration and vary $\lambda_8 :\{0.001-0.3\}$ and $\lambda_9 :\{0.01-2.0\}$ independently. For showing the annihilations $SS \to \Delta\D$, we highlight the region $M_S \gtrsim M_{\Delta} $. In the top panel of Fig.~\ref{fig:Omega-M2}, we show $\Omega h ^2$ versus $M_s$ for different choices of $\lambda_8:~ \{0.001-0.01\} ({\rm Blue}),~ \{0.01-0.05\} ({\rm Green}),~ \{0.05-0.1\} ({\rm Purple}),~\{0.1-0.2\} ({\rm Orange})$ while $\lambda_9 :\{0.01-2.0\}$ varies. Each range of chosen $\lambda_8$ actually shows a larger spread as has been pointed out in the bottom left panel, for example with $\lambda_8:~ \{0.001-0.01\}$. Evidently this is due to the large variety of $\lambda_9$ as chosen in the scan. Hence $\Omega h ^2$ versus $M_s$ with $\lambda_8:~ \{0.001-0.1\} $ and different choices of $\lambda_9:~ \{0.001-0.1\} ({\rm Purple}),~ \{0.1-0.5\} ({\rm Green}),~ \{0.5-1\} ({\rm Orange}),~\{1-2\} ({\rm Brown})$ are pointed out in the bottom right panel of Fig.~\ref{fig:Omega-M2}. The correct density is indicated through the red horizontal lines.

\begin{figure}[htb!]
$$
 \includegraphics[height=5.5cm]{M1-del500-l8.png}
 \includegraphics[height=5.5cm]{M1-del500-l8l9.png}
 $$
\caption{Left: Correct relic density region in $M_s-\lambda_8$ plane when $\lambda_8 :\{0.001-0.3\}$ and $\lambda_9 :\{0.01-2.0\}$ are varied independently. Right: We point out different ranges of $\lambda_9:\{0.1-0.2\} ({\rm Purple}),~\{0.2-0.3\} ({\rm Orange}),~\{0.3-0.4\} ({\rm Green}),~\{0.4-2.0\} ({\rm Brown}) $ in producing correct density in $M_s-\lambda_8$ plane. We choose the $M_{\Delta}=500~{\rm GeV}$ for illustration and focus in $M_S \gtrsim M_{\Delta} $.}
\label{fig:M1-l2}
\end{figure}

Now we turn to relic density allowed parameter space in $M_s-\lambda_8$ plane when $\lambda_8 :\{0.001-0.3\}$ and $\lambda_9 :\{0.01-2.0\}$ are varied independently and $M_{\Delta}=500~{\rm GeV}$ is chosen for illustration as shown in Fig.~\ref{fig:M1-l2}. As is clear from the relic density plot, with a large variation in the $SS -\Delta\D$ coupling, the allowed plane becomes much larger with even a very small value of $\lambda_8 \sim 10^{-3}$. Different contributions of $\lambda_9$ is shown in the right panel of Fig.~\ref{fig:M1-l2}. The white blanks in between is due to the coarseness of the scanning done in the numerical analysis and do not contain any physics. Similar feature will be observed with any other possible choice of $M_{\Delta}$.

\begin{figure}[htb!]
$$
 \includegraphics[height=5.8cm]{DDS.png}
  $$
 \caption{Feynman diagram for direct detection of the DM.}
 \label{fig:direct}
\end{figure}

Next we turn to direct search constraint for the model. Here the direct detection occurs through Higgs mediation as usual to Higgs portal DM as shown in Fig.~\ref{fig:direct} and the spin independent DM-nucleon cross-section reads:
\begin{eqnarray}
\sigma_{SI}^n=\frac{\alpha_n^2 \mu_{n}^2}{4\pi m_{S}^2}
\end{eqnarray}
 where the reduced mass $\mu_n=\frac{m_n m_{S}}{m_n+m_{S}}$ is expressed in terms of  nucleon mass $m_n$, and the nucleon form factors are given by:
  \begin{eqnarray}
  \alpha_n &=& m_n\sum_{u,d,s} f_{T_q}^{(n)} \frac{\alpha_q}{m_q} + \frac{2}{27} f_{T_g}^{(n)} \sum_{q=c,t,b}\frac{\alpha_q}{m_q} \nonumber \\
  &=& m_n\sum_{u,d,s} f_{T_q}^{(n)} \frac{\alpha_q}{m_q} + \frac{2}{27}(1-\sum_{u,d,s} f_{T_q}^{(n)})\sum_{q=c,t,b}\frac{\alpha_q}{m_q} \nonumber \\
  &=& \frac{m_n \lambda_i}{m_h^2}[(f_{T_u}^{(n)}+f_{T_d}^{(n)}+f_{T_s}^{(n)})+\frac{2}{9}(f_{T_u}^{(n)}+f_{T_d}^{(n)}+f_{T_s}^{(n)})]
  \end{eqnarray}

Here $n$ stands for both proton and neutron. For proton we choose : $f_{T_u}^p=0.0153$ , $f_{T_d}^p=0.0191$ , $f_{T_s}^p=0.0447$ as the default values in micrOMEGAs.

 It is of great importance to see what parameter space of the relic density allowed DM region of the scalar DM model is allowed by the spin independent direct search constraints by XENON100 \cite{Aprile:2012nq}, and LUX \cite{Akerib:2013tjd} data. This is what is presented in Fig.~\ref{fig:dd1} first with the simplified case of $\lambda_8=\lambda_9$.

 \begin{figure}[htb!]
$$
 \includegraphics[height=5.5cm]{dd-Del400.png}
 \includegraphics[height=5.5cm]{dd-Del700.png}
 $$
 $$
 \includegraphics[height=5.5cm]{DDscalar.png}
 $$
\caption{Spin-independent DM-nucleon effective cross-section in terms of DM mass for the case of a fixed $M_{\Delta}=400$ GeV (top left) and 700 GeV (top right) for points with correct relic density with $\lambda_8=\lambda_9$. XENON100, LUX updated constraints as well as XENON 1T prediction is shown in the figure. Bottom: The case for single component and two-component scalar singlet dark matter with Higgs portal interaction is shown.}
\label{fig:dd1}
\end{figure}

The most important conclusion from the analysis comes out from the direct search results. From Fig.~\ref{fig:dd1} we clearly see that the relic density allowed points of this model is allowed by the direct search constraint from updated LUX data when $M_s>M_{\Delta}$,  with $M_s$ starting from 400 GeV on the left panel and 700 GeV on the right panel Fig.~\ref{fig:dd1}. This is simply because with additional annihilation cross-section the required coupling to SM drops as shown in Fig~\ref{fig:M1-l} and hence the direct search cross-section also drops to keep the model alive. The smaller the $M_{\Delta}$, the larger is the allowed region. However, we have to abide by the constraints on the quadruplet mass from collider search bounds and we cant keep it as low as we wish.

\begin{figure}[htb!]
$$
 \includegraphics[height=6.5cm]{dd-del500.png}
 $$
 \caption{Spin-independent DM-nucleon effective cross-section in terms of DM mass for the case of a fixed $M_{\Delta}=500$ GeV when $\lambda_8 :\{0.001-0.3\}$ and $\lambda_9 :\{0.01-2.0\}$ are varied independently for points satisfying relic density. XENON100, LUX updated constraints as well as XENON 1T prediction is shown in the figure.}
\label{fig:dd2}
\end{figure}

The situation is even better for the case when $\lambda_8$ and $\lambda_9$ are varied as uncorrelated parameters as shown in Fig.~\ref{fig:dd2} as illustration for the case with $M_{\Delta}=500$ GeV. The green points which satisfy relic density for the particular case as indicated in Fig.~\ref{fig:M1-l2}, has a reduced SM Higgs-portal coupling owing to the annihilations to $\Delta\D$ that do not eventually contribute to direct search. Hence for $M_s>M_{\Delta}$, as the figure indicate, the direct search may get delayed till XENON1T. This is what makes the model alive for future direct search discoveries.

In summary the model is allowed by DM constraints in a large region of parameter space with relatively higher values of DM mass, larger than the quadruplet mass.

\section{Constraints from LHC}
\label{sec:collider}
This model provides an interesting avenue to test the neutrino mass generation mechanism at the LHC. The presence of the isospin $3/2$ scalar multiplet $\Delta$ and a pair of vector-like fermions $\Sigma$ can give rise to rich phenomenology at the LHC. The detailed study of collider signatures has already been studied in early literature \cite{Bambhaniya:2013yca}. From the dark matter perspective, which is main motivation of this work, the most important aspect is the limit on the mass of $\Delta$ from the latest LHC experimental results.

At the LHC, $\Delta^{\pm\pm\pm} \Delta^{\mp\mp\mp}$, $\Delta^{\pm\pm} \Delta^{\mp\mp}$ and $\Delta^{\pm} \Delta^{\mp}$ are pair produced via the s-channel $\gamma$ and Z exchanges, while $\Delta^{\pm\pm\pm} \Delta^{\mp\mp}$, $\Delta^{\pm\pm} \Delta^{\mp}$ and $\Delta^{\pm} \Delta^{0}$ are pair produced via s-channel $W^{\pm}$ exchange as shown in the Fig.~\ref{fig:12}.

 \begin{figure}[htb!]
$$
 \includegraphics[height=3cm]{1.png}
 \includegraphics[height=3cm]{2.png}
 $$
 \caption{Feynman diagrams for production of doubly and triply charged scalars at LHC.}
\label{fig:12}
\end{figure}

\begin{figure}[htb!]
$$
 \includegraphics[height=8cm]{3a.png}
 $$
 \caption{Feynman diagrams for decay of doubly- charged scalars. }
\label{fig:34}
\end{figure}

The production cross section of the doubly- and triply- charged scalars at the LHC for center-of-mass energy $\sqrt{s}$ = 13 TeV is shown in Fig.~\ref{fig:5} as a function of the mass parameter. For simplicity, we have taken the masses of the quadruplets to be the same \footnote{Constraints from the $\rho$ parameter dictates the splitting to be $< 38 ~{\rm GeV}$, and can be even smaller depending on the values of $\lambda_4$}.

The model has been implemented in CalcHEP package \cite{Belyaev:2012qa}. For the production cross-section, we used parton distribution function CT10 \cite{Dulat:2013kqa} from LHAPDF library \cite{Whalley:2005nh} with the renormalization and factorization scales being chosen to be the invariant mass of the constituent sub-process.

\begin{figure}[htb!]
$$
 \includegraphics[height=7cm]{5n.png}
 $$
\caption{Production cross-sections for triply and doubly charged scalars at the LHC ($\sqrt{S}$ = 13 TeV) as a function of  $M_{\Delta}$.}
\label{fig:5}
\end{figure}

\begin{figure}[htb]
$$
 \includegraphics[height=4.5cm]{n1.png}
 \includegraphics[height=4.5cm]{n2.png}
 $$
 $$
 \includegraphics[height=4.5cm]{n3.png}
 \includegraphics[height=4.5cm]{8n.png}
 $$
 \caption{ Branching ratio (BR) for different decay modes of $\Delta^{\pm\pm}$ as a function of $M_{\Delta}$ for three different values of $v_{\Delta}$ : $10^{-7}$ GeV (Top Left), $10^{-6}$ GeV (Top Right) and $10^{-5}$ GeV (Bottom Left). Bottom Right: Branching ratio (BR) for different decay modes of $\Delta^{\pm\pm}$ as a function of $v_{\Delta}$ considering neutral scalar mass 500 GeV and mass difference between two successive scalars $\Delta M=$ 1.6 GeV.}
\label{fig:66}
\end{figure}

As mentioned before, our main objective in this collider analysis is to check what constraint the LHC experiments put on the allowed mass of $\Delta$ which satisfy the dark matter constraints. In our model, $\Delta^{\pm\pm}$, which is a doubly charged scalar,  has three major decay modes : $\Delta^{\pm\pm} \rightarrow l^{\pm}l^{\pm}, W^{\pm}W^{\pm}, \Delta^{\pm}\pi^{\pm}$. All of these modes can give rise same sign dilepton in the final state for $\Delta^{\pm\pm}$ production at LHC. \footnote {A fourth decay mode, to $\Delta^{\pm}W^{{\star}^{\pm}}$, or $\Delta^{\pm\pm\pm} W^{{\star}^{\pm}}$ depending on whether the $\Delta^{\pm\pm\pm}$ is the heaviest or lightest in the quadruplet is possible, but its width is much smaller compared to the other three.}

The ATLAS Collaboration has recently searched \cite{ATLAS:2016pbt} for the doubly-charged Higgs boson in the same-sign di-electrons invariant mass spectrum with luminosity 13.9 $fb^{-1}$ at $\sqrt{S}$ = 13 TeV. Their observed lower mass limit for this doubly charged Higgs, assuming a 100$\%$ branching ratio to di-electrons, is 570 GeV, while the observed lower mass limit, assuming a 50$\%$ branching ratio to di-electrons, is 530 GeV. Our model also must comply with non-observation of excess in same sign dilepton search.

Hence, it is important to parametrise the branching fraction of the doubly charged scalar to be able to compute the same-sign dilepton final state signature arising from the pair production of $\Delta^{\pm\pm}$ at LHC. We note that the decay width for the decay mode $\Delta^{\pm\pm} \rightarrow l^{\pm}l^{\pm}$ is proportional to 1/$v_{\Delta}^{2}$, the decay width to $W^{\pm}W^{\pm}$ final state is proportional to $v_{\Delta}^{2}$, while the one to $\Delta^{\pm}\pi^{\pm}$ is independent of $v_{\Delta}$, and proportional to $(\Delta M)^{3}$. In  Fig.~\ref{fig:66}, we plot the relative branching ratios of $\Delta^{\pm\pm}$ as a function of $M_{\Delta}$ and $v_{\Delta}$. As expected, for a very small
$v_{\Delta}$, the decay to $l^{\pm}l^{\pm}$ dominate, whereas for higher values of $v_{\Delta}$, the mode $\Delta^{\pm}\pi^{\pm}$ dominate.

We use the dedicated search by the ATLAS Collaboration \cite{ATLAS:2016pbt} with luminosity 13.9 $fb^{-1}$ at $\sqrt{S}$ = 13 TeV, for the doubly charged scalar  di-electron resonance. Here, to put the most conservative bound on mass $M_{\Delta}$, we assume small value of $v_{\Delta}$ ($\approx 10^{-7}$ GeV) so that $\Delta^{\pm\pm}$ mostly decays into $l^{\pm}l^{\pm}$ with branching ratio nearly equal to one. From our  calculated productions cross sections for the doubly charged scalar, and the branching ratios as shown in Fig~\ref{fig:67}, we obtain a most conservative lower limit of $324$ GeV for the mass of $\Delta^{\pm\pm}$, assuming a 100$\%$ (33.33$\%$) branching ratio to $l^{\pm}l^{\pm}$($e^{\pm}e^{\pm}$). Thus a large mass range of  $\Delta$  satisfying the dark matter constraints are allowed by the latest LHC experimental search. We also note here, that with the assumption of $M_S> M_\Delta$, the additional charged scalars do not decay to DM and hence the branching fractions do not get altered.

\begin{figure}[htb!]
$$
 \includegraphics[height=6.8cm]{9.png}
 $$
 \caption{The model predicted, observed and expected 95$\%$ C.L. upper limits of the cross-section $\sigma$ (pp $\rightarrow \Delta^{\pm\pm} \rightarrow e^{\pm}e^{\pm})$ as a function of $\Delta^{\pm\pm}$ mass at $\sqrt{S}$ = 13 TeV. The limit is derived under the assumption that Br($\Delta^{\pm\pm} \rightarrow e^{\pm}e^{\pm}$) = 33.33$\%$. The red line under the light blue shaded region indicates the allowed mass region.}
\label{fig:67}
\end{figure}

\section{Summary and Conclusions}
\label{sec:conclusions}

In this work, we have analysed the possibility of having a viable DM as a EW scalar singlet $S$, that connects to SM via Higgs portal coupling in a framework that also cater to neutrino mass generation through the presence of a EW quadruplet, $\Delta$, and two EW triplet leptons $\Sigma$ and $\bar{\Sigma}$. An unbroken $Z_2$, under which $S$ is odd, makes the DM stable, while the rest are even. The DM phenomenology is shown to be crucially dictated by the interaction of $S$ with $\Delta$ on top of the Higgs portal coupling that it posses. This is because of the additional annihilation channels for the DM, $S$ to $\Delta$'s on top of SM particles whenever $M_S > M_\Delta$. In order to satisfy relic density constraints, the coupling of $S$ to the SM particles (Higgs-portal) is reduced compared to the case when annihilations only to SM has to contribute to the whole DM relic density. As a result, direct detection cross sections obtained for $M_S > M_\Delta$ satisfy the current LUX 2015 data for a large range of the DM parameter space. There is some possibility that the level of direct detection cross section predicted by the model might be within the reach of the XENON1T experiment, while it may go beyond as well.

This phenomena can be generalised by postulating whenever there is an additional channel for the DM to annihilate, but that doesn't give a contribution to direct search experiments for that DM, one can satisfy the stringent direct search bounds coming from non-observations of DM in terrestrial experiments. The particular example we have presented, has an additional motivation of explaining neutrino masses and the additional DM interaction naturally fits into the model framework with minimal assumptions.

The validity of the DM parameter space crucially depends on the mass required for the quadruplet $\Delta$. Therefore we attempted to evaluate the current bound on the mass of $\Delta$ from the latest LHC data. This methodology depends on the non-observation of a doubly charge scalar resonance decaying to same sign di-electron. We did not attempt to evaluate the bounds on the DM from LHC search as the one from collider is weaker than the one obtained from direct search.

To summarise the model is well motivated and has a rich phenomenology. It can be distinguished from the usual Higgs portal DM models from direct search prospect and with reduced SM coupling can easily survive LHC bounds. The model predicts signatures in leptonic final states at LHC through the productions of charged scalars. If signals in dilepton, four lepton and signals with higher leptonic multiplicity is seen with higher luminosity data, will indicate towards the existence of such framework.

\begin{acknowledgments}

The work of SN was supported in part by the United State Department of Energy Grant Number DE-SC 0016013. The work of SB is supported by DST-INSPIRE grant no. PHY/P/SUB/01 at IIT Guwahati.

\end{acknowledgments}

\section{References}