In [9]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import numpy.random as random
from IPython.display import Image
# seaborn is under active development and throws some scary looking warnings
import warnings 
# this will allow us to use the code in peace :) 
warnings.filterwarnings("ignore")

Lecture 17:

  • Learn how to use the seaborn package to produce beautiful plots

  • Learn about kernel density estimates

  • Learn appropriate ways of representing different types of data

seaborn

In this lecture, we will learn about seaborn. seaborn is a package with many tools for data visualization. It allows you to make pretty plots. Almost anything in seaborn can be done using matplotlib, but with seaborn's built-in functions you can reduce a lot of matplotlib code down to a single line. seaborn isn't just a pretty face. Its real power is in statistical data analysis. It has a lot of functions built in for visualizing the distribution of your data, for example.

Let's take a look at some of the plots we can make with this package. We can import it using:

In [2]:
import seaborn as sns

Unusual distributions, kernel density estimates and jointplots

In some cases, we have distributions of data that don't look like a simple (e.g., normal) distribution, for example, the data could be bimodal or have skewed shaped distributions.

Let's create some synthetic bimodal data by drawing from two separate normal/lognormal distributions and combining the two into two bimodal data sets. We do this by drawing from random.normal( ) twice for each of the two: $x$ and $y$.

In [3]:
xdata1=random.normal(20,25,5000) # first x draw
xdata2=random.normal(100,25,5000) # second x draw
ydata1=random.lognormal(2,0.1,8000) # first y draw
ydata2=random.lognormal(3,0.1,2000) # second y draw
xdata=np.append(xdata1,xdata2) # combine the two x data sets
ydata=np.append(ydata1,ydata2) # combine the two y data sets

When we plot our xdata as a histogram, we can see that we have a broadly bimodal distribution. For fun, let's also plot the mean of the distribution as a red line.

In [4]:
plt.hist(xdata,bins=50)
plt.axvline(np.mean(xdata),color='red',linewidth=3);

We can see that our mean lies right between the twin peaks. Describing this distribution with statistics meant for normal distributions (mean or standard deviation) is just plain wrong.

One of the ways of representing the distribution of a set of datapoints is known as a 'kernel density estimate' (kde). This is a useful way of showing the distribution of data. It places a 'kernel' (an assumed distribion at the data point level - usually a normal distribution) at each data point and then sums up the contributions from all the data points. This avoids the awkwardness of choice of bin size associated with histograms, for example. (We just picked 50 in the plot above - why 50?).

Here are some data represented on a bar plot in the lefthand plot. And on the right, we illustrate the idea behind kernal denisty estimates. The black lines are the locations of individual datapoints and the red dashed lines are the kernels at each point. The heavy blue line is the kernel density estimate (the sum of all the red dashed lines).

In [5]:
Image('Figures/KDE.png')
Out[5]:

Happily, we can plot kernel density estimates using the sns.kdeplot( ) function. The shade argument allows us to shade the area underneath the curve. In matplotlib, the same thing can be achieved using the function plt.fill_between

In [8]:
sns.kdeplot(xdata,shade=True);

We can also plot the density estimate and the histogram on top of one another using the sns.distplot( ) function.

In [11]:
sns.distplot(xdata);

As you can see, this is a lot quicker than how we were plotting our distribution in the previous lecture!

With our $ydata$ we can see that we also have a bimodal distribution, but there are far fewer data points in the wider mode (we only used 2000 of our 10000 points for this mode).

In [12]:
sns.distplot(ydata);

What if our data had both $x$ and $y$ components. For example, measurements of length and width from a set of shark teeth with two species in it. How would we visualize it? Let's just try plotting the $x$ data on the x axis and $y$ data on the y axis as dots on a regular matplotlib plot.

In [14]:
plt.plot(xdata,ydata,'.');

It's much harder to see the bimodal distribution in the $x$ data in this case, and we can't really see that there are so much fewer data points in the $y$ data mode than the more tightly clustered one. Also, this is not a plot you would expect from length and width dimensions of teeth from a single species - wouldn't that be a linear plot?

When we have a lot of datapoints, this type of plot gives us no easy way to estimate data density. Fortunately, seaborn includes a cool plot called 'sbn.jointplot( ). This combines the histograms and the scatter plot into a single graph.

In [15]:
sns.jointplot(xdata,ydata);

Jointplot can also do kernel density estimates! The kind argument on this plot gives us many other options plotting data as well:

In [11]:
help(sns.jointplot)
Help on function jointplot in module seaborn.axisgrid:

jointplot(x, y, data=None, kind='scatter', stat_func=None, color=None, height=6, ratio=5, space=0.2, dropna=True, xlim=None, ylim=None, joint_kws=None, marginal_kws=None, annot_kws=None, **kwargs)
    Draw a plot of two variables with bivariate and univariate graphs.
    
    This function provides a convenient interface to the :class:`JointGrid`
    class, with several canned plot kinds. This is intended to be a fairly
    lightweight wrapper; if you need more flexibility, you should use
    :class:`JointGrid` directly.
    
    Parameters
    ----------
    x, y : strings or vectors
        Data or names of variables in ``data``.
    data : DataFrame, optional
        DataFrame when ``x`` and ``y`` are variable names.
    kind : { "scatter" | "reg" | "resid" | "kde" | "hex" }, optional
        Kind of plot to draw.
    stat_func : callable or None, optional
        *Deprecated*
    color : matplotlib color, optional
        Color used for the plot elements.
    height : numeric, optional
        Size of the figure (it will be square).
    ratio : numeric, optional
        Ratio of joint axes height to marginal axes height.
    space : numeric, optional
        Space between the joint and marginal axes
    dropna : bool, optional
        If True, remove observations that are missing from ``x`` and ``y``.
    {x, y}lim : two-tuples, optional
        Axis limits to set before plotting.
    {joint, marginal, annot}_kws : dicts, optional
        Additional keyword arguments for the plot components.
    kwargs : key, value pairings
        Additional keyword arguments are passed to the function used to
        draw the plot on the joint Axes, superseding items in the
        ``joint_kws`` dictionary.
    
    Returns
    -------
    grid : :class:`JointGrid`
        :class:`JointGrid` object with the plot on it.
    
    See Also
    --------
    JointGrid : The Grid class used for drawing this plot. Use it directly if
                you need more flexibility.
    
    Examples
    --------
    
    Draw a scatterplot with marginal histograms:
    
    .. plot::
        :context: close-figs
    
        >>> import numpy as np, pandas as pd; np.random.seed(0)
        >>> import seaborn as sns; sns.set(style="white", color_codes=True)
        >>> tips = sns.load_dataset("tips")
        >>> g = sns.jointplot(x="total_bill", y="tip", data=tips)
    
    Add regression and kernel density fits:
    
    .. plot::
        :context: close-figs
    
        >>> g = sns.jointplot("total_bill", "tip", data=tips, kind="reg")
    
    Replace the scatterplot with a joint histogram using hexagonal bins:
    
    .. plot::
        :context: close-figs
    
        >>> g = sns.jointplot("total_bill", "tip", data=tips, kind="hex")
    
    Replace the scatterplots and histograms with density estimates and align
    the marginal Axes tightly with the joint Axes:
    
    .. plot::
        :context: close-figs
    
        >>> iris = sns.load_dataset("iris")
        >>> g = sns.jointplot("sepal_width", "petal_length", data=iris,
        ...                   kind="kde", space=0, color="g")
    
    Draw a scatterplot, then add a joint density estimate:
    
    .. plot::
        :context: close-figs
    
        >>> g = (sns.jointplot("sepal_length", "sepal_width",
        ...                    data=iris, color="k")
        ...         .plot_joint(sns.kdeplot, zorder=0, n_levels=6))
    
    Pass vectors in directly without using Pandas, then name the axes:
    
    .. plot::
        :context: close-figs
    
        >>> x, y = np.random.randn(2, 300)
        >>> g = (sns.jointplot(x, y, kind="hex")
        ...         .set_axis_labels("x", "y"))
    
    Draw a smaller figure with more space devoted to the marginal plots:
    
    .. plot::
        :context: close-figs
    
        >>> g = sns.jointplot("total_bill", "tip", data=tips,
        ...                   height=5, ratio=3, color="g")
    
    Pass keyword arguments down to the underlying plots:
    
    .. plot::
        :context: close-figs
    
        >>> g = sns.jointplot("petal_length", "sepal_length", data=iris,
        ...                   marginal_kws=dict(bins=15, rug=True),
        ...                   annot_kws=dict(stat="r"),
        ...                   s=40, edgecolor="w", linewidth=1)

Let's try the kind=kde option. (Warning, doing a 2d kde plot like this is a lot of work so this cell might take a while to run).

In [12]:
sns.jointplot(xdata,ydata,kind='kde');

Another type of plot that looks nice is the hexbin plot (kind of like a hexagonal 2d histogram)

In [13]:
sns.jointplot(xdata,ydata,kind='hexbin');

Multi dimensional data and pairplots

What if you have data with more than two dimensions? A good example of this is with isotopic data from Ocean Island Basalts. Isotopic systems are used to "finger-print" different sources of melt and in the mantle. It is used to characterize what is deep in the Earth using what gets brought up to form the ocean islands. By now there are data available for many different isotopic systems. Here we take a look at a small sample of what is in the GeoRoc database for ocean island basalts.

In [17]:
MantleArray=pd.read_csv('DataSets/GeoRoc/MantleArray_OIB.csv')
MantleArray.head()
Out[17]:
CITATION EPSILON_HF EPSILON_ND HF176_HF177 LAND/SEA (SAMPLING) LATITUDE (MAX.) LATITUDE (MIN.) LOCATION LOCATION COMMENT LONGITUDE (MAX.) ... PB208_PB204 ROCK NAME ROCK TYPE SAMPLE NAME SAMPLING TECHNIQUE SR87_SR86 TECTONIC SETTING TYPE OF MATERIAL UNIQUE_ID Year
0 [60] STILLE P. (1986) 14.534010 7.607708 0.283196 SUBAERIAL 19.83 19.83 HAWAIIAN ISLANDS NaN -155.42 ... 38.017 ANKARAMITE VOLCANIC ROCK samp. 79MK1 OUTCROP 0.70347 OCEAN ISLAND WHOLE ROCK 107-79MK1 1986
1 [60] STILLE P. (1986) 10.538041 6.378770 0.283083 SUBAERIAL 22.00 22.00 HAWAIIAN ISLANDS NAPALI MEMBER -159.50 ... 37.803 THOLEIITE VOLCANIC ROCK samp. KAU-1 OUTCROP 0.70384 OCEAN ISLAND WHOLE ROCK 5-KAU-1 1986
2 [60] STILLE P. (1986) 11.033117 6.281235 0.283097 SUBAERIAL 22.00 22.00 HAWAIIAN ISLANDS NAPALI MEMBER -159.50 ... 37.962 THOLEIITE VOLCANIC ROCK samp. 1D872-2 OUTCROP 0.70364 OCEAN ISLAND WHOLE ROCK NaN 1986
3 [60] STILLE P. (1986) 12.164719 5.696027 0.283129 SUBAERIAL 21.15 21.15 HAWAIIAN ISLANDS NaN -156.97 ... 37.751 THOLEIITE VOLCANIC ROCK samp. 71WMOL-1 OUTCROP 0.70378 OCEAN ISLAND WHOLE ROCK 782-71WMOL-1 1986
4 [60] STILLE P. (1986) 15.135173 5.891097 0.283213 SUBAERIAL 21.16 21.16 HAWAIIAN ISLANDS NaN -157.23 ... 37.754 THOLEIITE VOLCANIC ROCK samp. 71WMOL-3 OUTCROP 0.70376 OCEAN ISLAND WHOLE ROCK 783-71WMOL-3 1986

5 rows × 24 columns

There are lots of different isotope ratios available. In this lecture, we focus on four different isotopic ratios: $^{87}$Sr/$^{86}$Sr , $^{206}$Pb/$^{204}$Pb , $ \varepsilon$ Nd and $ \varepsilon$ Hf. (The ratios with an $ \varepsilon$ are the isotope ratios relative to a standard value in parts per 10000, this allows you to see the variation better as it is normally very small).

So, how do you plot multi-dimensional data? Fortunately, seaborn has a function for multi-dimensional data known as the sns.pairplot( ). This makes a plot that takes data with many dimensions and plots each possible two dimensional combination against one another in a grid where each of the rows and columns represents a dimension. When a particular dimension is plotted against itself, however, it gives the kernal density estimates for the different data types (set by the keyword argument, hue). In this case we will be using $LOCATION$ to group the data).

In our isotopic example, the top row plots $ \varepsilon$ Nd against all of the other isotope ratios in turn with the first plot on the left being the kernal density plot for each location with each in a different color.

Let's see it in action (warning, cell may take a while to run due to large number of locations):

In [26]:
sns.pairplot(MantleArray,\
             vars=['EPSILON_ND','SR87_SR86','PB206_PB204','EPSILON_HF'],\
             hue='LOCATION');

sns.pairplot( ) gives us an interesting way of seeing trends in multidimensional space. Some ocean islands have quite different isotopic compositions from others! It's important to note that some projections of these data (e.g. $^{206}$Pb/$^{204}$Pb vs $ \varepsilon$ Nd) show this effect more prominently whereas others, like $ \varepsilon$ Nd vs $ \varepsilon$ Hf, are more subtle. In this plot (second from the bottom on the left), the Austral-Cook Islands are all the way on the top to the right, and Pitcairn is way at the bottom to the left.

It might be surprising to some of you who know a little geophysics that there is this heterogeneity between different ocean islands. You might have thought that the mantle is very well mixed by convection. Many geochemists believe that the source regions for the plumes that form these basalts have different compositions because subduction brings down crustal material with different compositions to the base of the mantle, and this material then upwells in plumes that form the ocean islands. According to that hypothesis, the variation in isotopic ratios we see here comes from mixing between these recycled crustal end members and a mantle composition.

We will look further at this dataset when we learn how to plot things in three dimensions.

Box and whisker, violin plots and swarm plots

The kernal density estimate plots in the element versus itself spots in the sns.pairplot( ) plot example above plot all of the isotopic data on one plot. This doesn't really allow us to appreciate fully the variability at a given location very well. There are two additional plots that might give us a better idea of that: the 'box and whisker' plot and the 'violin' plot.

seaborn has functions (actually classes) for both of these, sns.boxplot( ) and sns.violinplot( ). In these plots, we can also place the data (and not just the distributions) and can use the handy sns.swarmplot( ) to spread out the data along the x-axis so that we can see them better. We'll go through each of these tricks in the following.

First, the boxplot:

In [18]:
help(sns.boxplot)
Help on function boxplot in module seaborn.categorical:

boxplot(x=None, y=None, hue=None, data=None, order=None, hue_order=None, orient=None, color=None, palette=None, saturation=0.75, width=0.8, dodge=True, fliersize=5, linewidth=None, whis=1.5, notch=False, ax=None, **kwargs)
    Draw a box plot to show distributions with respect to categories.
    
    A box plot (or box-and-whisker plot) shows the distribution of quantitative
    data in a way that facilitates comparisons between variables or across
    levels of a categorical variable. The box shows the quartiles of the
    dataset while the whiskers extend to show the rest of the distribution,
    except for points that are determined to be "outliers" using a method
    that is a function of the inter-quartile range.
    
    
    Input data can be passed in a variety of formats, including:
    
    - Vectors of data represented as lists, numpy arrays, or pandas Series
      objects passed directly to the ``x``, ``y``, and/or ``hue`` parameters.
    - A "long-form" DataFrame, in which case the ``x``, ``y``, and ``hue``
      variables will determine how the data are plotted.
    - A "wide-form" DataFrame, such that each numeric column will be plotted.
    - An array or list of vectors.
    
    In most cases, it is possible to use numpy or Python objects, but pandas
    objects are preferable because the associated names will be used to
    annotate the axes. Additionally, you can use Categorical types for the
    grouping variables to control the order of plot elements.    
    
    This function always treats one of the variables as categorical and
    draws data at ordinal positions (0, 1, ... n) on the relevant axis, even
    when the data has a numeric or date type.
    
    See the :ref:`tutorial <categorical_tutorial>` for more information.    
    
    Parameters
    ----------
    x, y, hue : names of variables in ``data`` or vector data, optional
        Inputs for plotting long-form data. See examples for interpretation.        
    data : DataFrame, array, or list of arrays, optional
        Dataset for plotting. If ``x`` and ``y`` are absent, this is
        interpreted as wide-form. Otherwise it is expected to be long-form.    
    order, hue_order : lists of strings, optional
        Order to plot the categorical levels in, otherwise the levels are
        inferred from the data objects.        
    orient : "v" | "h", optional
        Orientation of the plot (vertical or horizontal). This is usually
        inferred from the dtype of the input variables, but can be used to
        specify when the "categorical" variable is a numeric or when plotting
        wide-form data.    
    color : matplotlib color, optional
        Color for all of the elements, or seed for a gradient palette.    
    palette : palette name, list, or dict, optional
        Colors to use for the different levels of the ``hue`` variable. Should
        be something that can be interpreted by :func:`color_palette`, or a
        dictionary mapping hue levels to matplotlib colors.    
    saturation : float, optional
        Proportion of the original saturation to draw colors at. Large patches
        often look better with slightly desaturated colors, but set this to
        ``1`` if you want the plot colors to perfectly match the input color
        spec.    
    width : float, optional
        Width of a full element when not using hue nesting, or width of all the
        elements for one level of the major grouping variable.    
    dodge : bool, optional
        When hue nesting is used, whether elements should be shifted along the
        categorical axis.    
    fliersize : float, optional
        Size of the markers used to indicate outlier observations.
    linewidth : float, optional
        Width of the gray lines that frame the plot elements.    
    whis : float, optional
        Proportion of the IQR past the low and high quartiles to extend the
        plot whiskers. Points outside this range will be identified as
        outliers.
    notch : boolean, optional
        Whether to "notch" the box to indicate a confidence interval for the
        median. There are several other parameters that can control how the
        notches are drawn; see the ``plt.boxplot`` help for more information
        on them.
    ax : matplotlib Axes, optional
        Axes object to draw the plot onto, otherwise uses the current Axes.    
    kwargs : key, value mappings
        Other keyword arguments are passed through to ``plt.boxplot`` at draw
        time.
    
    Returns
    -------
    ax : matplotlib Axes
        Returns the Axes object with the plot drawn onto it.    
    
    See Also
    --------
    violinplot : A combination of boxplot and kernel density estimation.    
    stripplot : A scatterplot where one variable is categorical. Can be used
                in conjunction with other plots to show each observation.    
    swarmplot : A categorical scatterplot where the points do not overlap. Can
                be used with other plots to show each observation.    
    
    Examples
    --------
    
    Draw a single horizontal boxplot:
    
    .. plot::
        :context: close-figs
    
        >>> import seaborn as sns
        >>> sns.set(style="whitegrid")
        >>> tips = sns.load_dataset("tips")
        >>> ax = sns.boxplot(x=tips["total_bill"])
    
    Draw a vertical boxplot grouped by a categorical variable:
    
    .. plot::
        :context: close-figs
    
        >>> ax = sns.boxplot(x="day", y="total_bill", data=tips)
    
    Draw a boxplot with nested grouping by two categorical variables:
    
    .. plot::
        :context: close-figs
    
        >>> ax = sns.boxplot(x="day", y="total_bill", hue="smoker",
        ...                  data=tips, palette="Set3")
    
    Draw a boxplot with nested grouping when some bins are empty:
    
    .. plot::
        :context: close-figs
    
        >>> ax = sns.boxplot(x="day", y="total_bill", hue="time",
        ...                  data=tips, linewidth=2.5)
    
    Control box order by passing an explicit order:
    
    .. plot::
        :context: close-figs
    
        >>> ax = sns.boxplot(x="time", y="tip", data=tips,
        ...                  order=["Dinner", "Lunch"])
    
    Draw a boxplot for each numeric variable in a DataFrame:
    
    .. plot::
        :context: close-figs
    
        >>> iris = sns.load_dataset("iris")
        >>> ax = sns.boxplot(data=iris, orient="h", palette="Set2")
    
    Use ``hue`` without changing box position or width:
    
    .. plot::
        :context: close-figs
    
        >>> tips["weekend"] = tips["day"].isin(["Sat", "Sun"])
        >>> ax = sns.boxplot(x="day", y="total_bill", hue="weekend",
        ...                  data=tips, dodge=False)
    
    Use :func:`swarmplot` to show the datapoints on top of the boxes:
    
    .. plot::
        :context: close-figs
    
        >>> ax = sns.boxplot(x="day", y="total_bill", data=tips)
        >>> ax = sns.swarmplot(x="day", y="total_bill", data=tips, color=".25")
    
    Use :func:`catplot` to combine a :func:`pointplot` and a
    :class:`FacetGrid`. This allows grouping within additional categorical
    variables. Using :func:`catplot` is safer than using :class:`FacetGrid`
    directly, as it ensures synchronization of variable order across facets:
    
    .. plot::
        :context: close-figs
    
        >>> g = sns.catplot(x="sex", y="total_bill",
        ...                 hue="smoker", col="time",
        ...                 data=tips, kind="box",
        ...                 height=4, aspect=.7);

Let's look at the $\epsilon$ Hf data as a function of location (similar to the lower right hand corner plot in the pairplot example above:

In [19]:
box=sns.boxplot(data=MantleArray,x='LOCATION',y='EPSILON_HF')

Oh my. Pretty plot, but the x axis tick labels overlap so bad we can't even see them. To rotate them, we us the figure object as we did in Lecture 16. We will call the figure object box in the following to get and set the x tick labels.

In [20]:
box=sns.boxplot(data=MantleArray,x='LOCATION',y='EPSILON_HF')
plt.setp(box.get_xticklabels(), rotation=90);

The x axis labels are better, but the plot still doesn't look great. The black diamonds are known as 'fliers' and could represent "outliers".

As we have some pretty large "outliers" in this dataset, we need to either limit the range on the y axis or by making the plot larger. These fliers are annoying. It is possible to replace them with a small circle and exclude them from the box and whisker plot, using the argument fliersize=0.

You need to be aware that this can be a problematic practice as you are effectively hiding data from the viewer. Personally, I wouldn't do it... But lot's of folks do, so here is how:

In [23]:
plt.figure(figsize=(10,8))
box=sns.boxplot(data=MantleArray,x='LOCATION',y='EPSILON_HF',fliersize=0)
plt.setp(box.get_xticklabels(), rotation=90);

See! This looks much better. But looks can be deceiving!

Alternatively, we can put the data back on, but use a white dot with a black edge, and spread them out so we can see them.
In order to see the data points that overlap one another, seaborn has an option to plot the data points in a cloud using the function sns.swarmplot( ). Here we do both of these (and also make the points white with a black edge). The argument alpha=0.5 is used here to make the points semi transparent.

In [24]:
plt.figure(figsize=(20,10))
box=sns.boxplot(data=MantleArray,x='LOCATION',y='EPSILON_HF',fliersize=0)
plt.setp(box.get_xticklabels(), rotation=90);
sns.swarmplot(data=MantleArray,x='LOCATION',y='EPSILON_HF',color='white',edgecolor='black',linewidth=1,alpha=0.5);