<!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" />
<title>MASH v FLASH detailed simulation study</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/cosmo.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/jqueryui-1.11.4/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<script src="site_libs/navigation-1.1/codefolding.js"></script>
<link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/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) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (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 -->
<style type="text/css">
.code-folding-btn { margin-bottom: 4px; }
</style>
<script>
$(document).ready(function () {
window.initializeCodeFolding("hide" === "show");
});
</script>
<script>
$(document).ready(function () {
// move toc-ignore selectors from section div to header
$('div.section.toc-ignore')
.removeClass('toc-ignore')
.children('h1,h2,h3,h4,h5').addClass('toc-ignore');
// establish options
var options = {
selectors: "h1,h2,h3",
theme: "bootstrap3",
context: '.toc-content',
hashGenerator: function (text) {
return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_').toLowerCase();
},
ignoreSelector: ".toc-ignore",
scrollTo: 0
};
options.showAndHide = true;
options.smoothScroll = true;
// tocify
var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
.toc-content {
padding-left: 30px;
padding-right: 40px;
}
div.main-container {
max-width: 1200px;
}
div.tocify {
width: 20%;
max-width: 260px;
max-height: 85%;
}
@media (min-width: 768px) and (max-width: 991px) {
div.tocify {
width: 25%;
}
}
@media (max-width: 767px) {
div.tocify {
width: 100%;
max-width: none;
}
}
.tocify ul, .tocify li {
line-height: 20px;
}
.tocify-subheader .tocify-item {
font-size: 0.90em;
padding-left: 25px;
text-indent: 0;
}
.tocify .list-group-item {
border-radius: 0px;
}
</style>
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row-fluid">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>
<div class="toc-content col-xs-12 col-sm-8 col-md-9">
<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">MASHvFLASH</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>
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="https://github.com/willwerscheid/MASHvFLASH">
<span class="fa fa-github"></span>
</a>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<!-- Add a small amount of space between sections. -->
<style type="text/css">
div.section {
padding-top: 12px;
}
</style>
<div class="fluid-row" id="header">
<div class="btn-group pull-right">
<button type="button" class="btn btn-default btn-xs dropdown-toggle" data-toggle="dropdown" aria-haspopup="true" aria-expanded="false"><span>Code</span> <span class="caret"></span></button>
<ul class="dropdown-menu" style="min-width: 50px;">
<li><a id="rmd-show-all-code" href="#">Show All Code</a></li>
<li><a id="rmd-hide-all-code" href="#">Hide All Code</a></li>
</ul>
</div>
<h1 class="title toc-ignore">MASH v FLASH detailed simulation study</h1>
</div>
<p><strong>Last updated:</strong> 2018-07-22</p>
<strong>workflowr checks:</strong> <small>(Click a bullet for more information)</small>
<ul>
<li>
<p><details> <summary> <strong style="color:blue;">✔</strong> <strong>R Markdown file:</strong> up-to-date </summary></p>
<p>Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.</p>
</details>
</li>
<li>
<p><details> <summary> <strong style="color:blue;">✔</strong> <strong>Environment:</strong> empty </summary></p>
<p>Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.</p>
</details>
</li>
<li>
<p><details> <summary> <strong style="color:blue;">✔</strong> <strong>Seed:</strong> <code>set.seed(20180609)</code> </summary></p>
<p>The command <code>set.seed(20180609)</code> was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.</p>
</details>
</li>
<li>
<p><details> <summary> <strong style="color:blue;">✔</strong> <strong>Session information:</strong> recorded </summary></p>
<p>Great job! Recording the operating system, R version, and package versions is critical for reproducibility.</p>
</details>
</li>
<li>
<p><details> <summary> <strong style="color:blue;">✔</strong> <strong>Repository version:</strong> <a href="https://github.com/willwerscheid/MASHvFLASH/tree/044fb97dcd6a939dd85448f777624479a5f72588" target="_blank">044fb97</a> </summary></p>
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated. <br><br> Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use <code>wflow_publish</code> or <code>wflow_git_commit</code>). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
<pre><code>
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: docs/.DS_Store
Ignored: docs/images/.DS_Store
Ignored: output/.DS_Store
Unstaged changes:
Modified: analysis/MASHvFLASHgtex2.Rmd
Modified: analysis/mashvflash.Rmd
Modified: code/sims2.R
Modified: output/README.md
Modified: output/sims2mse.rds
Modified: output/sims2prash.rds
Modified: output/sims2prm.rds
Modified: output/sims2prpn.rds
</code></pre>
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes. </details>
</li>
</ul>
<details> <summary> <small><strong>Expand here to see past versions:</strong></small> </summary>
<ul>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
File
</th>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
<th style="text-align:left;">
Message
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/044fb97dcd6a939dd85448f777624479a5f72588/analysis/MASHvFLASHsims2.Rmd" target="_blank">044fb97</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-07-22
</td>
<td style="text-align:left;">
wflow_publish(“analysis/MASHvFLASHsims2.Rmd”)
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/willwerscheid/MASHvFLASH/8720e9c0eec4613a756a8113d69999cb1622cf7f/docs/MASHvFLASHsims2.html" target="_blank">8720e9c</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-07-01
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/1874a0d6a3e899671c92b176fb197918b9102aa5/analysis/MASHvFLASHsims2.Rmd" target="_blank">1874a0d</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-07-01
</td>
<td style="text-align:left;">
wflow_publish(c(“analysis/MASHvFLASHgtex2.Rmd”,
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/willwerscheid/MASHvFLASH/00ac712185e95359920e5107965a56d6968b29dd/docs/MASHvFLASHsims2.html" target="_blank">00ac712</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/81adb8ce99944c3f18a6e1a300b50b3d66f276b4/analysis/MASHvFLASHsims2.Rmd" target="_blank">81adb8c</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
<td style="text-align:left;">
wflow_publish(“analysis/MASHvFLASHsims2.Rmd”)
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/willwerscheid/MASHvFLASH/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/MASHvFLASHsims2.html" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/5d6259193e4d6246c25c85b88d9e8f50330ced16/analysis/MASHvFLASHsims2.Rmd" target="_blank">5d62591</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
<td style="text-align:left;">
wflow_publish(“analysis/MASHvFLASHsims2.Rmd”)
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/willwerscheid/MASHvFLASH/e875bf9fadfaa94260e2519e6baa0bb629f08a39/docs/MASHvFLASHsims2.html" target="_blank">e875bf9</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/7e4bbe002116cc56a1fb3861550bc93f86a5dc10/analysis/MASHvFLASHsims2.Rmd" target="_blank">7e4bbe0</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
<td style="text-align:left;">
wflow_publish(c(“analysis/MASHvFLASHsims2.Rmd”))
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/willwerscheid/MASHvFLASH/97fa87ceb554d2998662e546e58eceebcea3ce5d/docs/MASHvFLASHsims2.html" target="_blank">97fa87c</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/42cd89c485cdf6fb255679adcd67c8e2b098018e/analysis/MASHvFLASHsims2.Rmd" target="_blank">42cd89c</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
<td style="text-align:left;">
wflow_publish(c(“analysis/MASHvFLASHsims2.Rmd”,
</td>
</tr>
<tr>
<td style="text-align:left;">
html
</td>
<td style="text-align:left;">
<a href="https://cdn.rawgit.com/willwerscheid/MASHvFLASH/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/MASHvFLASHsims2.html" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
<td style="text-align:left;">
Build site.
</td>
</tr>
<tr>
<td style="text-align:left;">
Rmd
</td>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/411c3b41d268eab99f8c2d2198008b2301d63bdb/analysis/MASHvFLASHsims2.Rmd" target="_blank">411c3b4</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
<td style="text-align:left;">
wflow_publish(“analysis/MASHvFLASHsims2.Rmd”)
</td>
</tr>
</tbody>
</table>
</ul>
<p></details></p>
<hr />
<div id="introduction" class="section level2">
<h2>Introduction</h2>
<p>Here I study in detail simulations from a MASH model that extends the “model with independent, unique, and shared effects” from my larger <a href="MASHvFLASHsims.html">simulation study</a>. I run 10 simulations, each of which simulates data for 20 conditions and 1200 tests. I use 17 different covariance structures, each of which are used to simulate 25 tests. The other 775 tests are null across all conditions.</p>
</div>
<div id="independent-effects" class="section level2">
<h2>Independent effects</h2>
<p>Effects were nonnull for all conditions and generated independently from a <span class="math inline">\(N(0, \sigma^2)\)</span> distribution. I simulated <span class="math inline">\((1)\)</span> small independent effects (<span class="math inline">\(\sigma^2 = 2^2\)</span>), <span class="math inline">\((2)\)</span> large independent effects (<span class="math inline">\(\sigma^2 = 5^2\)</span>), and <span class="math inline">\((3)\)</span> independent effects of varying sizes (with <span class="math inline">\(\sigma^2\)</span> ranging from <span class="math inline">\(1^2\)</span> to <span class="math inline">\(5^2\)</span>).</p>
<p>Notice that cases 1 and 2 are covered by “canonical” covariance matrices in MASH, but 3 is not.</p>
</div>
<div id="identical-effects" class="section level2">
<h2>Identical effects</h2>
<p>Effects were nonnull for all conditions, with an effect size that was identical across conditions. The unique effect size was generated from a <span class="math inline">\(N(0, \sigma^2)\)</span> distribution. Similar to the above, I simulated <span class="math inline">\((4)\)</span> small identical effects (<span class="math inline">\(\sigma^2 = 2^2\)</span>) and <span class="math inline">\((5)\)</span> large identical effects (<span class="math inline">\(\sigma^2 = 5^2\)</span>).</p>
<p>Both of these cases are covered by canonical covariance matrices in MASH.</p>
</div>
<div id="rank-one-effects" class="section level2">
<h2>Rank-one effects</h2>
<p>These are similar to “identical effects” in that the covariance matrix has rank one (so that conditions 2-20 are always fixed multiples of condition 1), but here, the effect sizes vary.</p>
<p><span class="math inline">\((6)\)</span> The effect size for condition 1 was generated from a <span class="math inline">\(N(0, 1)\)</span> distribution, and conditions 1-20 were multiples evenly spaced between 1 and 5.</p>
<p>Unlike identical effects, rank-one effects are not directly modeled by canonical covariance matrices.</p>
</div>
<div id="unique-effects" class="section level2">
<h2>Unique effects</h2>
<p>Effects were nonnull in one condition only, with the nonnull effect simulated from a <span class="math inline">\(N(0, \sigma^2)\)</span> distribution. I simulated <span class="math inline">\((7)\)</span> small effects (<span class="math inline">\(\sigma^2 = 3^2\)</span>) unique to condition 1 and <span class="math inline">\((8)\)</span> large effects (<span class="math inline">\(\sigma^2 = 8^2\)</span>) unique to condition 2.</p>
<p>These effects are directly modeled by canonical covariance matrices.</p>
</div>
<div id="shared-effects" class="section level2">
<h2>Shared effects</h2>
<p>Effects were nonnull in several conditions. The nonnull effects were either identical across conditions or fixed multiples of one another. I included <span class="math inline">\((9)\)</span> medium effects (<span class="math inline">\(\sigma^2 = 3^2\)</span>) identical across 3 conditions (3-5), <span class="math inline">\((10)\)</span> small effects (<span class="math inline">\(\sigma^2 = 2^2\)</span>) identical across 10 conditions (1-10), and <span class="math inline">\((11)\)</span> effects of differing sizes over 5 conditions (6-10), with variances ranging from <span class="math inline">\(2^2\)</span> to <span class="math inline">\(4^2\)</span>.</p>
<p>None of these types of effects are modeled by canonical covariance matrices.</p>
</div>
<div id="random-covariance" class="section level2">
<h2>Random covariance</h2>
<p>I included three random covariance structures in which effects were nonnull across all conditions. In each case the random covariance matrix <span class="math inline">\(A^T A\)</span> was generated by sampling the entries of <span class="math inline">\(A\)</span> independently from a <span class="math inline">\(N(0, 2^2)\)</span> distribution. I included <span class="math inline">\((12)\)</span> a rank-5 random covariance matrix (with <span class="math inline">\(A \in \mathbb{R}^{5 \times 44}\)</span>), <span class="math inline">\((13)\)</span> a rank-10 random covariance matrix <span class="math inline">\((A \in \mathbb{R}^{10 \times 44})\)</span>, and <span class="math inline">\((14)\)</span> a full-rank random covariance matrix <span class="math inline">\((A \in \mathbb{R}^{44 \times 44})\)</span>.</p>
</div>
<div id="combinations-of-independent-identical-and-unique-effects" class="section level2">
<h2>Combinations of independent, identical, and unique effects</h2>
<p>Finally, I included several combinations of the above types of effects. In particular, I simulated <span class="math inline">\((15)\)</span> small identical effects (covariance type 4) plus large independent effects (type 2), <span class="math inline">\((16)\)</span> small independent effects (type 1) plus a large unique effect (type 8), and <span class="math inline">\((17)\)</span> small identical effects (type 4) plus a large unique effect (type 8).</p>
</div>
<div id="fitting-methods" class="section level2">
<h2>Fitting methods</h2>
<p>For each simulation, I fitted a MASH model and two FLASH models using the “one-hots last” method described in my larger <a href="MASHvFLASHsims.html">simulation study</a>. One FLASH model used the point-normal approach to solve the ebnm problem (<code>ebnm_fn = ebnm_pn</code>); the other used <code>ashr</code> (<code>ebnm_fn = ebnm_ash</code>).</p>
<p>First I load results.</p>
<pre class="r"><code>all.rrmses <- sqrt(readRDS("./output/sims2mse.rds"))
pr.pn <- readRDS("./output/sims2prpn.rds")
pr.ash <- readRDS("./output/sims2prash.rds")
pr.m <- readRDS("./output/sims2prm.rds")
# Nulls
nullidx <- c(7, 9, 11, 13, 15, 23)
null.names <- c("UnSm(7)", "UnLg(8)", "Sh3(9)", "Sh10(10)", "Sh5(11)", "Null")
# Effects that are independent or identical across conditions
indididx <- 1:6
indid.names <- c("IndSm(1)", "IndLg(2)", "IndDif(3)", "IdSm(4)", "IdLg(5)", "Rnk1(6)")
# Effects that are unique to one or several conditions
unshidx <- c(8, 10, 12, 14, 16)
unsh.names <- c("UnSm(7)", "UnLg(8)", "Sh3(9)", "Sh10(10)", "Sh5(11)")
# Other
otheridx <- c(17:22)
other.names <- c("Rn5(12)", "Rn10(13)", "Rnd(14)", "IdIn(15)", "InUn(16)", "IdUn(17)")</code></pre>
</div>
<div id="results-rrmse" class="section level2">
<h2>Results: RRMSE</h2>
<p>Here I give a detailed breakdown of the relative root mean-squared errors (that is, the RMSE for each fit object, divided by the RMSE that would be obtained by simply using the observed data <span class="math inline">\(Y\)</span> to estimate the “true effects” <span class="math inline">\(X\)</span>). In addition to calculating the RRMSE separately for each covariance type, I also separately consider null effects and nonnull effects.</p>
<p>MASH does much better in shrinking null effects towards zero, particularly for tests that are null across all conditions or that are unique to a single condition (covariance types 7-8). MASH also does consistently better in estimating nonnull effects when effects are independent across conditions (types 1-3). (In such cases, FLASH does worse than the naive estimate <span class="math inline">\(\hat{X} = Y\)</span>.)</p>
<p>Results for other covariance types are more varied: MASH appears to do better on identical effects (types 4-5, whose covariance structures are included as “canonical”), but slightly worse on rank-one effects (type 6, which is not accounted for by canonical covariance matrices). FLASH seems to outperform MASH on shared effects (types 9-11), but (somewhat surprisingly) MASH does better on random covariance matrices and on combinations involving independent effects (types 12-16), with an RRMSE near 1 (whereas FLASH is again outperformed by the naive estimate <span class="math inline">\(\hat{X} = Y\)</span>). Finally, FLASH does much better on a combination of identical and unique effects (type 17).</p>
<p>In general, the ash fits perform very similarly to the point-normal fits.</p>
<p>RRMSE for null effects is as follows (where “null” below refers to tests that are null across <em>all</em> conditions):</p>
<pre class="r"><code>par(mar = c(5,4,4,6))
legend.names <- c("FL-pn", "FL-ash", "MASH")
legend.args <- list(x="right", bty="n", inset=c(-0.25,0), xpd=T)
barplot(all.rrmses[,nullidx], names.arg = null.names, beside=T,
ylim=c(0, 1), ylab="RRMSE", main="RRMSE for null effects",
legend.text=legend.names, args.legend=legend.args)</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/mse_null-1.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of mse_null-1.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/mse_null-1.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/mse_null-1.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
<p>And RRMSE for nonnull effects is as follows.</p>
<pre class="r"><code>par(mar = c(5,4,4,6))
barplot(all.rrmses[,indididx], names.arg = indid.names, beside=T,
ylim=c(0, 2), ylab="RRMSE", main="RRMSE for nonnull effects",
legend.text=legend.names, args.legend=legend.args)</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/mse_nonnull-1.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of mse_nonnull-1.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/mse_nonnull-1.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/mse_nonnull-1.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>barplot(all.rrmses[,unshidx], names.arg = unsh.names, beside=T,
ylim=c(0, 2), ylab="RRMSE", main="",
legend.text=legend.names, args.legend=legend.args)</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/mse_nonnull-2.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of mse_nonnull-2.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/mse_nonnull-2.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/mse_nonnull-2.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>barplot(all.rrmses[,otheridx], names.arg = other.names, beside=T,
ylim=c(0, 2), ylab="RRMSE", main="",
legend.text=legend.names, args.legend=legend.args)</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/mse_nonnull-3.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of mse_nonnull-3.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/mse_nonnull-3.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/mse_nonnull-3.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
</div>
<div id="results-fprtpr" class="section level2">
<h2>Results: FPR/TPR</h2>
<p>As in the larger <a href="MASHvFLASHsims.html">simulation study</a>, I evaluate true and false positive rates using the built-in function <code>get_lfsr()</code> for MASH and by simulating from the posterior for FLASH. For each covariance structure, I plot the true positive rate for a given covariance structure against the <em>overall</em> false positive rate.</p>
<p>An examination of the ROC curves leads to somewhat different conclusions to the above. As above, MASH does much better with respect to independent effects (types 1-3), random covariance matrices (types 12-14), and a combination of identical and independent effects (type 15), while FLASH outperforms MASH on a combination of identical and unique effects (type 17). Interestingly, point-normal priors outperform ash priors in all of these cases.</p>
<p>MASH does slightly better with respect to both identical effects (types 4-5) and rank-one effects (type 6), even though FLASH did better on the latter in terms of RRMSE. All methods perform similarly for both unique and shared effects (types 7-11), with ash priors doing somewhat better than point-normal priors for unique effects. Both MASH and FLASH-pn do well on a combination of independent and unique effects, while FLASH-ash performs poorly.</p>
<p>In general, it appears that one can safely opt for the faster point-normal priors over the slower but more flexible ash priors.</p>
<pre class="r"><code>get_fpr <- function(pr) {
nullidx <- c(7, 9, 11, 13, 15)
fp <- 25 * rowSums(pr[, nullidx]) + 775 * (pr[, 23])
fp / (25 * length(nullidx) + 775)
}
plot_fprvtpr <- function(idx, typename) {
plot(get_fpr(pr.pn), pr.pn[, idx], type='l', col="red1", lty=2,
xlab="FPR", ylab="TPR", ylim=c(0, 1), main=typename)
lines(get_fpr(pr.ash), pr.ash[, idx], col="red4", lty=2)
lines(get_fpr(pr.m), pr.m[, idx], col="blue")
legend("bottomright", legend=c("Fl-pn", "Fl-ash", "MASH"),
col=c("red1", "red4", "blue"), lty=c(2, 2, 1))
}
par(mfrow=c(1, 3))
plot_fprvtpr(1, "Small independent (1)")
plot_fprvtpr(2, "Large independent (2)")
plot_fprvtpr(3, "Independent of varying size (3)")</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/tpr-1.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of tpr-1.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/tpr-1.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/tpr-1.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>plot_fprvtpr(4, "Small identical (4)")
plot_fprvtpr(5, "Large identical (5)")
plot_fprvtpr(6, "Rank-one (6)")</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/tpr-2.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of tpr-2.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/tpr-2.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/tpr-2.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>par(mfrow=c(1, 2))
plot_fprvtpr(8, "Small unique (7)")
plot_fprvtpr(10, "Large unique (8)")</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/tpr-3.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of tpr-3.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/tpr-3.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/17ae3f327922f59fc22e50d5b188f66ace7a38d3/docs/figure/MASHvFLASHsims2.Rmd/tpr-3.png" target="_blank">17ae3f3</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-24
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>par(mfrow=c(1, 3))
plot_fprvtpr(12, "Shared (3 conditions) (9)")
plot_fprvtpr(14, "Shared (10 conditions) (10)")
plot_fprvtpr(16, "Shared (varying sizes) (11)")</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/tpr-4.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of tpr-4.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/tpr-4.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>plot_fprvtpr(17, "Random rank-5 (12)")
plot_fprvtpr(18, "Random rank-10 (13)")
plot_fprvtpr(19, "Random full-rank (14)")</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/tpr-5.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of tpr-5.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/tpr-5.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
</tbody>
</table>
<p></details></p>
<pre class="r"><code>plot_fprvtpr(20, "Identical plus independent (15)")
plot_fprvtpr(21, "Independent plus unique (16)")
plot_fprvtpr(22, "Identical plus unique (17)")</code></pre>
<p><img src="figure/MASHvFLASHsims2.Rmd/tpr-6.png" width="672" style="display: block; margin: auto;" /></p>
<details> <summary><em>Expand here to see past versions of tpr-6.png:</em></summary>
<table style="border-collapse:separate; border-spacing:5px;">
<thead>
<tr>
<th style="text-align:left;">
Version
</th>
<th style="text-align:left;">
Author
</th>
<th style="text-align:left;">
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
<a href="https://github.com/willwerscheid/MASHvFLASH/blob/ce6e790072d0954a174fb1b1ab4a4c0228824e58/docs/figure/MASHvFLASHsims2.Rmd/tpr-6.png" target="_blank">ce6e790</a>
</td>
<td style="text-align:left;">
Jason Willwerscheid
</td>
<td style="text-align:left;">
2018-06-26
</td>
</tr>
</tbody>
</table>
<p></details></p>
</div>
<div id="code" class="section level2">
<h2>Code</h2>
<p>Click “Code” to view the code used to obtain the above results.</p>
<pre class="r"><code>devtools::load_all("/Users/willwerscheid/GitHub/flashr/")
library(mashr)
source("./code/fits.R")
source("./code/sims.R")
source("./code/utils.R")
set.seed(1)
n <- 20
p <- 1200
nsims <- 10
nsamp <- 200 # for sampling lfsr (FLASH fits)
ncol <- 25 # number of columns that exhibit each variance type
t <- 0.05 # "significance" threshold
Sigma <- list()
# Independent (small)
Sigma[[1]] <- diag(2^2, n)
# Independent (large)
Sigma[[2]] <- diag(5^2, n)
# Independent (different sizes)
sizes <- seq(1, 5, length.out=n)
Sigma[[3]] <- diag(sizes^2)
# Identical (small)
Sigma[[4]] <- matrix(2^2, nrow=n, ncol=n)
# Identical (large)
Sigma[[5]] <- matrix(5^2, nrow=n, ncol=n)
# Rank-one
Sigma[[6]] <- outer(sizes, sizes)
zeros <- matrix(0, nrow=n, ncol=n)
for (j in 7:11) {
Sigma[[j]] <- zeros
}
# Unique (small)
uniqsmidx <- 1
Sigma[[7]][uniqsmidx, uniqsmidx] <- 3^2
# Unique (large)
uniqlgidx <- 2
Sigma[[8]][uniqlgidx, uniqlgidx] <- 8^2
# Shared (3 conditions)
shar3idx <- 3:5
Sigma[[9]][shar3idx, shar3idx] <- matrix(3^2, nrow=3, ncol=3)
# Shared (10 conditions)
shar10idx <- 1:10
Sigma[[10]][shar10idx, shar10idx] <- matrix(2^2, nrow=10, ncol=10)
# Shared (5 conditions, different sizes)
shar5idx <- 6:10
sizes <- seq(2, 4, length.out=5)
Sigma[[11]][shar5idx, shar5idx] <- outer(sizes, sizes)
# Rank-5
A <- matrix(rnorm(n*5, 0, 2), nrow=5, ncol=n)
Sigma[[12]] <- t(A) %*% A
# Rank-10
A <- matrix(rnorm(n*10, 0, 2), nrow=10, ncol=n)
Sigma[[13]] <- t(A) %*% A
# Random
A <- matrix(rnorm(n*n, 0, 2), nrow=n, ncol=n)
Sigma[[14]] <- t(A) %*% A
# Large independent plus small identical
Sigma[[15]] <- Sigma[[2]] + Sigma[[4]]
# Small independent plus large unique
Sigma[[16]] <- Sigma[[1]] + Sigma[[8]]
# Small identical plus large unique
Sigma[[17]] <- Sigma[[4]] + Sigma[[8]]
ntypes <- 17
partnames <- c("IndSm", "IndLg", "IndDiff",
"IdentSm", "IdentLg", "Rank1",
"UniqSmNull", "UniqSmNonnull",
"UniqLgNull", "UniqLgNonnull",
"Shar3Null", "Shar3Nonnull",
"Shar10Null", "Shar10Nonnull",
"Shar5Null", "Shar5Nonnull",
"Rank5", "Rank10", "Random",
"IndIdent", "IndUniq", "IdentUniq",
"Null")
partxidx <- list(1:n, 1:n, 1:n, 1:n, 1:n, 1:n,
setdiff(1:n, uniqsmidx), uniqsmidx,
setdiff(1:n, uniqlgidx), uniqlgidx,
setdiff(1:n, shar3idx), shar3idx,
setdiff(1:n, shar10idx), shar10idx,
setdiff(1:n, shar5idx), shar5idx,
1:n, 1:n, 1:n, 1:n, 1:n, 1:n, 1:n)
partyidx <- list(1:ncol, ncol + 1:ncol, 2*ncol + 1:ncol,
3*ncol + 1:ncol, 4*ncol + 1:ncol, 5*ncol + 1:ncol,
6*ncol + 1:ncol, 6*ncol + 1:ncol,
7*ncol + 1:ncol, 7*ncol + 1:ncol,
8*ncol + 1:ncol, 8*ncol + 1:ncol,
9*ncol + 1:ncol, 9*ncol + 1:ncol,
10*ncol + 1:ncol, 10*ncol + 1:ncol,
11*ncol + 1:ncol, 12*ncol + 1:ncol,
13*ncol + 1:ncol, 14*ncol + 1:ncol,
15*ncol + 1:ncol, 16*ncol + 1:ncol,
(17*ncol + 1):p)
nparts <- length(partnames)
partition_by_type <- function(X) {
ret <- rep(0, nparts)
for (i in 1:nparts) {
ret[i] <- mean(X[partxidx[[i]], partyidx[[i]]])
}
names(ret) <- partnames
ret
}
mses <- matrix(0, nrow=3, ncol=nparts)
ts <- c(seq(.005, .05, by=.005), seq(.06, .1, by=.01), seq(.2, 1.0, by=.1))
pr.pn <- matrix(0, nrow=length(ts), ncol=nparts)
pr.ash <- matrix(0, nrow=length(ts), ncol=nparts)
pr.m <- matrix(0, nrow=length(ts), ncol=nparts)
for (i in 1:nsims) {
X <- matrix(0, nrow=n, ncol=p)
for (j in 1:ntypes) {
start_col = 1 + ncol*(j-1)
end_col = ncol*j
X[, start_col:end_col] <- t(MASS::mvrnorm(ncol, rep(0, n), Sigma[[j]]))
}
Y <- X + rnorm(n*p)
fl.pn <- fit_flash(Y, Kmax=30, methods=3, ebnm_fn="ebnm_pn") # OHL
fl.ash <- fit_flash(Y, Kmax=30, methods=3, ebnm_fn="ebnm_ash")
m <- fit_mash(Y)
base.se <- (Y - X)^2
base.mse <- partition_by_type(base.se)
fl.pn.se <- (flash_get_fitted_values(fl.pn$fits$OHL) - X)^2
fl.pn.mse <- partition_by_type(fl.pn.se) / base.mse
fl.ash.se <- (flash_get_fitted_values(fl.ash$fits$OHL) - X)^2
fl.ash.mse <- partition_by_type(fl.ash.se) / base.mse
m.se <- (t(get_pm(m$m)) - X)^2
m.mse <- partition_by_type(m.se) / base.mse
mses[1,] <- mses[1,] + fl.pn.mse
mses[2,] <- mses[2,] + fl.ash.mse
mses[3,] <- mses[3,] + m.mse
fl.pn.sampler <- flash_sampler(Y, fl.pn$fits$OHL, fixed="loadings")
fl.pn.samp <- fl.pn.sampler(nsamp)
fl.pn.lfsr <- flash_lfsr(fl.pn.samp)
fl.ash.sampler <- flash_sampler(Y, fl.ash$fits$OHL, fixed="loadings")
fl.ash.samp <- fl.ash.sampler(nsamp)
fl.ash.lfsr <- flash_lfsr(fl.ash.samp)
m.lfsr <- t(get_lfsr(m$m))
for (j in 1:length(ts)) {
fl.pn.signif <- fl.pn.lfsr <= ts[j]
fl.pn.pr <- partition_by_type(fl.pn.signif)
fl.ash.signif <- fl.ash.lfsr <= ts[j]
fl.ash.pr <- partition_by_type(fl.ash.signif)
m.signif <- m.lfsr <= ts[j]
m.pr <- partition_by_type(m.signif)
pr.pn[j,] <- pr.pn[j,] + fl.pn.pr
pr.ash[j,] <- pr.ash[j,] + fl.ash.pr
pr.m[j,] <- pr.m[j,] + m.pr
}
}
mses <- mses / nsims
pr.pn <- pr.pn / nsims
pr.ash <- pr.ash / nsims
pr.m <- pr.m / nsims
saveRDS(mses, "./output/sims2mse.rds")
saveRDS(pr.pn, "./output/sims2prpn.rds")
saveRDS(pr.ash, "./output/sims2prash.rds")
saveRDS(pr.m, "./output/sims2prm.rds")</code></pre>
</div>
<div id="session-information" class="section level2">
<h2>Session information</h2>
<pre class="r"><code>sessionInfo()</code></pre>
<pre><code>R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] workflowr_1.0.1 Rcpp_0.12.17 digest_0.6.15
[4] rprojroot_1.3-2 R.methodsS3_1.7.1 backports_1.1.2
[7] git2r_0.21.0 magrittr_1.5 evaluate_0.10.1
[10] stringi_1.1.6 whisker_0.3-2 R.oo_1.21.0
[13] R.utils_2.6.0 rmarkdown_1.8 tools_3.4.3
[16] stringr_1.3.0 yaml_2.1.17 compiler_3.4.3
[19] htmltools_0.3.6 knitr_1.20 </code></pre>
</div>
<!-- Adjust MathJax settings so that all math formulae are shown using
TeX fonts only; see
http://docs.mathjax.org/en/latest/configuration.html. This will make
the presentation more consistent at the cost of the webpage sometimes
taking slightly longer to load. Note that this only works because the
footer is added to webpages before the MathJax javascript. -->
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
"HTML-CSS": { availableFonts: ["TeX"] }
});
</script>
<hr>
<p>
This reproducible <a href="http://rmarkdown.rstudio.com">R Markdown</a>
analysis was created with
<a href="https://github.com/jdblischak/workflowr">workflowr</a> 1.0.1
</p>
<hr>
</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://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>