<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta charset="utf-8">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />


<meta name="author" content="Lei Sun" />


<title>Truncated ASH</title>

<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/bootstrap.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-1.1/default.css" rel="stylesheet" />
<script src="site_libs/highlightjs-1.1/highlight.js"></script>
<link href="site_libs/font-awesome-4.5.0/css/font-awesome.min.css" rel="stylesheet" />

<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
  pre:not([class]) {
    background-color: white;
  }
</style>
<script type="text/javascript">
if (window.hljs && document.readyState && document.readyState === "complete") {
   window.setTimeout(function() {
      hljs.initHighlighting();
   }, 0);
}
</script>



<style type="text/css">
h1 {
  font-size: 34px;
}
h1.title {
  font-size: 38px;
}
h2 {
  font-size: 30px;
}
h3 {
  font-size: 24px;
}
h4 {
  font-size: 18px;
}
h5 {
  font-size: 16px;
}
h6 {
  font-size: 12px;
}
.table th:not([align]) {
  text-align: left;
}
</style>


</head>

<body>

<style type = "text/css">
.main-container {
  max-width: 940px;
  margin-left: auto;
  margin-right: auto;
}
code {
  color: inherit;
  background-color: rgba(0, 0, 0, 0.04);
}
img {
  max-width:100%;
  height: auto;
}
.tabbed-pane {
  padding-top: 12px;
}
button.code-folding-btn:focus {
  outline: none;
}
</style>


<style type="text/css">
/* padding for bootstrap navbar */
body {
  padding-top: 51px;
  padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar)  */
.section h1 {
  padding-top: 56px;
  margin-top: -56px;
}

.section h2 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h3 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h4 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h5 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h6 {
  padding-top: 56px;
  margin-top: -56px;
}
</style>

<script>
// manage active state of menu based on current page
$(document).ready(function () {
  // active menu anchor
  href = window.location.pathname
  href = href.substr(href.lastIndexOf('/') + 1)
  if (href === "")
    href = "index.html";
  var menuAnchor = $('a[href="' + href + '"]');

  // mark it active
  menuAnchor.parent().addClass('active');

  // if it's got a parent navbar menu mark it active as well
  menuAnchor.closest('li.dropdown').addClass('active');
});
</script>


<div class="container-fluid main-container">

<!-- tabsets -->
<script>
$(document).ready(function () {
  window.buildTabsets("TOC");
});
</script>

<!-- code folding -->






<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">truncash</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="about.html">About</a>
</li>
<li>
  <a href="license.html">License</a>
</li>
      </ul>
      <ul class="nav navbar-nav navbar-right">
        <li>
  <a href="https://github.com/LSun/truncash">
    <span class="fa fa-github"></span>
     
  </a>
</li>
      </ul>
    </div><!--/.nav-collapse -->
  </div><!--/.container -->
</div><!--/.navbar -->

<div class="fluid-row" id="header">



<h1 class="title toc-ignore">Truncated ASH</h1>
<h4 class="author"><em>Lei Sun</em></h4>
<h4 class="date"><em>1/29/2017</em></h4>

</div>


<div id="introduction" class="section level2">
<h2>Introduction</h2>
<p>Suppose we have <span class="math inline">\(m + n\)</span> observations of <span class="math inline">\((\hat\beta_j, \hat s_j)\)</span> in two groups such that, with a pre-specified <span class="math inline">\(t\)</span>.</p>
<p><span class="math display">\[
\text{Group 1: }(\hat\beta_1, \hat s_1), \ldots, (\hat\beta_m, \hat s_m), \text{with } |\hat\beta_j/\hat s_j| \leq t
\]</span></p>
<p><span class="math display">\[
\text{Group 2: }(\hat\beta_{m+1}, \hat s_{m+1}), \ldots, (\hat\beta_{m+n}, \hat s_{m+n}), \text{with } |\hat\beta_j/\hat s_j| &gt; t
\]</span></p>
<p>For Group 1, we’ll only use the information that for each one, <span class="math inline">\(|\hat\beta_j/\hat s_j| \leq t\)</span>; that is, they are moderate observations. For Group 2, we’ll use the full observation <span class="math inline">\((\hat\beta_j, \hat s_j)\)</span>.</p>
<p>Now we have usual ASH assumptions that for both groups,</p>
<p><span class="math display">\[
\begin{array}{c}
\hat\beta_j | \hat s_j, \beta_j \sim N(\beta_j, \hat s_j^2)\\
\beta_j \sim \sum_k\pi_k N(0, \sigma_k^2)
\end{array}
\]</span></p>
<p>Now the estimate of <span class="math inline">\(\pi_k\)</span> from the marginal probability of the data in both groups is</p>
<p><span class="math display">\[
\max_{\{\pi_k\}}\sum\limits_{j=1}^m\log(\sum_k\pi_k(2\Phi(\frac{t\hat s_j}{\sqrt{\sigma_k^2+\hat s_j^2}}) - 1))
+
\sum\limits_{l=1}^{n}\log(\sum_k\pi_k N(\hat\beta_{m + l}; 0, \sigma_k^2 + \hat s_{m + l}^2))
\]</span></p>
<p>where <span class="math inline">\(\Phi(\cdot)\)</span> is the cdf of <span class="math inline">\(N(0, 1)\)</span>, and <span class="math inline">\(N(\cdot; 0, \sigma_k^2 + \hat s_{m + l}^2)\)</span> is the pdf of <span class="math inline">\(N(0, \sigma_k^2+\hat s_j^2)\)</span>.</p>
<p>Note that</p>
<ol style="list-style-type: decimal">
<li><p>The second part of the objective function is the standard objective function for the regular ASH. The first part of the objective function is new, but it would be as if we are using a different likelihood, so it wouldn’t make any difference numerically.</p></li>
<li><p>The first part of the objective,</p></li>
</ol>
<p><span class="math display">\[
\max_{\{\pi_k\}}\sum\limits_{j=1}^m\log(\sum_k\pi_k(2\Phi(\frac{t\hat s_j}{\sqrt{\sigma_k^2+\hat s_j^2}}) - 1))
\]</span></p>
<p>and its special case, when we set <span class="math inline">\(\hat s_j \equiv 1\)</span>,</p>
<p><span class="math display">\[
m\log(\sum_k\pi_k(2\Phi(\frac{t}{\sqrt{\sigma_k^2+1}}) - 1))
\]</span></p>
<p>will move more weight to <span class="math inline">\(\hat\pi_1\)</span>. Customarily, we assume <span class="math inline">\(\sigma_1 \leq \cdots \leq \sigma_K\)</span>, so to maximize this first part alone will be essentially to maximize <span class="math inline">\(\hat\pi_1\)</span>. In other words, if we only have Group 1 and no “extreme” observations in Group 2, the estimated prior of <span class="math inline">\(\beta\)</span> will be a point mass at <span class="math inline">\(0\)</span>. This is true no matter whether we are using <span class="math inline">\(z\)</span>-scores (thus making <span class="math inline">\(\hat s_j \equiv 1\)</span> effectively) or not, although using <span class="math inline">\(z\)</span>-scores will make the case more obviously.</p>
</div>
<div id="implementation" class="section level2">
<h2>Implementation</h2>
<p>Here is the main <code>truncash</code> function, utilizing <code>autoselect.mixsd</code>, <code>mixIP</code>, <code>normalmix</code> from <code>ashr</code>.</p>
<p><a href="https://github.com/LSun/truncash/blob/master/code/truncash.R"><code>truncash</code></a></p>
</div>
<div id="simulation" class="section level2">
<h2>Simulation</h2>
<div id="the-simplest-case-g-0.9delta_0-0.1n0-1-independent-observations" class="section level4">
<h4>The simplest case: <span class="math inline">\(g = 0.9\delta_0 + 0.1N(0, 1)\)</span>, independent observations</h4>
<p><span class="math inline">\(1000\)</span> pairs of <span class="math inline">\((\hat\beta, \hat s \equiv 1)\)</span> are simulated, so that <span class="math inline">\(900\)</span> of them come from <span class="math inline">\(N(0, 1)\)</span> and another 100 from <span class="math inline">\(N(\beta, 1)\)</span> with <span class="math inline">\(\beta \sim N(0, 1)\)</span>. The <span class="math inline">\(1000\)</span> pairs of observations are fed to both <code>ash</code> and <code>truncash</code>; <span class="math inline">\(\hat\pi_0\)</span> and MSE are obtained. Results from 100 runs are plotted below.</p>
<pre class="r"><code>source(&quot;../code/truncash.R&quot;)

set.seed(777)

sebetahat = rep(1, 1000)
t = qnorm(0.975)

pihat0_ash = pihat0_truncash = mse_ash = mse_truncash = mse_mle = lfsr0.05_ash = lfsr0.05_truncash = c()


for (i in 1:100) {
beta1 = rep(0, 900)
betahat1 = rnorm(length(beta1), beta1, 1)
beta2 = rnorm(100, 0, 1)
betahat2 = rnorm(length(beta2), beta2, 1)
beta = c(beta1, beta2)
betahat = c(betahat1, betahat2)

res.ash = ash(betahat, 1, mixcompdist = &quot;normal&quot;, method = &quot;fdr&quot;)
res.truncash = truncash(betahat, sebetahat, t)

pihat0_ash[i] = get_fitted_g(res.ash)$pi[1]
pihat0_truncash[i] = get_fitted_g(res.truncash)$pi[1]

mse_ash[i] = mean((get_pm(res.ash) - beta)^2)
mse_truncash[i] = mean((get_pm(res.truncash) - beta)^2)
mse_mle[i] = mean((betahat - beta)^2)

lfsr0.05_ash = mean(get_lfsr(res.ash) &lt;= 0.05)
lfsr0.05_truncash = mean(get_lfsr(res.truncash) &lt;= 0.05)
}

plot(pihat0_ash, pihat0_truncash, xlab = &quot;ash&quot;, ylab = &quot;truncash&quot;, main = expression(hat(pi)[0]),
     xlim = c(min(c(pihat0_ash, pihat0_truncash)), max(c(pihat0_ash, pihat0_truncash))),
     ylim = c(min(c(pihat0_ash, pihat0_truncash)), max(c(pihat0_ash, pihat0_truncash)))
     )
abline(0, 1, lty = 3, col = &quot;blue&quot;)
abline(v = 0.9, lty = 3, col = &quot;blue&quot;)
abline(h = 0.9, lty = 3, col = &quot;blue&quot;)
points(0.9, 0.9, col = &quot;red&quot;, cex = 1.5, pch = 19)</code></pre>
<p><img src="truncash_files/figure-html/unnamed-chunk-1-1.png" width="672" /></p>
<pre class="r"><code>plot(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9), xlab = &quot;ash&quot;, ylab = &quot;truncash&quot;,
     main = expression(paste(&quot;|&quot;, hat(pi)[0] - 0.9, &quot;|&quot;)), 
     xlim = c(0, max(c(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9)))),
     ylim = c(0, max(c(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9))))
     )
abline(0, 1, lty = 3, col = &quot;blue&quot;)
points(0.9, 0.9, col = &quot;red&quot;, cex = 1.5, pch = 19)</code></pre>
<p><img src="truncash_files/figure-html/unnamed-chunk-1-2.png" width="672" /></p>
<pre class="r"><code>plot(mse_ash, mse_truncash, xlab = &quot;ash&quot;, ylab = &quot;truncash&quot;, main = &quot;MSE&quot;,
     xlim = c(min(c(mse_ash, mse_truncash)), max(c(mse_ash, mse_truncash))),
     ylim = c(min(c(mse_ash, mse_truncash)), max(c(mse_ash, mse_truncash)))
     )
abline(0, 1, lty = 3, col = &quot;blue&quot;)</code></pre>
<p><img src="truncash_files/figure-html/unnamed-chunk-1-3.png" width="672" /></p>
</div>
</div>




</div>

<script>

// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
  $('tr.header').parent('thead').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
  bootstrapStylePandocTables();
});


</script>

<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>