Revolutions

Targeted Learning R Packages for Causal Inference and Machine Learning

by Sherri Rose
Assistant Professor of Health Care Policy
Harvard Medical School

Targeted learning methods build machine-learning-based estimators of parameters defined as features of the probability distribution of the data, while also providing influence-curve or bootstrap-based confidence internals. The theory offers a general template for creating targeted maximum likelihood estimators for a data structure, nonparametric or semiparametric statistical model, and parameter mapping.  These estimators of causal inference parameters are double robust and have a variety of other desirable statistical properties. 

Targeted maximum likelihood estimation built on the loss-based “super learning” system such that lower-dimensional parameters could be targeted (e.g., a marginal causal effect); the remaining bias for the (low-dimensional) target feature of the probability distribution was removed. Targeted learning for effect estimation and causal inference allows for the complete integration of machine learning advances in prediction while providing statistical inference for the target parameter(s) of interest. Further details about these methods can be found in the many targeted learning papers as well as the 2011 targeted learning book.

Practical tools for the implementation of targeted learning methods for effect estimation and causal inference have developed alongside the theoretical and methodological advances. While some work has been done to develop computational tools for targeted learning in proprietary programming languages, such as SAS, the majority of the code has been built in R.

Of key importance are the two R packages SuperLearner and tmle. Ensembling with SuperLearner allows us to use many algorithms to generate an ideal prediction function that is a weighted average of all the algorithms considered. The SuperLearner package, authored by Eric Polley (NCI), is flexible, allowing for the integration of dozens of prespecified potential algorithms found in other packages as well as a system of wrappers that provide the user with the ability to design their own algorithms, or include newer algorithms not yet added to the package. The package returns multiple useful objects, including the cross-validated predicted values, final predicted values, vector of weights, and fitted objects for each of the included algorithms, among others.

Below is sample code with the ensembling prediction package SuperLearner using a small simulated data set. 

library(SuperLearner)
##Generate simulated data##
set.seed(27)
n<-500
data <- data.frame(W1=runif(n, min = .5, max = 1),
W2=runif(n, min = 0, max = 1),
W3=runif(n, min = .25, max = .75),
W4=runif(n, min = 0, max = 1))
data <- transform(data, #add W5 dependent on W2, W3
W5=rbinom(n, 1, 1/(1+exp(1.5*W2-W3))))
data <- transform(data, #add Y dependent on W1, W2, W4, W5
Y=rbinom(n, 1,1/(1+exp(-(-.2*W5-2*W1+4*W5*W1-1.5*W2+sin(W4))))))
summary(data)
 
##Specify a library of algorithms##
SL.library <- c("SL.nnet", "SL.glm", "SL.randomForest")
 
##Run the super learner to obtain predicted values for the super learner as well as CV risk for algorithms in the library##
fit.data.SL<-SuperLearner(Y=data[,6],X=data[,1:5],SL.library=SL.library, family=binomial(),method="method.NNLS", verbose=TRUE)
 
##Run the cross-validated super learner to obtain its CV risk##
fitSL.data.CV <- CV.SuperLearner(Y=data[,6],X=data[,1:5], V=10, SL.library=SL.library,verbose = TRUE, method = "method.NNLS", family = binomial())
 
##Cross validated risks##
mean((data[,6]-fitSL.data.CV$SL.predict)^2) #CV risk for super learner
fit.data.SL #CV risks for algorithms in the library

The final lines of code return the cross-validated risks for the super learner as well as each algorithm considered within the super learner. While a trivial example with a small data set and few covariates, these results demonstrate that the super learner, which takes a weighted average of the algorithms in the library, has the smallest cross-validated risk and outperforms each individual algorithm. 

The tmle package, authored by Susan Gruber (Reagan-Udall Foundation), allows for the estimation of both average treatment effects and parameters defined by a marginal structural model in cross-sectional data with a binary intervention. This package also includes the ability to incorporate missingness in the outcome and the intervention, use SuperLearner to estimate the relevant components of the likelihood, and use data with a mediating variable. Additionally, TMLE and collaborative TMLE R code specifically tailored to answer quantitative trait loci mapping questions, such as those discussed in Wang et al 2011, is available in the supplementary material of that paper.

The multiPIM package, authored by Stephan Ritter (Omicia, Inc.), is designed specifically for variable importance analysis, and estimates an attributable-risk-type parameter using TMLE. This package also allows the use of SuperLearner to estimate nuisance parameters and produces additional estimates using estimating-equation-based estimators and g-computation. The package includes its own internal bootstrapping function to calculate standard errors if this is preferred over the use of influence curves, or influence curves are not valid for the chosen estimator.

Four additional prediction-focused packages are casecontrolSL, cvAUC, subsemble, and h2oEnsemble, all primarily authored by Erin LeDell (Berkeley). The casecontrolSL package relies on SuperLearner and performs subsampling in a case-control design with inverse-probability-of-censoring-weighting, which may be particularly useful in settings with rare outcomes. The cvAUC package is a tool kit to evaluate area under the ROC curve estimators when using cross-validation. The subsemble package was developed based on a new approach to ensembling that fits each algorithm on a subset of the data and combines these fits using cross-validation. This technique can be used in data sets of all size, but has been demonstrated to be particularly useful in smaller data sets. A new implementation of super learner can be found in the Java-based h2oEnsemble package, which was designed for big data. The package uses the H2O R interface to run super learning in R with a selection of prespecified algorithms.

Another TMLE package is ltmle, primarily authored by Joshua Schwab (Berkeley). This package mainly focuses on parameters in longitudinal data structures, including the treatment-specific mean outcome and parameters defined by a marginal structural model. The package returns estimates for TMLE, g-computation, and estimating-equation-based estimators.

The text above is a modified excerpt from the chapter "Targeted Learning for Variable Importance" by Sherri Rose in the forthcoming Handbook of Big Data (2015) edited by Peter Buhlmann, Petros Drineas, Michael John Kane, and Mark Van Der Laan to be published by CRC Press.

The most common R error messages

R has something of a reputation for generating, shall we say, obscure error messages like this:

Error in model.frame.default(formula = y ~ female + DNC + SE_region +  : could not find function "function (object, ...) \nobject"

One tip for dealing with error messages is to ignore everything between "Error in" and the colon: unless you are running a function that you wrote yourself, only the error message at the end is likely to be useful. If you're still stuck, another tip is to ask for help on Stackoverflow.com using the [r] tag, where you'll find more than 20,000 questions about R error messages.

Noam Ross has analyzed these questions to find the most commonly asked-about R error messages. Naturally, he used the stackr R package to interrogate the StackOverflow API, and downloaded around 10,000 error messages. He then used a regular expression to break the questions down into trigrams (sequences of 3 works) to be able to count which were the most common. On that basis, the most common types of error messages were:

1. "could not find function" errors, usually caused by typos or not loading a required package
2. "Error in if" errors, caused by non-logical data or missing values passed to R's "if" conditional statement
3. "Error in eval" errors, caused by references to objects that don't exist
4. "cannot open" errors, caused by attempts to read a file that doesn't exist or can't be accessed
5. "no applicable method" errors, caused by using an object-oriented function on a data type it doesn't support
6. "subscript out of bounds" errors, caused by trying to access an element or dimension that doesn't exist
7. package errors caused by being unable to install, compile or load a package.

Noam's full analysis is at the link below. In addition to providing insights about R's error messages, the trigram method he uses will be useful to anyone who needs to do frequency analysis on unstructured data.

Noam Ross (github): Common errors in R: An Empirical Investigation

Because it's Friday: Dominoes and Dance

There are plenty of great domino cascades on YouTube, but this is one of the best. There are a couple of obvious cuts, but using multiple takes is more than made up for by the perfect coordination with the A-Trak tune. Enjoy:

 

That's all for this week. See you again on Monday!

Replay: Reproducible data analysis with the checkpoint package

Thanks to all who attended my webinar earlier this week, Reproducibility with Revolution R Open and the Checkpoint Package. If you missed the live session, you can catch up with the slides and video replay which I've embedded below.

If you just want to check out the demo of the checkpoint package, it starts at 18:30 in the video below. If you want to follow along at home, you can download the demo script here.

  

Revolution Analytics webinars: Reproducibility with Revolution R Open and the Checkpoint Package

Review of "Hands-On Programming with R"

by Joseph Rickert

There have been well over a hundred books on R published within the last ten years. Most of these texts with titles like “Introduction Statistics with R” or “Time Series with R” offer the reader a way to jump right in and perform some concrete statistical analysis using R’s myriad built-in functions and extensive visualization features. And, while it is true that some R books appear to be little more than a rehash of basic documentation, there are nevertheless scores of carefully written texts from experts that not only illuminate some area of statistics but also demonstrate some good R programming as well. In no small way, I believe these works have contributed to the R’s popularity and growth by providing quality application level documentation.

Comparatively few books, however, are focused on teaching R programming itself. So it was a pleasant surprise when a copy of  Garrett Grolemund’s “Hands-On Programming with R: Write Your Own Functions and Simulations” (O’Reilly 2015) came my way. This is a superb book: well conceived, unusual in the choice of material and sufficiently streamlined (185 pages not including the appendices) to make it a non-stop beginning-to-end read.

At the very beginning Garrett says:

I want to help you become a data scientist, as well as a computer scientist, so this book will focus on programming skills that are most related to data science.

These skills have to do with solving what Garrett refers to as the "logistical problems" of data science. In the context of the R language, they include acquiring data, manipulating R objects, constructing custom functions, negotiating the R environment and above all, writing vectorized code.

Given the ambitious agenda, "Hands-On Programming with R" starts surprisingly slowly with arithmetic, assignment, useful R functions and basic housekeeping chores: getting help and looking for packages. Then, still slowly and deliberately the text discusses R objects, atomic vectors, data types and data structures. 48 pages in and Garrett is still lingering on attributes. But this discussion is more sophisticated than most authors attempt. The presentation of type, attributes and class, in particular the insight that the concept of class follows directly from attributes, is meant to cultivate a programmer's mindset.

Around page 65 when Garrett gets excited about subsetting the pace really picks up. If you are hooked and still reading like I was by page 112 you will have acquired a working knowledge of scoping rules and environments and be ready for the beguilingly lucid discussion of the S3 class system that begins on page 139. Even if you are an experienced R programmer you may want to borrow a copy of the book and read this. If you really know your stuff, you may not learn anything new, but I bet you will be hard pressed to do a better job of explaining S3 classes to someone else.

After S3 the text moves to considering loops as a prelude to its presentation of vectorized code. This section, which is really the final destination of the book, is exceptionally well done. First, vectorized code is characterized as code that takes advantage of three great features of the R language: fast logical tests, powerful subsetting operations and a multitude of built-in functions that permit element-wise execution. Then the text demonstrates how to put these ideas into practice.

As you can gather, I was impressed by the conceptual formulation of the material. However, the real strength of the text is its sharp presentation of essential elements of the  R language through a well-crafted, extended example that forms the spine of the book. "Hands-On Programming with R" is indeed a “hands-on” text that guides and challenges the reader to write good R code. A reader / coder who makes it to the end will have worked through several refinements of a small collection of functions that implement a fairly complex slot machine simulation. This example significantly raises the bar for selecting code examples in any R book. The simulation is rich enough to illustrate all of the R features presented in the text while allowing for refinement and polishing as the final form of the slot machine takes shape. The whole presentation is very tight. Garrett tells a pretty good story. During the final vectorized-code chapter I found myself reading with the delight of anticipation: “Just how is he going to make this code better?”.

I should also mention that the book is notable for what it does not include. This might be the first R book I have encountered that doesn’t develop any statistical models. Not a single regression is fit and there are no plots to speak of (3 histograms and a scatter plot). Certainly, this is the only R book I have come across that mentions data science in the preface that is not replete with Random Forest models and the like. Presumably, all of this will show up in the follow up book that Garrett promises in the preface.

"Hands-On Programming with R" presents but one carefully thought-through trajectory of many possible R language excursions. It is not to be compared with Hadley Wickham’s encyclopedic "Advanced R" and it contains only a fraction of the material you can find in Norm Matloff’s "The Art of R Programming". But, having worked through "Hands-On Programming with R" both of these texts should be accessible.

Garrett's book is a good read: a technical story with a plot and a few surprises that could help anyone starting out with the R language learn to write some pretty slick code.

Participate in the 2015 Rexer Data Mining Survey

Since 2007, Karl Rexer has been collecting data on the tools, skills and practices of statisticians and data miners. Over the years, his semi-annual Data Miner Survey has expanded in scope, and now includes research on the topics of data science, big data, and analytics applications on business. Unlike research conducted by big-name analyst firms, Karl shares the results of his surveys with everybody. You can read our review of the results of the 2013 Data Miner Survey here, which among many other things showed the rapid growth of R:

The 2015 Data Mining includes several new questions this year. In addition to being able to share your preferred tools and analytical methods, you'll also be able to rank the key skills of a data scientist, and even the preferred titles used by today's data professionals.

To be heard in this year's survey, follow the link below and use code AL79R5 . Survey results will be unveiled in late September at the 2015 Boston Predictive Analytics World event.

Rexer Analytics: Data Miner Survey 2015

R Package 'smbinning': Optimal Binning for Scoring Modeling

by Herman Jopia

What is Binning?

Binning is the term used in scoring modeling for what is also known in Machine Learning as Discretization, the process of transforming a continuous characteristic into a finite number of intervals (the bins), which allows for a better understanding of its distribution and its relationship with a binary variable. The bins generated by the this process will eventually become the attributes of a predictive characteristic, the key component of a Scorecard.

Why Binning?

Though there are some reticence to it [1], the benefits of binning are pretty straight forward:

• It allows missing data and other special calculations (e.g. divided by zero) to be included in the model.
• It controls or mitigates the impact of outliers over the model.
• It solves the issue of having different scales among the characteristics, making the weights of the coefficients in the final model comparable.

Unsupervised Discretization
Unsupervised Discretization divides a continuous feature into groups (bins) without taking into account any other information. It is basically a partiton with two options: equal length intervals and equal frequency intervals.

Equal length intervals

• Objective: Understand the distribution of a variable. 
• Example: The classic histogram, whose bins have equal length that can be calculated using different rules (Sturges, Rice, and others).
• Disadvantage: The number of records in a bin may be too small to allow for a valid calculation, as shown in Table 1.

Table 1. Time on Books and Credit Performance. Bin 6 has no bads, producing indeterminate metrics.

Equal frequency intervals

• Objective: Analyze the relationship with a binary target variable through metrics like bad rate.
• Example: Quartlies or Percentiles.
• Disadvantage: The cutpoints selected may not maximize the difference between bins when mapped to a target variable, as shown in Table 2

Table 2. Time on Books and Credit Performance. Different cutpoints may improve the Information Value (0.4969).

Supervised Discretization

Supervised Discretization divides a continuous feature into groups (bins) mapped to a target variable. The central idea is to find those cutpoints that maximize the difference between the groups.
In the past, analysts used to iteratively move from Fine Binning to Coarse Binning, a very time consuming process of finding manually and visually the right cutpoints (if ever). Nowadays with algorithms like ChiMerge or Recursive Partitioning, two out of several techniques available [2], analysts can quickly find the optimal cutpoints in seconds and evaluate the relationship with the target variable using metrics such as Weight of Evidence and Information Value.

An Example With 'smbinning'

Using the 'smbinning' package and its data (chileancredit), whose documentation can be found on its supporting website, the characteristic Time on Books is grouped into bins taking into account the Credit Performance (Good/Bad) to establish the optimal cutpoints to get meaningful and statistically different groups. The R code below, Table 3, and Figure 1 show the result of this application, which clearly surpass the previous methods with the highest Information Value (0.5353).

# Load package and its data 
library(smbinning) 
data(chileancredit) 
# Training and testing samples 
chileancredit.train=subset(chileancredit,FlagSample==1) 
chileancredit.test=subset(chileancredit,FlagSample==0) 
# Run and save results 
result=smbinning(df=chileancredit.train,y="FlagGB",x="TOB",p=0.05) 
result$ivtable
 
# Relevant plots (2x2 Page) 
par(mfrow=c(2,2)) 
boxplot(chileancredit.train$TOB~chileancredit.train$FlagGB, 
horizontal=T, frame=F, col="lightgray",main="Distribution") 
mtext("Time on Books (Months)",3) 
smbinning.plot(result,option="dist",sub="Time on Books (Months)") 
smbinning.plot(result,option="badrate",sub="Time on Books (Months)") 
smbinning.plot(result,option="WoE",sub="Time on Books (Months)")

 

Table 3. Time on Books cutpoints mapped to Credit Performance.

Figure 1. Plots generated by the package.

In the middle of the "data era", it is critical to speed up the development of scoring models. Binning, and more specifically, automated binning helps to reduce significantly the time consuming process of generating predictive characteristics, reason why companies like SAS and FICO have developed their own proprietary algorithms to implement this functionality on their respective software. For analysts who do not have these specific tools or modules, the R package 'smbinning' offers an statistically robust alternative to run their analysis faster.

 

For more information about binning, the package's documentation available on CRAN lists some references related to the algorithm behind it and its supporting website some references for scoring modeling development. 

References
[1] Dinero, T. (1996) Seven Reasons Why You Should Not Categorize Continuous Data. Journal of health & social policy 8(1) 63-72 (1996).
[2] Garcia, S. et al (2013) A Survey of Discretization Techniques: Taxonomy and Empirical Analysis in Supervised Learning. IEEE Transactions on Knowledge and Data Engineering, Vol. 25, No. 4, April 2013.

Tomorrow, 10AMPT: Live webinar on "checkpoint" package

A quick heads-up that tomorrow (Tuesday) at 10AM Pacific Time I'll be giving a live (and free) webinar: Reproducibility with Revolution R Open and the Checkpoint Package. I'll show you how to use the latest version of the checkpoint package to write R scripts without having to worry about unexpected package updates on CRAN. Here's the full abstract:

Good data analysis is reproducible. If someone else can’t independently replicate your results from your data, the consequences can be severe. With R, a challenge for reproducibility is the ever-changing package ecosystem: it's all too easy to develop an R script using packages, only to find collaborators will download later versions of those packages when they attempt to reproduce your results, making outcome unpredictable.

In this talk David will introduce the Reproducible R Toolkit, and the "checkpoint" package, included with Revolution R Open, and describe some best practices for writing reliable, reproducible R code with packages.

I'll also be including a demo of the checkpoint package in action. To sign up for the webinar (which also gets you an email when the replay and slides are available for download), follow the link below.

Revolution Analytics webinars: Reproducibility with Revolution R Open and the Checkpoint Package

Because it's Friday: Deep into the Mandelbrot Set

Think of a number. Now multiply it with itself. Now add the number you first thought of. Multiply the result with itself. Again add the number you first thought of. Keep repeating that process. If the number you have keeps getting larger and larger, the number you first thought of is not part of the Mandelbrot Set. (For example, starting with 0.5 you get 0.75 -> 1.0625 -> 1.628 -> 3.153 ... and so on off to infinity, so 0.5 is not part of the Mandelbrot Set.) 

You can also let the number you think of be a complex number of the form (a, b), where "a" is called the "real part" and b is called the "imaginary part". The only difference is that (a, b) multiplied by itself is (a²-b², 2ab), so the next step in the sequence is (a²-b²+a, 2ab+b), but otherwise the game is much the same. (This Jonathan Coulton song explains the rules in lyrical form.)

The trick is that you can use (a, b) as the X and Y coordinates on a 2-dimensional chart. Color the chart black if (a, b) is in the Mandelbrot Set, and otherwise color it according to the number of steps it takes to get the complex number at those coordinates to get above size 2 (5 steps in our example for 0.5 above). Then, you get the beautiful and familiar representation of the Mandelbrot Set. (You can easily write a program to do this in the R language, since it support complex arithmetic natively.)

One of the amazing things about the Mandelbrot set is that it's infinitely complex: as far as you zoom into its crenellated edges, it never smooths out. This video enlarges the Mandelbrot Set as much as expanding a proton to the size of the Universe -- and then billions and billions of times more. 

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

NYC is a city that does sleep, a bit

The On Broadway project collected more than 600,000 photographs taken near Broadway in New York City during a six-month period in 2014. If you're in New York, you can explore the images in an interactive installation at the New York Public Library though the end of this year. You can also explore them in your browser using this online app. 

In addition to being part of the team that created the exhibit, artist Lev Manovich used R to analyze some of the data behind the project. For example, the timestamps of the photographs can be used to measure relative activity of the denizens of New York throughout the day (at least, those prone to post photos to social media). Here's a chart of activity measured by number of photos posted throughout the 24 hours of the day (clearly done as a polar coordinates chart with ggplot2):

From this data, it's clear there's not a lot happening between 3AM and 5AM. Of course, activity levels aren't constant across all of Manhattan: here's a comparison of "Tourist New York" (Broadway south of 110th St) on the left, and uptown Manhattan (above 110th street) on the right:

Wherever you look, though, it seems New York does sleep sometimes! You can explore more analysis of the data in Lev Manovich's blog post linked below.

Software Studies Initiative: "A City That Never Sleeps?" - new data and analysis from "On Broadway" project