<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/reproducible-research-when-your-results-cant-be-reproduced/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://blog.appsilondatascience.com/rstats/2017/03/28/reproducible-research-when-your-results-cant-be-reproduced.html"> Appsilon Data Science - R language</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Even most sophisticated machine learning methods, most beautiful visualisations and perfect datasets can be useless if you do your research carelessly and ignore the context of your code execution. In this article, you will read about a few situations that can destroy your reproducibility and learn how to resolve and avoid them in the future.</p>
<h4>3 danger zones</h4>
<p>There are three places where things can go wrong if you don’t pay attention:</p>
<ol>
<li>R session context</li>
<li>Operating System (OS) context</li>
<li>Data versioning</li>
</ol>
<p>In this blogpost, I will focus on the first two. They are quite similar because both are related to software context of a research execution. While more and more data scientists pay attention to what happens in their R session, not many think about OS context, which can also significantly influence their research.</p>
<p>Data versioning is a huge topic itself and I will cover it in a separate article in the future. For now, let’s see some examples where things can go wrong, when you have different R session context or various Operating System configurations.</p>
<h4>R session context</h4>
<p>The basic elements of an R session context, that data scientists should be aware of are:</p>
<ul>
<li>R version</li>
<li>Packages versions</li>
<li>Using <code>set.seed()</code> for a reproducible randomisation</li>
<li>Floating point accuracy</li>
</ul>
<p><strong>R version</strong></p>
<p>Not many monitor <a target="_blank" rel="nofollow" href="https://cran.r-project.org/src/base/NEWS">this R language change log</a> on a regular basis, but sometimes a crucial changes are made and can lead to a failure in your code. One of such situation is described <a target="_blank" rel="nofollow" href="http://stackoverflow.com/questions/32234120/different-results-for-auto-arima-with-d-1-in-r-version-3-2-0-and-r-version-3-2-2/32234317#32234317">here</a>. Change in R(3.2.1) for “nested arima model has higher log-likelihood”, led to inconsistency of results for arima calculations with d >= 1.</p>
<p><strong>Packages versions</strong></p>
<p>Now, let’s see an example with a changing package version. I am going to use <code>dplyr</code> in versions <code>0.4.0</code> and <code>0.5.0</code>:</p>
<p><strong><em>dplyr 0.4.0</em></strong></p>
<figure>
<pre><code data-lang="r"><span>df</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>1</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>1</span><span>:</span><span>5</span><span>)</span><span>
</span><span>result</span><span> </span><span>&lt;-</span><span> </span><span>df</span><span> </span><span>%>%</span><span> </span><span>dplyr</span><span>::</span><span>distinct</span><span>(</span><span>x</span><span>)</span><span> </span><span>%>%</span><span> </span><span>sum</span><span>
</span><span>#</span><span> </span><span>result</span><span> </span><span>==</span><span> </span><span>3</span></code></pre>
</figure>
<p><strong><em>dplyr 0.5.0</em></strong></p>
<figure>
<pre><code data-lang="r"><span>df</span><span> </span><span>&lt;-</span><span> </span><span>data_frame</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>1</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>2</span><span>,</span><span> </span><span>2</span><span>),</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>1</span><span>:</span><span>5</span><span>)</span><span>
</span><span>result</span><span> </span><span>&lt;-</span><span> </span><span>df</span><span> </span><span>%>%</span><span> </span><span>dplyr</span><span>::</span><span>distinct</span><span>(</span><span>x</span><span>)</span><span> </span><span>%>%</span><span> </span><span>sum</span><span>
</span><span>#</span><span> </span><span>result</span><span> </span><span>==</span><span> </span><span>8</span></code></pre>
</figure>
<p><img alt="" src="https://i0.wp.com/blog.appsilondatascience.com/assets/article_images/2017-03-28-reproducible-research/boom.png?w=456" data-recalc-dims="1"> whoa! Version 0.5.0 introduced functionality breaking changes and result is now completely different.</p>
<p><strong>set.seed()</strong></p>
<p>When your code uses a randomisation and you want to have reproducible results, it is important to remember about setting a seed.</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span><span> </span><span>5</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>6</span><span> </span><span>7</span><span> </span><span>2</span><span> </span><span>5</span><span> </span><span>9</span><span>
</span><span>></span><span> </span><span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span><span> </span><span>5</span><span>)</span><span> </span><span># Different result
</span><span>[</span><span>1</span><span>]</span><span> </span><span>5</span><span> </span><span>7</span><span> </span><span>8</span><span> </span><span>2</span><span> </span><span>3</span></code></pre>
</figure>
<p>…and with setting seed:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>set.seed</span><span>(</span><span>42</span><span>)</span><span>
</span><span>></span><span> </span><span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span><span> </span><span>5</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>10</span><span> </span><span>9</span><span> </span><span>3</span><span> </span><span>6</span><span> </span><span>4</span><span>
</span><span>></span><span> </span><span>set.seed</span><span>(</span><span>42</span><span>)</span><span>
</span><span>></span><span> </span><span>sample</span><span>(</span><span>1</span><span>:</span><span>10</span><span>,</span><span> </span><span>5</span><span>)</span><span> </span><span># Result is the same
</span><span>[</span><span>1</span><span>]</span><span> </span><span>10</span><span> </span><span>9</span><span> </span><span>3</span><span> </span><span>6</span><span> </span><span>4</span></code></pre>
</figure>
<p>Random numbers returned by generator are in fact a sequence of numbers that can be reproduced. Setting a seed informs a generator where to start.</p>
<p><strong>Floating point accuracy</strong></p>
<p>Pay attention how expressions on floating point numbers are implemented. Let’s assume you want to implement an algorithm described in a recently published scientific paper. You get unexpected results in some cases and now try to debug it. And you eventually get this:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>x</span><span> </span><span>&lt;-</span><span> </span><span>3</span><span>
</span><span>></span><span> </span><span>y</span><span> </span><span>&lt;-</span><span> </span><span>2.9</span><span>
</span><span>></span><span> </span><span>x</span><span> </span><span>-</span><span> </span><span>y</span><span> </span><span>==</span><span> </span><span>0.1</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>FALSE</span></code></pre>
</figure>
<p><img alt="" src="https://i0.wp.com/blog.appsilondatascience.com/assets/article_images/2017-03-28-reproducible-research/boom.png?w=456" data-recalc-dims="1"> what? It can be solved by testing for “near equality”:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>x</span><span> </span><span>&lt;-</span><span> </span><span>3</span><span>
</span><span>></span><span> </span><span>y</span><span> </span><span>&lt;-</span><span> </span><span>2.9</span><span>
</span><span>></span><span> </span><span>all.equal</span><span>(</span><span>x</span><span> </span><span>-</span><span> </span><span>y</span><span>,</span><span> </span><span>0.1</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>TRUE</span></code></pre>
</figure>
<p>More about a floating point arithmetic can be <a target="_blank" rel="nofollow" href="https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html">found here</a>.</p>
<h4>Operating system context</h4>
<p>When your research is shared with others and you can’t predict which operating system someone is going to use, many things can go wrong. Any difference in one of the following elements:</p>
<ul>
<li>System packages versions</li>
<li>System locale</li>
<li>Environment variables</li>
</ul>
<p>can result in different execution results.</p>
<p>Many R packages depend on system packages. You can view it in a package DESCRIPTION file under the parameter <code>SystemRequirements</code>. Let’s look at the example with <code>png</code> package and running it on two different operating systems.</p>
<p>Package <code>png</code> DESCRIPTION file:</p>
<figure>
<pre><code data-lang="bash">Package: png
Version: 0.1-7
...
SystemRequirements: libpng
...</code></pre>
</figure>
<p>Let’s assume we wrote custom function <code>savePlotToPNG</code> that saves a plot to a file. We want to write an unit test for it and the idea is to compare produced result with a previously generated reference image. We generate reference plots on Ubuntu 16.04:</p>
<figure>
<pre><code data-lang="r"><span>reference.png.path</span><span> </span><span>&lt;-</span><span> </span><span>file.path</span><span>(</span><span>&quot;referencePlots&quot;</span><span>,</span><span> </span><span>&quot;ref0001.png&quot;</span><span>)</span><span>
</span><span>savePlotToPNG</span><span>(</span><span>reference.png.path</span><span>)</span></code></pre>
</figure>
<p>…and then we run tests on Debian 8 Jessie. We use exactly the same version of <code>png</code> package:</p>
<figure>
<pre><code data-lang="r"><span>reference.png.path</span><span> </span><span>&lt;-</span><span> </span><span>file.path</span><span>(</span><span>&quot;referencePlots&quot;</span><span>,</span><span> </span><span>&quot;ref0001.png&quot;</span><span>)</span><span>
</span><span>test.plot.path</span><span> </span><span>&lt;-</span><span> </span><span>tempfile</span><span>(</span><span>fileext</span><span> </span><span>=</span><span> </span><span>&quot;.png&quot;</span><span>)</span><span>
</span><span>savePlotToPNG</span><span>(</span><span>test.plot.path</span><span>)</span><span>
</span><span>identical</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>png</span><span>::</span><span>readPNG</span><span>(</span><span>test.plot.path</span><span>),</span><span>
</span><span>y</span><span> </span><span>=</span><span> </span><span>png</span><span>::</span><span>readPNG</span><span>(</span><span>reference.png.path</span><span>))</span><span>
</span><span>#</span><span> </span><span>Error</span><span>!</span><span> </span><span>libpng</span><span> </span><span>system</span><span> </span><span>package</span><span> </span><span>is</span><span> </span><span>different</span><span> </span><span>on</span><span> </span><span>Ubuntu</span><span> </span><span>16.04</span></code></pre>
</figure>
<p><img alt="" src="https://i0.wp.com/blog.appsilondatascience.com/assets/article_images/2017-03-28-reproducible-research/boom.png?w=456" data-recalc-dims="1"> …plots don’t match! This is caused by different <code>libpng</code> system package versions.</p>
<p><strong>Locale and system variables</strong></p>
<p>Ubuntu:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>x</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>&quot;-110&quot;</span><span>,</span><span> </span><span>&quot;-1&quot;</span><span>,</span><span> </span><span>&quot;-10&quot;</span><span>,</span><span> </span><span>&quot;10&quot;</span><span>,</span><span> </span><span>&quot;0&quot;</span><span>,</span><span> </span><span>&quot;110&quot;</span><span>,</span><span> </span><span>&quot;1&quot;</span><span>)</span><span>
</span><span>></span><span> </span><span>sort</span><span>(</span><span>x</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>&quot;0&quot;</span><span> </span><span>&quot;1&quot;</span><span> </span><span>&quot;-1&quot;</span><span> </span><span>&quot;10&quot;</span><span> </span><span>&quot;-10&quot;</span><span> </span><span>&quot;110&quot;</span><span> </span><span>&quot;-110&quot;</span></code></pre>
</figure>
<p>Windows:</p>
<figure>
<pre><code data-lang="r"><span>></span><span> </span><span>x</span><span> </span><span>&lt;-</span><span> </span><span>c</span><span>(</span><span>&quot;-110&quot;</span><span>,</span><span> </span><span>&quot;-1&quot;</span><span>,</span><span> </span><span>&quot;-10&quot;</span><span>,</span><span> </span><span>&quot;10&quot;</span><span>,</span><span> </span><span>&quot;0&quot;</span><span>,</span><span> </span><span>&quot;110&quot;</span><span>,</span><span> </span><span>&quot;1&quot;</span><span>)</span><span>
</span><span>></span><span> </span><span>sort</span><span>(</span><span>x</span><span>)</span><span>
</span><span>[</span><span>1</span><span>]</span><span> </span><span>&quot;-1&quot;</span><span> </span><span>&quot;-10&quot;</span><span> </span><span>&quot;-110&quot;</span><span> </span><span>&quot;0&quot;</span><span> </span><span>&quot;1&quot;</span><span> </span><span>&quot;10&quot;</span><span> </span><span>&quot;110&quot;</span></code></pre>
</figure>
<p><img alt="" src="https://i0.wp.com/blog.appsilondatascience.com/assets/article_images/2017-03-28-reproducible-research/boom.png?w=456" data-recalc-dims="1"> …okayyy <img alt="🙂" class="wp-smiley" src="https://s.w.org/images/core/emoji/2/72x72/1f642.png"> String sorting in R is based on the locale which is different for Windows and Linux systems. Read more about setting locales <a target="_blank" rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/base/html/locales.html">here</a>.</p>
<h4>Tools for R session level reproducibility</h4>
<p>Data scientist should be aware of the traps described above. These are the foundations of the reproducible research. Below there is a list of tools that support managing packages versions:</p>
<p><strong><a target="_blank" rel="nofollow" href="https://rstudio.github.io/packrat/">Packrat</a></strong> – popular package that allows to create and manage a private library attached to your project.</p>
<p><strong><a target="_blank" rel="nofollow" href="https://github.com/gmbecker/switchr">Switchr</a></strong> – handy abstraction over <code>.libPaths()</code>, <code>lib.loc</code> and other lower level library settings. Helpful in creating and managing libraries.</p>
<p><strong><a target="_blank" rel="nofollow" href="https://mran.microsoft.com/documents/rro/reproducibility/">Checkpoint</a></strong> – helpful when you need to find older versions of a specific package.</p>
<h4>Tools for OS level reproducibility</h4>
<p>The best approach is to do your research in an isolated environment, that can be shared as an image. Inside of this environment you can also use any tool for R session level reproducibility, like Packrat, Switchr or Checkpoint.</p>
<p>To create an isolated environment you can use virtualization or conteinarization. There are many tools that can be used for that: <a target="_blank" rel="nofollow" href="https://www.vagrantup.com/">Vagrant</a>, <a target="_blank" rel="nofollow" href="https://www.docker.com/">Docker</a>, <a target="_blank" rel="nofollow" href="https://coreos.com/rkt/">rkt</a> and many more. I recommend using Docker because it’s easy to use and has a good community support.</p>
<p>Example structure of your <code>~/project</code>:</p>
<figure>
<pre><code data-lang="r"><span>Dockerfile</span><span>
</span><span>install_libraries.R</span><span>
</span><span>R</span><span>/</span><span> </span><span>&lt;--</span><span> </span><span>contains</span><span> </span><span>your</span><span> </span><span>research</span><span> </span><span>R</span><span> </span><span>scripts</span></code></pre>
</figure>
<p>Let’s look at a sample <code>Dockerfile</code> that describes the steps necessary for recreation of an isolated environment:</p>
<figure>
<pre><code data-lang="bash">FROM rocker/rstudio:3.3.2 

MAINTAINER Paweł Przytuła <span>&quot;join@appsilondatascience.com&quot;</span>

RUN <span>echo</span> <span>&quot;en_US.UTF-8 UTF-8&quot;</span> >> /etc/locale.gen <span>\</span>
  <span>&amp;&amp;</span> locale-gen en_US.utf8 <span>\</span>
  <span>&amp;&amp;</span> /usr/sbin/update-locale <span>LANG</span><span>=</span>en_US.UTF-8

ENV LC_ALL en_US.UTF-8
ENV LANG en_US.UTF-8

<span># Set the CRAN mirror</span>
RUN <span>echo</span> <span>&quot;local({</span><span>\n</span><span> r &lt;- getOption(</span><span>\&quot;</span><span>repos</span><span>\&quot;</span><span>)</span><span>\n\</span><span>
r[</span><span>\&quot;</span><span>CRAN</span><span>\&quot;</span><span>] &lt;- </span><span>\</span><span>
</span><span>\&quot;</span><span>https://cran.r-project.org</span><span>\&quot;\n\</span><span>
options(repos = r)</span><span>\n\</span><span>
})</span><span>\n</span><span>&quot;</span> >> /etc/R/Rprofile.site

COPY install_libraries.R /code/install_libraries.R

RUN R -f /code/install_libraries.R

CMD <span>[</span><span>&quot;/init&quot;</span><span>]</span></code></pre>
</figure>
<p><code>install_libraries.R</code> contains code that reproduces an R environment. It can be either some Packrat code or something like this:</p>
<figure>
<pre><code data-lang="r"><span>install.packages</span><span>(</span><span>&quot;devtools&quot;</span><span>)</span><span>
</span><span>devtools</span><span>::</span><span>install_version</span><span>(</span><span>&quot;dplyr&quot;</span><span>,</span><span> </span><span>version</span><span> </span><span>=</span><span> </span><span>&quot;0.5.0&quot;</span><span>)</span><span>
</span><span>devtools</span><span>::</span><span>install_version</span><span>(</span><span>&quot;stringr&quot;</span><span>,</span><span> </span><span>version</span><span> </span><span>=</span><span> </span><span>&quot;1.1.0&quot;</span><span>)</span><span>
</span><span>devtools</span><span>::</span><span>install_version</span><span>(</span><span>&quot;forecast&quot;</span><span>,</span><span> </span><span>version</span><span> </span><span>=</span><span> </span><span>&quot;7.3&quot;</span><span>)</span></code></pre>
</figure>
<p>After that it’s enough to run a following command to build an image:</p>
<figure>
<pre><code data-lang="bash">docker build -t myimagename:1.0.0 .</code></pre>
</figure>
<p>When an image is ready you should either save it as a file using <code>docker save -o myimagename.img myimagename:1.0.0</code> or upload it to <a target="_blank" rel="nofollow" href="https://hub.docker.com/">Docker Hub</a> – <a target="_blank" rel="nofollow" href="https://docs.docker.com/engine/getstarted/step_six/">steps described here</a>.</p>
<p>Using a previously saved Docker image you can start a container with the following:</p>
<figure>
<pre><code data-lang="bash">docker run -d -p 8787:8787 -v ~/project/R:/R myimagename:1.0.0</code></pre>
</figure>
<p>Because image is based on <code>rocker/rstudio</code>, you can access an RStudio session via <code>http://localhost:8787</code> (login <code>rstudio</code>, password <code>rstudio</code>). Your <strong>R/</strong> directory will be available inside the running container via <code>/R</code> path. Awesome, isn’t it?!</p>
<p>With a Docker image like this, you have a control over whole software context needed for a reproducible research. You can easily extend it with more packages and share it without a fear of losing reproducibility. I recommend this approach in any data science project, which results must be reproduced on more than one machine and also when working in a bigger data science team.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://blog.appsilondatascience.com/rstats/2017/03/28/reproducible-research-when-your-results-cant-be-reproduced.html"> Appsilon Data Science - R language</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/YoMs1BVFB50"></a></div>