« June 2014 | Main | August 2014 »
R is a very powerful language for creating custom data visualizations, but during the development process sometimes you make a mistake and things go horribly wrong. But sometime serendipity intervenes, and the (unintended) result can be quite beautiful. Accidental aRt, if you will. Curated by Kara Woo and Erika Mudrak, this fantastic Tumblr captures beautiful but unintended examples from R and other data visualization tools. Here are a couple of recent favourites (click for the originals):
I could definitely see some of these hanging on my wall. (Hmm ... I have an idea for this impressive-looking gadget.) Check out the full gallery at the link below.
Tumblr: Accidental aRt
by Andrie de Vries
In my previous post I wrote about how to identify and visualize package dependencies. Within hours, Duncan Murdoch (member of R-core) identified some discrepancies between my list of dependencies and the visualisation. Since then, I fixed the dispecrancies. In this blog post I attempt to clarify the issues involved in listing package dependencies.
In miniCRAN I expose two functions that provides information about dependencies:
Take a look at some examples. I illustrate with the the package chron, because chron neatly illustrates the different roles of Imports, Suggests and Enhances:
The function pkgDep() exposes not only these dependencies, but also also all recursive dependencies. In other words, it answers the question which packages need to be installed to satsify all dependencies of dependencies.
This means that the algorithm is as follows:
The resulting list of packages should then contain the complete list necessary to satisfy all dependencies. In code:
> library(miniCRAN)
> tags <- "chron" > pkgDep(tags, suggests=FALSE, enhances=FALSE, includeBasePkgs = TRUE) [1] "chron" "graphics" "stats"
> pkgDep(tags, suggests = TRUE, enhances=FALSE) [1] "chron" "RColorBrewer" "dichromat" "munsell" "plyr" "labeling" [7] "colorspace" "Rcpp" "digest" "gtable" "reshape2" "scales" [13] "proto" "MASS" "stringr" "ggplot2"
> pkgDep(tags, suggests = TRUE, enhances=TRUE) [1] "chron" "RColorBrewer" "dichromat" "munsell" "plyr" "labeling" [7] "colorspace" "Rcpp" "digest" "gtable" "reshape2" "scales" [13] "proto" "MASS" "stringr" "lattice" "ggplot2" "zoo"
Created by Pretty R at inside-R.org
To create an igraph plot of the dependencies, you can use the function makeDepGraph() and plot the results:
Created by Pretty R at inside-R.org
Note how the dependencies expand to zoo (enhanced), scales and ggplot (suggested) and then recursively from there to get all the Imports and LinkingTo dependencies.
In my previous post I tried to plot the most popular package tags on StackOverflow. Using the updated functionality in the miniCRAN functions, it is easier to understand the structure of the dependencies:
> tags <- c("ggplot2", "data.table", "plyr", "knitr", + "shiny", "xts", "lattice") > pkgDep(tags, suggests = TRUE, enhances=FALSE) [1] "ggplot2" "data.table" "plyr" "knitr" "shiny" "xts" [7] "lattice" "digest" "gtable" "reshape2" "scales" "proto" [13] "MASS" "Rcpp" "stringr" "RColorBrewer" "dichromat" "munsell" [19] "labeling" "colorspace" "evaluate" "formatR" "highr" "markdown" [25] "mime" "httpuv" "caTools" "RJSONIO" "xtable" "htmltools" [31] "bitops" "zoo" "SparseM" "survival" "Formula" "latticeExtra" [37] "cluster" "maps" "sp" "foreign" "mvtnorm" "TH.data" [43] "sandwich" "nlme" "Matrix" "bit" "codetools" "iterators" [49] "timeDate" "quadprog" "Hmisc" "BH" "quantreg" "mapproj" [55] "hexbin" "maptools" "multcomp" "testthat" "mgcv" "chron" [61] "reshape" "fastmatch" "bit64" "abind" "foreach" "doMC" [67] "itertools" "testit" "rgl" "XML" "RCurl" "Cairo" [73] "timeSeries" "tseries" "its" "fts" "tis" "KernSmooth" > set.seed(1) > plot(makeDepGraph(tags, includeBasePkgs=FALSE, suggests=TRUE, enhances=TRUE), + legendPosEdge = c(-1, -1), legendPosVertex = c(1, -1), vertex.size=10, cex=0.5)
Created by Pretty R at inside-R.org
After my previous post, Duncan Murdoch pointed out that the package rgl, suggested by knitr, appeared in the list, but not in the plot. This new version of the function fixes this bug, which was introduced because I retrieved the suggests dependencies incorrectly:
EDIT:
A few hours ago the miniCRAN went live on CRAN. Find miniCRAN at http://cran.r-project.org/web/packages/miniCRAN/index.html
As global warming causes sea levels to rise, the risk of flooding for coastal settlements also rises over time. A recent analysis by Reuters find that incidents of coastal flooding along the Eastern seaboard of the United States have surged in recent years as the sea level steadily rises.
Flood levels have been exceeded in six eastern cities by an average of 20 days or more since 2001. The analysis is based on more than 25 million hourly tide-gauge readings compared to flood thresholds set by the National Oceanic and Atmospheric Administration (NOAA).
As noted by reporter Ryan McNeill on Twitter, the analysis was done in R, after using the Ruby's RinRuby library to download the data from NOAA:
@dancow @znmeb I used Ruby to download and process the data, then I fed it into R using this gem: https://t.co/GGgnU5d04b
— Ryan McNeill (@McNeill_Tweets) July 10, 2014
Reuters plans to publish a broad examination of rising sea levels later this year. You can read the July 10 article at the link below.
Reuters: Exclusive: Coastal flooding has surged in U.S., Reuters finds (via Sharon Machlis)
I was honoured to be invited earlier this month to the Directions of Statistical Computing meeting in Brixen, Italy. DSC is one of two meetings run by the R Project and unlike the useR! conference, DSC is a much smaller and intimate meeting (DSC 2014 had about 30 participants). If you haven't come across DSC meeting before (quite possible, given that it had last been held in 2009), R Core Group member Martyn Plummer has a nice overview of DSC.
A focus of the first day of the conference was on the performance of R computation engine. The organizers invited representatives from all of the "alternative" R engine implementations, and I believe it marked the first time that developers involved with pqR, Renjin, FastR, and Riposte and TERR were gathered in the same place. (The CXXR project was unfortunately not represented.) Jan Vitek [slides] presented a fascinating comparison of the various projects, based on his interviews with the developers.
It was interesting to see the commonalities in many of the approaches. Three projects, Renjin [slides], FastR [slides] and Riposte [slides] use just-in-time compilation and an optimized bytecode engine. All have achieved impressive performance gains, but have struggled with compatibility (and especially being able to run the 6000+ CRAN packages). But it's clear that their work is having an influence on R itself: Thomas Kalibera [slides] (who previously worked on the FastR project) is working with Luke Tierney and Jan Vitek to improve the performance of R's bytecode interpreter.
Other approaches are also being pursued to improve the performance of the R engine. Luke Tierney [slides] described new improvements in R 3.1 to streamline the reference counting system, and noted that several of the performance improvements implemented by Radford Neal [slides] in pqR have already been incorporated into the R engine. And Helena Kotthaus [slides] has done some very exciting work to profile the performance of the R engine which has already led to performance improvements when virtual memory is being used.
Overall, it was exciting to see collaboration and research into R as a language, and especially the attention from the computer science community to the implementation of R. As Robert Gentleman (co-creator of R and conference lead) noted, R now has a new community beyond statisticians and data scientists: computer scientists. It's exciting to see how R is incorporating learning and innovation from this new community.
For more on DSC 2014, see the reports from Martyn Plummer on Day 1 and Day 2 of the conference. The full program, with links to download the slide presentation, is at the link below.
DSC 2014: Schedule (and slide downloads)
by Joseph Rickert
If I had to pick just one application to be the “killer app” for the digital computer I would probably choose Agent Based Modeling (ABM). Imagine creating a world populated with hundreds, or even thousands of agents, interacting with each other and with the environment according to their own simple rules. What kinds of patterns and behaviors would emerge if you just let the simulation run? Could you guess a set of rules that would mimic some part of the real world? This dream is probably much older than the digital computer, but according to Jan Thiele’s brief account of the history of ABMs that begins his recent paper, R Marries NetLogo: Introduction to the RNetLogo Package in the Journal of Statistical Software, academic work with ABMs didn’t really take off until the late 1990s.
Now, people are using ABMs for serious studies in economics, sociology, ecology, socio-psychology, anthropology, marketing and many other fields. No less of a complexity scientist than Doyne Farmer (of Dynamic Systems and Prediction Company fame) has argued in Nature for using ABMs to model the complexity of the US economy, and has published on using ABMs to drive investment models. in the following clip of a 2006 interview, Doyne talks about building ABMs to explain the role of subprime mortgages on the Housing Crisis. (Note that when asked about how one would calibrate such a model Doyne explains the need to collect massive amounts of data on individuals.)
Fortunately, the tools for building ABMs seem to be keeping pace with the ambition of the modelers. There are now dozens of platforms for building ABMs, and it is somewhat surprising that NetLogo, a tool with some whimsical terminology (e.g. agents are called turtles) that was designed for teaching children, has apparently become a defacto standard. NetLogo is Java based, has an intuitive GUI, ships with dozens of useful sample models, is easy to program, and is available under the GPL 2 license.
As you might expect, R is a perfect complement for NetLogo. Doing serious simulation work requires a considerable amount of statistics for calibrating models, designing experiments, performing sensitivity analyses, reducing data, exploring the results of simulation runs and much more. The recent JASS paper Facilitating Parameter Estimation and Sensitivity Analysis of Agent-Based Models: a Cookbook Using NetLogo and R by Thiele and his collaborators describe the R / NetLogo relationship in great detail and points to a decade’s worth of reading. But the real fun is that Thiele’s RNetLogo package lets you jump in and start analyzing NetLogo models in a matter of minutes.
Here is part of an extended example from Thiele's JSS paper that shows R interacting with the Fire model that ships with NetLogo. Using some very simple logic, Fire models the progress of a forest fire.
Snippet of NetLogo Code that drives the Fire model
to go if not any? turtles ;; either fires or embers [ stop ] ask fires [ ask neighbors4 with [pcolor = green] [ ignite ] set breed embers ] fade-embers tick end ;; creates the fire turtles to ignite ;; patch procedure sprout-fires 1 [ set color red ] set pcolor black set burned-trees burned-trees + 1 end
The general idea is that turtles represent the frontier of the fire run through a grid of randomly placed trees. Not shown in the above snippet is the logic that shows that the entire model is controlled by a single parameter representing the density of the trees.
This next bit of R code shows how to launch the Fire model from R, set the density parameter, and run the model.
# Launch RNetLogo and control an initial run of the # NetLogo Fire Model library(RNetLogo) nlDir <- "C:/Program Files (x86)/NetLogo 5.0.5" setwd(nlDir) nl.path <- getwd() NLStart(nl.path) model.path <- file.path("models", "Sample Models", "Earth Science","Fire.nlogo") NLLoadModel(file.path(nl.path, model.path)) NLCommand("set density 70") # set density value NLCommand("setup") # call the setup routine NLCommand("go") # launch the model from R
Here we see the Fire model running in the NetLogo GUI after it was launched from RStudio.
This next bit of code tracks the progression of the fire as a function of time (model "ticks"), returns results to R and plots them. The plot shows the non-linear behavior of the system.
# Investigate percentage of forest burned as simulation proceeds and plot library(ggplot2) NLCommand("set density 60") NLCommand("setup") burned <- NLDoReportWhile("any? turtles", "go", c("ticks", "(burned-trees / initial-trees) * 100"), as.data.frame = TRUE, df.col.names = c("tick", "percent.burned")) # Plot with ggplot2 p <- ggplot(burned,aes(x=tick,y=percent.burned)) p + geom_line() + ggtitle("Non-linear forest fire progression with density = 60")
As with many dynamical systems, the Fire model displays a phase transition. Setting the density lower than 55 will not result in the complete destruction of the forest, while setting density above 75 will very likely result in complete destruction. The following plot shows this behavior.
RNetLogo makes it very easy to programatically run multiple simulations and capture the results for analysis in R. The following two lines of code runs the Fire model twenty times for each value of density between 55 and 65, the region surrounding the pahse transition.
d <- seq(55, 65, 1) # vector of densities to examine res <- rep.sim(d, 20) # Run the simulation
The plot below shows the variability of the percent of trees burned as a function of density in the transition region.
My code to generate plots is available in the file: Download NelLogo_blog while all of the code from Thiele's JSS paper is available from the journal website.
Finally, here are a few more interesting links related to ABMs.
R is a functional language, which means that your code often contains a lot of ( parentheses ). And complex code often means nesting those parentheses together, which make code hard to read and understand. But there's a very handy R package — magrittr, by Stefan Milton Bache — which lets you transform nested function calls into a simple pipeline of operations that's easier to write and understand.
Hadley Wickham's dplyr package benefits from the %>% pipeline operator provided by magrittr. Hadley showed at useR! 2014 an example of a data transformation operation using traditional R function calls:
Here's the same code, but rather than nesting one function call inside the next, data is passed from one function to the next using the %>% operator:
You can read this version aloud to easily get a sense of what it does: the flights data frame is filtered (to remove missing values of the dep_delay variable), grouped by hours within days, the mean delay is calculated withn groups, and returns the mean delay for those hours with more than 10 flights.
You can use the %>% operator with standard R functions — and even your own functions — too. The rules are simple: the object on the left hand side is passed as the first argument to the function on the right hand side. So:
It's even possible to pass in data to something other than the first argument of the function using a . (dot) operator to mark the place where the object goes — see the magrittr vignette for details.
This new "pipelining" operation is a really useful addition to the R language, and R developers are starting to use it to make their code simpler to write and maintain. Hadley Wickham's newest R package, tidyr, makes it easy to clean up data sets for analysis by stringing together operations like "gather" and "spread" using the %>% operator.
And speaking of pipelining, you may have been wondering where the name "magrittr" comes from. Here's the answer:
The only other question is: will Stefan be making this coffee mug available?
magrittr vignette: Ceci n'est pas un pipe
Follow David on Twitter: @revodavid
