Revolutions

Statistics vs Data Science vs BI

As someone who trained as a statistician, I've always struggled with that title. I love the rigor and insight that Statistics brings to data analysis, but let's face it: Statistics — the name — has always had a bit of a branding problem. Telling someone I was a statistician was more likely to conjure up images of me counting runs at a baseball (or cricket) game than pursuing serious science. And the image of what Statistics ideally is about — collaborative, interactive, applied, fun — was too often subsumed by the stereotype image — isolated, actuarial, ivory tower, report driven. (And hey, even actuaries can be fun sometimes.)

That's why I'm a fan of the term "data scientist" — it embodies everything that Statistics always should be, without the baggage and tradition of the term "statistician". So I enjoyed participating in yesterday's Kalido webinar "Data Scientist: Your Must-Have Business Investment Now" where I could make the following contrast between the images of Statisticians and Data Scientists:

(A quick aside on the "Data Size" row above: while the unstructured or unaggregated data source data that data scientist work with can be in the terabytes range or even large, by the time it's cleaned and prepared for statistical modeling, a file in the gigabytes range is even more typical — even at "Big Data" companies like Facebook. This is a topic I cover in more detail in my recent Strata talk on real-time predictive analytics.)

So bottom line: while I am a statistician, and I love Statistics dearly, I do prefer to call myself a Data Scientist today, because it better represents to me what Statistics really is to me (if that makes sense). And that's certainly not to diminish the achievements of those who do call themselves Statistician. In particular, I want to recognize George Box: a true hero of mine, coiner of the idiom "all models are wrong, but some are useful", and one of the nicest people I ever met, who sadly passed away in March.

On the other hand, I have no qualms about making a competitive comparison between Data Science and Business Intelligence:

You can get the details of how I differentiate Statistics and Data Science and BI, and hear other perspectives on Data Science from fellow data scientists Carla Gentry and Gregory Piatetsky in the slides and replay of the webinar provided by Kalido at the link below.

 

Kalido: Data Scientist: Your Must-Have Business Investment NOW

Trevor Hastie presents glmnet: lasso and elastic-net regularization in R

by Joseph Rickert

Even a casual glance at the R Community Calendar shows an impressive amount of R user group activity throughout the world: 45 events in April and 31 scheduled so far for May. New groups formed last month in Knoxville, Tennessee (The Knoxville R User Group: KRUG) and Sheffield in the UK (The Sheffield R Users). An this activity seems to be cumulative. This month, the Bay Area R User’s Group (BARUG) expects to hold its 52nd and 53rd meet ups while the Sydney Users of R Forum (SURF) will hold its 50th. Everywhere R user groups are sponsoring high quality presentations and making them available online, but the Orange County R User Group is pushing the envelope with respect to sophistication and reach. Last Friday, I attended a webinar organized by this group where Professor Trevor Hastie of Stanford University presented Sparse Linear Models with demonstrations using GLMNET. This was a world-class presentation and quite a coup for Orange County to have Professor Hastie present.

The glmnet package written Jerome Friedman, Trevor Hastie and Rob Tibshirani contains very efficient procedures for fitting lasso or elastic-net regularization paths for generalized linear models. So far the glmnet function can fit gaussian and multiresponse gaussian models, logistic regression, poisson regression, multinomial and grouped multinomial models and the Cox model. The efficiency of the glmnet algorithm comes from using cyclical coordinate descent in the optimization process and from Jerome Friedman's underlying Fortran code.

Although Professor Hastie’s presentation was primarily concerned with fitting models for the wide problem (the number of explanatory variables is much larger than the number of observations) the lasso and elastic-net algorithms are just as applicable to data sets with large numbers of observations. It is likely that in the future we will see glmnet implementations for variable selection on datasets with thousands of variables and hundreds of millions of observations. The following graph shows the regularization paths for the coefficients of a model fit the HIV data from one Professor Hastie’s examples.

Each curve represents a coefficient in the model. The x axis is a function of lambda, the regularization penalty parameter. The y axis gives the value of the coefficient. The graph shows how the coefficients “enter the model” (become non-zero) as lambda changes. The following code, based on an example from the webinar, produces the plot and also shows how easy it is to perform cross-validation.

library(glmnet)	    # load the package	
load("hiv.rda")	    # HIV data
class(hiv.train)    # The data are stored as a list 
names(hiv.train)    # The names of the list elements are x and y   
dim(hiv.train$x)    # The explanatory data consists of 704 observations of
		    # 208 binary mutation variables
head(hiv.train[[1]])	# Look at the explanatory data
head(hiv.train[[2]])    # Look at the response data: changes in susceptibility to antiviral drugs
fit=glmnet(hiv.train$x,hiv.train$y)	# fit the model
plot(fit,xvar="lambda", main="HIV model coefficient paths") # Plot the paths for the fit
fit		    # look at the fit for each coefficient
#
cv.fit=cv.glmnet(hiv.train$x,hiv.train$y)	# Perform cross validation on the fited model
plot(cv.fit)	# Plot the mean sq error for the cross validated fit as a function
                # of lambda the shrinkage parameter
		# First vertical line indicates minimal mse
		# Second vertical line is one sd from mse: indicates a smaller model
		# is "almost as good" as the minimal mse model
tpred=predict(fit,hiv.test$x)	# Predictions on the test data
mte=apply((tpred-hiv.test$y)^2,2,mean)	# Compute mse for the predictions
points(log(fit$lambda),mte,col="blue",pch="*")	# overlay the mse predictions on the plot
legend("topleft",legend=c("10 fold CV","Test"),pch="*",col=c("red","blue"))

Created by Pretty R at inside-R.org

Don’t be content with this partial example. Professor Hastie and The Orange County R User Group have graciously made the slides, code and data available at this link. The webinar is well worth watching in its entirety.

 

As you might expect, Professor Hastie gives a masterful presentation: lucid, clear and succinct. This is inspite of the fact that Professor Hastie begins the presentation by commenting that it was his first webinar ever and that he was a little uncomfortable talking to his screen. (I think anyone who has ever given a webinar can relate to this: you talk to the screen and no energy from the audience comes back. Nothing is more disruptive to efforts to be enthusiastic than silence.) Nevertheless, Professor Hastie presents a difficult topic with a clarity that carries his audience along, and he is completely unphased by the inevitable glitch. Watch how he handles the upside down slide. You can download his slides, R scripts and data from the link below.

Trevor Hastie: Sparse Linear Models, with demonstrations using GLMNET

SAS Big Data Analytics Benchmark (Part Two)

by Thomas Dinsmore

On April 26, SAS published on its website an undated Technical Paper entitled Big Data Analytics: Benchmarking SAS, R and Mahout.  In the paper, the authors (Allison J. Ames, Ralph Abbey and Wayne Thompson) describe a recent project to compare model quality, product completeness and ease of use for two SAS products together with open source R and Apache Mahout.

This is the second post in a two-part series.  Last week, I covered simple mistakes and errors. this week's post will cover the authors' methodology and findings.

I have contacted the authors and asked if they plan to release their data to the community.  To date, I have received no response. 

Methodology

The methodological problems in the SAS benchmark paper are too numerous to cite comprehensively, but I'll review the most glaring issues.

Choice of Products to Test:  Although the authors compiled the list of algorithms to be tested through a survey of SAS users, no single SAS product supports all three.  R supports all three, as does Mahout (with a workaround). 

The choice of products for comparison is strangely "apples and oranges".  Mahout is an early-stage analytics library; R is a mature analytic programming environment; SAS Rapid Predictive Modeler is a semi-automated ensemble modeling framework, and SAS High Performance Analytics Server is a high-end analytics platform designed for use in an MPP computing environment.  Comparing the two open source products with two high-end SAS applications on “ease of use” is a straw man comparison.  One would expect the two high-end SAS applications to be easier to use; that’s why SAS charges a premium for them.    

Oddly enough, SAS chose to omit its own product SAS/STAT from the benchmark. As an analytic programming language, SAS/STAT offers a more apt comparison to R, and it is far more widely used than either of the two high-end SAS products.

Analytic Workflow:  In the real world, data resides in external sources, and must be loaded into each of the respective platforms; analysts have choices when marshalling and loading data, and can do so in a way that best suits the application. The authors used SAS to prepare the data for use in all four products, starting with structured data sets loaded into Hadoop. This biases the "ease of use" assessment toward SAS.  

Structured data stored in Hadoop is not "native" to Hadoop but structured externally and loaded into the Hadoop file system; the authors did not assess the effort needed to do this, which biases the "ease-of-use" assessment towards SAS High Performance Analytics Server, which cannot operate directly on raw data. It also biases the assessment against Mahout, since the principal advantage of Mahout is its ability to work with data stored natively in Hadoop.  

The authors' decision to split the data into training and validation sets in SAS rather than using native splitting for each tool is puzzling. The stated reason for using this approach — that random seeds may differ among the tools — is fallacious; random samples of sufficient size are generalizable to one another even when produced with different random seeds.

Computing Environments: The computing environments used in this benchmark are vastly different. This is attributable in part to the widely different requirements of the products compared, but nevertheless there are some unexplained anomalies. There is no justification for using machines with very different memory profiles to host open source R and SAS Enterprise Miner; the two machines can and should be identical. By the same token, a five node Hadoop cluster is not equivalent to a sixty node MPP Greenplum appliance.  

Assumptions About User Skill: The authors tested the open source tools and the SAS high-end applications using a naïve "hands-off" approach that depends on "out of the box" functionality.  This is a trivial comparison that simply proves the obvious. When Deep Blue beats Kasparov at chess, it's news; when Deep Blue beats a three-year-old at chess, it's not news.  

Findings

The authors propose five main findings:

(1) The types and depth of classification models that could be run in SAS products and R outnumbered the options for Mahout. 

(2) The effort required by individual modelers to prepare a representative table of results was much greater in both R and Mahout than it was in SAS High-Performance Analytics Server and SAS Rapid Predictive Modeler for SAS Enterprise Miner.

(3) The object-oriented programming in R led to memory-management problems that did not occur in SAS products and Mahout; thus, the size of customer data that could be analyzed in R was considerably limited. 

(4) Overall accuracy was "comparable" across software in "most instances". 

(5) When the customer modeling scenarios were evaluated based on event precision and rank ordering, SAS HighPerformance Analytics Server models were more accurate. 

The first two points are true but trivial. Anyone with a smattering of knowledge understands that Mahout is an early stage project with limited functionality. The two SAS tools are high-end applications with a number of automated features. That's fine, there's nothing wrong with that for users willing to pay the license fees.  A reasonably competent R user can produce a comparable table of results with little effort; moreover, the table can be easily customized, which is not so easy to do with "out-of-the-box" tools.

The third point is silly. The authors attribute memory management issues encountered with R to its object-oriented nature. A more obvious explanation is that the authors deployed R on a machine with much less memory than the machine they used for SAS.

Points four and five relate to model quality, and deserve some analysis.

According to the authors, the overall misclassification rate (the inverse of overall accuracy) achieved by the four products is "comparable in most instances." They do not disclose the overall accuracy statistics, which leaves the reader in the dark about what they mean by "comparable" and which instances the overall accuracy was not comparable.

However, the SAS products produced models with higher precision and lift in the first decile. The authors choose to highlight this finding, reporting only the precision and lift statistics in the summary tables, while failing to report overall accuracy statistics.

If the models produced by SAS are more precise but not more accurate, it means that they were produced using techniques that minimize false positives at the expense of false negatives (since the overall error rate is the same). There are situations where this is a good thing — for example, in marketing campaigns where the goal is to optimize campaign ROI, the analytic goal is to target a list as precisely as possible. On the other hand, there are circumstances where the analyst must also be concerned about false negatives — for example, in a medical diagnostic situation where the cost of a false positive is "patient receives unnecessary dose" and the cost of a false negative is "patient dies."

While it is generally best to evaluate models using a range of measures and over the entire range of possible predictions (through a gains chart or ROC analysis), in the absence of information about the cost of errors, accuracy is the best single measure. And, as the authors admit, all of the products produced models of comparable accuracy.

But let's stipulate that increased precision is a good thing where models are equally accurate. How did the SAS tools produce the improved precision? As noted in the Discussion section on page 11, the SAS products accomplished this through use of well-known techniques:

• Oversampling the target event;
• Modification of model priors;
• Missing value imputation;
• Variable selection.

All of these methods are available in R; the authors simply chose not to use them. Which brings us back to a basic flaw in the authors' approach: it's news if SAS' automated tooling can beat a trained and reasonably competent analyst using R; it's not news if SAS' automated tooling can beat a passive and untrained analyst.

The reader should take note that SAS' automated tooling would also outperform a passive and untrained analyst using SAS/STAT. While all of the techniques detailed above can be implemented in SAS/STAT through the SAS Programming Language — a point recognized by the authors — they are not "automatic."

While there is a role for semi-automated model building tools in analytics, evidence from the marketplace suggests that working analysts in the real world prefer the flexibility and control offered by an analytic programming environment over "black-boxy" applications. This is true in the SAS user community as well as the broader community of analysts; users of SAS/STAT and the SAS Programming Language far outnumber users of the high end SAS applications included in this "benchmark."

Derek Norton, Andrie de Vries, Joseph Rickert, Bill Jacobs, Mario Inchiosa, Lee Edlefsen and David Smith all contributed to this post. 

Learn how to analyze data with R with Coursera's "Data Analysis" videos

If you didn't manage to catch Coursera's Data Analysis course, don't despair. Instructor Jeff Leek has made the course videos available on YouTube, which you can review at your leisure to learn how to plan, carry out, and communicate analyses of real data sets with R. (The course assumes you already have familiarity with R, so if you're new to the R language start with the videos from the Computing for Data Analysis course.) Begin with the Course Introduction for "Data Analysis":

 

The lecture videos are organized by week (the first is at the bottom of the list, so work from bottom to top):

• Week 1
• Week 2
• Week 3
• Week 4
• Week 5
• Weeks 6-7

By the way, "Data Analysis" was an extremely popular course: just over 100,000 students enrolled in the course, of whom more than half watched videos. 20,000 students completed quizzes, and 5,500 completed the data analysis assignments. Jeff Leek offers further reflections on the course in this Simply Statistics podcast.

Coursera hasn't yet scheduled another run of the course, but keep an eye out for the next one if you want to get the certificate. In the meantime, check out the videos linked below.

Jeff Leek (YouTube): Coursera Data Analysis

Stepwise Regression for Big Data with RevoScaleR

by Joseph Rickert

In a recent blog post, Revolution's Thomas Dinsmore announced stepwise regression for big data as a new feature of Revolution R Enterprise 6.2 that is scheduled for general availability later this month. Today, I would like to provide a simple example of doing stepwise regression with rxLinMod() (the RevoScaleR analog of lm()),  using a 100,000 row subset of the Million Song data set. (Look here for a description of the variables contained in the file.) This file is big enough to make RevoScaleR functionality interesting, but small enough so that it can also be processed with lm() and step().

The following code runs a stepwise regression with RevoScaleR's rxLinMod() and rxStepControl() functions.

# STEPWISE REGRESSION EXAMPLE WITH rxLinMod and rxStepControl
#
# Access the data
MSdir <- "C:/Users/Joseph/Documents/DATA/Million Song/MillionSong _XDF"
fileName <- "MS_ABVZ.xdf"
songs <- file.path(MSdir,fileName)
# Look at a summary of the data
rxGetInfo(songs,getVarInfo=TRUE)
# Specify the linear model
form <- formula(duration ~ artist.hotttnesss +
             artist.familiarity + track.7digitalid + release.7digitalid  +
	     year + mode + mode.confidence + key + key.confidence + time.signature + 
	     time.signature.confidence + start.of.fade.out + end.of.fade.in + tempo + loudness)
# Run the linear model				
system.time(mod <- rxLinMod(formula = form,
	                     data = songs,
			     blocksPerRead=10000))
 
summary(mod)             # Look at a summary of the model
# Specify the scope for the stepwise regression
scope <- list(
    lower = ~ loudness,
    upper = ~ artist.hotttnesss + artist.familiarity + track.7digitalid + 
	          release.7digitalid  + year + mode + mode.confidence + key + 
			  key.confidence + time.signature + time.signature.confidence + 
			  start.of.fade.out + end.of.fade.in + tempo + loudness)
 
# Set up the variable selection parameter
varsel <- rxStepControl(method = "stepwise", scope = scope)
# Run the stepwise regression
system.time(rxlm.step <- rxLinMod(form, data = songs,
 	              blocksPerRead=100000,
	              variableSelection = varsel,
                      verbose = 1, 
		      dropMain = FALSE, 
		      coefLabelStyle = "R"))

Created by Pretty R at inside-R.org

Notice that the output from the function rxStepControl() is used to set the variableSelection parameter of rxLinMod().

The RevoScaleR code is very similar to code one would write using lm() and step():

# Code to turn file into a data frame and run with lm and step
# Read the data from a .xdf file into a data frame
MSdf <- rxXdfToDataFrame(songs, maxRowsByCols=4000000)
dim(MSdf)
# 
system.time(rlm.mod <- lm(form, data = MSdf))
summary(rlm.mod)
 
system.time(rlm.step <- step(rlm.mod, direction = "both", scope = scope, trace = 1))
#user  system elapsed 
  38.56    4.12   17.89 

Created by Pretty R at inside-R.org

The output from the RevoScaleR stepwise regression is included in the file Output (download Output) and is also similar to what is produced by lm() and step(). Notice, however, that it took step() nearly 18 seconds to run while the entire stepwise regression only took 0.16 seconds to run with rxLinMod(). We expect that, in general, computation time for rxLinMod()with rxStepControl() increase linearly with the number of observations.

A few lists for data scientists and statisticians

Looking for more resources on the web or people to follow on Twitter? Here are some lists you may find useful:

• 100 Savvy Sites on Statistics, which includes 17 sites that focus on R Programming
• Kalido offers this list of 30 influential data scientists on Twitter
• Big Data Republic is taking votes for this list of 100 influential tweeters on Big Data (voting closes on April 15).

 

Lots of data != "Big Data"

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:

1. Lots of data does not necessarily equate to “Big Data”
2. For exploratory modeling you want to work in an environment that allows for the rapid prototyping and provides the statistical tools for model evaluation and visualizations. There is no better environment that R for this kind of work, and the Revolution’s distribution of R offers the ability to work with very large samples.
3. The ability draw random samples from large data sets is the way to balance accuracy against computational constraints.

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.

R's Garden of Probability Distributions

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")       
        }

 

 

Flowchart: How to learn survey analysis with R

In a recent talk to the DC R User Group, Anthony Damico presented the following handy flowchart for learning to do survey analysis with R (actually, it's a pretty good flowchart for learning R for any application):

Since they're not clickable above, here are the resource links:

• Learn R by watching two‐minute videos on http://twotorials.com
• Read the “Getting Started with R” Guide on http://flowingdata.com
• Enrol in the free "Computing for Data Analysis" course at http://coursera.com
• Analyze Survey Data with the scripts on http://asdfree.com

Chart creator Anthony Damico is a Statistical Analyst at the Henry J. Kaiser Family Foundation, and also the author of the R SAScii package for importing public data. The first 14 minutes of his presentation, which is a great introduction to survey analysis with R, can be viewed at the link below.

Civil Statistician: DC R Meetup: “Analyze US Government Survey Data with R”

Votamatic predicted the election with R

While Nate Silver got a lot of the attention for correctly forecasting the US presidential election, other forecasters were just as succesful. Drew Linzer used the R language to build the statistical model behind votamatic.org, and was able to predict the outcome of the election months before most pundits.

Drew's model initially relied mostly on fundamental quantities: the president’s net approval-disapproval rating in June, the percent change in GDP from Q1 to Q2 of 2012, and whether the incumbent party has held the presidency for two or more terms. On that basis, Drew forecast on June 23 that the outcome (in electoral college votes) would be Obama 332 votes, Romney 206. Over time, the model used Bayesian statistics to gradually incorporate real-time polling data, and used smoothing methods to account for the fact that many state polls were sporadic. Nonetheless, the forecast never changed much, and remained around 332:206 right up to election day:

The final election result? Obama 332, Romney 206. 

Drew described his methodology in an interview with the LA Times:

On Nov. 6, I predicted that Obama would win 332 electoral votes, with 206 for Romney. But I also predicted the exact same outcome on June 23, and the prediction barely budged through election day.

How is this possible? Statistics. I did it by systematically combining information from long-term historical factors — economic growth, presidential popularity and incumbency status — with the results of state-level public opinion polls. The political and economic "fundamentals" of the race indicated at the outset that Obama was on track to win reelection. The polls never contradicted this, even after the drop in support for Obama following the first presidential debate. In fact, state-level voter preferences were remarkably stable this year; varying by no more than 2 or 3 percentage points over the entire campaign (as compared to the 5% to 10% swings in 2008).

The actual mechanics of my forecasts were performed using a statistical model that I developed and posted on my website, votamatic.org. While quantitative election forecasting is still an emerging area, many analysts were able to predict the result on the day of the election by aggregating the polls. The challenge remains to improve estimates of the outcome early in the race, and use this information to better understand what campaigns can accomplish and how voters make up their minds.

He also shared with me that he used R (and his polCA package for latent class analysis) to create the entire forecast:

Everything's done in R — data processing and graphics — and the model is fitted using WinBUGS. The website is just a WordPress blog that automatically pulls the R image files from a public dropbox.

(That's a neat trick for sharing graphics from R, by the way: write a script that writes R images to a local DropBox folder, and let DropBox take care of the web publishing by simply linking to the online DropBox file.)

You can find more details on the votamatic.org methodology at the link below. Also, if you're near San Francisco, Drew will be giving a talk to the Bay Area R User Group on February 12.

votamatic.org: How it Works