Revolutions

Because it's Friday: Levitation with Sound Waves

No camera trickery here: engineers at the University of Tokyo have demonstrated how to levitate and move polystyrene particles using only sound waves (via Gizmodo).

 

The basic principle is fairly straightforward: pairs of speakers generate a standing wave (at ultrasonic frequencies) at a focus point, and the motion of the air molecules is enough to keep the particles suspended. With some clever engineering, by varying the phase of the sounds, the can move the focus point (and also the particles). Sadly, it's had to see how this could be scaled up to macro-sized objects though.

That's all for this week! See you back here on the blog on Monday.

Princeton's guide to linear modeling and logistic regression with R

If you're new to the R language but keen to get started with linear modeling or logistic regression in the language, take a look at this "Introduction to R" PDF, by Princeton's Germán Rodríguez. (There's also a browsable HTML version.)

In a crisp 35 pages it begins by taking you through the basics of R: simple objects, importing data, and graphics. Then, it works through several examples of linear models (formula basics, fitting a model, model diagnostics, analysis of variance and even regression spones). Finally, there's a section on Generalized Linear Models, with a focus on logistic regression. The document doesn't attempt to explain all of the capabilities of R, but instead works through a series of examples to teach by demonstration. All of the datasets used in the guide are available online, so it's easy to follow along from home.

Also by the same author: a guide to R for Stata users.

Germán Rodríguez: Introducing R

A First Look at rxDForest()

by Joseph RIckert

Last July, I blogged about rxDTree() the RevoScaleR function for building classification and regression trees on very large data sets. As I explaned then, this function is an implementation of the algorithm introduced by Ben-Haim and Yom-Tov in their 2010 paper that builds trees on histograms of data and not on the raw data itself. This algorithm is designed for parallel and distributed computing. Consequently, rxDTree() provides the best performance when it is running on a cluster: either an Microsoft HPC cluster or a Linux LSF cluster.

rxDForest() (new with Revolution R Enterprise 7.0) uses rxDTree() to take the next logical step and implement a random forest type algorithm for building both classification and regression forests. Each tree of the ensemble constructed by rxDForest() is built with a bootstrap sample that uses about 2/3 of the original data. The data not used in builting a particular tree is used to make predictions with that tree. Each point of the original data set is fed through all of the trees that were built without it. The decision forest prediction for that data point is the statistical mode of the individual tree predictions. (For classification problems the prediction is a majority vote, for regression problems the prediction is the mean of the predictions.) 

Only a couple of parameters need to be set to fit a decision forest. nTree specifies the number of trees to grow and  mTry spedifies the number of variables to sample as split candidates at each tree node. Of course, many more parameters can be set to control the algorithm, including the parameters that control the underlying rxDTree() algorithm.

The following is a small example of the rxDForest() fucntion using the mortgage default dataset that can be downloaded from Revolution Analytic's website. Here are the first three lines of data.

  creditScore houseAge yearsEmploy ccDebt year default
1    615        10        5          2818 2000    0
2    780        34        5          3575 2000    0 
3    735        12        1          3184 2000    0 

The idea is to see if the variables creditScore, houseAge etc. are useful in predicting a default. The RevoScaleR R code in the file Download RxDForest reads in the mortgage data, splits the data into a training file and a test file, uses rxDTree() to build a single tree (just to see what one looks like for this file) and plots the tree. Then rxDForest() is run against the training file to to build an ensemble model and this model run against the test file to make predictions. Finally, the code plots the ROC curve for the decision forest ensemble model.

Here is what the first few nodes of the tree looks like. (The full tree is printed at the bottom of the code in the file above.)

Call:
rxDTree(formula = form1, data = "mdTrain", maxDepth = 5)
File: C:\Users\Joe.Rickert\Documents\Revolution\RevoScaleR\mdTrain.xdf
Number of valid observations: 8000290
Number of missing observations: 0

Tree representation:
n= 8000290

node), split, n, deviance, yval
* denotes terminal node

1) root 8000290 39472.30000 4.958445e-03
2) ccDebt< 9085.5 7840182 21402.25000 2.737309e-03
4) ccDebt< 7844 7384170 8809.46500 1.194447e-03 

He is a plot of the right part of the tree drawn with RevoScaleR's creatTreeView() function that enables plot() to put the graph in your browser.

 

 

And, finally, here is the ROC curve for the decision Forest model. (The text output describing the model is also in the file containing the code.)

 

I plan to try rxDForest() out on a cluster with a bigger data set. When I do, I will let you know. 

NYT's 4th Down Bot gives the SuperBowl edge to the Broncos

Who will win the SuperBowl this Sunday: Seattle or Denver? As pundits around the country weigh in with their predictions, you might want to check out the analysis from the New York Times' 4th Down Bot, which compares the coaches' calls on fourth down plays with what historical statistics and a point-forecasting model indicate would have been the ideal play.

According to the Bot, the Seahawks made of 6 of 8 4th down calls correctly in their championship-winning game (one of those was "complicated" according to the Bot, so let's give it to Seattle). In the Bronco's championship game, Denver made 5 of 6 calls correctly. So at least according to the Bot, Denver has the edge over Seattle in 4th down calls. (I'm going to get run out of Seattle for saying that!)

As you're watching the game on Sunday, you can see the Bots calls in real-time on Twitter (follow NYT4thDownBot), and see if you agree with the model or with the coaches calls.

By the way, the New York Times makes extensive use of R in its sports analysis, from fantasy football to draft picks to baseball.

New York Times: 4th Down Bot

Quantitative Finance Applications in R - 3: Plotting xts Time Series

by Daniel Hanson, QA Data Scientist, Revolution Analytics

Introduction and Data Setup

Last time, we included a couple of examples of plotting a single xts time series using the plot(.) function (ie, said function included in the xts package).  Today, we’ll look at some quick and easy methods for plotting overlays of multiple xts time series in a single graph.  As this information is not explicitly covered in the examples provided with xts and base R, this discussion may save you a bit of time.

To start, let’s look at five sets of cumulative returns for the following ETF’s:

SPY SPDR S&P 500 ETF Trust
QQQ PowerShares NASDAQ QQQ Trust
GDX Market Vectors Gold Miners ETF
DBO PowerShares DB Oil Fund (ETF)
VWO Vanguard FTSE Emerging Markets ETF

We first obtain the data using quantmod, going back to January 2007:

library(quantmod)
tckrs <- c("SPY", "QQQ", "GDX", "DBO", "VWO")
getSymbols(tckrs, from = "2007-01-01")

Then, extract just the closing prices from each set:

SPY.Close <- SPY[,4]
QQQ.Close <- QQQ[,4]
GDX.Close <- GDX[,4]
DBO.Close <- DBO[,4]
VWO.Close <- VWO[,4]

What we want is the set of cumulative returns for each, in the sense of the cumulative value of $1 over time.  To do this, it is simply a case of dividing each daily price in the series by the price on the first day of the series.  As SPY.Close[1], for example, is itself an xts object, we need to coerce it to numeric in order to carry out the division:

SPY1 <- as.numeric(SPY.Close[1])
QQQ1 <- as.numeric(QQQ.Close[1])
GDX1 <- as.numeric(GDX.Close[1])
DBO1 <- as.numeric(DBO.Close[1])
VWO1 <- as.numeric(VWO.Close[1])

Then, it’s a case of dividing each series by the price on the first day, just as one would divide an R vector by a scalar.  For convenience of notation, we’ll just save these results back into the original ETF ticker names and overwrite the original objects:

SPY <- SPY.Close/SPY1
QQQ <- QQQ.Close/QQQ1
GDX <- GDX.Close/GDX1
DBO <- DBO.Close/DBO1
VWO <- VWO.Close/VWO1

We then merge all of these xts time series into a single xts object (à la a matrix):

basket <- cbind(SPY, QQQ, GDX, DBO, VWO)

Note that is.xts(basket)returns TRUE. We can also have a look at the data and its structure:

> head(basket)
           SPY.Close QQQ.Close GDX.Close DBO.Close VWO.Close
2007-01-03 1.0000000  1.000000 1.0000000        NA 1.0000000
2007-01-04 1.0021221  1.018964 0.9815249        NA 0.9890886
2007-01-05 0.9941289  1.014107 0.9682540 1.0000000 0.9614891
2007-01-08 0.9987267  1.014801 0.9705959 1.0024722 0.9720154
2007-01-09 0.9978779  1.019889 0.9640906 0.9929955 0.9487805
2007-01-10 1.0012025  1.031915 0.9526412 0.9517923 0.9460847

> tail(basket)
           SPY.Close QQQ.Close GDX.Close DBO.Close VWO.Close
2014-01-10  1.302539        NA 0.5727296  1.082406 0.5118100
2014-01-13  1.285209  1.989130 0.5893833  1.068809 0.5053915
2014-01-14  1.299215  2.027058 0.5750716  1.074166 0.5110398
2014-01-15  1.306218  2.043710 0.5826177  1.092707 0.5109114
2014-01-16  1.304520  2.043941 0.5886027  1.089411 0.5080873
2014-01-17  1.299003  2.032377 0.6070778  1.090647 0.5062901

Note that we have a few NA values here.  This will not be of any significant consequence for demonstrating plotting functions, however.

We will now look how we can plot all five series, overlayed on a single graph.  In particular, we will look at the plot(.) functions in both the zoo and xts packages.

Using plot(.) in the zoo package

The xts package is an extension of the zoo package, so coercing our xts object basket to a zoo object is a simple task:

 zoo.basket <- as.zoo(basket)

Looking at head(zoo.basket) and tail(zoo.basket), we will get output that looks the same as what we got for the original xts basket object, as shown above; the date to data mapping is preserved. The plot(.) function provided in zoo is very simple to use, as we can use the whole zoo.basket object as input, and the plot(.) function will overlay the time series and scale the vertical axis for us with the help of a single parameter setting, namely the screens parameter.

Let’s now look at the code and the resulting plot in the following example, and then explain what’s going on:

# Set a color scheme:
tsRainbow <- rainbow(ncol(zoo.basket))
# Plot the overlayed series
plot(x = zoo.basket, ylab = "Cumulative Return", main = "Cumulative Returns",
        col = tsRainbow, screens = 1)
# Set a legend in the upper left hand corner to match color to return series
legend(x = "topleft", legend = c("SPY", "QQQ", "GDX", "DBO", "VWO"),
       lty = 1,col = tsRainbow)

 

We started by setting a color scheme, using the rainbow(.) command that is included in the base R installation.  It is convenient as R will take in an arbitrary positive integer value and select a sequence of distinct colors up to the number specified.  This is a nice feature for the impatient or lazy among us (yes, guilty as charged) who don’t want to be bothered with picking out colors and just want to see the result right away.

Next, in the plot(.) command, we assign to x our “matrix” of time series in the zoo.basket object, labels for the horizontal and vertical axes (xlab, ylab), a title for the graph (main), the the colors (col). Last, but crucial, is the parameter setting screens = 1, which tells the plot command to overlay each series in a single graph.

Finally, we include the legend(.) command to place a color legend at the upper left hand corner of the graph.  The position (x) may be chosen from the list of keywords "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center"; in our case, we chose "topleft".  The legend parameter is simply the list of ticker names.  The lty parameter refers to “line type”, and by setting it to 1, the lines in the legend are shown as solid lines, and as in the plot(.) function, the same color scheme is assigned to the parameter col.

Back to the color scheme, we may at some point need to show our results to a manager or a client, so in that case, we probably will want to choose colors that are easier on the eye.  In this case, one can just store the colors into a vector, and then use it as an input parameter.  For example, set

myColors <- c("red", "darkgreen", "goldenrod", "darkblue", "darkviolet")

Then, just replace col = tsRainbow with col = myColors in the plot and legend commands:

plot(x = zoo.basket, xlab = "Time", ylab = "Cumulative Return",
        main = "Cumulative Returns", col = myColors, screens = 1)
legend(x = "topleft", legend = c("SPY", "QQQ", "GDX", "DBO", "VWO"),
       lty = 1, col = myColors)

We then get a plot that looks like this:

 

Using plot(.) in the xts package

While the plot(.) function in zoo gave us a quick and convenient way of plotting multiple time series, it didn’t give us much control over the scale used along the horizontal axis.  Using plot(.) in xts remedies this; however, it involves doing more work.  In particular, we can no longer input the entire “matrix” object; we must add each series separately in order to layer the plots.  We also need to specify the scale along the vertical axis, as in the xts case, the function will not do this on the fly as it did for us in the zoo case.

We will use individual columns from our original xts object, basket.  By using basket rather than basket.zoo, this tells R to use the xts version of the function rather than the zoo version (à la an overloaded function in traditional object oriented programming).  Let’s again look at an example and the resulting plot, and then discuss how it works:

plot(x = basket[,"SPY.Close"], xlab = "Time", ylab = "Cumulative Return",
main = "Cumulative Returns", ylim = c(0.0, 2.5), major.ticks= "years",
        minor.ticks = FALSE, col = "red")
lines(x = basket[,"QQQ.Close"], col = "darkgreen")
lines(x = basket[,"GDX.Close"], col = "goldenrod")
lines(x = basket[,"DBO.Close"], col = "darkblue")
lines(x = basket[,"VWO.Close"], col = "darkviolet")
legend(x = 'topleft', legend = c("SPY", "QQQ", "GDX", "DBO", "VWO"),
      lty = 1, col = myColors)

As mentioned, we need to add each time series separately in this case in order to get the desired overlays.  If one were to try x = basket in the plot function, the graph would only display the first series (SPY), and a warning message would be returned to the R session.  So, we first use the SPY series as input to the plot(.) function, and then add the remaining series with the lines(.) command.  The color for each series is also included at each step (the same colors in our myColors vector).

As for the remaining arguments in the plot command, we use the same axis and title settings in xlab, ylab, and main.  We set the scale of the vertical axis with the ylim parameter; noting from our previous example that VWO hovered near zero at the low end, and that DBO reached almost as high as 2.5, we set this range from 0.0 to 2.5.  Two new arguments here are the major.ticks and minor.ticks settings. The major.ticks argument represents the periods in which we wish to chop up the horizontal axis; it is chosen from the set

{"years", "months", "weeks", "days", "hours", "minutes", "seconds"}

In the example above, we chose "years".  The minor.ticks parameter can take values of TRUE/FALSE, and as we don’t need this for the graph, we choose FALSE.  The same legend command that we used in the zoo case can be used here as well (using myColors to indicate the color of each time series plot). Just to compare, let’s change the major.ticks parameter to "months" in the previous example. The result is as follows:

 

Wrap-up

A new package, called xtsExtra, includes a new plot(.) function that provides added functionality, including a legend generator.  However, while it is available on R-Forge, it has not yet made it into the official CRAN repository.  More sophisticated time series plotting capability can also be found in the quantmod and ggplot2 packages, and we will look at the ggplot2 case in an upcoming post.  However, for plotting xts objects quickly and with minimal fuss, the plot(.) function in the zoo package fills the bill, and with a little more effort, we can refine the scale along the horizontal axis using the xts version of plot(.).  R help files for each of these can be found by selecting plot.zoo and plot.xts respectively in help searches.

John Chambers recounts the history of S and R

"R has had a revolutionary effect on the way statistics are communicated." So says John Chambers: one of the members of the R-core team overseeing R; and co-inventor of the S language. In this interview with Trevor Hastie (his co-author on Statistical Models in S), John Chambers recounts his involvement in the birth of the S language in 1976, and how it evolved over the years to become the inspiration for the R language.

 

(via Siamak Faridani.) One interesting tid-bit from the video: John Chambers owns the original CD-ROM (serial number #1) of R 1.0.0, released on February 29 2001 2000, and signed by all the members of R-core.

In 1998, John Chambers won the ACM Software Systems award for the S language, which the ACM heralded as having "forever altered how people analyze, visualize, and manipulate data". He generously used the prize money to create the John M Chambers Statistical Software Award, which recognizes innovative statistical software written by an undergraduate or graduate student. If you know a deserving recipient, the deadline for nominations for the 2014 award is February 18.

Because it's Friday: Fly, fly away

In 2012, filmmakers Sophie and Liberty took a canoe trip on the Shannon River in Ireland. You can see their approach to a fog-enshrouded island in the still frames at the start of this video, and then they encountered sometime breathtaking in motion:

(Via Vimeo Staff Picks.) It still amazes me that simple rules of individual agents (all one bird knows is to avoid predators and not collide with other birds) can create such complex (and beautiful) emergent behaviour. On a similar note, tracing the paths of indvidual birds create their own distinct patterns:

 

(via FlowingData). We'll be back on Monday, but if you're out and about this weekend, take a closer look at the birds in flight. You might see them in a different light.

Demo this Wednesday: Drag-and-drop to create R-based workflows

Want to see how you can use a drag-and-drop user interface to run and share R code? Check out our webinar next Wednesday January 29 (hosted by Alteryx and Revolution Analytics): Creating Value That Scales with Revolution Analytics & Alteryx. 

In the webinar, Dan Putler (Alteryx's Data Artisan in Residence) will demonstrate the Alteryx drag-and-drop Alteryx GUI, which provides access to many of R functions (and which you can extend with your own). And because it connects with Revolution R Enterprise, you also benefit from its performance and big-data capabilities. During the demo, Dan will show:

• How to create a complete data analysis workflow with the Alteryx drag-and-drop user interface;
• How to run R functions as workflow nodes, and create your own custom R-based nodes;
• How to analyze big data files without memory limitations, thanks to the Revolution R Enterprise XDF out-of-memory data object;
• How to publish an Alteryx workflow to a gallery, with an easy-to-use web-based interface that anyone can use to control its operation.

If you're an R user looking for an easier way to interact with R, or looking for a way to share access to your R code running on big data in an easy-to-use web-based interface, you'll want to register for the webinar at the link below.

Webinar Jan 29: Creating Value That Scales with Alteryx and Revolution Analytics

Book review: "Doing Data Science" by Rachel Schutt and Cathy O'Neil

by Joseph Rickert

Every once in a while a single book comes to crystallize a new discipline. If books still have this power in the era of electronic media, "Doing Data Science, Straight Talk from the Frontline" by Rachel Schutt and Cathy O’Neil: O'Reilly, 2013 might just be the book that defines data science. "Doing Data Science", which is based on a course that Rachel taught at Columbia University and to which Cathy contributed, is ambitious and multidimensional. It presents data science in all of its messiness as an open-ended practice that is coalescing around an expanding class of problems; problems which are yielding to an interdisciplinary approach that includes ideas and techniques from statistics, computer science, machine learning, social science and other disciplines.

The book is neither a statistics nor a machine learning text, but there are plenty of examples of statistical models and machine learning algorithms. There is enough R code in the text to get a beginner started on real problems with tools that are immediately useful. There is Python code, a bash shell script, mention of JSON and a down to earth discussion of Hadoop and MapReduce that many should find valuable. My favorite code example is the bash script (p 105) that fetches an Enron spam file and performs some basic word count calculations. Its almost casual insertion into the text, without fanfare and little explanation, provides a low key example of the kinds of baseline IT/ programmer skills that a newly minted statistician must acquire in order to work effectively as a data scientist.

"Doing Data Science" is fairly well balanced in its fusion of the statistics and machine learning world views, but Rachel’s underlying bias as a PhD statistician comes through when it counts. The grounding in linear models and the inclusion of time series models establish the required inferential skills. The discussion of causality shows how statistical inference is essential to obtaining a deep understanding of how things really work, and the chapter on epidemiology provides a glimpse into just how deep and difficult are the problems that statisticians have been wrestling with for generations. (I found the inclusion of this chapter in a data science book to be a delightful surprise.)

It is not only the selection of material, however, that betrays the book's statistical bias. When the authors take on the big questions their language indicates a statistical mindset. For example, in the discussion following "In what sense does data science deserve the word “science” in its name?" (p114) the authors write: “Every design choice you make can be formulated as an hypothesis, against which you will use rigorous testing and experimentation to validate or refute”. This is the language of a Neyman/Pearson trained statistician trying to pin down the truth. It stands in stark contrast with the machine learning viewpoint espoused in a quote by Kaggle’s Jeremy Howard who, when asked “Can you see any downside to the data-driven, black-box approach that dominates on Kaggle?”, replies:

Some people take the view that you don’t end up with a richer understanding of the problem. But that’s just not true: The algorithms tell you what’s important and what’s not. You might ask why those things are important, but I think that’s less interesting. You end up with a predictive model that works. There is not too much to argue about there.

So, whether you are doing science or not might just be in your intentions and point of view. Schutt and O’Neil do a marvelous job of exploring the tension between the quest for understanding and and the blunt success of just getting something that works.

An unusual aspect of the book is its attempt to understand data science as a cultural phenomenon and to place the technology in a historical and social context. Most textbooks in mathematics, statistics and science make no mention of how things came to be. Their authors are just under too much pressure to get on with presenting the material to stop and and discuss “just what were those guys thinking?”. But Schutt and O’Neill take the time, and the book is richer for it. Mike Driscoll and Drew Conway, two practitioners who early on recognized that data science is something new, are quoted along with other contemporary data scientists who are shaping the discipline both through their work and how they talk about it.

A great strength of the book is its collection of the real-world, big-league examples contributed by the guest lecturers to Rachel’s course.  Doug Perlson of Real Direct, Jake Hofman of Microsoft Research, Brian Dalessandro and Claudia Perlich both of Media6Degrees, Kyle Teague of GetGlue, William Cukierski of Kaggle, David Huffaker of Google, Matt Gattis of Hutch.com, Mark Hansen of Columbia University, Ian Wong of Square, John Kelley of Morningside Analytics and David Madigan, Chair of the Columbia’s Statistics Department, all bring thoughtful presentations of difficult problems with which they have struggled. The perspective and insight of these practicing data scientists and statisticians is invaluable. Claudia Perlich’s discussion of data leakage alone is probably worth the price of the book.

A minor fault of the book is the occasional lapse into the hip vulgar. Someone being “pissed off” and talking about a model “that would totally suck” are probably innocuous enough phrases, but describing a vector as “huge ass” doesn’t really contribute to clarity. In a book that stresses communication, language counts. Nevertheless, "Doing Data Science" is a really “good read”. The authors have done a remarkable job of integrating class notes, their respective blogs, and the presentations of the guest speakers into a single, engaging voice that mostly speaks clearly to the reader. I think this book will appeal to a wide audience. Beginners asking the question “How do I get into data science?” will find the book to be a guide that will take them a long way. Accomplished data scientists will find a perspective on their profession that they should appreciate as being both provocative and valuable. "Doing Data Science" argues eloquently for a technology that respects humanist ideals and ethical considerations. We should all be asking "What problems should I be working on?", "Am I doing science or not?", and "What are the social and ethical implications of my work?". Finally, technical managers charged with assembling a data science team, and other interested outsiders, should find the book helpful in getting beyond the hype and and having a look at what it really takes to squeeze insight from data.

Fast and easy data munging, with dplyr

RStudio's Hadley Wickham has just introduced a new package for filtering, selecting, restructuring and aggregating tabular data in R: the dplyr package. It's similar in concept to Hadley's original plyr package from 2009, but with several key improvements:

• It works exclusively with data in R data frames;
• It can process data in remote databases (with the transformations done in-database — only the result is returned to R);
• It introduces a "grammar of data manipulation", allowing you to string together operations with the %.% operator;
• And it's much, much faster than plyr or standard R operations (most processing is done in parallel in C++).

For example, here's how you'd calculate the Top 5 all-time batters based in number of runs scored (based on 2012 MLB data, courtesy of the Lahman package)

Batting %.%
  group_by(playerID) %.%
  summarise(total = sum(R)) %.%
  arrange(desc(total)) %.%
  head(5)

   playerID total
1 henderi01  2295
2  cobbty01  2246
3 bondsba01  2227
4  ruthba01  2174
5 aaronha01  2174

There are plenty of ways you could solve the same problem in standard R, but here's one way:

totals <- aggregate(. ~ playerID, data=Batting[,c("playerID","R")], sum)
ranks <- sort.list(-totals$R)
totals[ranks[1:5],]

The dplyr version took 0.036 seconds on my MacBook Air, compared to 0.266 seconds for the "standard R" version — about a 7x speedup. More importantly, it'd been so long since I used R's aggregate function that it took me about 5 minutes and several debugging attempts before I got it right. (To be fair, I copy and pasted the dplyr example from Hadley, but still: it looks a lot cleaner and less error prone.) There are certainly faster and easier ways to do this in standard R, but that's kind of the point: a standard "grammar of data manipulation" provides consistency, and with a speed boost to boot!

The dplyr package is available on CRAN now. You can read more about dplyr in the introductory vignette and in Hadley's blog post linked below.

RStudio blog: Introducing dplyr