title: Showing data and model predictions in one plot using Seaborn and Statsmodels
body: <p>Seaborn is a great package for doing some high-level plotting with pretty outputs. However, I'm struggling a little with using Seaborn to overlay both data and model predictions from an externally-fit model. In this example I am fitting models in Statsmodels that are too complex for Seaborn to do out-of-the-box, but I think the problem is more general (i.e. if I have model predictions and want to visualise both them and data using Seaborn).</p>
  
  <p>Let's start with imports and a dataset:</p>
  
  <pre><code>import numpy as np
  import pandas as pd
  import seaborn as sns
  import statsmodels.formula.api as smf
  import patsy
  import itertools
  import matplotlib.pyplot as plt
  
  np.random.seed(12345)
  
  # make a data frame with one continuous and two categorical variables:
  df = pd.DataFrame({'x1': np.random.normal(size=100),
                       'x2': np.tile(np.array(['a', 'b']), 50),
                       'x3': np.repeat(np.array(['c', 'd']), 50)})
  
  # create a design matrix using patsy:
  X = patsy.dmatrix('x1 * x2 * x3', df)
  
  # some random beta weights:
  betas = np.random.normal(size=X.shape[1])
  
  # create the response variable as the noisy linear combination of predictors:
  df['y'] = np.inner(X, betas) + np.random.normal(size=100)
  </code></pre>
  
  <p>We fit a model in statsmodels containing all predictor variables and their interactions:</p>
  
  <pre><code># fit a model with all interactions
  fit = smf.ols('y ~ x1 * x2 * x3', df).fit()
  print(fit.summary())
  </code></pre>
  
  <p>Since in this case we have all combinations of variables specified, and our model predictions are linear, it would suffice for plotting to add a new "predictions" column to the dataframe containing the model predictions. However, that's not very general (imagine our model is nonlinear and so we want our plots to show smooth curves), so instead I make a new dataframe with all combinations of predictors, then generate predictions:</p>
  
  <pre><code># create a new dataframe of predictions, using pandas' expand grid:
  def expand_grid(data_dict):
      """ A port of R's expand.grid function for use with Pandas dataframes.
  
      from http://pandas.pydata.org/pandas-docs/stable/cookbook.html?highlight=expand%20grid
  
      """
      rows = itertools.product(*data_dict.values())
      return pd.DataFrame.from_records(rows, columns=data_dict.keys())
  
  
  # build a new matrix with expand grid:
  
  preds = expand_grid(
                  {'x1': np.linspace(df['x1'].min(), df['x1'].max(), 2),
                   'x2': ['a', 'b'],
                   'x3': ['c', 'd']})
  preds['yhat'] = fit.predict(preds)
  </code></pre>
  
  <p>The <code>preds</code> dataframe looks like this:</p>
  
  <pre><code>  x3        x1 x2      yhat
  0  c -2.370232  a -1.555902
  1  c -2.370232  b -2.307295
  2  c  3.248944  a -1.555902
  3  c  3.248944  b -2.307295
  4  d -2.370232  a -1.609652
  5  d -2.370232  b -2.837075
  6  d  3.248944  a -1.609652
  7  d  3.248944  b -2.837075
  </code></pre>
  
  <p>Since Seaborn plot commands (unlike <code>ggplot2</code> commands in R) appear to accept one <em>and only one</em> dataframe, we need to merge our predictions into the raw data:</p>
  
  <pre><code># append to df:
  merged = df.append(preds)
  </code></pre>
  
  <p>We can now plot the model predictions along with the data, with our continuous variable <code>x1</code> as the x-axis:</p>
  
  <pre><code># plot using seaborn:
  sns.set_style('white')
  sns.set_context('talk')
  g = sns.FacetGrid(merged, hue='x2', col='x3', size=5)
  # use the `map` method to add stuff to the facetgrid axes:
  g.map(plt.plot, "x1", "yhat")
  g.map(plt.scatter, "x1", "y")
  g.add_legend()
  g.fig.subplots_adjust(wspace=0.3)
  sns.despine(offset=10);
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/ZlxSX.png" alt="enter image description here"></p>
  
  <p>So far so good. Now imagine that we didn't measure the continuous variable <code>x1</code>, and we only know about the other two (categorical) variables (i.e., we have a 2x2 factorial design). <strong>How can we plot the model predictions against data in this case?</strong></p>
  
  <pre><code>fit = smf.ols('y ~ x2 * x3', df).fit()
  print(fit.summary())
  
  preds = expand_grid(
                  {'x2': ['a', 'b'],
                   'x3': ['c', 'd']})
  preds['yhat'] = fit.predict(preds)
  print(preds)
  
  # append to df:
  merged = df.append(preds)
  </code></pre>
  
  <p>Well, we can plot the model predictions using <code>sns.pointplot</code> or similar, like so:</p>
  
  <pre><code># plot using seaborn:
  g = sns.FacetGrid(merged, hue='x3', size=4)
  g.map(sns.pointplot, 'x2', 'yhat')
  g.add_legend();
  sns.despine(offset=10);
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/l6ung.png" alt="enter image description here"></p>
  
  <p>Or the data using <code>sns.factorplot</code> like so:</p>
  
  <pre><code>g = sns.factorplot('x2', 'y', hue='x3', kind='point', data=merged)
  sns.despine(offset=10);
  g.savefig('tmp3.png')
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/dOpq1.png" alt="enter image description here"></p>
  
  <p>But I do not see how to produce a plot similar to the first one (i.e. lines for model predictions using <code>plt.plot</code>, a scatter of points for data using <code>plt.scatter</code>). The reason is that the <code>x2</code> variable I'm trying to use as the x-axis is a string / object, so the pyplot commands don't know what to do with them.</p>