Stata and R in a jupyter notebook¶

The jupyter notebook project is now designed to be a 'language agnostic' web-application front-end for any one of many possible software language kernels. We've been mostly using python but there are in fact several dozen other language kernels that can be made to work with it including Julia, R, Matlab, C, Go, Fortran and Stata.

The ecosystem of libraries and packages for scientific computing with python is huge and constantly growing but there are still many statistics and econometrics applications that are available as built-in or user-written modules in Stata that have not yet been ported to python or are just simply easier to use in Stata. On the other hand there are some libraries such as python pandas and different visualization libraries such as seaborn or matplotlib that give features that are not available in Stata.

Fortunately you don't have to choose between using Stata or python, you can use them both together, to get the best of both worlds.

Jupyter running an R kernel¶

R is a powerful open source software environment for statistical computing. R has R markdown which allows you to create R-markdown notebooks similar in concept to jupyter notebooks. But you can also run R inside a jupyter notebook (indeed the name 'Jupyter' is from Julia, iPython and R).

See my notebook with notes on Research Discontinuity Design for an example of a jupyter notebook running R. To install an R kernel see the IRkernel project.

Jupyter with a Stata Kernel¶

Kyle Barron has created a stata_kernel that offers several useful features including code-autocompletion, inline graphics, and generally fast responses.

For this to work you must have a working licensed copy of Stata version 14 or greater on your machine.

Python and Stata combined in the same notebook¶

Sometimes it may be useful to combine python and Stata in the same notebook. Ties de Kok has written a nice python library called ipystata that allows one to execute Stata code in codeblocks inside an ipython notebook when preceded by a %%stata magic command.

This workflow allows you to pass data between python and Stata sessions and to display Stata plots inline. Compared to the stata_kernel option the response times are not quite as fast.

The remainder of this notebook illustrates the use of ipystata.

A sample ipystata session¶

For more details see the example notebook and documentation on the ipystata repository.

In [1]:

%matplotlib inline
import seaborn as sns
import pandas as pd
import statsmodels.formula.api as smf
import ipystata

The following opens a Stata session where we load a dataset and summarize the data. The -o flag following the `%%Stata``` magic instructs it to output or return the dataset in Stata memory as a pandas dataframe in python.

In [2]:

%%stata -o life_df
sysuse lifeexp.dta
summarize

(Life expectancy, 1998)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
      region |         68         1.5    .7431277          1          3
     country |          0
   popgrowth |         68    .9720588    .9311918        -.5          3
        lexp |         68    72.27941    4.715315         54         79
       gnppc |         63    8674.857    10634.68        370      39980
-------------+---------------------------------------------------------
   safewater |         40        76.1    17.89112         28        100

Let's confirm the data was returned as a pandas dataframe:

In [3]:

life_df.head(3)

Out[3]:

	region	country	popgrowth	lexp	gnppc	safewater
0	Eur & C.Asia	Albania	1.2	72	810.0	76.0
1	Eur & C.Asia	Armenia	1.1	74	460.0	NaN
2	Eur & C.Asia	Austria	0.4	79	26830.0	NaN

A simple generate variable command and ols regression in Stata:

In [4]:

%%stata -o life_df
gen lngnppc = ln(gnppc)
regress lexp lngnppc

(5 missing values generated)

      Source |       SS           df       MS      Number of obs   =        63
-------------+----------------------------------   F(1, 61)        =     97.09
       Model |  873.264865         1  873.264865   Prob > F        =    0.0000
    Residual |  548.671643        61  8.99461709   R-squared       =    0.6141
-------------+----------------------------------   Adj R-squared   =    0.6078
       Total |  1421.93651        62  22.9344598   Root MSE        =    2.9991

------------------------------------------------------------------------------
        lexp |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     lngnppc |   2.768349   .2809566     9.85   0.000     2.206542    3.330157
       _cons |   49.41502   2.348494    21.04   0.000     44.71892    54.11113
------------------------------------------------------------------------------

And the same regression using statsmodels and pandas:

In [5]:

model = smf.ols(formula = 'lexp ~ lngnppc',
                data = life_df)
results = model.fit()
print(results.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   lexp   R-squared:                       0.614
Model:                            OLS   Adj. R-squared:                  0.608
Method:                 Least Squares   F-statistic:                     97.09
Date:                Tue, 08 May 2018   Prob (F-statistic):           3.13e-14
Time:                        14:25:48   Log-Likelihood:                -157.57
No. Observations:                  63   AIC:                             319.1
Df Residuals:                      61   BIC:                             323.4
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     49.4150      2.348     21.041      0.000      44.719      54.111
lngnppc        2.7683      0.281      9.853      0.000       2.207       3.330
==============================================================================
Omnibus:                       21.505   Durbin-Watson:                   1.518
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               46.140
Skew:                          -1.038   Prob(JB):                     9.57e-11
Kurtosis:                       6.643   Cond. No.                         52.7
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Back to ipython¶

Let's change one of the variables in the dataframe in python:

In [6]:

life_df.popgrowth = life_df.popgrowth * 100

In [7]:

life_df.popgrowth.mean()

Out[7]:

97.20587921142578

And now let's push the modified dataframe into the Stata dataset with the -d flag:

In [8]:

%%stata -d life_df
summarize

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
       index |         68        33.5    19.77372          0         67
      region |         68          .5    .7431277          0          2
     country |          0
   popgrowth |         68    97.20588    93.11918        -50        300
        lexp |         68    72.27941    4.715315         54         79
-------------+---------------------------------------------------------
       gnppc |         63    8674.857    10634.68        370      39980
   safewater |         40        76.1    17.89112         28        100
     lngnppc |         63    8.250023    1.355677   5.913503   10.59613

A Stata plot:

In [9]:

%%stata -d life_df --graph 
graph twoway (scatter lexp lngnppc) (lfit lexp lngnppc)

Now on the python side use lmplot from the seaborn library to graph a similar scatter and fitted line but by region.

In [10]:

sns.set_style("whitegrid")
g=sns.lmplot(y='lexp', x='lngnppc', col='region', data=life_df,col_wrap=2)