Writing a Knitr Engine to run code through an interpreter of your choice!

Preamble

Recently I took a job with the Department of Fisheries and Oceans translating some statistical software from python into R. The package uses kernel smoothing to estimate the probability distribution of a variable of interest relative to one or more covariates. I’m employed for the next month or so, and hopefully by the end we’ll be submitting the package to CRAN. I’ll write a little more about that as it gets closer to happening.

I’ve gotten to the point in the development process where I need to think hard about optimization and benchmarking. I’m in the habit of keeping my bosses/colleagues up to date on my work through html files I make using R markdown (the same way I write this blog). R markdown has drastically changed the way I work, allowing me to weave narrative and code together to make a more compelling exploration of my progress. I was hoping to benchmark my code against the original python code to see how they stack up speedwise and thought it would be great to have the python benchmarks right in there with my R code. I knew knitr provided facillities for including non-R code in documents so I set to work figuring our how to do it.

Back to basics

For those of you who don’t know the fantastic one-two punch that is R markdown and Knitr let me tell you a little bit about it. R, since at least the early 2000’s, has provided facillities for (s)weaving code and documents together. The aim was to mesh analysis and reporting into one cohesive action in the name of reproducibility. Though there is still much progress to be made moving toward this ideal the tools have gotten better.

In 2012, Yihui Xie introduced the world to knitr, a neater knitter of code and documentation. It has (almost?) completely supplanted sweave as the primary tool to acheive this goal. The benefits of knitr are many, it allows you to create html and pdfs, include code, figures, and LaTeX all in the same document, and rerun the analysis any time the included code changes (but only if you want it to). Needless to say I’m a big fan.

The next major innovation in the reproducible research game was R Markdown. Markdown is family of markup languages introduced by John Gruber in 2004 designed to facillitate writing html documents in plain text. The goal was to have readable plain text documents where the html tags could be generated on demand. R markdown is variant designed to work primarily with R, created by the folks at R Studio. If you use R Studio, the integration of R markdown is seamless, you write an R markdown formatted file, the documents is converted to a format knitr understands, and then knitr makes you a beautiful document with all of your work ready to be shared with the world. I highly recommend using these tools for just about any authoring job short of manuscript preparation.

Now that I’ve given you a brief introduction, lets dive into my oddessy making a knitr engine.

The problem

As I mentioned in the preamble, I wanted to run R and python code side-by-side in my document so that both sets of benchmarks could be readily compared. Knitr comes pre-loaded with a quite a few engines for other languages:

library(knitr)
names(knit_engines$get())

##  [1] "awk"       "bash"      "coffee"    "gawk"      "haskell"  
##  [6] "node"      "perl"      "python"    "Rscript"   "ruby"     
## [11] "sas"       "scala"     "sed"       "sh"        "zsh"      
## [16] "highlight" "Rcpp"      "tikz"      "dot"       "c"        
## [21] "fortran"   "asy"       "cat"       "asis"

A quick scan of the available engines shows that knitr should be perfectly capable of executing python code. So I tried:

```{r, engine = “python”}
#My python code
```

Which sadly did not work as I had hoped. The python interpreter used by knitr in that situation was the stock python interpreter. It failed to properly import the modules I needed and was untenable for my objectives. “Alas, if only there was an ipython engine” I thought to myself, and then promptly decided I should probably build one.

Taking apart the engine

To figure out how to build my own engine, I decided to have a look at the engine provided for me. A convenience of R that I think more people should take advantage of is the ability to get code for any function you like. Provided it’s written in R in many cases you can figure out how it works. You can skim over this code, I’ll get back to the important bits.

knit_engines$get("python")

## function (options) 
## {
##     engine = options$engine
##     code = if (engine %in% c("highlight", "Rscript", "sas", "haskell")) {
##         f = basename(tempfile(engine, ".", switch(engine, sas = ".sas", 
##             Rscript = ".R", ".txt")))
##         writeLines(c(switch(engine, sas = "OPTIONS NONUMBER NODATE PAGESIZE = MAX FORMCHAR = '|----|+|---+=|-/<>*' FORMDLIM=' ';", 
##             haskell = ":set +m"), options$code), f)
##         on.exit(unlink(f))
##         switch(engine, sas = {
##             saslst = sub("[.]sas$", ".lst", f)
##             on.exit(unlink(c(saslst, sub("[.]sas$", ".log", f))), 
##                 add = TRUE)
##             f
##         }, haskell = paste("-e", shQuote(paste(":script", f))), 
##             f)
##     }
##     else paste(switch(engine, bash = "-c", coffee = "-e", node = "-e", 
##         perl = "-e", python = "-c", ruby = "-e", scala = "-e", 
##         sh = "-c", zsh = "-c", NULL), shQuote(paste(options$code, 
##         collapse = "\n")))
##     code = if (engine %in% c("awk", "gawk", "sed", "sas")) 
##         paste(code, options$engine.opts)
##     else paste(options$engine.opts, code)
##     cmd = options$engine.path %n% engine
##     out = if (options$eval) {
##         message("running: ", cmd, " ", code)
##         system2(cmd, code, stdout = TRUE, stderr = TRUE)
##     }
##     else ""
##     if (!options$error && !is.null(attr(out, "status"))) 
##         stop(paste(out, collapse = "\n"))
##     if (options$eval && engine == "sas" && file.exists(saslst)) 
##         out = c(readLines(saslst), out)
##     engine_output(options, options$code, out)
## }
## <environment: namespace:knitr>

First good sign, all the code is R with no references to byte-code and it’s short enough to not take all day to figure out. First thing to note about the code is it’s a function, and not only that it’s a function of one argument, options. The first line of code takes a named element of options called engine (so options is a named list containing things like the name of the engine). My intuition was that options was a list of every component of chunk, and the rest of the function largely confirmed that intuition. Other named elements of options used in the function include:

• code
• engine.opts
  
• engine.path
• error

The shape of the code the begins to come into focus. The function primarily takes your code out of the chunk (element code), converts it into shell command that runs engine with some flags and a quoted argument containing your code. This gets executed by the function system2 which is used to execute shell commands and retrieve output. This output is then passed to a function engine_output along with the code and options, the results of which are used to build your document. So from this point it seamed obvious what needed to happen.

1. Concatenate engine.path and engine to create the command that needs to be run
2. Convert code into a string appropriate for use in shell commands
3. Paste flags onto code so that the engine knows to expect a string of code
4. Run the command and the modified code through the engine with system2, saving the output
5. Pass everything to engine_output

Pretty straight forward really, and except for some simple housekeeping that’s all I needed to write.

Writing the new engine

Now I began writing my own engine (and it didn’t work the first time, some of the housekeeping I didn’t figure out right away, but lets pretend it did). I called my engine ipythonKnitEngine but you can call your engine whatever you want.

ipythonKnitEngine <- function(options){
  
  enginePath <- options$engine.path
  if(grepl("/$", enginePath)) enginePath <- paste0(enginePath, "/")
  
  engine <- paste0(enginePath, options$engine)
  
  code <- paste("-c", shQuote(paste(options$code,
                                    collapse = "\n")))
  
  out <- system2(engine, code, stdout = TRUE, stderr = TRUE)
  
  engine_output(options, options$code, out)
  } 

As of writing this, the engine requires you specify the path to your ipython executable, its a kink I’ll probably iron out in the near future, but this was a quick fix to finish a document, so it still needs some polish.

For those who aren’t familiar shQuote converts a character vector to a string useable in your shell. The other potentially obtuse bit is the collapse = "\n" to paste. When paste takes a character vector of length n, by default it makes a vector of length n, when collapse is not NULL, all of those elements are concatenated together separated by the collapse argument. Otherwise pretty standard stuff.

The last necessary little bit of funny business is engine_output isn’t actually available from the knitr package to users, but you can force your engine to use knitr’s namespace by setting your engine’s environment:

environment(ipythonKnitEngine) <- environment(knit_engines$set)

Now the environment of the function is set to the knitr namespace and it will be able to find engine_output.

And finally we can set the new engine and begin using it!

knit_engines$set(ipython = ipythonKnitEngine)

Using the engine

To use the engine in your R Markdown / Knitr documents just run library(knitr) in your first code chunk, and include the function declaration for your engine along with the environment and knit_engines$set statements above and then you can use your engine in your document.

Now we can run

```{r, engine = “ipython”, engine.path = “path/to/ipython/”}
#Python code to execute
```
with impunity and it works like a charm!.

Wrap-Up

I hope you liked this little tutorial. If you like to say “nuts to writing my own, I’ll just use yours” you can find all this code in a convenient format on github.

Chris

99th ESA Meeting in Sacramento

Coming to you live from Sacramento California, it’s Datamancy on the Road!

I’m posting this with extreme pre-talk jitters, I will be presenting tomorrow morning in conference room 314 at the Sacramento conference centre. I’m hoping for the best.
If you’re reading this after the fact and have come to have another look at my slides, thanks for your interest! I really love what I do and enjoy sharing it just as much, skip to the bottom to get to the slides.

Conference Recap

Environmental Sequencing Workshop

The ESA 2014 annual meeting has been a whirlwind of wonderful science and people. I was lucky enough to get to sit in on a wonderful workshop lead by Holy Bik, Kelley Thomas, Dorota Porazinska, and Simon Creer on environmental sequencing approaches. It was a blast to hear from such passionate people in such an interesting field.
I went in to the session hoping to learn a little more about the state of the field in order to hopefully answer a irksome problem my supervisor has been grappling with. Our model species (Splachnum ampullaceum and Splachnum pensylvanicum) are poorly studied genetically so efficient markers don’t exist for detecting their spores in the environment. I wanted to tap into the thoughts of people who deal with this kind of problem daily.
I’ve been inspired to revist the prior work that has been done on developing markers for our species, and perhaps employ some more modern techniques for identifying them. So kudos to the leaders for putting the workshop together.

Talks, Talks, Posters, and More Talks

I really enjoyed getting a whirlwind tour of what’s going on in ecology research. I focused mostly on catching community ecology and modelling talks. I saw many other young scientists put themselves on display for what must have been an intimidating crowd.

I really enjoyed talks from:

Aaron David on microbe innoculation sources for new plants. Interesting work in teasing out whether new plant roots gain their fungal endophytes from the parent plant, conspecific neighbours, or the surrounding soil.

Yue “Max” Li from the Chesson lab, looked at coexistence using a long-standing field experiment.
Jiaqi Tan looked at coexistence mechanisms amongst protists. He used a cool strategy to identify dietary strategy amongst the protist, then used dietary preference as a proxy for niche to estimate niche differentiation (amongst other things).

Dan Scholes from Ken Paige’s lab in Illinois looking at the genetic basis of plant overcompensation in response to herbivory. I had the pleasure of having dinner with Dan and Ken, fascinating work that reminded me of the Arabidopsis days of my undergrad.

Those are just a few of my highlights, I’d continue but it is exhausting looking for personal websites for young scientists.

Slides!

Examining Coexistence between Splachnaceae Mosses with Individual-Based Modelling from Chris Hammill

Regular Expressions in R: A Primer

Preamble

Hi readers,
Sorry for the delay in posting, I’ve been hard at work getting my talk ready for the ESA annual meeting in Sacramento next week. I’ll definitely post a recap and link to my slides after the conference.

I’m on a workacation on Lake of Bays in Ontario (my home province) at the moment, enjoying the lake-side life on my breaks.

I’m continuing work on my Build Your Bar project, dabbling in natural language processing. I find myself using regular expressions (regex) frequently enough to think it warrants a post. The R help on Regex ?regex is long and technical and makes no attempt to convince you of it’s importance. So here is a quick primer on using regex in R, the focus will be on perl-like regular expressions.

History

Regular Expressions have been a part of computing since the early days of home computing. When the primary input-output interface to a computer was text, an efficient means of finding the exact text you were interested in was necessary. “Simple” shorthands for complicated series of characters were developed. The idea for regular expressions originated in the work of Stephen Kleene(1956) on regular sets. Just over a decade later, Ken Thompson imported the ideas into the Unix operating system via the text editing software QED. For a more detailed history of regex I suggest reading Staffan Noteberg’s brief review.

R’s Native Regex

R supports two, or imprecisely three, regex schemes. The two true regex schemes are Extended which follows the POSIX 1003.2 standard and Perl-Like my personal favourite that adapts the incredibly powerful regex syntax of the perl programming language. The last scheme is Fixed which is less of a regex style than a lack thereof, fixed mode will take characters absolutely literally, except for R’s own escape characters.

By default R uses extended regex, the unfortunate consequence is if you choose to follow my lead and use perl-style regex, you have to be a bit more verbose in your code and specify perl = TRUE each time you use a regex function¹. I however feel it’s a small price to pay for getting to use perl-like regex.

More about Perl-Like

The help file for regex offers the advice to head to pcre web page, but you will likely want to go straight to the source: http://perldoc.perl.org/perlre.html. The perl documentation page on the current state of their regex syntax is a handy reference for anyone who plans on using perl-like regex.

Unfortunately for us R users, many of perl’s really handy regex features aren’t available to us. Native perl regex allow for concurrent search and accession, meaning you can search for multiple tokens within a string and access the returned tokens without storing the matches in an intermediary object first. More simply, in perl, searching and editing a string is essentially the same action. In R, searching and accession are handled by completely seperate functions, and regular expressions aren’t tokenized to allow multiple operations on the same search criterion.

These differences are due to the fact that perl was designed for text manipulation, it is self-described as the “swiss army knife” of scripting languages allowing you to functionality from multiple languages with perl as the middle-man. The lingua franca of all programming languages is the source code, which perl can elegantly cut, trim, and glue together. R was built for statistics, and so string manipulation is more of an acquired skill. However, let it never be said you can’t do powerful string manipulation in R, it’s just not true. String manipulation in R is clunky, maybe, but powerful enough to do anything you need.

Finding That Special String

I love string manipulation because it appeals to my love of puzzle solving. Finding just the right regex to get what you’re looking for really gives me the thrill of the hunt.

• So you want to find a number in the middle of a word, but only after an m and only in the last word of string: REGEX
• So you want to trim whitespace from the end of the string but only if there are 3 or more spaces: REGEX
• Want to find brand names with between one and three capitalized words in a list of drink ingredients REGEX

Hopefully your interest is piqued and you now want to learn all about using regex like a wizard. If so, great, because I’m going to show you how! But first, a key R idiosyncracy that needs to be addressed.

R Idiosyncracy

R by default recognizes backslashes “\” as the beginning of an escape sequence, also called a metacharacter. Perl also uses a backslash to indicate a metacharacter so there is a clash. In native perl “\d” is a metacharacter meaning all digits, but when R sees “\d” it thinks you’re invoking the metacharacter in an inappropriate way. To get the behaviour you’d expect, you need to escape the behaviour of the backslash to have it pass the un-escaped regex. You do this by adding an extra backslash “\\d” behaves exactly as you’d expect “\d” to.
This situation gets particularly comical when you want to use a literal backslash in a regular expression. In native perl “\\” gets you a literal backslash, but in R you need to escape BOTH backslashes to get the desired behaviour “\\\\” gets you a literal backslash (as I am typing this in R markdown which by default follows R escaping rules, I just typed 8 backslashes to show you four, ugh).

First Steps (Wild-Cards and Quantifiers)

The first things you’re going to need to learn to become a Perl-Like regex master is to use wild-cards and quantifiers. You may have run into wild-cards in any number of different places, google used to support them for example. Wild-cards are regex speak for “just give me anything”.

Unless you’ve worked with regex before you’ve probably never worked closely with quantifiers. Quantifiers are regex speak for “give me some number of these”. If you’ve seen the wild-cards “*” and “?” from DOS, “?” is a wild-card meaning any one character, “*” is a wild-card with an indefinite quantifier meaning give me any number of any character.

Character	Meaning
.	Match any 1 character
*	Match ZERO or more of the previous character
+	Match ONE or more of the previous character
?	Match ZERO or ONE of the previous character
{N, M}	Match atleast N but fewer than M of the previous character Leaving M empty ({N,}) means N or more

With these tools at your disposal you can begin your odyssey into the world of pattern matching. Let’s see some examples

Regex	Sample Matches
"Al.*"	"Alfred", "Albert", "Alphonso", "Alex likes going to the movies"
"5? ?Apples"	"5 Apples", "5Apples", " Apples", "Apples"
"(Na )+Batman!"	"Na Na Na Na Na Na Na Na Batman!" but not "Batman!"

As you may have noticed, the special quantifiers (“*“,”+“, and”?“) are all just convenient shorthands for (”{0,}“,”{1,}“, and”{0,1}“), you’ll find you use the special quantifiers more frequently.

Brackets, Groups, and Alternation

I’ve hinted at this next group of behaviours in the previous set of examples. The “(Na )+” I used to match the batman theme song is an example of grouping. Using parantheses groups a series of characters such that they can be quantified like a single character. The above regex states that I want one or more sets of the sequence “Na”.

Suppose you weren’t satisfied that you found all the ways people transcribe the batman theme, for example you have a friend use the phoneme “Da” to represent the notes in the song. Never fear, bracketed character classes to the rescue. You can allow the regular expression to choose one of several characters. Bracketed character classes are denoted with square braces “[]”.

Lastly, maybe the vowel is all wrong too. Maybe you think people might be using phonemes like “Duh” or “Nah”. You can specify these specific combinations within a group, and let the regex match any group with an alternation. Alternations are specified with “|”.

Regex	Sample Matches
"(Na )+Batman!"	"Na Na Na Na Na Na Na Na Batman!"
"([ND]a )+Batman!"	"Na Na Na Na Na Na Na Na Batman!" and "Da Na Na Na Da Na Na Na Batman!"
"[ND]((a)|(uh)|(ah))"	"Duh", "Na", "Dah"

Now you should be able to get just about any set of characters, and in any amount as you could possibly want.

Positioning and Quantifier Modification

Another useful tool is the ability to specificy where in a string you’d like to match. Two handy positioning modifiers are allowed in a regular experession. You can specify the beginning of a string “^” and the end of a string “$”

Regex	Sample Matches
"^Hi"	"Hi there reader" but not "You there reader, Hi"
"Goodbye$"	"Dear reader, Goodbye" but not "Goodbye dear reader"

Quantifier modification is a topic I’ve touched on in a previous post where I demonstrated how to process html to extract tag values. The question mark has a double meaning depending on where it is found in a regular expression. Normally it is a treated as a quantifier meaning {0,1}, however when it appears immediately after another quantifier, it instructs that quantifier to behave non-greedily.

Greedy is a reference to a quantifier’s habit of trying to capture as many characters as possible. If given a choice between matching one or seven characters, it will choose seven. We can modify this behaviour to just take one.
Take for example the body of a very simple html document:

<body>
<p> Generic things on the <b>Internet<b> have grown <b> tiresome</b></p>
</body>

If I was curious about all the things the author thought worthy of bolding, I might first try the regex
"<b>.+</b>"

But I would get one match <b>Internet<b> have grown <b> tiresome</b>
My indefinite quantifier would match all the way to the very last </b> in the document, obviously not what I wanted. So instead I should tell the “+” quantifier to not be so greedy
"<b>.+?</b>"

Then it would correctly match "<b>Internet</b>" and "<b>Tiresome</b>".

Special Characters

The last type of metacharacter I’d like to show you are the special character classes, useful in bracketed character classes and on there own or combined with quantifiers. Some useful character class metacharacters are

Character	Meaning
\w	any word character - all numbers, letters, and underscores
\W	non-word character - anything except the above
\d	any single digit number
\s	any space character (regular space or tab)

And with that you should be able to match just about any pattern in a string that your heart desires, here are a few interesting examples.
Stringr comes pre-loaded with a method to trim whitespace from a string, with regex we can do that by saying match some number of whitespaces at the beginning or end of the string like so:

gsub(pattern = \"(^ +)|( +$)\", replacement = "", string)²

Or as I mentioned earlier, maybe you want to find brand names in a set of ingredients, the ingredient set was comma delimited, and brands of note were tagged with ® the html code for the tm symbol. Brand names also have each word capitalized and tend to be between 1-4 words and may include hyphens. To pull out each brand in the ingredient list you could do this:

matchPull(pattern = "([A-Z][\\w-']*? ){0,3}[A-Z][\\w-']+?®")³

The one thing you may not recognize is the bracketed character class "[A-Z]" which means exactly what it looks, all the capitals from A through Z.

Outro

I hope from this post you’ve learned a little bit about the application of regular expressions in R and can see some opportunities to use them in your own work. The ability to precisely and effectively specify exactly the string you’re interested makes performing any text mining tasks orders of magnitude easier. Also if you’re just interested in organizing and curating data sets for thesis work for example the ability to accurately edit information stored in strings is a massive advantage over doing manual edits in software such as excel.

As always, thanks for reading!

-Chris

---

1. I recently discovered the wonderful package stringr another of Hadley Wickham’s packages. The syntax for specifying perl-like mode with stringr functions is different, you just wrap your regex string in a call to the perl function like so: stringr_foo("string", pattern = perl("pattern")). Stringr is a dependency of several of Wickham’s other packages so you likely already have it installed.↩
2. The function gsub is R’s equivalent of find and replace all↩
3. The matchPull function is my convenience function for finding and extracting a matched text from a string. Similar in idea to str_extract from stringr, but I wrote matchPull before I knew about stringr, and I’m too lazy to learn someone else’s convenience function when I have my own. Plus as a hold over from my days coding in java I think underscores in function names are yucky. Code for matchPull is available in my post here and also comes in my friendlyShiny package.↩

Advanced ggplot2: How to add overlays, underlays, and make multi-axis plots

Intro

Hi Readers,

Previously I mentioned I’m working on a contract for the Newfoundland government. For this contract I’ve had to make a number of figures for a report that’s coming out in the near future on the state of Caribou in the province. When I began the contract I was new to ggplot2, a nice graphics package by R guru Hadley Wickham which has become de rigueur in the circles I find myself in. Now I use it all the time when I want to make high quality figures.

For its many virtues there are drawbacks that come with producing figures with ggplot2. The primary difficulty is the syntax is inobvious without having read Wickham’s book (or the Wilkinson book that inspired it). Additionally there is an unfortunate dearth of documentation on the finer points of package.
Like many attempts to simplify tricky topics, the simplification adds its own baggage that confounds its initial objective. However, ggplot2 is immensely powerful and for most applications it is the best option available, I recommend it, but be sure to devote some time to learning it.

In this post I want to walk you through how to make a watermarked plot. Sometimes a bit of background information can keep a viewer in tune with the message of the figure.

A Cut to the Chase

One key frustration I experienced was the
overly-pedanticomission of dual axis plotting. Granted there are many work arounds to this problem and with a little digging you’re likely to find Wickham’s work around, which exposes the tender underbelly of the package. This workaround showed me just enough of the inner workings of the package to be able to meet my client’s demands for a three-panelled plot with a watermark (things that aren’t available through the surface user interface of ggplot2.) I’d like to show you how you to do it.

The Setup

First lets get the packages we need installed and loaded:

install.packages(c("ggplot2", "grid", "gtable"))

library(ggplot2)
library(grid)
library(gtable)

The second and third packages are probably new to you so I’ll spend a paragraph to explain what they’re all about.

Wickham’s ggplot2 package was built on the pre-existing infrastructure of grid graphics. Grid is a low level plotting interface that suspends the user above the gory details of specifying an image tuple¹ by tuple, but not so far as to lose any power. Grid-like plotting may be familiar to you if you have ever programmed graphics from scratch in a lower level language than R.

The grid objects (grobs) produced by ggplot2 are hopelessly complicated and innacessible for direct manipulation. However, the package gtable exists to help us bridge the divide. So we use gtable to convert a ggplot2 image into a table of grobs, from this table we extract just the peices we want so that we can manually add them to our plots later using the grid package directly.

Back to the task at hand. We’re going to need some data for our plots. Lets make a quartic and look at its first three derivatives. We’ll use the original function as a watermark to compare the derivatives

x <- -10:10

y4 <- x^4 + 3*x^3 + 5*x^2 - 10*x
y3 <- 4*x^3 + 9*x^2 + 10*x - 10
y2 <- 12*x^2 + 18*x + 10
y1 <- 24*x + 18

polyData <- data.frame(x, y4, y3, y2, y1)

Getting the Watermark

To pull out our watermark line, we first have to put it in a plot

quartPlot <- ggplot(data = polyData, aes(x = x, y = y4)) + 

                    geom_line(colour = "grey80", alpha = .75, size = 1)

Now we build our gtable from quartPlot. The structure of a gtable is hopelessly complicated so ignore this if you’re not in the mood for some head scratching. What I didn’t borrow straight from Wickham I gleaned through agonizing trial and error.

By searching through the “grob” element of your gtable for the element that matches up with the “name” element “panel”, we can find the plotting panel where our line lives. The panel grob has 3 children, some gridlines, backgrounds, and our line. The line turned out to be child number 2.

grobTab <- ggplot_gtable(ggplot_build(quartPlot))
quartLine <- grobTab$grob[[which(grobTab$layout$name == "panel")]]$children[[2]]

So now we have our watermark line ready to be added to our future plots.

The Other Plots

The other plots are going to need special themes to ensure the watermark is visible. Themes are how ggplot2 controls the final look of a plot. By specifying the background of both the panel and plot windows as “element_blank()” you ensure that your plots to cover your watermark (this took me far too long to figure out, I had the watermark superimposed but highly transparent because I had given up on having it in the background.) I also removed the grid lines because they’re obtrusive².

basePlot <- ggplot(data = polyData, aes(x = x)) +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        panel.grid = element_blank())

cubicPlot <- basePlot + geom_line(aes(y = y3), size = 1.5) + ylab("Cubic") + 
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank())

quadPlot <- basePlot + geom_line(aes(y = y2), size = 1.5) + ylab("Quadratic") + 
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank())

linePlot <- basePlot + geom_line(aes(y = y1), size = 1.5) + ylab("Line")

Note that for the quadratic and cubic plots I also removed all traces of the x axis, so the plots can visually share the x-axis (although in reality they do not).

Setting up Viewports

A useful grid-ism to understand for working with ggplot2 is the concept of a viewport. Viewports are regions of the screen where things can be plotted. They’re specified as proportions of the screen for their width and height, with an x and y coordinates for the middle of the viewport. I could explain further but it will be quicker to show them in action

vp1 <- viewport(width = 1, height = .33, y = 1/6, x = 1/2)
vp2 <- viewport(width = 1, height = .33, y = 3/6, x = 1/2)
vp3 <- viewport(width = 1, height = .33, y = 5/6, x = 1/2)

This code creates three viewports each occupying a third of the screen vertically.
We also need a viewport for our watermark but this takes some trial and error to fit well

quartLine$vp <- viewport(width = .8, height = .91, x = .55, y = .55)

The Final Product

grid.newpage()
grid.draw(quartLine)

print(linePlot, vp = vp1)
print(quadPlot, vp = vp2)
print(cubicPlot, vp = vp3)

And voila, there you have a plot with a watermark to show off, there are lots of ways to customize these points, for my contract I included axis lines like so

grid.newpage()
grid.draw(quartLine)

print(linePlot + theme(axis.line = element_line(colour = "black")), vp = vp1)
print(quadPlot + theme(axis.line = element_line(colour = "black")), vp = vp2)
print(cubicPlot + theme(axis.line = element_line(colour = "black")), vp = vp3)

For highly technical work it is important to note that the panels will not necessarily be in alignment. For example, because the longest axis label on the cubic plot is a digit longer than the other plots it’s drawing panel (where the plotting happens) is one digit narrower.

This can probably be circumvented (though I’ve never tried) by manually specifying the size of the label box in your themes.

Outro

So there is how you can get around some of the limitations of ggplot2 to make really nifty figures for your own work and analyses.

If you’re really intrigued, maybe investigate grid more deeply and shed the overhead of ggplot2, or you can delve deeper into how to combine the two to make a tricked-out personalized ggplot2 workflow to make some astonishing figures.

I hope you enjoyed this little journey through the innards of ggplot2.

-Chris

---

1. If you’re not familiar with the term tuple it’s a very useful term for vectors of coordinates or quantities. I wish I had learned it earlier in life. From high-school onward we’re used to specifying cartesian coordinates like (x,y). This bracketed construct of coordinates form a couple that represents a location in 2D space. When drawing an image we might want to tell our computers to draw a point at (x,y) and give it colour c, we can specify this as a triple (x,y,c). Adding two more quantities would give us a quintuple. By now you’ve noticed the pattern. This extends to as many quantities as you’d like, the genotype of an individual at many loci could be treated as an n-tuple.↩
2. I haven’t gone back to check yet, but I think if you set all the panel elements to “element_blank()” in the plot we used for the watermark you can ensure that you only have one child in your panel. This would make life easier because you wouldn’t have to search for which child is the line, you could just take element one every time.↩

Introducing FriendlyShiny: A package to make interactivity with R easy

I currently am working on a contract for the Government of Newfoundland and Labradour’s Department of Environment & Conservation, and they have me making graphics for an upcoming report. While working on these graphics I repeatedly find myself fixing graphic parameters to try to improve the aesthetic appeal. It’s a tedious process including a lot more “guess-and-check” than I’d really like to admit.

If you’ve read a few of my other posts you’d have caught that I recently learned shiny in order to bring you the Drink Name Generator. I thought to myself surely there must be a way to leverage the power of shiny to perform parameter fine tuning in a way that doesn’t create more overhead than it’s worth.

By the time I had that thought I was already lost. I spent most of the past two days putting together a small package that gives a (in my opinion) intuitive interface to shiny. The package allows for applets to be created from a pre-existing R code with only a few minor changes.

Without further ado, the README from my GitHub

FriendlyShiny is my attempt at making the wonderful reactive code abilities of R’s shiny package more accessible to novice users, and folks who want interactive code quickly. FriendlyShiny provides a simple syntax for specifying reactive elements in code chunk without the overhead of designing the user interface and coding the applet by hand.

To allow interactivity for a code chunk, it just needs to be wrapped in an interact function call.

Consider trying to teach a math class about line specifications. You’ve taught them a line can be specified by y = mx + b but you’d like to show them how the line changes with its parameters.

You could write demostration code like this to plot a line, and then you could bring the plot up in front of your students

  slope <- 1
  intercept <- 0
  
  plot(0,
       ylim = c(0, 15),
       xlim = c(0, 10),
       xlab = "x",
       ylab = "y",
       col = 0)
  
  abline(intercept, slope)

Now you could go back to your code and change the values of slope and intercept by hand, OR you can make it interactive very simply like so:

interact(
{
  slope <- sI("slope", 1, min = 1/3, max = 3)
  intercept <- nI("intercept", 0)

  plot(0,
       ylim = c(0, 15),
       xlim = c(0, 10),
       xlab = "x",
       ylab = "y",
       col = 0)
  
  abline(intercept, slope)
  
}, outputType = "plot")

All you need to do is wrap your code in curly braces and a call to interact, then you need to specify you want a slider to help you choose the slope, and a numerical input box to specify the intercept (you probably want a slider for that too, but I’m sure you can see how to fix it).

The last argument after the curly braces tells the interact function whether to render your output as text or as a plot. Text output is useful for fine tuning parameters for modelling parameters and plot output is useful for getting graphing parameters juuust right.

Supported Widgets

Currently friendlyShiny supports:

• SliderInput: sI(name, start, min, max, step = NULL)
• NumericInput: nI(name, start)
• TextInput: tI(name, start = "")
• SelectInput (i.e. Drop down box): dI(name, type, start, ...) where … are your choices (comma delimited) and type can be “character” or “numeric”
• RadioButtons: rI(name, type, start, …) exactly like dI but a different choice mode
• CheckboxInput: cI(name, start = FALSE) toggles a logical between true and false

How to install

If you’re on here you probably know how to do this better than I do, but you can install this package to play with and improve (if you’re up to the task) by running the following commands in R:

install.packages("devtools") #if you don't already have devtools installed
library(devtools)

install_git("cfhammill/friendlyShiny")

Outro

I hope you like my package! I sincerely hope someone takes the reigns on this project, anyone can take my code and turn it into something amazing. I don’t intend to polish it much more than this, I just thought shiny should be accessible to everyone. If you notice any glaring errors please let me know.

-Chris

Build Your Bar Project: What's In A Drink Name / Drink Name Generator

This post is a short-attention-span version of my previous post, plus a little gadget I wrote to generate drink names.

Recap

I have been looking at a data set containing names and recipes for ~17000 cocktails. The original goal of the project was to use publicly available recipe data to inform a decision of what drinks to buy for a small home bar. The project since spun out of control because I realized the wealth of interesting trivia in the data.

How Do People Name Their Cocktails

Conceptually it might be worthwhile to separate out drinks into several naming classes based on their popularity. There are popular drinks that can be found at nearly any bar worth its salt, there are specialty drinks which are often variants of a popular drink with ingredient tweaks, and there are unique drinks that are unlikely to be found at any bar and will require a recipe if you want one.

Each class is generally named differently. Popular drinks have readily recognizable names. Specialty drinks often have parts of a popular drink name with additional words indicating the variations involved. The unique drinks do not follow a pattern as closely as the other two, reflecting the heterogeneity of people creating these unique drinks.

In order to learn a bit about how words are distributed within drink names I collated all the drink names in my data set and performed some simple analyses. After clean-up, I made a word cloud containing the top 100 words used in drink names, sized according to their frequency. This analysis does not consider the popularity of the drink itself, but the data does contain some indication of popularity (e.g. more popular drinks are more likely to have multiple recipes with the same name).

Booze Cloud!

 

Click to embignify

An interesting follow up to this analysis would be to re-weight the word frequencies according to how popular the drinks they are found in are. This will make the data look more like what you might see on a cocktail menu.

Drink Name Generator

Part of the data exploration I calculated the frequencies of the top 500 words, as well as the distribution of the number of words in a cocktail name. The largest name in my data set was 8 words, with two being the most frequent.

I constrained my drink name generator to produce drink names between 2 and 5 words. The generator generates a random number between 0 and 1, it compares that number to a look-up table corresponding to the cumulative distribution function (cdf) of word number frequencies. That number of words is generated using the same method but from the cdf of word name frequencies.
The generator then returns your new drink name. The app was written using shiny, a way to write small apps using R, which I will probably talk about in a future post.

Here’s the generator!

Make sure to try a few, some are pretty funny, a few gems so far are "Monkey Nut Punch" and "Flaming Panty Sweat".
Post your best in the comments!

Click here to open the app in a new window

Chris

Build Your Bar Project: Synthesis and Exploratory Data Analysis

Hello readers, this post will be a good one I promise. I’ve started using R Markdown which seems like it will greatly increase the speed with which I can give you analyses.

Diving Back In

We left off last time after having downloaded a ton of html files containing drink recipe data. The first thing to do is to have a look inside the html file. This step is critical for looking at how to pull out the information that we want

Exerpts

Found the drink name

<title> 73 Bus #2 recipe</title> 

Found some key words

<meta content="73 bus #2, 73, bus, #2, gin,<br/>
  triple sec, lime juice, cranberry juice, drink recipe, drink, recipe,<br/>
  alcoholic drink recipe, cocktail recipe, cocktail, mixed drink, martini"<br/> 
  name="keywords">

Found the hierarchy (drink class)

<div class="pm" style="margin-top:20px;"><a href="/cat/1/
  ">Cocktails</a>
  > <a href="/cat/14/">Short drinks</a> 
  > <a href="/cat/141/">by base-ingredient</a> 
  > <a href="/cat/40/">gin-based</a></div>

Pulling out the information

So at this point I’ve identified seven variables I’d like to track for each drink

1. Drink name
2. Drink hierarchy(class)
3. Ingredients
4. Keywords
5. Number of ratings
6. Average rating (out of 10)
7. URL

In order to get at these values I need to design a regular expression that will only capture the tag of interest. I discussed regular expressions briefly in my previous post, we will rely heavily on the non-greedy quantifier “.+?” I discussed here.

Data Extraction

In order to pull the data out I used R (surprise surprise). I wrote two small accessory functions. One to make string manipulation easier, and one to remove html tags.

#Use either regexpr(default) or grexepr to match elements of interest
#Extract and return them using regmatches
matchPull <- function(pattern, text, invert = FALSE, global = FALSE, ...){
  if(global){
    match <- gregexpr(pattern, text, ...)
  } else {
    match <- regexpr(pattern, text, ...)
  }
  
  pulled <- regmatches(text, match, invert)
  if(length(pulled) == 0) pulled <- NA
  
  pulled
}

# Remove html tags, note the use of the .*? quantifier, 
# a cousin of .+? that can match 0 characters
# Where .+? matches 1 of more.
stripTags <- function(text){
  gsub("<.*?>", "", text, perl = TRUE)
}

First step lets bring in one of the many html files we downloaded and try to extract all of the important data. We’ll use the examples noted above to practice.

fileCon <- file(siteName, blocking = FALSE)
site <- paste0(readLines(con = fileCon), collapse = "\n")
close(fileCon)

# Pull out the whole title, remove the tages
# and remove the word recipe which flanks each recipe name
name <-  matchPull("<title>.*?</title>", site,
                   ignore.case = TRUE, perl = TRUE)
name <- stripTags(name)
name <- sub(" recipe$", "", name, perl = TRUE)
  
# Pull out just the meta tag with the name keywords,
# the pull out the contents, and remove quotes
keywords <- matchPull("<meta content=.*?name=\"keywords\"", site, 
                      ignore.case = TRUE, perl = TRUE)
keywords <- matchPull("\".*?\"", keywords, perl = TRUE)
keywords <- gsub("\"", "", keywords)

# Pull out the division of class "pm" style "yadda-yadda"
# and remove all tags
hierarchy <- matchPull("<div class=\"pm\" style=\"margin-top:20px;\">.*?</div>", 
                       site, ignore.case = TRUE, perl = TRUE)
hierarchy <- stripTags(hierarchy)

Once we’ve figured out how to get all the useful data out of one file, we can encase it in a function that returns one row of data, and apply that to all the files we downloaded (after testing it on a much smaller subset). After that we’ll have a data frame containing all the juicy data, which is much easier to work with. Suppose we’ve encased our processing in {processSite <- function(siteName)} we can apply it many sites all at once by wrapping it in another function

processSites <- function(siteList){
  frameSeed <- processSite(siteList[1])
  drinkFrame <- frameSeed[rep(1, length(siteList)),]
  
  #This sapply structure is basically just a for loop
  sapply(1:length(siteList), function(i){
    drinkFrame[i,] <<- processSite(siteList[i]) 
  })
  
  drinkFrame
}

fileNames <- list.files()#Make sure you've set your working directory before this
fileNames <- fileNames[grepl("\\.html", fileNames)] #grab just .html files

#Try it out on the first 6 files
practiceNames <- head(fileNames)
practiceData <- processSites(practiceNames)

#After inspecting practice data for quality, process them all
drinkData <- processSites(fileNames)

Processing ~17000 sites took R around 6 minutes (wow) on my computer, producing a 10Mb data.frame, which I saved so that I never have to run this code again. Now that the data is in, we can begin with the fun parts. I’ll some skip the quality control steps to get right to the meaty stuff.

Exploratory Data Analysis (EDA)

The cornerstone of any data related project is poking and proding the data to figure out what’s in there. Make some histograms, correlation matrices, and any other simple data visualizations you think might be informative. This step is probably my favourite because it

1. Helps recognize general patterns
2. Identifies issues with data quality
3. It’s fun to watch your data set begin to tell its first story

I’ll present one EDA that I thought was fun. I began getting interested in what people liked to name their drinks, surely there would be some cool patterns in that. I decided I wanted to try my hand at making a word cloud with some of the most common words in drink names. To make a word cloud using R I used the wordcloud package and the tm (text mining) package.

library(wordcloud)
library(tm)

load("drinkData.rda") #Bring in our drinkData from the last step

# Filter out drinks that are neither cocktails nor shots 
# by looking in their hierarchy
drinksFrame <- drinksTable[grepl("(cocktails)|(shots)", 
                                 drinksTable$hierarchy, 
                                 ignore.case = TRUE, 
                                 perl = TRUE),]

# I had a few cases of multiple duplicates, 
# this loop keeps tacks on Alt to duplicated names
# Repeats until no duplicates are found
while(anyDuplicated(drinksFrame$name) != 0){
  dupeNamed <- duplicated(drinksFrame$name)
  drinksFrame$name[dupeNamed] <- 
    paste(drinksFrame$name[dupeNamed],"Alt",sep=" ")
}

#Clean up now empty dupeNamed vector
rm(dupeNamed)

#Convert from character vector to one long string of words
nameVector <- paste0(tolower(drinksFrame$name), collapse = " ") 

# Use tm's built in functions to remove stopwords, see below for a note,
# Also remove alt (because I put it there)
# As well as punctuation. I removed numbers because the site named 
# duplicates with sequential numbers
# And "2" was one of the most popular words 
nameVectorCleaned <- removeWords(nameVector, c(stopwords("english"), "alt"))
nameVectorCleaned <- removePunctuation(nameVectorCleaned)
nameVectorCleaned <- removeNumbers(nameVectorCleaned)

# Split the cleaned string back into a vector of words (tokens really)
# Separated by white space
nameVector <- unlist(strsplit(nameVectorCleaned, " "))

#Then use table to count instances of each word
#Remove number one which was an empty string
#(an unfortunate consequence of our splitting algorithm)
nameFreqs <- table(nameVector)[-1]
nameWords <- as.character(names(nameFreqs))
namesFreqs <- as.numeric(nameFreqs)

freqOrder <- order(nameFreqs, decreasing = TRUE) #Create an ordering vector
top100 <- head(freqOrder, 100) #Indices of the 100 most popular words

#Make a 10 inch by 10 inch pdf to hold the wordcloud  
pdf("boozeNameCloud.pdf", width = 10, height = 10)

#Plot the words, ordered and sized by frequency
wordcloud(nameWords[top100], namesFreqs[top100], 
          scale = c(12,1), random.order = FALSE)

#Close up shop and admire our work
dev.off()

I re-ran the code to make the wordcloud pdf multiple times, because there is something stochastic in the word placement. After a few tries the words aligned and suited my aesthetic tastes. And so I give you

Booze Cloud!

Click to embignify