« February 2013 | Main | April 2013 »
On Tuesday and Wednesday this week, the US Supreme Court heard two landmark cases related to same-sex marriage rights in this country. Simultaneously, to show solidarity for the cause of civil rights for gay members of society, thousands of Facebook users changed their profile pictures to a red equality symbol at the suggestion of the Human Rights Campaign. Most profile pictures looked like the one shown at right, but there were also many variants used as well. (By the way, here's R code use can use to generate the red equality logo, with thanks to Josef Fruehwald.)
The change in Facebook profiles provided some interesting information about the support of advances in gay rights across the USA. Some clever data science from Facebook's Eytan Bakshy reveals that there was indeed a significant increase in the rate of profile picture changes during the oral hearings:
(The chart above is the result of a time-series analysis using the R language, after extracting the normal day-night cycle of profile updates.) Moreover, by looking at the rate of changes of Facebook profiles across every county in the USA in a choropleth, we see those regions where support for gay rights is the strongest:
While it's unlikely that the Justices will consult Facebook profiles as part of their deliberations, the proceedings did seem to be favourable to the advance of gay rights. (Especially for the second case regarding DOMA, which if overturned will be a hugely positive outcome for thousands of couples already married, including me and my husband.) But for those of us for whom this cause means so much personally, it was truly heartwarming to see that sea of red from allies all across the county. What strikes me so much about the map above is that support is coming from all across the country — not one county is a totally unsupportive white. And that's something that all of America can be proud of.
Facebook (Eytan Bakshy): Showing Support for Marriage Equality on Facebook (via Sean J. Taylor)
by Joseph Rickert
When talking with data scientists and analysts — who are working with large scale data analytics platforms such as Hadoop — about the best way to do some sophisticated modeling task it is not uncommon for someone to say, "We have all of the data. Why not just use it all?" This sort of comment often initially sounds pragmatic and reasonable to almost everyone. After all, wouldn’t a model based on all of the data be better than a model based on a subsample? Well, maybe not — it depends, of course, on the problem at hand as well as time and computational constraints. To illustrate the kinds of challenges that large data sets present, let’s just look at something very simple using the airlines data set from the 2009 ASA challenge.
Here are some of the results for a regression of ArrDelay on CRSDepTime with a random sample of 12,283 records drawn from that data set:
# Coefficients:
# Estimate
Std. Error t value Pr(>|t|)
# (Intercept)
-0.85885 0.80224 -1.071
0.284
# CRSDepTime 0.56199
0.05564 10.100 2.22e-16
# Multiple
R-squared: 0.008238
# Adjusted R-squared: 0.008157
And here are some results from the same model using 120,947,440 records:
#Coefficients:
# Estimate
Std. Error t value Pr(>|t|)
# (Intercept)
-2.4021635 0.0083532 -287.6 2.22e-16 ***
# CRSDepTime 0.6990404
0.0005826 1199.9 2.22e-16 ***
# Multiple
R-squared: 0.01176
# Adjusted R-squared: 0.01176
More data data didn’t yield an obvously better model! I don’t think anyone would really find this to be much of a surprise. We are dealing with a not very good model to begin with. Nevertheless, the example does provide the opportunity to investigate how estimates of the coefficients change with sample size. This next graph shows the coeffients of the slope plotted against sample size with sample sizes ranging from 12,283 to 12,094,709 records. Each regression was done on a random sample that includes about 12,000 points more than the previous one. The graph also shows the standard estimate for the confidence interval for the coefficient at each point in red. Notice that after some initial instability, the coefficient estimates settle down to something close to the value of beta obtained using all of the data.
The rapid approach to the full-data-set value of the coefficient is even
more apparent in the following graph that shows the difference between the
estimated value of the beta coefficient at each sample and the value obtained using
all of the data. The maximum difference from the fourth sample on is 0.07. This is pretty close indeed. In cases like this, if you
believed that your samples were representative of the entire data set, working
with all of the data to evaluate possible models would be a waste of time an possibly counterproductive.
I am certainly not arguing that one never wants to use all of the data. For one thing, when scoring a model or making predictions the goal is to do something with all of the records. Moreover, in more realistic modeling situations where there are thousands of predictor variables 120M observations might not be enough data to conclude anything. A large model can digest degrees of freedom very quickly and severely limit the ability to make any kind of statistical inference. I do want to argue, however, that with large data sets the ability to work with random samples of the data confers the freedom to examine several models quickly with considerable confidence that results would be decent estimates of what would be obtained in using the full data set.
I did the random sampling and regressions in my little example using functions from Revolution Analytics RevoScaleR package. Initially, all of the data was read from the csv files that comprise the FAA data set into the binary .xdf file format that is used by the RevoScaleR package. Then the random samples were selected by using the rxDataStep function of RevoScaleR. This function was designed to quickly manipulate large data sets. The code below reads a record, draws a random number with a value between 1 and 9999 and assigns it to the variable urns.
rxDataStep(inData = working.file, outFile = working.file, transforms=list(urns = as.integer(runif(.rxNumRows,1,10000))), overwrite=TRUE)
Random samples for each regression were drawn by looping throught the appropriate values of the variables. Notice how the call to R’s runif() function happens within the transforms parameter of rxDataStep. It took about 33 seconds to do the full regression on my laptop which made it feasible to undertake the extravagent number of calculations necessary to do the 1,000 regressions in a few hours after dinner.
I think there are three main take aways from this exercise:
To my way of thinking, the single most important capability to implement in any large scale data platform that is going to support sophisticated analytics is the ability to quickly construct, high quality random samples.
You may think your simulations take a long time to run, but this one's been running for 200 years: (click to enlarge)
The chart above, by Tim Granger, shows just 14 years of data from the Oxford / Cambridge "Bumps" races, a unique tournament designed to find the best rowing crew in the Oxbridge system. Because the race rivers are so narrow and winding, a traditional race with multiple competitors in parallel lanes is out of the question. Instead, boats line up in single file, and if during the course of the race a boat manage to bump into (or come alongside) the one ahead, those two boats will swap positions for the next race. The boat at the head of the pack at the end of the multi-race tournament is declared the winner ... and the same ordering is held for the next tournament in six months time. This has been going on since eighteen-fracking-fifteen!
Sure it's really more of a bubble sort than a simulation, and teams do change ability over time so the sorting isn't stable, but still: that's a really neat way to determine the best Oxbridge rowing team, given the constraints. This Quora post by Ferenc Huszár explains more of the details.
That's all for this week — enjoy your weekend!
If you're laying down a friendly bet on the March Madness games or just tweaking your fantasy roster, this NCAA Data Visualizer by Rodrigo Zamith will be a boon. Just choose two teams to compare head-to-head, choose an attribute to compare them on. You can look at more than a dozen invividual player attributes (e.g. points scored, assists, 3-point shots made) or over 20 team attributes (e.g. fouls, wins, and turnovers), and compare against other top players or NCAA team averages.
For example, here's how the top players match up for individual points scored for tomorrow's Oregon/Saint Louis game:
All of the data used was scraped from the NCAA website for the 2012-2103 season. Rodrigo created the visualizations using the R programming language, which means he can used some advanced techniques like the boxplots above. Rather than just showing the median or highest/lowest scores for each player, he can show all of each player's games, and highlight the bulk (the middle 50%) of their games using the box. But you can still see, for example, that Saint Louis's Kwamain Mitchell had a couple of hot streaks last season and was in fact their highest scorer last season — a fact that you might miss just looking at averages.
You can play with the NCAA data visualizer and choose your own teams and stats to compare at the link below.
Rodrigo Zamith: Visualizing Season Performance by NCAA Tournament Teams (via Myles Harrison)
by Joseph Rickert
If you type ?Distributions at the R console you get a list of the 21 probability distributions included in the stats package that ships with base R. The same list appears in the Introduction to R Manual on CRAN and in most of the many fine introductory books available for the R language. These are indeed fundamental distributions, sufficient for most elementary work in probability and statistics. The fact that the R functions implementing these distributions all follow same syntax greatly eases a beginner's task of trying to get some useful work done with a minimum of memorization.
The following figure shows plots of the cumulative distribution pgamma()and probability density function dgamma() along with the histogram of random draws from a gamma distribution rgamma(2,2)with shape and scale parameters both set to 2.
However, if a person isn’t familiar with how information about R is organized on CRAN, he or she might conclude: “that’s it” or most of it anyway, with respect to R and probability distributions. Imagine the surprise then of a person with such modest expectations about R’s probability distributions accidently stumbling into the overgrown garden of R’s Probability Distributions Task View. I think my first reaction was kind of glazed over inability to take it all in.
However, if you just let your eyes relax and pick out a flower with which you are familiar, binomial for example, you can see that the chief gardener Christophe Dutang, listed as the maintainer of the Task View, and the eight individuals whom acknowledges have done a remarkable job of organizing the distributions according to their genus (discrete or continuous), species (binomial in this case) and variety (truncated binomial and zero inflated binomial). I can’t imagine the number of volunteer hours took to assemble this page, and keeping it up to date can’t be easy either. I spent a half hour or so just trying to count the distributions. Not counting copulas, random matrices and other exotica I came up with 31 discrete, 133 continuous and 9 mixture distributions. Others may count more or less depending on how they group things together. It seems as if few people outside of the folks at Wikipedia have given much thought to the taxonomy of probability distributions and only Mathematica 9 which includes 130 probability distributions comes close to cultivating so many distributions in one coherent system. (To be fair, the online documentation for SAS, Matlab and SPSS is so distributed that it is difficult to determine how many probability distrbutions have ben implemented in these software packages.)
While the Probability Distributions Task view may be the place to start for information about probability distributions, the complete R documentation is itself an open ended, organic system that depends on the communication style of package authors and the experiences of everyone who leaves a record of their attempts to work with probability distributions.
The entire ecosystem of R documentation for a probability distribution function starts with the command line help ( e.g. ?pgamma) and the package pdf on CRAN that includes the function, but may also include, vignettes, external web pages, blog posts and questions and discussions on help bulletin boards such as the R mailing lists and StackOverflow. For some typical examples, consider that the actuar package from Vincet Goulet et al. which provides a number of distributions of interest to acturies has six vignettes, while Thomas Yee's VGAM package for Vector Generalized Linear and Additive Models, a source for many R probability distributions, has a web page as well as a vignette.
John D. Cook’s clickable diagram for elementary probability distributions is hosted on his private website while and the paper by Delignette-Muller et al. on fitting distributions with R’s fitdistrplus package is hosted on an academic website. Mage's post from December 2011 on fitting distributions in R is an example of the many blog posts that deserve a second look.
As a final example of how the community comes to play a part of the extended documentation for R, consider my attempt get a handle on the Cauchy distribution. Here I ran the below and got four very different looking plots. This is not unexpected given that I’m working with random draws from a probability distribution for which both the mean and variance are not defined. But why only two bins for the histograms?
Well, I wasn’t the first person to pause for a moment over this. Someone recently asked
this question on StackOverflow and received some good advice.
Hats off and thank you to everyone involved in cultivating R’s garden of probability distributions
# Cauchy plots n <- 10000 location <- -1 scale <- 4 par(mfrow=c(2,2)) # Make four plots for(i in 1:4){ y <- rcauchy(n, location, scale) hist(y, freq = FALSE, col = rainbow(6), main="random draw from rcauchy(-1,4)") fd <- function(y)dcauchy(y,shape,scale) curve(fd, col = "black", add = TRUE,lwd=2) rug(y,col="grey") }
When I lived in Paris years ago, I worked near Gare du Nord, but my friend Jenny lived near République. If we wanted to meet up after work, we'd just meet halfway along the Orange Métro line, around Gare de l'Est. Easy. Since that's within walking distance we wouldn't actually take the Métro, but Métro stations are useful waypoints in Paris even if you're travelling by foot.
But when we wanted to meet with our other friend Eric who would travel from Porte de Versailles, it was tricky to find a mutually-agreeable place to meet, so that we'd all be taking an approximately equally long Métro ride. If only we'd had Jean-Robert's Metro Meeting Point app, it would have made the decision an easy one.
The app works by representing the Métro network as a connected graph, and finding the station that minimises the maximum number of stops from each of 2 or 3 or more starting points. This is possible thanks to the fact that the RATP (the Paris transit agency) has published open data on Métro stops and their geographic location. All it took was for Jean-Robert to write a script in the R language to import the data and create the network as an R object, and then use the Djikstra algorithm and minimax optimization to find the best meeting point. The web interface for the app was done using RStudio's Shiny.
For complete details on how Jean-Robert created this cool and useful app from open data, check out his blog post linked below. As for me, I'll be daydreaming about Spring evenings in a Paris café with a kir royale in hand.
Jean-Robert blog: Optimal Meeting Point on the Paris Metro
Follow David on Twitter: @revodavid
