<!-- Generated by pkgdown: do not edit by hand -->
<!DOCTYPE html>
<html>
  <head>
  <meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">

<title>Fit variable selection model using variational approximation methods. — varbvs • varbvs</title>

<!-- jquery -->
<script src="https://code.jquery.com/jquery-3.1.0.min.js" integrity="sha384-nrOSfDHtoPMzJHjVTdCopGqIqeYETSXhZDFyniQ8ZHcVy08QesyHcnOUpMpqnmWq" crossorigin="anonymous"></script>
<!-- Bootstrap -->

<link href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-BVYiiSIFeK1dGmJRAkycuHAHRg32OmUcww7on3RYdg4Va+PmSTsz/K68vbdEjh4u" crossorigin="anonymous">
<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js" integrity="sha384-Tc5IQib027qvyjSMfHjOMaLkfuWVxZxUPnCJA7l2mCWNIpG9mGCD8wGNIcPD7Txa" crossorigin="anonymous"></script>

<!-- Font Awesome icons -->
<link href="https://maxcdn.bootstrapcdn.com/font-awesome/4.6.3/css/font-awesome.min.css" rel="stylesheet" integrity="sha384-T8Gy5hrqNKT+hzMclPo118YTQO6cYprQmhrYwIiQ/3axmI1hQomh7Ud2hPOy8SP1" crossorigin="anonymous">


<!-- pkgdown -->
<link href="../pkgdown.css" rel="stylesheet">
<script src="../jquery.sticky-kit.min.js"></script>
<script src="../pkgdown.js"></script>

<!-- mathjax -->
<script src='https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>

<!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->


  </head>

  <body>
    <div class="container template-reference-topic">
      <header>
      <div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <a class="navbar-brand" href="../index.html">varbvs</a>
    </div>
    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
        <li>
  <a href="../index.html">Home</a>
</li>
<li>
  <a href="../articles/index.html">Vignettes</a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
      </ul>
      
      <ul class="nav navbar-nav navbar-right">
        <li>
  <a href="http://github.com/pcarbo/varbvs">
    <span class="fa fa-github fa-lg"></span>
     
  </a>
</li>
      </ul>
    </div><!--/.nav-collapse -->
  </div><!--/.container -->
</div><!--/.navbar -->

      
      </header>

      <div class="row">
  <div class="col-md-9 contents">
    <div class="page-header">
    <h1>Fit variable selection model using variational approximation methods.</h1>
    </div>

    
    <p>Compute fully-factorized variational approximation for
  Bayesian variable selection in linear (family = gaussian) or logistic
  regression (family = binomial). More precisely, find the "best"
  fully-factorized approximation to the posterior distribution of the
  coefficients, with spike-and-slab priors on the coefficients. By
  "best", we mean the approximating distribution that locally minimizes
  the Kullback-Leibler divergence between the approximating distribution
  and the exact posterior.</p>
    

    <pre><span class='fu'>varbvs</span> (<span class='no'>X</span>, <span class='no'>Z</span>, <span class='no'>y</span>, <span class='kw'>family</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='st'>"gaussian"</span>,<span class='st'>"binomial"</span>), <span class='no'>sigma</span>, <span class='no'>sa</span>,
          <span class='no'>logodds</span>, <span class='no'>alpha</span>, <span class='no'>mu</span>, <span class='no'>eta</span>, <span class='no'>update.sigma</span>, <span class='no'>update.sa</span>,
          <span class='no'>optimize.eta</span>, <span class='no'>initialize.params</span>, <span class='kw'>nr</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>sa0</span> <span class='kw'>=</span> <span class='fl'>1</span>, <span class='kw'>n0</span> <span class='kw'>=</span> <span class='fl'>10</span>,
          <span class='kw'>tol</span> <span class='kw'>=</span> <span class='fl'>1e-4</span>, <span class='kw'>maxiter</span> <span class='kw'>=</span> <span class='fl'>1e4</span>, <span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)</pre>
    
    <h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a> Arguments</h2>
    <table class="ref-arguments">
    <colgroup><col class="name" /><col class="desc" /></colgroup>
    <tr>
      <th>X</th>
      <td><p>n x p input matrix, where n is the number of samples,
           and p is the number of variables. X cannot be sparse,
           and cannot have any missing values (NA).</p></td>
    </tr>
    <tr>
      <th>Z</th>
      <td><p>n x m covariate data matrix, where m is the number of
           covariates. Do not supply an intercept as a covariate (i.e.,
           a column of ones), because an intercept is automatically
           included in the regression model. For no covariates, set
           <code>Z = NULL</code>. The covariates are assigned an improper,
           uniform prior. Although improper priors are generally not
           advisable because they can result in improper posteriors and
           Bayes factors, this choice allows us to easily integrate out
           these covariates.</p></td>
    </tr>
    <tr>
      <th>y</th>
      <td><p>Vector of length n containing observations of binary
          (<code>family = "binomial"</code>) or continuous (<code>family =
          "gaussian"</code>) outcome. For a binary outcome, all entries
          of y must be 0 or 1.</p></td>
    </tr>
    <tr>
      <th>family</th>
      <td><p>"gaussian" for linear regression model, or "binomial"
                for logistic regression model.</p></td>
    </tr>
    <tr>
      <th>sigma</th>
      <td><p>Candidate settings for the residual variance
  parameter. Must be of the same length as inputs sa and logodds (or
  have length equal to the number of columns of logodds). Only used for
  linear regression, and will generate an error if <code>family =
  "binomial"</code>. If missing, residual variance parameter is automatically
  fitted to data by computing approximate maximum-likelihood (ML)
  estimate.</p></td>
    </tr>
    <tr>
      <th>sa</th>
      <td><p>Hyperparameter sa is the prior variance of regression
    coefficients for variables that are included in the model. This
    prior variance is always scaled by sigma (for logistic regression,
    we take sigma = 1). Scaling the variance of the coefficients in this
    way is necessary to ensure that this prior is invariant to
    measurement scale (e.g., switching from grams to kilograms). This
    input specifies the candidate settings for sa, of the same length as
    inputs sigma and logodds (or have length equal to the number of
    columns of logodds). If missing, prior variance is automatically
    fitted to data by compute approximate maximum (ML) estimates, or
    maximum a posteriori estimates when <code>n0 &gt; 0</code> and <code>sa0 &gt;
    0</code>.</p></td>
    </tr>
    <tr>
      <th>logodds</th>
      <td><p>Hyperparameter logodds is the prior log-odds that a
    variable is included in the regression model; it is defined as
    \(logodds = log10(q/(1-q)),\) where q is the prior probability
    that a variable is included in the regression model. Note that we
    use the base-10 logarithm instead of the natural logarithm because
    it is usually more natural to specify prior log-odds settings in
    this way. The prior log-odds may also be specified separately for
    each variable, which is useful is there is prior information about
    which variables are most relevant to the outcome. This is
    accomplished by setting logodds to a p x ns matrix, where p is the
    number of variables, and ns is the number of hyperparameter
    settings. Note it is not possible to fit the logodds parameter; if
    logodds input is not provided as input, then it is set to the
    default value when sa and sigma are missing, and otherwise an error
    is generated.</p></td>
    </tr>
    <tr>
      <th>alpha</th>
      <td><p>Good initial estimate for the variational parameter alpha
    for each hyperparameter setting. Either missing, or a p x ns matrix,
    where p is the number of variables, and ns is the number of
    hyperparameter settings.</p></td>
    </tr>
    <tr>
      <th>mu</th>
      <td><p>Good initial estimate for the variational parameter mu for
    each hyperparameter setting. Either missing, or a p x ns matrix,
    where p is the number of variables, and ns is the number of
    hyperparameter settings.</p></td>
    </tr>
    <tr>
      <th>eta</th>
      <td><p>Good initial estimate of the additional free parameters
    specifying the variational approximation to the logistic regression
    factors. Either missing, or an n x ns matrix, where n is the number
    of samples, and ns is the number of hyperparameter settings.</p></td>
    </tr>
    <tr>
      <th>update.sigma</th>
      <td><p>Setting this to TRUE ensures that sigma is always
    fitted to data, in which case input vector sigma is used to provide
    initial estimates.</p></td>
    </tr>
    <tr>
      <th>update.sa</th>
      <td><p>Setting this to TRUE ensures that sa is always fitted
    to data, in which case input vector sa is used to provide initial
    estimates.</p></td>
    </tr>
    <tr>
      <th>optimize.eta</th>
      <td><p>When <code>optimize.eta = TRUE</code>, eta is fitted to
    the data during the inner loop coordinate ascent updates, even when
    good estimates of eta are provided as input.</p></td>
    </tr>
    <tr>
      <th>initialize.params</th>
      <td><p>If FALSE, the initialization stage of the
    variational inference algorithm (see below) will be skipped, which
    saves computation time for large data sets.</p></td>
    </tr>
    <tr>
      <th>nr</th>
      <td><p>Number of samples of "model PVE" to draw from posterior.</p></td>
    </tr>
    <tr>
      <th>sa0</th>
      <td><p>Scale parameter for a scaled inverse chi-square prior on
    hyperparameter sa. Must be &gt;= 0.</p></td>
    </tr>
    <tr>
      <th>n0</th>
      <td><p>Number of degrees of freedom for a scaled inverse chi-square
    prior on hyperparameter sa. Must be &gt;= 0. Large settings of
    <code>n0</code> provide greater stability of the parameter estimates for
    cases when the model is "sparse"; that is, when few variables are
    included in the model.</p></td>
    </tr>
    <tr>
      <th>tol</th>
      <td><p>Convergence tolerance for inner loop.</p></td>
    </tr>
    <tr>
      <th>maxiter</th>
      <td><p>Maximum number of inner loop iterations.</p></td>
    </tr>
    <tr>
      <th>verbose</th>
      <td><p>If <code>verbose = TRUE</code>, print progress of algorithm
    to console.</p></td>
    </tr>
    </table>
    
    <h2 class="hasAnchor" id="regression-models"><a class="anchor" href="#regression-models"></a>Regression models</h2>

    <p>Two types of outcomes (y) are modeled: (1) a continuous outcome,
    also a "quantitative trait" in the genetics literature; or (2) a
    binary outcome with possible values 0 and 1. For the former, set
    <code>family = "gaussian"</code>, in which case, the outcome is
    i.i.d. normal with mean \(u0 + Z*u + X*b\) and variance sigma, in
    which u and b are vectors of regresion coefficients, and u0 is the
    intercept. In the second case, we use logistic regression to model
    the outcome, in which the probability that y = 1 is equal to
    \(sigmoid(u0 + Z*u + X*b).\) See <code>help(sigmoid)</code> for a
    description of the sigmoid function. Note that the regression always
    includes an intercept term (u0).</p>
    
    <h2 class="hasAnchor" id="co-ordinate-ascent-optimization-procedure"><a class="anchor" href="#co-ordinate-ascent-optimization-procedure"></a>Co-ordinate ascent optimization procedure</h2>

    <p>For both regression models, the fitting procedure consists of an inner
    loop and an outer loop. The outer loop iterates over each of the
    hyperparameter settings (sa, sigma and logodds). Given a setting of
    the hyperparameters, the inner loop cycles through coordinate ascent
    updates to tighten the lower bound on the marginal likelihood,
    \(Pr(Y | X, sigma, sa, logodds)\). The inner loop coordinate
    ascent updates terminate when either (1) the maximum number of inner
    loop iterations is reached, as specified by <code>maxiter</code>, or (2)
    the maximum difference between the estimated posterior inclusion
    probabilities is less than <code>tol</code>.</p>
<p>To provide a more accurate variational approximation of the posterior
    distribution, by default the fitting procedure has two stages. In the
    first stage, the entire fitting procedure is run to completion, and the
    variational parameters (alpha, mu, s, eta) corresponding to the maximum
    lower bound are then used to initialize the coordinate ascent updates
    in a second stage. Although this has the effect of doubling the
    computation time (in the worst case), the final posterior estimates
    tend to be more accurate with this two-stage fitting procedure.</p>
    
    <h2 class="hasAnchor" id="variational-approximation"><a class="anchor" href="#variational-approximation"></a>Variational approximation</h2>

    <p>Outputs alpha, mu and s specify the approximate posterior distribution
    of the regression coefficients. Each of these outputs is a p x ns
    matrix. For the ith hyperparameter setting, alpha[,i] is the
    variational estimate of the posterior inclusion probability (PIP)
    for each variable; mu[,i] is the variational estimate of the
    posterior mean coefficient given that it is included in the model;
    and s[,i] is the estimated posterior variance of the coefficient
    given that it is included in the model. These are also the
    quantities that are optimized as part of the inner loop coordinate
    ascent updates. An additional free parameter, eta, is needed for
    fast computation with the logistic regression model <code>(family =
    "binomial")</code>. The fitted value of eta is returned as an n x ns
    matrix.</p>
<p>The variational estimates should be interpreted carefully,
    especially when variables are strongly correlated. For example,
    consider the simple scenario in which 2 candidate variables are
    closely correlated, and at least one of them explains the outcome
    with probability close to 1. Under the correct posterior
    distribution, we would expect that each variable is included with
    probability ~0.5. However, the variational approximation, due to the
    conditional independence assumption, will typically get this wrong,
    and concentrate most of the posterior weight on one variable (the
    actual variable that is chosen will depend on the starting
    conditions of the optimization). Although the individual PIPs are
    incorrect, a statistic summarizing the variable selection for both
    correlated variables (e.g., the total number of variables included
    in the model) should be reasonably accurate.</p>
<p>More generally, if variables can be reasonably grouped together
    based on their correlations, we recommend interpreting the variable
    selection results at a group level. For example, in genome-wide
    association studies (see the vignettes) ,a SNP with a high PIP
    indicates that this SNP is probably associated with the trait, and
    one or more nearby SNPs within a chromosomal region, or ``locus,''
    may be associated as well. Therefore, we interpreted the GWAS
    variable selection results at the level of loci, rather than at the
    level of individual SNPs.
    Also note that special care is required for interpreting the results
    of the variational approximation with the logistic regression
    model. In particular, interpretation of the individual estimates of
    the regression coefficients (e.g., the posterior mean estimates
    <code>fit$mu</code>) is not straightforward due to the additional
    approximation introduced on the individual nonlinear factors in the
    likelihood. As a general guideline, only the relative magnitudes of
    the coefficients are meaningful.</p>
    
    <h2 class="hasAnchor" id="averaging-over-hyperparameter-settings"><a class="anchor" href="#averaging-over-hyperparameter-settings"></a>Averaging over hyperparameter settings</h2>

    <p>In many settings, it is good practice to account for uncertainty in
    the hyperparameters when reporting final posterior quantities. For
    example, hyperparameter sa is often estimated with a high degree of
    uncertainty when only a few variables are included in the model.
    Provided that (1) the hyperparameter settings sigma, sa and logodds
    adequately represent the space of possible hyperparameter settings
    with high posterior mass, (2) the hyperparameter settings are drawn
    from the same distribution as the prior, and (3) the
    fully-factorized variational approximation closely approximates the
    true posterior distribution, then final posterior quantities can be
    calculated by using logw as (unnormalized) log-marginal
    probabilities.</p>
<p>Even when conditions (1), (2) and/or (3) are not satisfied, this can
    approach can still often yield reasonable estimates of averaged
    posterior quantities. The examples below demonstrate how final
    posterior quantities are reported by function
    <code><a href='summary.varbvs.html'>summary.varbvs</a></code> (see <code>help(summary.varbvs)</code> for
    more details). To account for discrepancies between the prior on
    (sigma,sa,logodds) and the sampling density used to draw candidate
    settings of the hyperparameters, adjust the log-probabilities by
    setting <code>fit$logw &lt;- fit$logw + logp/logq</code>, where logp is the
    log-density of the prior distribution, and logq is the log-density
    of the sampling distribution. (This is importance sampling; see, for
    example, R. M. Neal, Annealed importance sampling, <em>Statistics
      and Computing</em>, 2001.)</p>
    
    <h2 class="hasAnchor" id="prior-on-proportion-of-variance-explained"><a class="anchor" href="#prior-on-proportion-of-variance-explained"></a>Prior on proportion of variance explained</h2>

    <p>Specifying the prior variance of the regression coefficients (sa) can
  be difficult, which is why we have included the option of fitting this
  hyperparameter to the data (see input argument <code>update.sa</code>
  above). However, in many settings, especially when a small number of
  variables are included in the regression model, it is preferrable to
  average over candidate settings of sa instead of fitting sa to the
  data. To choose a set of candidate settings for sa, we have advocated
  for setting sa indirectly through a prior estimate of the proportion
  of variance in the outcome explained by the variables (abbreviated as
  PVE), since it is often more natural to specify the PVE rather than
  the prior variance (see references below). This is technically only
  suitable or the linear regression model (<code>family = "gaussian"</code>),
  but could potentially be used for the linear regression model in an
  approximate way.</p>
<p>For example, one could approximate a uniform prior on the PVE by
  drawing the PVE uniformly between 0 and 1, additionally specifying
  candidate settings for the prior log-odds, then computing the prior
  variance (sa) as follows:</p><pre>sx &lt;- sum(var1.cols(X))
sa &lt;- PVE/(1-PVE)/(sigmoid(log(10)*logodds)*sx)</pre><p>Note that this calculation will yield <code>sa = 0</code> when <code>PVE =
  0</code>, and <code>sa = Inf</code> when <code>PVE = 1</code>.
    Also, bear in mind that if there are additional covariates (Z)
  included in the linear regression model that explain variance in Y,
  then it will usually make more sense to first remove the linear
  effects of these covariates before performing this calculation. The
  PVE would then represent the prior proportion of variance in the
  residuals of Y that are explained by the candidate variables. For an
  example of how to do this, see the code for function <code>varbvs</code>, in
  file varbvs.R, under "preprocessing steps". Alternatively, one could
  include the matrix Z in the calculation above, taking care to ensure
  that the covariates are included in the model with probability 1.</p>
    
    <h2 class="hasAnchor" id="memory-requirements"><a class="anchor" href="#memory-requirements"></a>Memory requirements</h2>

    <p>Finally, we point out that the optimization procedures were
    carefully designed so that they can be applied to very large data
    sets; to date, this code has been tested on data sets with &gt;500,000
    variables and &gt;10,000 samples. An important limiting factor is the
    ability to store the data matrix X in memory. To reduce memory
    requirements, in the MATLAB interface we require that X be single
    precision, but this option is not available in R. Additionally, we
    mostly avoid generating intermediate products that are of the same
    size as X. Only one such intermediate product is generated when
    <code>family = "gaussian"</code>, and none for <code>family = "binomial"</code>.</p>
    
    <h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2>

    <p>An object with S3 class <code>c("varbvs","list")</code>.</p>
<p></p>
<dt>family</dt><dd><p>Either "gaussian" or "binomial".</p></dd>

  <dt>n</dt><dd><p>Number of data samples used to fit model.</p></dd>

  <dt>sigma</dt><dd><p>Settings for sigma (family = "gaussian" only).</p></dd>

  <dt>sa</dt><dd><p>Settings for prior variance parameter.</p></dd>

  <dt>logodds</dt><dd><p>Prior log-odds settings.</p></dd>

  <dt>prior.same</dt><dd><p>TRUE if prior is identical for all variables. When
        logodds is a p x ns matrix, <code>prior.same = FALSE</code>.</p></dd>

  <dt>sa0</dt><dd><p>Scale parameter for prior on hyperparameter sa.</p></dd>

  <dt>n0</dt><dd><p>Degrees of freedom for prior on hyperparameter sa.</p></dd>

  <dt>update.sigma</dt><dd><p>If TRUE, sigma was fit to data for each setting of
    prior logodds (<code>family = "gaussian"</code> only).</p></dd>

  <dt>update.sa</dt><dd><p>If TRUE, sa was fit to data for each setting of prior
    logodds.</p></dd>

  <dt>logw</dt><dd><p>An array with ns elements, in which <code>logw[i]</code> is the
    variational lower bound on the marginal log-likelihood for setting i
    of the hyperparameters. These provide approximate values of the
    marginal log-likelihood for each hyperparameter setting.</p></dd>

  <dt>w</dt><dd><p>Normalized weights (posterior probabilities) for each of the
    hyperparameter settings computed from <code>logw</code> using
    <code>normalizelogweights</code>.</p></dd>
  
  <dt>alpha</dt><dd><p>Variational estimates of posterior inclusion
    probabilities for each hyperparameter setting.</p></dd>

  <dt>mu</dt><dd><p>Variational estimates of posterior mean coefficients for
    each hyperparameter setting.</p></dd>

  <dt>s</dt><dd><p>Variational estimates of posterior variances for each
    hyperparameter setting.</p></dd>

  <dt>pip</dt><dd><p>"Averaged" posterior inclusion probabilities computed as a
    weighted average of the individual PIPs (<code>alpha</code>), with weights
    given by <code>w</code>.</p></dd>

  <dt>beta</dt><dd><p>"Averaged" posterior mean regression coefficients.</p></dd>
  
  <dt>mu.cov</dt><dd><p>Posterior mean regression coefficients for covariates,
    including intercept, for each hyperparameter setting.</p></dd>
  
  <dt>eta</dt><dd><p>Additional variational parameters for <code>family =
    "binomial"</code> only.</p></dd>
  
  <dt>optimize.eta</dt><dd><p>If TRUE, eta was fit to data (<code>family =
  "binomial"</code> only).</p></dd>

  <dt>pve</dt><dd><p>For each hyperparameter setting, and for each variable,
    mean estimate of the proportion of variance in outcome explained
    conditioned on variable being included in the model (<code>family =
    "gaussian"</code> only).</p></dd>
    
  <dt>model.pve</dt><dd><p>Samples drawn from posterior distribution giving
  estimates of proportion of variance in outcome (y) explained by fitted
  variable selection model. This is for <code>family = "gaussian"</code>
  only). To obtain the posterior mean estimate of the proportion of
  variance explained (PVE), for example, simply type
  <code>mean(fit$model.pve)</code>.</p></dd>

    
    <h2 class="hasAnchor" id="references"><a class="anchor" href="#references"></a>References</h2>

    <p>P. Carbonetto and M. Stephens (2012). Scalable variational inference
  for Bayesian variable selection in regression, and its accuracy in
  genetic association studies. <em>Bayesian Analysis</em> <b>7</b>,
  73--108.</p>
<p>Y. Guan and M. Stephens (2011). Bayesian variable selection regression
  for genome-wide association studies and other large-scale
  problems. <em>Annals of Applied Statistics</em> <b>5</b>, 1780--1815.</p>
<p>X. Zhou, P. Carbonetto and M. Stephens (2013). Polygenic modeling with
  Bayesian sparse linear mixed models. <em>PLoS Genetics</em> <b>9</b>,
  e1003264.</p>
    
    <h2 class="hasAnchor" id="see-also"><a class="anchor" href="#see-also"></a>See also</h2>

    <p><code><a href='summary.varbvs.html'>summary.varbvs</a></code>, <code><a href='varbvscoefcred.html'>varbvscoefcred</a></code>,
  <code><a href='varbvs-internal.html'>varbvspve</a></code>, <code><a href='varbvs-internal.html'>varbvsnorm</a></code>,
  <code><a href='varbvs-internal.html'>varbvsbin</a></code>, <code><a href='varbvs-internal.html'>varbvsbinz</a></code>,
  <code><a href='normalizelogweights.html'>normalizelogweights</a></code>, <a href='varbvs-internal.html'>varbvs-internal</a></p>
    

    <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
    <pre class="examples">
<span class='co'># LINEAR REGRESSION EXAMPLE</span>
<span class='co'># -------------------------</span>
<span class='co'># Data are 200 uncorrelated ("unlinked") single nucleotide polymorphisms</span>
<span class='co'># (SNPs) with simulated genotypes, in which the first 20 of them have an</span>
<span class='co'># effect on the outcome. Also generate data for 3 covariates.</span>
<span class='no'>maf</span> <span class='kw'>&lt;-</span> <span class='fl'>0.05</span> + <span class='fl'>0.45</span>*<span class='fu'>runif</span>(<span class='fl'>200</span>)
<span class='no'>X</span>   <span class='kw'>&lt;-</span> (<span class='fu'>runif</span>(<span class='fl'>400</span>*<span class='fl'>200</span>) <span class='kw'>&lt;</span> <span class='no'>maf</span>) + (<span class='fu'>runif</span>(<span class='fl'>400</span>*<span class='fl'>200</span>) <span class='kw'>&lt;</span> <span class='no'>maf</span>)
<span class='no'>X</span>   <span class='kw'>&lt;-</span> <span class='fu'>matrix</span>(<span class='fu'>as.double</span>(<span class='no'>X</span>),<span class='fl'>400</span>,<span class='fl'>200</span>,<span class='kw'>byrow</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
<span class='no'>Z</span>   <span class='kw'>&lt;-</span> <span class='fu'><a href='rand.html'>randn</a></span>(<span class='fl'>400</span>,<span class='fl'>3</span>)

<span class='co'># Generate the ground-truth regression coefficients for the variables</span>
<span class='co'># (X) and additional 3 covariates (Z). Adjust the QTL effects so that</span>
<span class='co'># the variables (SNPs) explain 50 percent of the variance in the</span>
<span class='co'># outcome.</span>
<span class='no'>u</span>    <span class='kw'>&lt;-</span> <span class='fu'>c</span>(-<span class='fl'>1</span>,<span class='fl'>2</span>,<span class='fl'>1</span>)
<span class='no'>beta</span> <span class='kw'>&lt;-</span> <span class='fu'>c</span>(<span class='fu'>rnorm</span>(<span class='fl'>20</span>),<span class='fu'>rep</span>(<span class='fl'>0</span>,<span class='fl'>180</span>))
<span class='no'>beta</span> <span class='kw'>&lt;-</span> <span class='fl'>1</span>/<span class='fu'>sd</span>(<span class='fu'>c</span>(<span class='no'>X</span> <span class='kw'>%*%</span> <span class='no'>beta</span>)) * <span class='no'>beta</span>

<span class='co'># Generate the quantitative trait measurements.</span>
<span class='no'>y</span> <span class='kw'>&lt;-</span> <span class='fu'>c</span>(-<span class='fl'>2</span> + <span class='no'>Z</span> <span class='kw'>%*%</span> <span class='no'>u</span> + <span class='no'>X</span> <span class='kw'>%*%</span> <span class='no'>beta</span> + <span class='fu'>rnorm</span>(<span class='fl'>400</span>))

<span class='co'># Fit the variable selection model.</span>
<span class='no'>fit</span> <span class='kw'>&lt;-</span> <span class='fu'>varbvs</span>(<span class='no'>X</span>,<span class='no'>Z</span>,<span class='no'>y</span>,<span class='kw'>logodds</span> <span class='kw'>=</span> <span class='fu'>seq</span>(-<span class='fl'>3</span>,-<span class='fl'>1</span>,<span class='fl'>0.1</span>))
<span class='fu'>print</span>(<span class='fu'>summary</span>(<span class='no'>fit</span>))

<span class='co'># Compute the posterior mean estimate of hyperparameter sa.</span>
<span class='no'>sa</span> <span class='kw'>&lt;-</span> <span class='fu'>with</span>(<span class='no'>fit</span>,<span class='fu'>sum</span>(<span class='no'>sa</span> * <span class='no'>w</span>))

<span class='co'># Compare estimated outcomes against observed outcomes.</span>
<span class='no'>y.fit</span> <span class='kw'>&lt;-</span> <span class='fu'>predict</span>(<span class='no'>fit</span>,<span class='no'>X</span>,<span class='no'>Z</span>)
<span class='fu'>print</span>(<span class='fu'>cor</span>(<span class='no'>y</span>,<span class='no'>y.fit</span>))

<span class='co'># LOGISTIC REGRESSION EXAMPLE</span>
<span class='co'># ---------------------------</span>
<span class='co'># Data are 100 uncorrelated ("unlinked") single nucleotide polymorphisms</span>
<span class='co'># (SNPs) with simulated genotypes, in which the first 10 of them have an</span>
<span class='co'># effect on the outcome. Also generate data for 2 covariates.</span>
<span class='no'>maf</span> <span class='kw'>&lt;-</span> <span class='fl'>0.05</span> + <span class='fl'>0.45</span>*<span class='fu'>runif</span>(<span class='fl'>100</span>)
<span class='no'>X</span>   <span class='kw'>&lt;-</span> (<span class='fu'>runif</span>(<span class='fl'>750</span>*<span class='fl'>100</span>) <span class='kw'>&lt;</span> <span class='no'>maf</span>) + (<span class='fu'>runif</span>(<span class='fl'>750</span>*<span class='fl'>100</span>) <span class='kw'>&lt;</span> <span class='no'>maf</span>)
<span class='no'>X</span>   <span class='kw'>&lt;-</span> <span class='fu'>matrix</span>(<span class='fu'>as.double</span>(<span class='no'>X</span>),<span class='fl'>750</span>,<span class='fl'>100</span>,<span class='kw'>byrow</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
<span class='no'>Z</span>   <span class='kw'>&lt;-</span> <span class='fu'><a href='rand.html'>randn</a></span>(<span class='fl'>750</span>,<span class='fl'>2</span>)

<span class='co'># Generate the ground-truth regression coefficients for the variables</span>
<span class='co'># (X) and additional 2 covariates (Z).</span>
<span class='no'>u</span>    <span class='kw'>&lt;-</span> <span class='fu'>c</span>(-<span class='fl'>1</span>,<span class='fl'>1</span>)
<span class='no'>beta</span> <span class='kw'>&lt;-</span> <span class='fu'>c</span>(<span class='fl'>0.5</span>*<span class='fu'>rnorm</span>(<span class='fl'>10</span>),<span class='fu'>rep</span>(<span class='fl'>0</span>,<span class='fl'>90</span>))

<span class='co'># Simulate the binary trait (case-control status) as a coin toss with</span>
<span class='co'># success rates given by the logistic regression.</span>
<span class='no'>y</span> <span class='kw'>&lt;-</span> <span class='fu'>as.double</span>(<span class='fu'>runif</span>(<span class='fl'>750</span>) <span class='kw'>&lt;</span> <span class='fu'><a href='varbvs-internal.html'>sigmoid</a></span>(-<span class='fl'>1</span> + <span class='no'>Z</span> <span class='kw'>%*%</span> <span class='no'>u</span> + <span class='no'>X</span> <span class='kw'>%*%</span> <span class='no'>beta</span>))

<span class='co'># Fit the variable selection model.</span>
<span class='no'>fit</span> <span class='kw'>&lt;-</span> <span class='fu'>varbvs</span>(<span class='no'>X</span>,<span class='no'>Z</span>,<span class='no'>y</span>,<span class='st'>"binomial"</span>,<span class='kw'>logodds</span> <span class='kw'>=</span> <span class='fu'>seq</span>(-<span class='fl'>2</span>,-<span class='fl'>0.5</span>,<span class='fl'>0.5</span>))
<span class='fu'>print</span>(<span class='fu'>summary</span>(<span class='no'>fit</span>))</pre>
  </div>
  <div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
    <h2>Contents</h2>
    <ul class="nav nav-pills nav-stacked">
      <li><a href="#arguments">Arguments</a></li>
      
      <li><a href="#regression-models">Regression models</a></li>

      <li><a href="#co-ordinate-ascent-optimization-procedure">Co-ordinate ascent optimization procedure</a></li>

      <li><a href="#variational-approximation">Variational approximation</a></li>

      <li><a href="#averaging-over-hyperparameter-settings">Averaging over hyperparameter settings</a></li>

      <li><a href="#prior-on-proportion-of-variance-explained">Prior on proportion of variance explained</a></li>

      <li><a href="#memory-requirements">Memory requirements</a></li>

      <li><a href="#value">Value</a></li>

      <li><a href="#references">References</a></li>

      <li><a href="#see-also">See also</a></li>
      
      <li><a href="#examples">Examples</a></li>
    </ul>

    <h2>Author</h2>
    Peter Carbonetto <a href='mailto:peter.carbonetto@gmail.com'>peter.carbonetto@gmail.com</a>
  </div>
</div>

      <footer>
      <div class="copyright">
  <p>Developed by Peter Carbonetto, Matthew Stephens, David Gerard.</p>
</div>

<div class="pkgdown">
  <p>Site built with <a href="http://hadley.github.io/pkgdown/">pkgdown</a>.</p>
</div>

      </footer>
   </div>

  </body>
</html>
