Revolutions

A first look at rxBTrees

by Joseph Rickert

The gradient boosting machine as developed by Friedman, Hastie, Tibshirani and others, has become an extremely successful algorithm for dealing with both classification and regression problems and is now an essential feature of any machine learning toolbox. R’s gbm() function (gbm package) is a particularly well crafted implementation of the gradient boosting machine that served as a model for the rxBTrees() function which was released last year as a feature of Revolution R Enterprise 7.3. You can think of rxBTRees(), as a scaled up version gbm() that is designed to work with massive data sets in distributed environments.

The basic idea underlying the gbm is that weak learners may be combined in an additive, iterative fashion to produce highly accurate classifiers that are resistant to overfitting. While in theory any weak learner can be incorporated into a gbm, classification and regression trees, or decision trees as they are called in RRE, are particularly convenient. As Kuhn and Johnson (p 205) point out:

1. They have the flexibility to be weak learners by simply restricting their depth
2. Separate trees can be easily added together
3. Trees can be generated very quickly

rxBTrees() scales gbm() by building the boosting algorithm around RRE’s rxDTree() decision tree implementation. As I described in a previous post, this algorithm scales to work with large data sets by building trees on histograms, a summary of the data, and not on the data itself. rxBTrees(), like rxDTree() is also a proper parallel external memory algorithm (PEMA) This means two things:

1. It uses a “chunking” strategy to work with a block of data at a time that has been read from disk. Since there is no need to keep all of the data in memory at one time, rxBtrees() can build a gbm on arbitrarily large files.
2. PEMA’s are parallel algorithms built on top of RRE’s distributed computing infrastructure. This allows them to take full advantage of the opportunities for parallelism implicit in the underlying hardware by setting up a “compute context”. In the example below, the compute context is set to local parallel to take advantage of the cores on my PC. RRE also allows establishing compute contexts that enable PEMAs to run in parallel on clusters.

The following function fits a rxBTrees() model to a training data subset of the a mortgage data set. rxBTrees() contains many input parameters. Some are common to all RevoScaleR PEMAs, others are inherited from rxDTree(), and still others such as nTree and learningRate have direct counterparts in R's gbm() function.  In this function call, I have explicitly called out parameters that pertain to fitting a gbm model

# Fit a Boosted Trees Mode on the training data set
form <- formula(default ~ houseAge + creditScore + yearsEmploy + ccDebt)
 
BT_model <- rxBTrees(formula = form, 
                     data= "mdTrain",
                     lossFunction = "bernoulli",
                     minSplit = NULL,    # Min num of obs at node before a split is attempted
                     minBucket = NULL,   # Min num of obs at a terminal node
                     maxDepth = 1,       # 2^maxDepth = number of leaves on tree
                     cp = 0,             # Complexity parameter
                     maxSurrogate = 0,   # Max surrogate splits retained, 0 improves perf
                     useSurrogate = 2,   # Determines how surrogates are used
                     surrogateStyle = 0, # Sets surrogate algorithm
                     nTree = 100,        # Num of trees = num of boosting iterations
                     mTry = NULL,        # Num vars to sample as a split candidate 
                     sampRate = NULL,    # % of observations to sample for each tree
                                         # default = 0.632
                     importance = TRUE,  # Importance of predictors to be assessed
                     learningRate = 0.1, # Also known as shrinkage
                     maxNumBins = NULL,  # Max num bins in histograms
                                         # default: min(1001,max(sqrt(num of obs)))
                     verbose = 1,        # Amount of status sent ot console
                     computeContext="RxLocalParallel",  # use all cores available in PC
                     overwrite = TRUE)   # Allow xdf output file to be overwritten

Some especially noteworthy parameters are maxDepth, maxNumBins and computeContext. maxDepth accomplishes the same purpose as  interaction.depth in the gbm() function. However, they are not the same parmeter. Setting interaction.depth = 1 produces stump trees with two leaves each because the number of leaves = interaction.depth + 1. In rxBTRees(), stumps are produced by setting maxDepth = 1. But, this is because 2^maxDepth gives the number of leaves.

MaxNumBins determines the number of bins to be used in making histograms of the data. If anyone ever wanted to compare gbm() directly with rxBTrees() on a small data set for example, setting MaxNumBins to the number of points in the data set would mimic gbm(). Note, however, that this would saddle rxBTrees() with all of the overhead of computing histograms while receiving none of the benefits.

The computeContext parameter is the mechanism for connecting a PEMA to the underlying hardware and that enables functions such as rxBTrees() to run in parallel on a various clusters using distributed data. Some of the details of how this happens have been described in a previous post, but I hope to be able to follow up this post with an example of rxBTrees() running on Hadoop in the not too distant future.

The data set used to build the model above is available on Revolution Anaytics' data set web page both as an xdf file and as multiple csv files. The following command shows what the data looks like. 

> rxGetInfo(mortData,getVarInfo=TRUE,numRows=3)

File name: C:\DATA\Mortgage Data\mortDefault\mortDefault2.xdf 
Number of observations: 1e+07 
Number of variables: 6 
Number of blocks: 20 
Compression type: none 
Variable information: 
Var 1: creditScore, Type: integer, Low/High: (432, 955)
Var 2: houseAge, Type: integer, Low/High: (0, 40)
Var 3: yearsEmploy, Type: integer, Low/High: (0, 15)
Var 4: ccDebt, Type: integer, Low/High: (0, 15566)
Var 5: year, Type: integer, Low/High: (2000, 2009)
Var 6: default, Type: integer, Low/High: (0, 1)
Data (3 rows starting with row 1):
  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

In building the model, I divided the mortgage default data randomly into a training and test data sets. Then, I fit an rxBTrees() model on the training data and used the test data to compute the ROC curve shown below.

 

Along the way, the model also produced the out of bag error versus number of trees plot which shows that for this data set one could probably get by with only building 70 or 80 trees.

The code I used to build the model is available here: Download Post_model. If you do not have Revolution R Enterprise available to you at your company, you can try out this model for yourself in the AWS cloud by taking advantage of Revolution Analytics' 14 day free trial.

Take the first step towards building a gbm on big data.

Updated checkpoint package: faster reproducibility with more feedback

A new version of the checkpoint package for R has just been released on CRAN. With the checkpoint package, you can easily:

• Write R scripts or projects using CRAN package versions from a specific point in time;
• Share R scripts with others that will automatically install the appropriate package versions (no need to manually install CRAN packages);
• Write R scripts that use older versions of packages, or packages that are no longer available on CRAN;
• Install packages (or package versions) visible only to a specific project, without affecting other R projects or R users on the same system;
• Manage multiple projects that use different package versions;
• Write and share code R whose results can be reproduced, even if new (and possibly incompatible) package versions are released later.

The biggest change with this new version of checkpoint is that it checks whether package versions have already been installed for your project: running a checkpoint'ed script for the second time has virtually no overhead now. There's also better feedback during the scanning phase while it determines which (if any) packages need to be installed. You can also now control where packages are installed on your system (the default is a .checkpoint folder in your home directory), and check the specific version of R in use when checkpoint is called.

If you're interested in learning more about how checkpoint works, you can read the new vignette Using checkpoint for reproducible research, or follow its development at the checkpoint Github page. You can install checkpoint from CRAN now by running install.packages("checkpoint") at the R command line.

MRAN: checkpoint: Install Packages from Snapshots on the Checkpoint Server for Reproducibility

Some thoughts on Vim

by Gary R. Moser
Director of Institutional Research and Planning
The California Maritime Academy

I recently contacted Joseph Rickert about inviting Vim guru Drew Niel (web: vimcasts.org, book: "Practical Vim: Edit Text at the Speed of Thought") to speak at the Bay Area R User Group group. Due to Drew's living in Great Britain that might not be easily achieved, so Joe generously extended an invitation for me to share a bit about why I like Vim here.

When it comes to text editors, there are a number of strong contenders to choose from (...lots and lots of them). if you've found a tool with a rich feature set that makes you more productive (and more importantly, is something you like) you should probably stick with that. Given that editing text of various types is such a big part of our lives, however, it's worth turning a critical eye toward how well our tools facilitate this specific activity.

Do you spend a lot of time editing text? Yeah, me too. Personally this includes R code, Markdown, Latex, as well as notes and outlines (including this article). When I rediscovered Vim fairly recently, it was from the perspective of being a proficient R user. Therefore, my values and beliefs about what constitutes good software are framed by my experience with R. That includes being OK with an initial learning curve.

I have been using RStudio as my primary IDE for R since it was offered.  It's great software; they did something not too long ago that I really appreciate - they added a stripped-down Vim editing mode. Vim pushes users to the documentation pretty quickly (usually as the result of accidentally deleting large chunks of work), and as I dug in and began to discover its full functionality, I came to realize how much I was missing out on by using the emulator. The ability to set and jump to marks in a document, or utilizing multiple registers for copy/paste are two good examples of essential but missing features in RStudio.

Vim has been described as "a language for text editing," which I think is a useful way to view it. At the risk of sounding snotty, I would compare the experiences of using Vim (or another good editor) versus a plain-jane text editor to that of playing chess versus checkers. That is, there's an element of strategic and intentional action compared to simply playing one of a limited set of moves over and over again.

One of the things that makes Vim so interesting and different from other editors stems from its origins. As the result of being developed in the context of severe constraints (slow networks, no digital displays, limited system resources, and no mouse), Vim - then "ed" - had to accomplish the greatest amount of work with the least number of keystrokes. This requirement led to the development of a vast number of very specific commands that can be combined in useful ways. Drew Neil artfully compares this to playing notes, chords, and melodies on a piano. It's also an appropriate comparison for setting one's expectations toward becoming a skilled Vim user! Michael Mrozekon's humorous plot cleverly suggests that, not unlike R, Vim doesn't hold your hand.

It also speaks to my point about specificity. Emacs, for example, can be extended to be a web client or music player, hence the rabbit-hole learning curve, but isn't that getting away from the primary task of text editing?

The fundamental way that Vim differs from most other text editors is that it is explicitly modal; all software is technically modal in certain ways (that is, the same keys do different things under different circumstances), but with Vim it is a central design feature. Essentially, what this means is that by switching modes, a different keyboard comes into existence under your fingers.  Because Vim has four modes, and a very rich and terse set of key-bindings, it's like having four+ keyboards in one. The keyboard cheat sheet is a useful reference, especially in the beginning.

Warning: after becoming familiar with Vim's basic functionality, going back to a typical text editor feels rather clumsy.

Vim as an interface to R using the Vim-R-plugin is mostly good for how I use it, but I expect to be dialing-in Vim for a long time before it's got all the features I want. I don't mind this, but I can see how someone else might. I encourage you consider your own tools and how well they facilitate your most frequent tasks. If you're an RStudio user, try giving Vim mode a go. A visit to www.vim.org will connect you to the resources you'll need.

The impact of open source software on the data science revolution

It's hard to overstate the role of open-source software in the data science revolution. Tools like Hadoop, Spark, R, and Python are essential parts of the modern data science toolkit. These tools are likewise part of the solutions built by the Consulting Services group at Revolution Analytics. Our VP of Professional Services, Neera Talbert, shares her view on the impact of open source software in an article for Inside BigData published today. Here's Neera on the typical day of a data scientist:

Unlike a pure statistician, a data scientist is also expected to write code and understand business. Data science is a multi-disciplinary practice requiring a broad range of knowledge and insight. It’s not unusual for a data scientist to explore a fresh set of data in the morning, create a model before lunch, run a series of analytics in the afternoon and brief a team of digital marketers before heading home at night.

You can read the entire article at the link below.

Inside BigData: Open Source Software Fuels a Revolution in Data Science

 

 

Because it's Pi Day: Pie and Paella

Today is a day of joy for math geeks in the USA, and a day of bemusement for maths geeks everywhere else: Pi Day. (If only April had 31 days!) In fact, today is a very special once-in-a-century Pi day: at 9:26:53.5898 this morning was a "pi instant" where the time and date read out pi to 10 significant digits.

I'm celebrating Pi day on a cooking weekend with friends on an island getaway in Puget Sound. We began with pi-ella:

and of course, finished with apple pi:

Happy pi day everyone, and we'll see you back here on Monday.

New Whitepaper: Connect R to other applications with DeployR

DeployR is a server-based framework that provides simple, secure R integration for application developers. It's available in two editions: DeployR Open, which is free and open-source; and Revolution R Enterprise DeployR, which adds a scalable grid framework and enterprise authentication features for production applications integrated with R. 

If you're looking for an overview of what DeployR is and how you can use it to access R from other applications, we've just released a new white paper, Using DeployR to Solve the R Integration Problem. This 14-page document steps you through the basic DeployR workflow: the R user creates an script and publishes it to the server; an application developer links to the R script via the DeployR API; and an application user sees some results of an R calculation in a desktop, mobile or web-based application (most likely without being aware that R was involved at all). 

You can download the whitepaper (free registration required) from the Whitepapers section of the Revolution Analytics website. For the technically-minded, you may also want to explore the DeployR Application Developer Documentation. You can also check out the list of our other most popular white papers and webinars on R, data science and big data analytics. 

Revolution Analytics White Papers: Using DeployR to Solve the R Integration Problem

 

A Monte Carlo Simulation for Pi Day

by Joseph Rickert

What will you be doing at 26 minutes and 53 seconds past 9 this coming Saturday morning? I will probably be running simulations. I have become obsessed with an astounding result from number theory and have been trying to devise Monte Carlo simulations to get at it. The result, well known to number theorists says: choose two integers at random; the probability that they will be coprime is 6/π² ! Here, π materializes out of thin air. Who could have possibly guessed this? - well, Leonhard Euler, apparently, and this sort of magic seems to be quite common in number theory. 

More formally, the theorem Euler proved goes something like this: Let P_N be the probability that two randomly chosen integers in {1, 2, ... N} are coprime. Then, as N goes to infinity P_N goes to 6/π² . Well, this seems to be a little different. You don't actually have to sample from an infinite set. So, I asked myself, would a person who is not familiar with this result, but who is allowed to do some Monte Carlo simulations, have a reasonable chance of guessing the answer? I know that going to infinity would be quite a trip, but I imagined that I could take a few steps and see something interesting. How about this: choose some boundary, N, and draw lots of pairs of random numbers. Count the number of coprime pairs M. Then sqrt( 6 / (M/N)) will give an estimate of π. As N gets bigger and bigger you should see the digits of the average of the estimates marching closer and closer to π - 3.1, 3.14, 3.141 etc. 

Before you go off and try this, let me warn you: there is no parade of digits. For a modest 100,000 draws with an N around 10,000,000 you should get an estimate of π close to 3.14. Then, even with letting N get up around 1e13 with a million draws, you won't do much better. The following code (compliments of a colleague who introduced me to mapply) performs 100 simulations, each with 100,000 draws as the range varies from +/- 1,000,000 to +/- 1e13. (The code runs pretty quickly on my laptop.)

# Monte Carlo estimate of pi
library( numbers )
library(ggplot2)
 
set.seed(123)
 
bigRange <- seq(1e6,1e13,by=1e11)
M <- length(bigRange)
 
draws <- 1e5
 
for(i in 1:M){
  maxRange <- bigRange[i]
  print(bigRange[i])
  min <- -maxRange
  max <- maxRange
 
  r1 <- round( runif( n = draws, min = min, max = max ) )
  r2 <- round( runif( n = draws, min = min, max = max ) )
 
  system.time( coprimeTests <- mapply( coprime, r1, r2 ) )
 
  prob[i] <- sum( coprimeTests ) / draws
  print(prob[i])
  piEst[i] <- sqrt( 6 / prob[i] )
  print(piEst[i])
}
 
 
piRes2 <- data.frame(bigRange,prob,piEst)
p2 <- ggplot(piRes2,aes(bigRange,piEst))
p2 + geom_line() +
  geom_point() +
  xlab("Half Range for Random Draws")+
  ylab("Estimate of Pi") + 
  ggtitle("Expanding Range Simulation") +
  geom_smooth(method = "lm", se=FALSE, color="black")

Here is the "time series" plot of the results.

The mean is 3.142959. The slight upward trend is most likely a random number induced illusion.

lm(formula = piRes$piEst ~ piRes$index)
 
Coefficients:
(Intercept)  piRes$index  
  3.142e+00    1.225e-05  

The problem as, I have framed it, is probably beyond the reach of a naive Monte Carlo approach. Nevertheless, on Saturday, when I have some simulation, time I will try running 100,000,000 draws. This should get me another digit of π since the accuracy of the mean increases as sqrt(N). 

The situation, however, is not as dismal as I have been making it out. Don't imagine that just because the problem is opaque to a brute force Monte Carlo effort that it is without the possibility of computational illumination. Euler's proof, mentioned above, turns on recognizing that the expression for the probability of two randomly chosen integers being coprime may be expressed as the mathematical function z(2). The poof yields π² = 6z(2), The wizards of R have reduced this calculation to the trivial. The Rmpfr package, which allows the use of arbitrarily precise numbers instead of R's double precision numbers, includes the function zeta(x)! So, here, splendidly arrayed, is π.

library(Rmpfr)
sqrt(6*zeta(2))
 
> sqrt(6*zeta(2))
1 'mpfr' number of precision  128   bits 
[1] 3.141592653589793238462643383279502884195

Happy Pi Day! 

In case you missed it: February 2015 roundup

In case you missed them, here are some articles from February of particular interest to R users.

The John M. Chambers Statistical Software Award announcement for 2015.

The new R package "distcomp" allows researchers to collaborate on data spread across multiple sites.

David Smith's interview with theCUBE on R, data science, and Microsoft's acquisition of Revolution Analytics.

R is used to measure impact of climate change, and other Strata keynote presentations. 

Some tricks for monitoring the progress on parallel R jobs using foreach. 

"Analytics Marketplaces" are all the rage today, but CRAN was there first. 

A preview of some of the major R-related conferences and events of 2015. 

The checkpoint package has been updated to make it even easier to run R scripts with fixed R package versions.

A tutorial to introduce R to users of Microsoft Excel.

R used to assess the "virality" of posts on new-media sites like Buzzfeed.

A review of the HP Workshop on Distributed Computing in R featuring Luke Tierney, Dirk Eddelbuettel, Martin Morgan, Simon Urbanek and other R luminaries.

R is a top-ranked language on GitHub, according to statistics like number of repositories, code pushes and number of forks. 

You can use the rcrunchbase package to access data on startup companies. 

The R package "syuzhet" applies sentiment analysis to novels to infer their dramatic arc.

The new "quickcheck" package provides assertion-based testing with random inputs for R.

A theorem for calculating an upper bound for the generalization error of a machine learning classifier.

Some practical advice for sharing Shiny applications with shinyapps.io.

A visualization of Paris's street orientations reveals the history of the city.

General interest stories (not related to R) in the past month included: how our brains trick us into seeing the wrong colors, an hilarious parody of cooking shows, a new ASA website to help journalists with Statistics, the statistical model behind the rules of cricket, and why rivers meander.

As always, thanks for the comments and please send any suggestions to me at [email protected]. Don't forget you can follow the blog using an RSS reader, via email using blogtrottr, or by following me on Twitter (I'm @revodavid). You can find roundups of previous months here.

Creating progress bars with foreach parallel processing

by Andrie de Vries

In my previous post, I demonstrated how to get some status of running jobs on a parallel back end. However, I stopped short of actually demonstrating progress bars. In this post I demonstrate how to do this.

The StackOverflow question How do you create a progress bar when using the “foreach()” function in R? () asks the question how to create progress bars with foreach.

This question attracted some helpful answers, but ultimately didn’t answer the question.  In particular, Steve Weston said this approach is flawed, since

It won't work well with doParallel because doParallel only calls the combine function after all of the results have been returned, since it is implemented by calling the parallel clusterApplyLB function. This technique only with works well with backends that call the combine function on-the-fly, like doRedis, doMPI, doNWS, and (defunct?) doSMP.

My experiments

However, this set me thinking.  Could it be that one could trigger a tcltk progress bar in each worker?

Thus I started to experiment, and came up with the following code.  On my windows machine, using doParallel, this seems to work reliably (but I could not get this approach to work on Linux using doMC).

When you run this code, you should see synchronous progress bars, one for each worker:

Conclusion

Although it takes some work, and there isn’t a single solution for all cases, it is possible to get progress bars for at least some parallel back ends.

In this post I demonstrated how to create progress bars using doParallel on Windows, but failed to replicate on Linux.

Can you do better?

R 3.1.3 now available

R 3.1.3, the final update in the R 3.1 series, has been released. As of this writing only the source distribution is currently available, but expect binary builds for Windows, Mac and various Linux platforms to appear soon on your local CRAN mirror.

As has become usual in March, this release is primarily for minor bugs and improvements in R, as the development switches focus to the next major release (R 3.2.0 is expected on April 16). Improvements include: ability to download files greater than 2GB on 32-bit builds; better handling of complex (imaginary) numbers for several functions; improved command completion in the Windows GUI; and improved performance when using S4 classes. (See the complete list here.)

r-announce mailing list: R 3.1.3 is released