Reproducible research

A minimal standard for data analysis and other scientific computations is the reproducibility. That is, the method, code and data are assembled in a way so that another group can recreate all of the results.

Markup and Markdown

The essential idea in a markup language is that it consists of ordinary text, plus signs which indicate how to change the formatting or meaning of the text. Markdown is a lightweight markup language with plain text formatting syntax. It is intended to be as easy-to-read and easy-to-write as is feasible. Markdown is designed so that it can be converted to HTML and dozens of output formats, like PDFs, Word files, slideshows, and more.

R Markdown

R Markdown is a file format for making dynamic documents with R. An R Markdown document is written in Markdown and contains chunks of embedded R code. R Markdown provides an unified authoring framework for data science, combining your code, its results, and your prose commentary.

https://kieranhealy.org/resources/

Prerequisites

You need knitr and rmarkdown packages, but you don’t need to explicitly install them or load them, as R Studio automatically does both when needed.

Paragraph Breaks and Forced Line Breaks

To insert a break between paragraphs, include a single completely blank line.

To force a line break, put two blank spaces at the end of a line.

To insert a break between paragraphs, include a single completely blank line.

To force a line break, put *two* blank
spaces at the end of a line.

Title, Author, Date, Output Format, Table of Contents

You caan also use the header to tell R Markdown whether you want it to render to HTML (the default), PDF, or something else. See http://rmarkdown.rstudio.com/lesson-9.html for more about the output format.

---
title: Your Title
author: Your Name
date: Today
output:
 pdf_format
  toc: yes
  toc_float: no
---

Code Chunks and Their Results

A code chunk is simply an off-set piece of code by itself. It is preceded by ```{r} on a line by itself, and ended by a line which just says ```. The code itself goes in between. Here, for instance, is some code which loads a data set from a library, and makes a scatter plot.

library(MASS)
data(cars)
plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)")

First, notice how the code is included, nicely formatted, in the document. Second, notice how the output of the code is also automatically included in the document. If your code outputs numbers or text, those can be included too:

with(cars, cor(speed, dist))

## [1] 0.8068949

library(ggplot2)

# What does the data look like?
DT::datatable(cars) # Interactive table

plot <- ggplot(cars, aes(x = speed, y = dist)) + geom_point() + #
  geom_smooth(method = "lm") + theme_bw() + # add linear model and make black and white
  xlab("Speed (mph)") + ylab("Stopping distance (ft)") # Add axis labels
plotly::ggplotly(plot) # Make the plot interactive

## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Inline Code

Code output can also be seamlessly incorporated into the text, using inline code. This is code not set off on a line by itself, but beginning with `r and ending with `. For example, the sample correlation coefficient between speed and dist is 0.8068949.

Changing Image Sizes and Alignments

There are a bunch of options for adjusting the placement of the figures which R produces. fig.align controls the horizontal alignment (left, right, or center). The options fig.height and fig.width let you specify the desired height or width of the figure.

plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)", las = 1)
lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")
title(main = "cars data")

plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)", las = 1, log = "xy")
title(main = "cars data (logarithmic scales)")
lines(lowess(cars$speed, cars$dist, f = 2/3, iter = 3), col = "red")

The figure size can be globally defined in the header of your R Markdown file.

---
output:
    fig_height: 9
    fig_width: 6
---

Tables

The default print-out of matrices, tables, etc. from R Markdown is frankly ugly.

fm1 <- lm(log(dist) ~ log(speed), data = cars)
coef(summary(fm1))

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -0.7296687  0.3758457 -1.941405 5.809204e-02
## log(speed)   1.6023912  0.1395376 11.483579 2.259356e-15

knitr::kable(summary(fm1)$coef, digits=4)

## Warning in summary(fm1)$coef: partial match of 'coef' to 'coefficients'

pander::panderOptions("digits", 4)
pander::pander(fm1)

## Warning in x$coeff: partial match of 'coeff' to 'coefficients'

Caching

If document rendering becomes time consuming due to long computations or plots that are expensive to generate, you can ask R Markdown to keep track of whether a chunk of code has changed, and only re-run it if it has. This is called caching the chunk. Let’s see the effect of caching by executing an inline function. The Fibonacci numbers (\(F_n\)) are defined by \(F_0 = 0\), \(F_1 = 1\), and \(F_n = F_{n-1} + F_{n-2}\) for \(n \ge 1\).

R_fibonacci = function(x){
  if (x == 0 || x == 1) return(x);
  return (R_fibonacci(x-1) + R_fibonacci(x-2))}

The below code will take several seconds to compute at the first run. On the second run, knitr will use the cached output without evaluating the code (unless you change it!)

R_fibonacci(30)

## [1] 832040

However, be careful with caching. One issue is that a chunk of code which hasn’t changed itself might call on results of earlier, modified chunks, and then we would want to re-run the downstream chunks. There are options for manually telling R Markdown “this chunk depends on this earlier chunk”, but it’s generally easier to let it take care of that, by setting the autodep=TRUE option.

1. If you load a package with the library() or require() commands, R Markdown isn’t smart enough to check whether the package has changed (or indeed been installed, if you were missing it). So that won’t trigger an automatic re-running of a cached code chunk.
2. To manually force re-running all code chunks, the easiest thing to do is to delete the directory R Markdown will create (named something like filename_cache) which it uses to store the state of all code chunks.

Rcpp Code Chunks

This section is optional. Those who have a C compiler installed can turn on the evaluation of three chunks below. If you don’t have the C compiler, please ignore this section. Sometimes R code just isn’t fast enough. You’ve used profiling to figure out where your bottlenecks are, and you’ve done everything you can in R, but your code still isn’t fast enough. You may want to consider the Rcpp package to accelerate computations, which permit to rewrite functions in C++.

#include <Rcpp.h>

// [[Rcpp::export]]
int Rcpp_fibonacci(const int x) {
    if (x == 0 || x == 1) return(x);
    return (Rcpp_fibonacci(x - 1)) + Rcpp_fibonacci(x - 2);
}

Rcpp_fibonacci(30L)

The improved performance by using the Rcpp function over the naive R function can be seen by benchmarking these two functions. To understand performance, you can also profile execution of two functions using Rprof or profvis.

library(microbenchmark)
op <- microbenchmark(R_fibonacci(10),Rcpp_fibonacci(10L),times=10L)
pander::pander(op)

If you want use an Rcpp function in R console, you must compile it or install it using a local R package. For more details, please see http://www.rcpp.org.

Setting Defaults for All Chunks

You can tell R to set some defaults to apply to all chunks where you don’t specifically over-ride them. Here are the ones I generally use:

# Need the knitr package to set chunk options
library(knitr)

# Set knitr options for knitting code into the report:
# Don't print out code (echo)
# Save results so that code blocks aren't re-run unless code changes (cache),
# or a relevant earlier code block changed (autodep), but don't re-run if the
# only thing that changed was the comments (cache.comments)
# Don't clutter R output with messages or warnings (message, warning)
# This will leave error messages showing up in the knitted report
opts_chunk$set(echo=TRUE,
               cache=FALSE, autodep=TRUE, cache.comments=FALSE,
               message=FALSE, warning=FALSE)

This sets some additional options beyond the ones I’ve discussed, like not re-running a chunk if only the comments have changed (cache.comments = FALSE), and leaving out messages and warnings. I would typically give this set-up chunk itself the option include=FALSE.

Individual chunk options will overide the global chuck options.

Tangle

You can extract the source code chunk out of the R Markdown document. This procedure is called “tangling”, which can be done by using knit() function with tangle=TRUE option in the R console. Then, you can execute the “tangled” code by source() to get the results without rendering the R Markdown document.

knitr::knit("yourdocument.rmd", tangle=TRUE)
source("yourdocument.R")

Prerequisites

If your output format is html, equations are displayed using the MathJax JavaScript library. Note that this library is loaded from the MathJax website so readers of your document must be online to see the rendered equations. On the other hand, if your output format is pdf, then you need to install a program called LaTeX. I recommend to install TeX Live for Linux or Window machines and MacTeX for OS X.

\begin{align*}
P(|X-\mu| > k) & = P(|X-\mu|^2 > k^2)\\
& \leq \frac{\mathbb{E}\left[|X-\mu|^2\right]}{k^2}\\
& \leq \frac{\mathrm{Var}[X]}{k^2}
\end{align*}

\[
\begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{pmatrix}
\]