Using jjplot to explore tipping behavior

In this post, I’ll show off some recent changes to jjplot that we think are really cool. To help motivate these changes, I’ll walk through them using the tips dataset included with the reshape package.

  • Improved faceting along multiple dimensions. This shows a scatter plot of how much males and females tip on each day of the week, along with a best-fit lines. The black, dashed line shows the best-fit across all data points. Points/lines are otherwise colored by day. I’ll leave it to you to guess why the slope is higher for men on Saturday, but lower on Sunday.

    jjplot(tip ~ (abline() : group(fit(), by = day: sex) +
    point(alpha = 0.5)) : color(day) +
    abline(lty = "dashed") : fit() + total_bill,
    data = tips,
    facet.y = day, facet.x = sex)

  • New stats/geoms such as area/density. Here we’ll make a density plot of the tip fraction, that is, the tip amount over the total bill. The black density shows the overall density, while each each overlaid density shows the density just for points in that panel.

    jjplot(~ area() : group(density(), by = day:sex) : color(day, alpha = 0.5) +
    area() : group(density(), by = day) +
    I(tip / total_bill),
    data = tips,
    facet.y = day, facet.x = sex,
    xlab = "tip fraction",
    ylab = "")

  • Custom geoms/stats. We want to make it easier for the community to augment the system. Right now, the syntax is still sort of opaque and we’re working on it, but you can already get a custom stat just be naming your function jjplot.stat.*. For example, below we define a new kmeans stat. We then cluster the points and draw a best-fit line for each cluster.

    jjplot.stat.kmeans <- function(state, K, use.y = FALSE) {
    if (use.y) {
    km <- kmeans(cbind(state$data$x, state$data$y), K)
    } else {
    km <- kmeans(state$data$x, K)
    state$data$cluster <- factor(km$cluster)
    jjplot(tip ~ point() +
    abline() : group(fit(), cluster) : kmeans(3) +
    data = tips)

  • Coloring on derived statistics. You may have noticed in the earlier examples that the color syntax has changed. We figured color should be kind of like sort — it’s a pseudo-statistic which can be inserted anywhere in a statistics stack. This means that it becomes easy to color based off of derived statistics. In this example, we make the previous plot much more useful by coloring the fits and points according to the assigned cluster.

    jjplot(tip ~ (point() +
    abline() : group(fit(), cluster)) : color(cluster) : kmeans(3) +
    data = tips)

Let us know what you think! P.S. A release on CRAN is coming very soon…


Filed under Uncategorized

ePluribus: Ethnicity on Social Networks

is the name of the paper I wrote with Lars, Itamar, and Cameron. It will appear at this year’s ICWSM. You may commence bating those breaths.


Filed under Uncategorized

A few jjplot updates

Eytan and I have been actively exploring lots of crazy new ideas in jjplot, a new plotting library for R. Here’s a quick rundown of recent changes. We’d love to hear what you guys think

  1. Formulae. The old way of expressing the series of geoms and stats that form the plot was cumbersome. Putting a series of commands in the … leads to annoying issues such as poorer error handling. More importantly, because it can only express a series of statements, it becomes unclear which stats affect which geoms, making it impossible to express more complicated combinations.

    We believe that formulae are a good solution to this. Layers are separated by ‘+’ operations. Interactions between stats and geoms are expressed via the interaction operator ‘:’. This allows us to gracefully express arbitrary trees of stats and geoms. An example of a jittered scatter plot:


    jjplot(x = x, y = y, data = data,


    jjplot(y ~ point() : jitter(xfactor = 1) + x, data = data)

    The leftmost and rightmost terms correspond to the y and x aesthetics. For a simple case such as this, formulae might not seem like much of an improvement. But consider a more complex example:

    jjplot( ~ line(lty="dashed", col = "red") : hist() +
    bar(width = 0.1) : hist() : jitter(xfactor = 1) +
    Sepal.Length, data = iris)

    Reading from the right, this says to take iris$Sepal.Length, jitter it, bin the data, and bar plot the result. This is cool because it’s immediately clear that you’re stacking stats, plotting a histogram of the jittered data. The first term does the same thing, except that it does a hist() statistic without the jitter, and draws this as a red line.

    By using parentheses, you can also apply a stat to multiple stats/geoms.

    jjplot( ~ (point(col = "blue", size=3) +
    line(col = "red", lty="dashed") +
    bar(width=0.25)) : hist() +
    Petal.Length, data = iris)

    Here we’re just plotting a histogram but with some extra geoms on top for some extra flair.

    We think this notation is a simple and elegant way of expressing what interacts with what.

  2. Facets This way of thinking about facets is somewhat controversial among us. Normally, facets conflate two concepts: how you compute statistics and how you plot them. This means that you compute statistics on facet subsets, then you plot each subset in a separate panel. Well, currently jjplot takes a different tack, treating facets as merely a command to plot different subsets of the data in different panels. To see what this implies, consider

    df <- data.frame(state = rownames(state.x77),
    region = state.region,
    jjplot(Murder ~ abline(lty = "dashed") : fit() +
    abline() : group(fit(), by = region) +
    point() + Income,
    data = df, color = region, facet = region)

    The first two terms simply do a scatter plot. The next line does lm fits on each subset. Note that you have to be explicit with the grouping. With old semantics, you’d have an implicit group by on the facet variable, but because we aren’t combining the grouping and the faceting anymore, you have to spell it out. The first line shows you the effect of leaving out the grouping operator: you get a fit over all the data that appears on all panels. This is something I’ve always wanted to do and it seems to also be persistent question on stack overflow (e.g., “how do I draw a line at the facet/global mean on each facet panel?”). Hopefully this formulation makes it obvious.

  3. Sorting Another persistent question is how to perform sorting on factor scales. Because of the ease of stacking stats in the formula formulation, we think it makes sense to add a few special stats/geoms. One of them is the sort stat. This performs an identity operation on the data frame but also appends some metadata about how to order things which is then intercepted when the scales are created. Here are some usage examples:

    df <- data.frame(name = factor(letters),
    value = rnorm(26 * 6),
    type = rep(factor([1:6]), each = 26))
    jjplot(name ~ point() + value,
    data = df, color = type, facet = type)

    The first plot is the data unsorted.

    jjplot(name ~ point() : sort(y = value) + value,
    data = df, color = type, facet = type)

    The second plot sorts according to the mean value associated with each factor across all facets (remember no grouping!). Like relevel, the sort statistic can take a function argument to specify how multiple points should be sorted.

    jjplot(name ~ point() : group(sort(y = value), by=type) + value,
    data = df, color = type, facet = type)

    The last plot wraps the sort in a group by, meaning that each facet panel has its own sorting order.

All of this awesomeness is available in the current svn repo. Check it out!


Filed under Uncategorized

R LDA package updated to version 1.2 and an ideal-point model for political blogs

I’ve been on a bit of a R tear lately. Today you should see a new version of the R lda package. This version has lots of fixes including a working mmsb demo with the latest version of ggplot2, corrected RTM code, improved likelihood reporting, better documentation, and much more. Grab it from CRAN today! Special thanks to the following people for bug reports/feature requests (sorry if I forgot anyone):

  • Edo Airoldi
  • Jordan Boyd-Graber
  • Khalid El-Arini
  • Roger Levy
  • Solomon Messing
  • Joerg Reichardt

One of the new features is a method to make sLDA predictions on response variables conditioned on documents. In the demo accompanying the package, I fit an sLDA model to a corpus of political blogs tagged as being either liberal or conservative. With this fitted model, I can now use the new predict method to predict the political bent of each of the blogs within a continuous space. The density plot of these predictions is given below, broken down by the the original conservative/liberal label (color of shading).

I like how there’s some bimodality for each contingency — a moderate group and a more extreme group. The model also predicts a heavy tail of super-conservative blogs. There is a real notable bump down by -3. I dunno if this represents reality; it’s probably worthwhile to do more extensive model checking.


Filed under Uncategorized

jjplot: Yet another plotting library for R

Those of you who follow this blog know that making (somewhat) pretty plots is an abiding interest of mine. Many of the plots I’ve made in the past were done using the great ggplot2 package. But recently Eytan Bakshy and I have been tinkering with our own plotting library, jjplot, as a playground for various ideas we’ve had. As the name indicates, it is heavily inspired by hadley’s library. Our library doesn’t do quite as much as ggplot2, and ours is liable to be much buggier. But it’s still fun to play with. Here are some examples of what jjplot can do:

  • Bar plots with fills controlled by the values.

    df <- data.frame(x = 1:50, y = rnorm(50))
    jjplot(x, y, data = df, fill = y, = "black"))

  • Boxplots.

    df <- data.frame(state = rownames(state.x77), region = state.region, state.x77)
    jjplot(region, Income, data = df, fill = region,, by = region),

  • Scatter plot, colored by factor, with alpha blending. This also demonstrates how statistics can be used to visualize different aspects of the data simultaneously.

    df <- data.frame(x = rnorm(10000) + (1:4) * 1, f = factor(c('A', 'B', 'C', 'D')))
    df$y <- c(-6, -2, 2, 4) * df$x + rnorm(10000)
    jjplot(x + 2, y, data = df, alpha = 0.10, color = f, jjplot.point(),, by = f), jjplot.abline(),, jjplot.hline(lty = "dashed"))

  • An example of log scales and the CCDF statistic.

    df <- data.frame(x=rlnorm(1000,2,2.5))
    jjplot(x, data = df, jjplot.ccdf(density=TRUE), jjplot.point(), log='xy')

Lots more demos and documentation are here. To install visit and install the downloaded package using

R CMD INSTALL jjplot_1.0.tar.gz

We’re eager to hear your feedback!


Filed under Uncategorized

Axl Rose by any other name…

In a post a while ago, I wondered how much info about the band one could glean just by looking at the name. I mean, shouldn’t it be obvious that a band named “Trauma” should be heavy metal?

This was the genesis of a collaboration between me and Matt Hoffman. We wanted to see if you could improve genre prediction using the names of the bands. Unfortunately, neither of us had enough time to really get this project going, but I thought I’d share what results we did get in hopes that someone else will pick up the torch.

To start off, we need a large training set of band/genre mappings. We opted for the DBPedia Infobox mine that you can find at infochimps. (For those who don’t know, they’ve done some awesome data mining to grab all the structured info from Wikipedia infoboxes). I did some cleaning up and have put up the list of artists and genres (the artist in each line of the first file is associated with the genres on the corresponding line of the second file).

You might have noticed that Wikipedia is pretty crazy when it comes to genre definitions (because god forbid we confuse Melodic Death Metal and Power Metal). This craziness makes it hard to map the artists to any canonicalized genre set (such as CAL-500). I tried a bunch of techniques to do this canonicalization (including doing my own crawl of Wikipedia with all sorts of heuristics). None of it worked very well for mapping genres to a canonicalized set, but it did let me make a really cool graph of connections between genres. Eventually, we came to the conclusion that we needed human judgments. We got mechanical turkers to label Wikipedia genres with CAL-500 genres. Those results are here.

With that training set in place, I decided to explore the data to see if there truly were correlations between substrings of artist names and genres. The plot below shows the prevalence in each genre of artists containing “death” (red) or “boyz” (blue) in their name. The green dots show the overall distribution of genres among artists in Wikipedia.

The graph shows that bands containing “death” in their name are much more likely to be Rock, Alternative, Metal/Hard Rock, or Alternative. Conversely, they are less likely to be Jazz, Hip-Hop, or Soul. In contrast, bands containing “boyz” in their name are overwhelmingly Hip-Hop. This confirmed my intuition and seemed promising to me, so we went ahead and developed a classifier for the CAL-500 data set. The techniques we tried were:

  • names (corrLDA)– the correlated topic model fit to the Wikipedia data. Predictions use only names.
  • names (NB) – naive Bayes fit to the Wikipedia data. Predictions use only names.
  • names (LR) – logistic regression fit to the Wikipedia data. Predictions use only names.
  • baseline – Predictions use the baseline frequency of genres on Wikipedia. Predictions do not use any information about the instances.
  • svm – SVM fit using MFCC features. Predictions use both names and audio.
  • svm + names (corrLDA) – SVM fit using MFCC features plus the results of names (corrLDA). Predictions use both names and audio.
  • svm + names (NB)– SVM fit using MFCC features plus the results of names (NB). Predictions use both names and audio.
  • svm + names (LR)– SVM fit using MFCC features plus the results of names (LR). Predictions use both names and audio.

The plot below shows the precision-recall for each of these techniques. As you can see, it’s not very promising. The SVM will outclass any technique which uses the name along; otherwise all of the name techniques look about the same. It looks like we might get a small bump by combining SVM with names (LR) but it’s hard to tell.

But precision-recall may not be the right metric. After all, pop and rock are so frequent that you will probably predict pop for every single item in the test set before you even make any other prediction. Something which is perhaps more meaningful is to look at the rank of the correct labels on a per-test-instance level; the lower the rank, the better the model is at making predictions. Boxplots of the ranks are given below.

We see slightly different patterns when we look at the ranks. Without using any audio data, the naive Bayes technique performs best and manages to get a non-trivial bump beyond the baseline. When audio is included, the names add something, but not much. Interestingly, the names (LR) technique which looked like it might help us at precision-recall actually does a bit worse when you look at the rank. On the other hand, SVM + names (corrLDA) has the same median as SVM, but manages to do a better job at some of the difficult-to-predict cases, leading to a smaller interquartile range.

In sum, names give us something — unfortunately, it’s not a whole lot.


Filed under Uncategorized

The cost of a sample

I once heard it on good authority that Gelman says you usually don’t need more than 12 samples. Well, as a result of a discussion with Sam Gershman (sorry Sam for not answering the actual question you asked!), I wondered if that was true; that is, if under reasonable assumptions it might be better to take a small number of samples. Caveat: there’s probably lots of work on this already, but where would the fun be in that?

Ok, let’s assume that your goal is to estimate \mathbb{E}_{z \sim p(z | x)}[f(z)], where p(z | x) represents some distribution on hidden variables over which you are trying to compute a function, f. For the usual reasons, it’s intractable to compute this exactly, so you’re going to use a sampler. Let’s assume

  • that your sampler has mixed and that you’re getting independent samples (that condition alone should give you fair warning that what I’m about to say is of little practical value);
  • f is bounded (say between 0 and 1);
  • to obtain n samples from the sampler costs some amount, say R(n).

More samples are usually better, because they’ll give you a better representation of the true distribution, i.e. \mathbb{E}_{z \sim \hat{p_n}(z | x)}[f(z)] \rightarrow \mathbb{E}_{z \sim p(z | x)}[f(z)], where \hat{p_n} is the distribution obtained by using n samples. Unfortunately, more samples come at a cost here, so you don’t want too many. How should you tradeoff then?

We can define a loss by \ell = R(n) + |\mathbb{E}_{z \sim \hat{p_n}(z | x)}[f(z)] - \mathbb{E}_{z \sim p(z | x)}[f(z)]|, that is, how far off our sampled estimate is from the truth, plus the cost of obtaining those samples. Using Hoeffding, we can bound the loss \ell < R(n) + \epsilon with probability 1 - 2 \exp( -2 n \epsilon^2). This expression gives you something to think about when you're trying to decide how many samples to take — more samples loosen the bound but increase its probability.

If your cost is linear, R(n) = a n, you might want to choose
something like n = \frac{\epsilon}{a}, which gives you a loss of \ell < 2 \epsilon with probability 1 - 2 \exp(-2 \epsilon^3 / a).

The plot below shows what might happen if you make such a choice. Here, I've let the posterior be an equiprobable binomial distribution. The function I'm computing is the identity f(z) = z. The curves show the loss, \ell for various choices of the cost parameter a as a function of the number of samples. The dots show the chosen values of n for each value of a; the horizontal lines show the 80% loss bound for these choices.

Turns out for some reasonable values, you really should stick to about 12 samples.


Filed under Uncategorized