Inference (with sklearn)¶

sklearn Basics¶

Most sklearn models have a standardized interface. An estimator will have a fit method that takes X and y in the form of arrays. Once you've estimated the model, you can also use the predict method that takes X values and returns predicted y values. Certain estimators, especially classifiers, will have a combined fit_predict method that does both in a single step.

This library wont do a lot of specification work for you and expects numeric arras (sometimes called "design matrices"). However, it does have a lot of preprocessing capabilities for numerical data. To get from something like a pandas DataFrame to numeric arrays, you'll need an intermediate library. Some options include patsy, formulaic, and my very own fastreg. In fact, formulaic is in some sense the evolution of patsy, so you might want to stick with that.

First let's generate some fake data to play with.

We can generate design matrices by creating a Formula object with an R (or "Wilkinson") style formula specification that references the columns of data. We then apply that to our DataFrame yielding a design matrix for the endogenous variable (y) and the exogenous variables (X)

Next let's try a setting with categorical variables. We'll generate a fairly large dataset here just to make sure things scale well.

The syntax for generating categoricals is putting a C() around the name like a function. The design matrix for a categorical variable isn't just a single column with the value, since the value very well may not be a number. Instead it's a "one-hot" encoded matrix, where each row is a unit vector with a one in the position corresponding to which label that row is. So if the categorical was US states and our label was "NY", it would a be a length 50 vector with 49 zeros and a one in the position corresponding to the "NY" label. Even in this case, the matrix would be mostly zeros, and for larger categoricals, the problem would be even worse. For this reason, we ask it to return a sparse matrix.

Make sure you understand why the dimensions of X are what they are. Here's some info on the values for the first row.

Note that the way we encode categoricals is important. Since we have included an explicit constant (Intercept) term, if we had an indicator for each possible value of a categorical (say id1=0 through id1=9), our X matrix would be rendered singular, as the sum of each of the id1 columns would be the same as the contant column (all ones). So we need to drop one value for each categorical variable. Typically this is the first value encountered, though most libraries allow you to choose which value to drop, or to drop the most common value encountered. And below is the original data for the first row.

Let's look at the speed dimension

If you want to use fastreg, the syntax is pretty similar.

Linear Regression¶

First let's do our simple linear regression. Not so bad for N=10. Remember to aware of that intercept, as the default is fit_intercept=True, but we included the intercept in our design matrix already.

Next up are the sparse results. Note that in both cases y is returned as a matrix, so we need to turn it into just a column.

Our coefficients are all returned in one big vector. But we can extract the estimated values for the second categorical id2. Here we see that we get a pretty good fit, though the estimating has some issue differentiating between the intercept and the base value of the categorical.

Generalized Linear Models (GLM)¶

Got integer or count data? Don't sacrifice your ideals! Run a Poisson regression. This is but one option in the broader class of GLM. These are models in which the outcome variable (y) follows a certain distribution (say a Poisson distribution), and the parameters of this distribution (say the Poisson rate) are linear functions (or non-linear functions of linear forms) of the exogenous (X) variables.

Let's use a handy routine from fastreg to generate some random Poisson data.

Next we'll use fastreg again to generate the design matrices. Note that we could also have used formulaic here.

Finally, let's run the Poisson estimation. Note that GLM estimation isn't quite as friendly and stable your typical OLS, as it's running some kind of maximum likelihood or moment routine in the background. For instance the below routine won't actually succeed unless you add in alpha=1e-12. The default of alpha=0 leads to erroneous results.

Here we can again check the fit on the larger categorical and see that it's pretty good, though noisier than the OLS results.

We can also run a full Poisson regression with fastreg alone, though it uses a slower maximum likelihood routine.

Here again the coefficients line up pretty well.

K-means Clustering¶

Let's do a brief interlude on classification before we get into the heavier inference models. With K-means clustering, for a given integer value of $K$, we wish to find the $K$ centroids in a given space that minimize the mean-squared distance overall points to their nearest centroid. This is by far the most common spatial clustering technique. What we get out at the end is a list of centroids $c_k$, and an assignment of each point $x_i$ to their nearest centroid. The objective function for K-means is then

The downside to K-means is that you have to pick $K$ in advance. You can fiddle around with different values, but there's not way to choose it automatically. Some other clustering algorithms go a bit in this direction, but there's no silver bullet.

Now let's generate some random clusters!

Running the estimation is pretty straightforward. There are number of runtime parameters that are listed in the sklearn documentation. Finally, we extract the estimated centroid locations and labels (which centroid each point is assigned to, integers from 0 to K-1).

Plotting the estimated labels, we can see that we did a pretty good job considering the level of overlap!

Advanced Estimators¶

Next we'll move on to some of the more advanced regression options. These start to tackle some of the issues of model selection. In the regular OLS world we just kind of choose our regressors and hope for the best. Handling things like non-linearities and interaction terms is up to our intuition. This more advanced set of estimators mostly uses decision trees and bootstrapping or model averaging to implicity do model selection and achieve much better prediction outcomes with minimal fine tuning. The cost is that they are much harder to interpret. You don't just get a regression coefficient that you can point to.

We need to make some data sufficiently complex to make the more advanced learning methods sweat a little bit. This would include things like non-linearities and interaction terms. Below, I'm putting in a cubic term in x1 and an interaction term whereby the effect of x2 is dependent on the value of the categorical id1.

Random Forest¶

Same basic routine for making design matrices. Note that we just put in the raw variables here, figuring out non-linearities and interactions is up to the estimator. Our main hyperparameter here is the max_depth option which determines the number of layers in the decision tree. Try playing around with different numbers to see how it does.

Now we can compare this to how a regular probably misspecified linear regression would do.

Okay, what if we magically knew the true underlying model? Well, then we can do really well!

Gradient Boosting¶

Finally, we have gradient boosting. This achieves similar goals to a random forest, but can often achieve better prediction accuracy. Because it operates using decision trees (like random forest) it's sort dual use as both a classified and regressor. We can use the same data as before and check the observed fit. Looks pretty similar to random forest in this case.

There's another well known gradient boosting package called xgboost that is widely used in industry. This conveniently has the same interface type as an sklearn estimator. It has a ton of different options you can use, but below you can see that for the default options it does extremely well and is quite fast.

Lasso¶

Forgot about Lasso. Actually it does pretty bad in this case, but that's probably because we don't have any non-linearity

Let's see what happens if we manually introduce all cubic interaction terms and let Lasso figure it out. We're doing some advanced fastreg here too. If we set $\alpha$ very low, we reduce to properly specified OLS, but in practice we would want to use a non-zero $\alpha$.