Download the Jupyter Notebook for this section: sklearn.ipynb

scikit-learn

[29]:
import pyhdfe
import numpy as np
from sklearn import datasets, linear_model

pyhdfe.__version__
[29]:
'0.1.0'

In this tutorial, we’ll use the boston data set from scikit-learn to demonstrate how pyhdfe can be used to absorb fixed effects before running regressions.

First, load the data set and create a matrix of fixed effect IDs. We’ll use a dummy for the Charles river and an index of accessibility to radial highways.

[30]:
boston = datasets.load_boston().data
ids = boston[:, [3, 8]]
ids
[30]:
array([[0., 1.],
       [0., 2.],
       [0., 2.],
       ...,
       [0., 1.],
       [0., 1.],
       [0., 1.]])

Next, choose our variables: per capita crime rate, proportion of residential land zoned for lots over 25,000 square feet, and proportion of non-retail business acres per town.

[31]:
variables = boston[:, :3]
variables
[31]:
array([[6.3200e-03, 1.8000e+01, 2.3100e+00],
       [2.7310e-02, 0.0000e+00, 7.0700e+00],
       [2.7290e-02, 0.0000e+00, 7.0700e+00],
       ...,
       [6.0760e-02, 0.0000e+00, 1.1930e+01],
       [1.0959e-01, 0.0000e+00, 1.1930e+01],
       [4.7410e-02, 0.0000e+00, 1.1930e+01]])

The create function initializes an Algorithm for fixed effect absorption that can residualize matrices with Algorithm.residualize. We’ll use the default algorithm. You may want to try other algorithms if it takes a long time to absorb fixed effects into your data.

[32]:
algorithm = pyhdfe.create(ids)
residualized = algorithm.residualize(variables)
residualized
[32]:
array([[-1.08723516e-01, -2.20167195e+01, -2.65583593e+00],
       [-5.59754167e-02, -2.04166667e+01, -2.56083333e+00],
       [-5.59954167e-02, -2.04166667e+01, -2.56083333e+00],
       ...,
       [-5.42835164e-02, -4.00167195e+01,  6.96416407e+00],
       [-5.45351644e-03, -4.00167195e+01,  6.96416407e+00],
       [-6.76335164e-02, -4.00167195e+01,  6.96416407e+00]])

We can now run a regression of per capita crime rate on the other two variables and our fixed effects.

[33]:
y = residualized[:, [0]]
X = residualized[:, 1:]
regression = linear_model.LinearRegression()
regression.fit(X, y)
regression.coef_
[33]:
array([[-6.97058632e-05,  5.53038164e-02]])