Hello everyone!

So far, we’ve covered 9 topics from Exploratory Data Analysis to Time Series Analysis in Python. Today we are going to have a look at one of the most popular and practical machine learning algorithms: gradient boosting. You can find more detailed math in the nbviewer format of the article.

Article Outline

1. Introduction and history of boosting
2. Gradient Boosting Machine algorithm
3. Loss functions
4. Assignment #10
5. Useful resources

1. Introduction and history of boosting

Almost everyone in machine learning has heard about gradient boosting. Many data scientists include this algorithm in their data scientist’s toolbox because of the good results it yields on any given (unknown) problem.

Furthermore, XGBoost is often the standard recipe for winning ML competitions. It is so popular that the idea of stacking XGBoosts has become a meme. Moreover, boosting is an important component in many recommender systems; sometimes, it is even considered a brand. Let’s look at the history and development of boosting.

Boosting was born out of the question: is it possible to get one strong model from a large amount of relatively weak and simple models?
By saying “weak models”, we do not mean simple basic models like decision trees but models with poor accuracy performance, where poor is a little bit better than random.

A positive mathematical answer to this question was identified, but it took a few years to develop fully functioning algorithms based on this solution e.g. AdaBoost. These algoritms take a greedy approach: first, they build a linear combination of simple models (basic algorithms) by re-weighing the input data. Then, the model (usually a decision tree) is built on earlier incorrectly predicted objects, which are now given larger weights.

Many machine learning courses study AdaBoost — the ancestor of GBM (Gradient Boosting Machine). However, since AdaBoost merged with GBM, it has become apparent that AdaBoost is just a particular variation of GBM.

The algorithm itself has a very clear visual interpretation and intuition for defining weights. Let’s have a look at the following toy classification problem where we are going to split the data between the trees of depth 1 (also known as ‘stumps’) on each iteration of AdaBoost. For the first two iterations, we have the following picture:

The size of point corresponds to its weight, which was assigned for an incorrect prediction. On each iteration, we can see that these weights are growing — the stumps cannot cope with this problem. Although, if we take a weighted vote for the stumps, we will get the correct classifications:

Here is a more detailed example of AdaBoost where, as we iterate, we can see the weights increase, especially on the border between classes.

AdaBoost works well, but the lack of explanation for why the algorithm is successful sewed the seeds of doubt. Some considered it a super-algorithm, a silver bullet, but others were skeptical and believed AdaBoost was just overfitting.

The overfitting problem did indeed exist, especially when data had strong outliers. Therefore, in those types of problems, AdaBoost was unstable. Fortunately, a few professors in the statistics department at Stanford, who had created Lasso, Elastic Net, and Random Forest, started researching the algorithm. In 1999, Jerome Friedman came up with the generalization of boosting algorithms development — Gradient Boosting (Machine), also known as GBM. With this work, Friedman set up the statistical foundation for many algorithms providing the general approach of boosting for optimization in the functional space.

CART, bootstrap, and many other algorithms have originated from Stanford’s statistics department. In doing so, the department has solidified their names in future textbooks. These algorithms are very practical, and some recent works have yet to be widely adopted. For example, check out glinternet.

Not many video recordings of Friedman are available. Although, there is a very interesting interview with him about the creation of CART and how they solved statistics problems (which is similar to data analysis and data science today) more than 40 years ago. There is also a great lecture from Hastie, a retrospective on data analysis from one of the creators of methods that we use everyday.

In general, there has been a transition from engineering and algorithmic research to a full-fledged approach to building and studying algorithms. From a mathematical perspective, this is not a big change — we are still adding (or boosting) weak algorithms and enlarging our ensemble with gradual improvements for parts of the data where the model was inaccurate. But, this time, the next simple model is not just built on re-weighted objects but improves its approximation of the gradient of overall objective function. This concept greatly opens up our algorithms for imagination and extensions.

History of GBM

It took more than 10 years after the introduction of GBM for it to become an essential part of the data science toolbox. GBM was extended to apply to different statistics problems: GLMboost and GAMboost for strengthening already existing GAM models, CoxBoost for survival curves, and RankBoost and LambdaMART for ranking.

Many realizations of GBM also appeared under different names and on different platforms: Stochastic GBM, GBDT (Gradient Boosted Decision Trees), GBRT (Gradient Boosted Regression Trees), MART (Multiple Additive Regression Trees), and more. In addition, the ML community was very segmented and dissociated, which made it hard to track just how widespread boosting had become.

At the same time, boosting had been actively used in search ranking. This problem was rewritten in terms of a loss function that penalizes errors in the output order, so it became convenient to simply insert it into GBM. AltaVista was one of the first companies who introduced boosting to ranking. Soon, the ideas spread to Yahoo, Yandex, Bing, etc. Once this happened, boosting became one of the main algorithms that was used not only in research but also in core technologies in industry.

ML competitions, especially Kaggle, played a major role in boosting’s popularization. Now, researchers had a common platform where they could compete in different data science problems with large number of participants from around the world. With Kaggle, one could test new algorithms on the real data, giving algoritms oppurtunity to “shine”, and provide full information in sharing model performance results across competition data sets. This is exactly what happened to boosting when it was used at Kaggle (check interviews with Kaggle winners starting from 2011 who mostly used boosting). The XGBoost library quickly gained popularity after its appearance. XGBoost is not a new, unique algorithm; it is just an extremely effective realization of classic GBM with additional heuristics.

This algorithm has gone through very typical path for ML algorithms today: mathematical problem and algorithmic crafts to successful practical applications and mass adoption years after its first appearance.

2. Gradient Boosting Machine algorithm

We are going to solve the problem of function approximation in a general supervised learning setting. We have a set of features X and target variables y which we use to restore the dependence y = f(x). We restore the dependence by approximating f(x) and by understanding which approximation is better when we use the loss function L(y,f) , which we want to minimize:

At this moment, we do not make any assumptions regarding the type of dependence f(x), the model of our approximation, or the distribution of the target variable. We only expect that the function L(y, f) is differentiable. Our formula is very general; let’s define it for a particular data set with a population mean. Our expression for minimizing the loss of the data is the following:

Unfortunately, the number of these functions is not just large, its functional space is infinite-dimensional. That is why it is acceptable for us to limit the search space by some family of functions. This simplifies the objective a lot because now we have a solvable optimization of parameter values.

Simple analytical solutions for finding the optimal parameters often do not exist, so the parameters are usually approximated iteratively. To start, we write down the empirical loss function that will allow us to evaluate our parameters using our data. Additionally, let’s write out our approximation for a number of M iterations as a sum:

Then, the only thing left is to find a suitable, iterative algorithm to minimize the last expression. Gradient descent is the simplest and most frequently used option. We define the gradient and add our iterative evaluations to it (since we are minimizing the loss, we add the minus sign). Our last step is to initialize our first approximation and choose the number of iterations M. Let’s review the steps for this inefficient and naive algorithm:

Functional gradient descent

Let’s imagine for a second that we can perform optimization in the function space and iteratively search for the approximations as functions themselves. We will express our approximation as a sum of incremental improvements, each being a function.

Nothing has happened yet; we have only decided that we will search for our approximation not as a big model with plenty of parameters (as an example, neural network), but as a sum of functions, pretending we move in functional space.

In order to accomplish this task, we need to limit our search by some function family. There are a few issues here — first of all, the sum of models can be more complicated than any model from this family; secondly, the general objective is still in functional space. Let’s note that, on every step, we will need to select an optimal coefficient. For step t, the problem is the following:

Here is where the magic happens. We have defined all of our objectives in general terms, as if we could have trained any kind of model for any type of loss functions. In practice, this is extremely difficult, but, fortunately, there is a simple way to solve this task.

Knowing the expression of loss function’s gradient, we can calculate its value on our data. So, let’s train the models such that our predictions will be more correlated with this gradient (with a minus sign). In other words, we will use least squares to correct the predictions with these residuals. For classification, regression, and ranking tasks, we will minimize the squared difference between pseudo-residuals r and our predictions. For step t, the final problem looks the following:

Friedman’s classic GBM algorithm

We can now define the classic GBM algorithm suggested by Jerome Friedman in 1999. It is a supervised algorithm that has the following components:

• a dataset;
• number of iterations M;
• choice of loss function with a defined gradient;
• choice of function family of base algorithms with the training procedure;
• additional hyperparameters (for example, in decision trees, the tree depth);

The only thing left is the initial approximation. For simplicity, for an initial approximation, a constant value gamma is used. The constant value, as well as the optimal coefficient, are identified via binary search or another line search algorithm over the initial loss function (not a gradient). So, we have our GBM algorithm described as follows:

Step-By-Step example: How GBM Works

Let’s see an example of how GBM works. In this toy example, we will restore a noisy function

This is a regression problem with a real-valued target, so we will choose to use the mean squared error loss function. We will generate 300 pairs of observations and approximate them with decision trees of depth 2. Let’s put together everything we need to use GBM:

• Toy data
• Number of iterations M = 3
• The mean squared error loss function L(y, f) = (y - f)²
• Gradient of L(y, f) is just residuals r = (y-f)
• Decision trees as base algorithms
• Hyperparameters of the decision trees: trees depth is equal to 2

We will run GBM and draw two types of graphs: the current approximation (blue graph) and every tree built on its pseudo-residuals (green graph). The graph’s number corresponds to the iteration number:

By the second iteration, our trees have recovered the basic form of the function. However, at the first iteration, we see that the algorithm has built only the “left branch” of the function. This was due to the fact that our trees simply did not have enough depth to build a symmetrical branch at once, and it focused on the left branch with the larger error. Therefore, the right branch appeared only after the second iteration.

The rest of the process goes as expected — on every step, our pseudo-residuals decreased, and GBM approximated the original function better and better with each iteration. However, by construction, trees cannot approximate a continuous function, which means that GBM is not ideal in this example. To play with GBM function approximations, you can use the awesome interactive demo in this blog called Brilliantly wrong:

3. Loss functions

If we want to solve a classification problem instead of regression, what would change? We only need to choose a suitable loss function L(y, f). This is the most important, high-level moment that determines exactly how we will optimize and what characteristics we can expect in the final model.

As a rule, we do not need to invent this ourselves — researchers have already done it for us. Today, we will explore loss functions for the two most common objectives: regression and binary classification.

Regression loss functions

Let’s start with a regression problem for y, a real number. In order to choose the appropriate loss function, we need to consider which of the properties of the conditional distribution (y|x) we want to restore. The most common options are:

• L(y, f) = (y — f)² a.k.a. L2 loss or Gaussian loss. It is the classical conditional mean, which is the simplest and most common case. If we do not have any additional information or requirements for a model to be robust, we can use the Gaussian loss.
• L(y, f) = |y — f| a.k.a. L1 loss or Laplacian loss. At the first glance, this function does not seem to be differentiable, but it actually defines the conditional median. Median, as we know, is robust to outliers, which is why this loss function is better in some cases. The penalty for big variations is not as heavy as it is in L2.
• Lq loss or Quantile loss. Instead of median, it uses quantiles. We can see that this function is asymmetric and penalizes the observations which are on the right side of the defined quantile.

Let’s use loss function L_q on our data. The goal is to restore the conditional 75%-quantile of cosine. Let us put everything together for GBM:

For our initial approximation, we will take the needed quantile of y. However, we do not know anything about optimal coefficients, so we’ll use standard line search. The results are the following:

The overall results of GBM with quantile loss function are the same as the results with quadratic loss function offset by ~ 0.135. But if we were to use the 90%-quantile, we would not have enough data due to the fact that classes would become unbalanced. We need to remember this when we deal with non-standard problems.

For regression tasks, many loss functions have been developed, some of them with extra properties. For example, they can be robust like in the Huber loss function. For a small number of outliers, the loss function works as L2, but after a defined threshold, the function changes to L1. This allows for decreasing the effect of outliers and focusing on the overall picture.

We can illustrate this with the following example. Data is generated from the function y = sin(x) / x with added noise, a mixture from normal and Bernoulli distributions. We show the functions on graphs A-D and the relevant GBM on F-H (graph E represents the initial function):

Original size.

In this example, we used splines as the base algorithm. See, it does not always have to be trees for boosting?

We can clearly see the difference between the functions L2, L1, and Huber loss. If we choose optimal parameters for the Huber loss, we can get the best possible approximation among all our options. The difference can be seen as well in the 10%, 50%, and 90%-quantiles.

Unfortunately, Huber loss function is supported only by very few popular libraries/packages; h2o supports it, but XGBoost does not. It is relevant to other things that are more exotic like conditional expectiles, but it may still be interesting knowledge.

Classification loss functions

Now, let’s look at the binary classification problem. We saw that GBM can even optimize non-differentiable loss functions. Technically, it is possible to solve this problem with a regression L2 loss, but it wouldn’t be correct.

The distribution of the target variable requires us to use log-likehood, so we need to have different loss functions for targets multiplied by their predictions. The most common choices would be the following:

Let’s generate some new toy data for our classification problem. As a basis, we will take our noisy cosine, and we will use the sign function for classes of the target variable. Our toy data looks like the following (jitter-noise is added for clarity):

We will use logistic loss to look for what we actually boost. So, again, we put together what we will use for GBM:

This time, the initialization of the algorithm is a little bit harder. First, our classes are imbalanced (63% versus 37%). Second, there is no known analytical formula for the initialization of our loss function, so we have to look for it via search:

Our optimal initial approximation is around -0.273. You could have guessed that it was negative because it is more profitable to predict everything as the most popular class, but there is no formula for the exact value. Now let’s finally start GBM, and look what actually happens under the hood:

The algorithm successfully restored the separation between our classes. You can see how the “lower” areas are separating because the trees are more confident in the correct prediction of the negative class and how the two steps of mixed classes are forming. It is clear that we have a lot of correctly classified observations and some amount of observations with large errors that appeared due to the noise in the data.

Weights

Sometimes, there is a situation where we want a more specific loss function for our problem. For example, in financial time series, we may want to give bigger weight to large movements in the time series; for churn prediction, it is more useful to predict the churn of clients with high LTV (or lifetime value: how much money a client will bring in the future).

The statistical warrior would invent their own loss function, write out the gradient for it (for more effective training, include the Hessian), and carefully check whether this function satisfies the required properties. However, there is a high probability of making a mistake somewhere, running up against computational difficulties, and spending an inordinate amount of time on research.

In lieu of this, a very simple instrument was invented (which is rarely remembered in practice): weighing observations and assigning weight functions. The simplest example of such weighting is the setting of weights for class balance. In general, if we know that some subset of data, both in the input variables and in the target variable, has greater importance for our model, then we just assign them a larger weight w(x,y). The main goal is to fulfill the general requirements for weights:

Weights can significantly reduce the time spent adjusting the loss function for the task we are solving and also encourages experiments with the target models’ properties. Assigning these weights is entirely a function of creativity. We simply add scalar weights:

It is clear that, for arbitrary weights, we do not know the statistical properties of our model. Often, linking the weights to the values y can be too complicated. For example, the usage of weights proportional to |y| in L1 loss function is not equivalent to L2 loss because the gradient will not take into account the values of the predictions themselves.

We mention all of this so that we can understand our possibilities better. Let’s create some very exotic weights for our toy data. We will define a strongly asymmetric weight function as follows:

With these weights, we expect to get two properties: less detailing for negative values of X and the form of the function, similar to the initial cosine. We take the other GBM’s tunings from our previous example with classification including the line search for optimal coefficients. Let’s look what we’ve got:

We achieved the result that we expected. First, we can see how strongly the pseudo-residuals differ; on the initial iteration, they look almost like the original cosine. Second, the left part of the function’s graph was often ignored in favor of the right one, which had larger weights. Third, the function that we got on the third iteration received enough attention and started looking similar to the original cosine (also started to slightly overfit).

Weights are a powerful but risky tool that we can use to control the properties of our model. If you want to optimize your loss function, it is worth trying to solve a more simple problem first but add weights to the observations at your discretion.

Conclusion

Today, we learned the theory behind gradient boosting. GBM is not just some specific algorithm but a common methodology for building ensembles of models. In addition, this methodology is sufficiently flexible and expandable — it is possible to train a large number of models, taking into consideration different loss-functions with a variety of weighting functions.

Practice and ML competitions show that, in standard problems (except for image, audio, and very sparse data), GBM is often the most effective algorithm (not to mention stacking and high-level ensembles, where GBM is almost always a part of them). Also, there are many adaptations of GBM for Reinforcement Learning (Minecraft, ICML 2016). By the way, the Viola-Jones algorithm, which is still used in computer vision, is based on AdaBoost.

In this article, we intentionally omitted questions concerning GBM’s regularization, stochasticity, and hyper-parameters. It was not accidental that we used a small number of iterations M = 3 throughout. If we used 30 trees instead of 3 and trained the GBM as described, the result would not be that predictable:

4. Assignment #10

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

5. Useful resources

• Original article about GBM by Jerome Friedman
• “Gradient boosting machines, a tutorial”, paper by Alexey Natekin, and Alois Knoll
• Chapter in Elements of Statistical Learning from Hastie, Tibshirani, Friedman (page 337)
• Wiki article about Gradient Boosting
• Frontiers tutorial article about GBM
• Video-lecture by Hastie about GBM at h2o.ai conference

Author: Alexey Natekin, OpenDataScience founder, Machine Learning Evangelist. Translated and edited by Olga Daykhovskaya, Anastasia Manokhina, Egor Polusmak, and Yuanyuan Pao.

Open Machine Learning Course. Topic 10. Gradient Boosting was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

Time series forecasting finds wide application in data analytics. These are only some of the conceivable predictions of future trends that might be useful:

• The number of servers that an online service will need next year;
• The demand for a grocery product at a supermarket on a given day;
• The tomorrow closing price of a tradable financial asset.

For another example, we can make a prediction of some team’s performance and then use it as a baseline: first to set goals for the team, and then to measure the actual team performance relative to the baseline.

There are quite a few different methods to predict future trends, for example, ARIMA, ARCH, regressive models, neural networks.

In this article, we will look at Prophet, a library for time series forecasting released by Facebook and open-sourced on February 23, 2017. We will also try it out in the problem of predicting the daily number of posts published on Medium.

Article outline

1. Introduction

2. The Prophet Forecasting Model

3. Practice with Prophet

• 3.1 Installation in Python
• 3.2 Dataset
• 3.3 Exploratory visual analysis
• 3.4 Making a forecast
• 3.5 Forecast quality evaluation
• 3.6 Visualization

4. Box-Cox Transformation

5. Summary

6. References

The following material is better viewed as a Jupyter notebook and can be reproduced locally with Jupyter if you clone the course repository.

1. Introduction

According to the article by Facebook Research, Prophet was initially developed for the purpose of creating high quality business forecasts. This library tries to address the following difficulties common to many business time series:

• Seasonal effects caused by human behavior: weekly, monthly and yearly cycles, dips and peaks on public holidays;
• Changes in trend due to new products and market events;
• Outliers.

The authors claim that, even with default settings, in many cases, their library produces forecasts as accurate as those delivered by experienced analysts.

Moreover, Prophet has a number of intuitive and easily interpretable customizations that allow gradually improving the quality of the forecasting model. What is especially important, these parameters are quite comprehensible even for non-experts in time series analysis, which is a field of data science requiring certain skill and experience.

By the way, the original article is called “Forecasting at Scale”, but it is not about the scale in the “usual” sense, that it’s addressing computational and infrastructure problems of a large number of working programs. According to the authors, Prophet should scale well in the following 3 areas:

• Accessibility to a wide audience of analysts, possibly without profound expertise in time series;
• Applicability to a wide range of distinct forecasting problems;
• Automated performance estimation of a large number of forecasts including flagging of potential problems for their subsequent inspection by the analyst.

2. The Prophet Forecasting Model

Now, let’s take a closer look at how Prophet works. In its essence, this library utilizes the additive regression model y(t) comprising the following components:

where:

• Trend g(t) models non-periodic changes;
• Seasonality s(t) represents periodic changes;
• Holidays component h(t) contributes information about holidays and events.

Below, we will consider some important properties of these model components.

Trend

The Prophet library implements two possible trend models for g(t).

The first one is called Nonlinear, Saturating Growth. It is represented in the form of the logistic growth model:

where:

• C is the carrying capacity (that is the curve’s maximum value).
• k is the growth rate (which represents “the steepness” of the curve).
• m is an offset parameter.

This logistic equation allows modelling non-linear growth with saturation, that is when the growth rate of a value decreases with its growth. One of the typical examples would be representing the growth of the audience of an application or a website.

Actually, CC and kk are not necessarily constants and may vary over time. Prophet supports both automatic and manual tuning of their variability. The library can itself choose optimal points of trend changes by fitting the supplied historical data.

Also, Prophet allows analysts to manually set changepoints of the growth rate and capacity values at different points in time. For instance, analysts may have insights about dates of past releases that prominently influenced some key product indicators.

The second trend model is a simple Piecewise Linear Model with a constant rate of growth. It is best suited for problems without saturating growth.

Seasonality

The seasonal component s(t) provides a flexible model of periodic changes due to weekly and yearly seasonality.

Weekly seasonal data is modeled with dummy variables. Six new variables are added: monday, tuesday, wednesday, thursday, friday, saturday, which take values 0 or 1 depending on the day of the week. The feature sunday is not added because it would be a linear combination of the other days of the week, and this fact would have an adverse effect on the model.

Yearly seasonality model in Prophet relies on Fourier series.

Since version 0.2 you can also use sub-daily time series and make sub-daily forecasts as well as employ the new daily seasonality feature.

Holidays and Events

The component h(t) represents predictable abnormal days of the year including those on irregular schedules, e.g., Black Fridays.

To utilize this feature, the analyst needs to provide a custom list of events.

Error

The error term ϵ(t) represents information that was not reflected in the model. Usually it is modeled as normally distributed noise.

Prophet Benchmarking

For a detailed description of the model and algorithms behind Prophet refer to the paper “Forecasting at scale” by Sean J. Taylor and Benjamin Letham.

The authors also compared their library with several other methods for time series forecasting. They used Mean Absolute Percentage Error (MAPE) as a measure of prediction accuracy. In this research, Prophet has shown substantially lower forecasting error than the other models.

Let’s look closer at how the forcasting quality was measured in the article. To do this, we will need the formula of Mean Absolute Percentage Error.

Let yi be the actual (historical) value and ŷi be the forecast value given by our model.

Then ei=yi−ŷi is the forecast error and pi=ei/yi is the relative forecast error.

We define

MAPE is widely used as a measure of prediction accuracy because it expresses error as a percentage and thus can be used in model evaluations on different datasets.

In addition, when evaluating a forecasting algorithm, it may prove useful to calculate MAE (Mean Absolute Error) in order to have a picture of errors in absolute numbers. Using previously defined components, its equation will be

A few words about the algorithms that Prophet was compared with. Most of them are quite simple and often are used as a baseline for other models:

• naive is a simplistic forecasting approach where we predict all future values relying solely on the observation at the last available point of time.
• snaive (seasonal naive) is a model that makes constant predictions taking into account information about seasonality. For instance, in the case of weekly seasonal data for each future Monday, we would predict the value from the last Monday, and for all future Tuesdays we would use the value from the last Tuesday and so on.
• mean uses the averaged value of data as a forecast.
• arima stands for Autoregressive Integrated Moving Average, see Wikipedia for details.
• ets stands for Exponential Smoothing, see Wikipedia for more.

3. Practice with Facebook Prophet

3.1 Installation in Python

First, you need to install the library. Prophet is available for Python and R. The choice will depend on your personal preferences and project requirements. Further in this article we will use Python.

In Python you can install Prophet using PyPI:

$ pip install fbprophet

In R you can find the corresponing CRAN package. Refer to the documentation for details.

Let’s import the modules that we will need, and initialize our environment:

3.2 Dataset

We will predict the daily number of posts published on Medium.

First, we load our dataset:

df = pd.read_csv('../../data/medium_posts.csv', sep='\t')

Next, we leave out all columns except published and url. The former corresonds to the time dimension while the latter uniquely identifies a post by its URL. Along the way we get rid of possible duplicates and missing values in the data:

df = df[['published', 'url']].dropna().drop_duplicates()

Next, we need to convert published to the datetime format because by default pandas treats this field as string-valued.

df['published'] = pd.to_datetime(df['published'])

Let’s sort the dataframe by time and take a look at what we’ve got:

df.sort_values(by=['published']).head(n=3)

Medium’s public release date was August 15, 2012. But, as you can see from the data above, there are at least several rows with much earlier publication dates. They have somehow turned up in our dataset, but they are hardly legitimate ones. We will just trim our time series to keep only those rows that fall onto the period from August 15, 2012 to June 25, 2017:

df = df[(df['published'] > '2012-08-15') & 
        (df['published'] < '2017-06-26')].
         sort_values(by=['published'])
df.head(n=3)

df.tail(n=3)

As we are going to predict the number of published posts, we will aggregate and count unique posts at each given point in time. We will name the corresponding new column posts:

aggr_df = df.groupby('published')[['url']].count()
aggr_df.columns = ['posts']

In this practice, we are interested in the number of posts a day. But at this moment all our data is divided into irregular time intervals that are less than a day. This is called a sub-daily time series. To see it, let’s print out the first 3 rows:

aggr_df.head(n=3)

To fix this, we need to aggregate the post counts by “bins” of a date size. In time series analysis, this process is referred to as resampling. And if we reduce the sampling rate of data it is often called downsampling.

Luckily, pandas has a built-in functionality for this task. We will resample our time index down to 1-day bins:

daily_df = aggr_df.resample('D').apply(sum)
daily_df.head(n=3)

3.3 Exploratory visual analysis

As always, it may be helpful and instructive to look at a graphical representation of your data.

We will create a time series plot for the whole time range. Displaying data over such a long period of time can give clues about seasonality and conspicuous abnormal deviations.

First, we import and initialize the Plotly library, which allows creating beautiful interactive plots:

We also define a helper function, which will plot our dataframes throughout the article:

Let’s try and plot our dataset as is:

plotly_df(daily_df, title='Posts on Medium (daily)')

High-frequency data can be rather difficult to analyze. Even with the ability to zoom in provided by Plotly, it is hard to infer anything meaningful from this chart apart from the prominent upward and accelerating trend.

To reduce the noise, we will resample the post counts down to weekly bins. Besides binning, other possible techniques of noise reduction include Moving-Average Smoothing and Exponential Smoothing, among others.

We save our downsampled dataframe in a separate variable because further in this practice we will work only with daily series:

weekly_df = daily_df.resample('W').apply(sum)

Finally, we plot the result:

plotly_df(weekly_df, title='Posts on Medium (weekly)')

This downsampled chart proves to be somewhat better for an analyst’s perception.

One of the most useful functions that Plotly provides is the ability to quickly dive into different periods of timeline in order to better understand the data and find visual clues about possbile trends, periodic and irregular effects.

For example, zooming-in on a couple of consecutive years shows us time points corresponding to Christmas holidays, which greatly influence human behaviors.

Now, we’re going to omit the first few years of observations, up to 2015. First, they won’t contribute much into the forecast quality in 2017. Second, these first years, having very low number of posts per day, are likely to increase noise in our predictions, as the model would be forced to fit this abnormal historical data along with more relevant and indicative data from the recent years.

daily_df = daily_df.loc[daily_df.index >= '2015-01-01']
daily_df.head(n=3)

To sum up, from visual analysis we can see that our dataset is non-stationary with a prominent growing trend. It also demonstrates weekly and yearly seasonality and a number of abnormal days in each year.

3.4 Making a forecast

Prophet’s API is very similar to the one you can find in sklearn. First we create a model, then call the method fit, and, finally, make a forecast. The input to the method fit is a DataFrame with two columns:

• ds (datestamp) must be of type date or datetime.
• y is a numeric value we want to predict.

To get started, we’ll import the library and mute unimportant diagnostic messages:

from fbprophet import Prophet
import logging

logging.getLogger().setLevel(logging.ERROR)

Let’s convert our dataframe to the format required by Prophet:

df = daily_df.reset_index()
df.columns = ['ds', 'y']
df.tail(n=3)

The authors of the library generally advise to make predictions based on at least several months, ideally, more than a year of historical data. Luckily, in our case we have more than a couple of years of data to fit the model.

To measure the quality of our forecast, we need to split our dataset into the historical part, which is the first and biggest slice of our data, and the prediction part, which will be located at the end of the timeline. We will remove the last month from the dataset in order to use it later as a prediction target:

prediction_size = 30
train_df = df[:-prediction_size]
train_df.tail(n=3)

Now we need to create a new Prophet object. Here we can pass the parameters of the model into the constructor. But in this article we will use the defaults. Then we train our model by invoking its fit method on our training dataset:

m = Prophet()
m.fit(train_df);

Using the helper method Prophet.make_future_dataframe, we create a dataframe which will contain all dates from the history and also extend into the future for those 30 days that we left out before.

future = m.make_future_dataframe(periods=prediction_size)
future.tail(n=3)

We predict values with Prophet by passing in the dates for which we want to create a forecast. If we also supply the historical dates (as in our case), then in addition to the prediction we will get an in-sample fit for the history. Let's call the model's predict method with our future dataframe as an input:

forecast = m.predict(future)
forecast.tail(n=3)

In the resulting dataframe you can see many columns characterizing the prediction, including trend and seasonality components as well as their confidence intervals. The forecast itself is stored in the yhat column.

The Prophet library has its own built-in tools for visualization that enable us to quickly evaluate the result.

First, there is a method called Prophet.plot that plots all the points from the forecast:

m.plot(forecast);

This chart doesn’t look very informative. The only definitive conclusion that we can draw here is that the model treated many of the data points as outliers.

The second function Prophet.plot_components might be much more useful in our case. It allows us to observe different components of the model separately: trend, yearly and weekly seasonality. In addition, if you supply information about holidays and events to your model, they will also be shown in this plot.

Let’s try it out:

m.plot_components(forecast);

As you can see from the trend graph, Prophet did a good job by fitting the accelerated growth of new posts at the end of 2016. The graph of weekly seasonality leads to the conclusion that usually there are less new posts on Saturdays and Sundays than on the other days of the week. In the yearly seasonality graph there is a prominent dip on Christmas Day.

3.5 Forecast quality evaluation

Let’s evaluate the quality of the algorithm by calculating the error metrics for the last 30 days that we predicted. For this, we will need the observations yi and the corresponding predicted values ŷi.

Let’s look into the object forecast that the library created for us:

print(', '.join(forecast.columns))

ds, trend, trend_lower, trend_upper, yhat_lower, yhat_upper, seasonal, seasonal_lower, seasonal_upper, seasonalities, seasonalities_lower, seasonalities_upper, weekly, weekly_lower, weekly_upper, yearly, yearly_lower, yearly_upper, yhat

We can see that this dataframe contains all the information we need except for the historical values. We need to join the forecast object with the actual values y from the original dataset df. For this we will define a helper function that we will reuse later:

Let’s apply this function to our last forecast:

cmp_df = make_comparison_dataframe(df, forecast)
cmp_df.tail(n=3)

We are also going to define a helper function that we will use to gauge the quality of our forecasting with MAPE and MAE error measures:

Let’s use our function:

for err_name, err_value in 
 calculate_forecast_errors(cmp_df, prediction_size).items():
    print(err_name, err_value)

MAPE 22.7243579814
MAE 70.452857085

As a result, the relative error of our forecast (MAPE) is about 22.72%, and on average our model is wrong by 70.45 posts (MAE).

3.6 Visualization

Let’s create our own visualization of the model built by Prophet. It will comprise the actual values, forecast and confidence intervals.

First, we will plot the data for a shorter period of time to make the data points easier to distinguish. Second, we will show the model performance only for the period that we predicted, that is the last 30 days. It seems that these two measures should give us a more legible plot.

Third, we will use Plotly to make our chart interactive, which is great for exploring.

We will define a custom helper function show_forecast and call it (for more on how it works please refer to the comments in the code and the documentation):

At first glance, the prediction of the mean values by our model seems to be sensible. The high value of MAPE that we got above may be explained by the fact that the model failed to catch on to increasing peak-to-peak amplitude of weakly seasonality.

Also, we can conclude from the graph above that many of the actual values lie outside the confidence interval. Prophet may not be suitable for time series with unstable variance, at least when the default settings are used. We will try to fix this by applying a transform to our data.

4. Box-Cox Transformation

So far we have used Prophet with the default settings and the original data. We will leave the parameters of the model alone. But despite this we still have some room for improvement. In this section, we will apply the Box–Cox transformation to our original series. Let’s see where it will lead us.

A few words about this transformation. This is a monotonic data transformation that can be used to stabilize variance. We will use the one-parameter Box–Cox transformation, which is defined by the following expression:

We will need to implement the inverse of this function in order to be able to restore the original data scale. It is easy to see that the inverse is defined as:

The corresponding function in Python is implemented as follows:

First, we prepare our dataset by setting its index:

train_df2 = train_df.copy().set_index('ds')

Then, we apply the function stats.boxcox from Scipy, which applies the Box–Cox transformation. In our case it will return two values. The first one is the transformed series and the second one is the found value of λ that is optimal in terms of the maximum log-likelihood:

train_df2['y'], lambda_prophet = stats.boxcox(train_df2['y'])
train_df2.reset_index(inplace=True)

We create a new Prophet model and repeat the fit-predict cycle that we have already done above:

m2 = Prophet()
m2.fit(train_df2)
future2 = m2.make_future_dataframe(periods=prediction_size)
forecast2 = m2.predict(future2)

At this point, we need to revert the Box–Cox transformation with our inverse function and the known value of λ:

for column in ['yhat', 'yhat_lower', 'yhat_upper']:
    forecast2[column] = inverse_boxcox(forecast2[column],
                                       lambda_prophet)

Here we will reuse our tools for making the comparison dataframe and calculating the errors:

cmp_df2 = make_comparison_dataframe(df, forecast2)
for err_name, err_value in calculate_forecast_errors(cmp_df2, prediction_size).items():
     print(err_name, err_value)

MAPE 11.5921879552
MAE 39.072031256

So, we can definitely state an increase in the quality of the model.

Finally, let’s plot our previous performance with the latest results side-by-side. Note that we use prediction_size for the third parameter in order to zoom in on the interval being predicted:

show_forecast(cmp_df, prediction_size, 100, 'No transformations')
show_forecast(cmp_df2, prediction_size, 100, 'Box–Cox transformation')

We see that the forecast of weekly changes in the second graph is much closer to the real values now.

5. Summary

We have taken a look at Prophet, an open-source forecasting library that is specifically targeted at business time series. We have also done some hands-on practice in time series prediction.

As we have seen, the Prophet library does not make wonders, and its predictions out-of-box are not ideal. It is still up to the data scientist to explore the forecast results, tune model parameters and transform data when necessary.

However, this library is user-friendly and easily customizable. The sole ability to take into account abnormal days that are known to the analyst beforehand might make a difference in some cases.

All in all, the Prophet library is worth being a part of your analytical toolbox.

6. References

• Official Prophet repository on GitHub.
• Official Prophet documentation.
• Sean J. Taylor, Benjamin Letham “Forecasting at scale” — scientific paper explaining the algorithm which lays the foundation of Prophet.
• Forecasting Website Traffic Using Facebook’s Prophet Library — Prophet overview with an example of website traffic forecasting.
• Rob J. Hyndman, George Athanasopoulos “Forecasting: principles and practice” — a very good online book about time series forecasting.

Author: Egor Polusmak, Head of Data Analytics at Mail.Ru Search.

Open Machine Learning Course. Topic 9. Part 2. Predicting the future with Facebook Prophet was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

Hi there!

We continue our open machine learning course with a new article on time series.

Let’s take a look at how to work with time series in Python, what methods and models we can use for prediction; what’s double and triple exponential smoothing; what to do if stationarity is not you favorite game; how to build SARIMA and stay alive; how to make predictions using xgboost. All of this will be applied to (harsh) real world example.

For more data science articles consider following my blog: DataWondering.com

Article outline

1. Introduction
    — Basic definitions
    — Quality metrics
2. Move, smoothe, evaluate
    — Rolling window estimations
    — Exponential smoothing, Holt-Winters model
    — Time-series cross validation, parameters selection
3. Econometric approach
    — Stationarity, unit root
    — Getting rid of non-stationarity
    — SARIMA intuition and model building
4. Linear (and not quite) models on time series
    — Feature extraction
    — Linear models, feature importance
    — Regularization, feature selection
    — XGBoost

The following content is better viewed and reproduced as a Jupyter-notebook

In my day to day job I encounter time series-connected tasks almost every day. The most frequent question is — what will happen with our metrics in the next day/week/month/etc. — how many players will install the app, how much time will they spend online, how many actions users will do, and so forth. We can approach prediction task using different methods, depending on the required quality of the prediction, length of the forecasted period, and, of course, time we have to choose features and tune parameters to achieve desired results.

Introduction

Small definition of time series:
Time series — is a series of data points indexed (or listed or graphed) in time order.

Therefore data is organized around relatively deterministic timestamps, and therefore, compared to random samples, may contain additional information that we will try to extract.

Let’s import some libraries. First and foremost we will need statsmodels library that has tons of statistical modeling functions, including time series. For R afficionados (that had to move to python) statsmodels will definitely look familiar as it supports model definitions like ‘Wage ~ Age + Education’.

As an example let’s use some real mobile game data on hourly ads watched by players and daily in-game currency spent:

Forecast quality metrics

Before actually forecasting, let’s understand how to measure the quality of predictions and have a look at the most common and widely used metrics

• R squared, coefficient of determination (in econometrics it can be interpreted as a percentage of variance explained by the model), (-inf, 1] sklearn.metrics.r2_score
• Mean Absolute Error, it is an interpretable metric because it has the same unit of measurement as the initial series, [0, +inf)
  sklearn.metrics.mean_absolute_error
• Median Absolute Error, again an interpretable metric, particularly interesting because it is robust to outliers, [0, +inf)
  sklearn.metrics.median_absolute_error
• Mean Squared Error, most commonly used, gives higher penalty to big mistakes and vise versa, [0, +inf)
  sklearn.metrics.mean_squared_error
• Mean Squared Logarithmic Error, practically the same as MSE but we initially take logarithm of the series, as a result we give attention to small mistakes as well, usually is used when data has exponential trends, [0, +inf)
  sklearn.metrics.mean_squared_log_error
• Mean Absolute Percentage Error, same as MAE but percentage, — very convenient when you want to explain the quality of the model to your management, [0, +inf), not implemented in sklearn

Excellent, now we know how to measure the quality of the forecasts, what metrics can we use and how to translate the results to the boss. Little thing is left — building the model.

Move, smoothe, evaluate

Let’s start with a naive hypothesis — “tomorrow will be the same as today”, but instead of a model like ŷ(t)=y(t−1) (which is actually a great baseline for any time series prediction problems and sometimes it’s impossible to beat it with any model) we’ll assume that the future value of the variable depends on the average n of its previous values and therefore we’ll use moving average.

Out: 116805.0

Unfortunately we can’t make this prediction long-term — to get one for the next step we need the previous value to be actually observed. But moving average has another use case — smoothing of the original time series to indicate trends. Pandas has an implementation available DataFrame.rolling(window).mean(). The wider the window - the smoother will be the trend. In the case of the very noisy data, which can be very often encountered in finance, this procedure can help to detect common patterns.

Smoothing by last 4 hours
plotMovingAverage(ads, 4)

Smoothing by last 12 hours
plotMovingAverage(ads, 12)

Smoothing by 24 hours — we get daily trend
plotMovingAverage(ads, 24)

As you can see, applying daily smoothing on hour data allowed us to clearly see the dynamics of ads watched. During the weekends the values are higher (weekends — time to play) and weekdays are generally lower.

We can also plot confidence intervals for our smoothed values
plotMovingAverage(ads, 4, plot_intervals=True)

And now let’s create a simple anomaly detection system with the help of the moving average. Unfortunately, in this particular series everything is more or less normal, so we’ll intentionally make one of the values abnormal in the dataframe ads_anomaly

Let’s see, if this simple method can catch the anomaly
plotMovingAverage(ads_anomaly, 4, plot_intervals=True, plot_anomalies=True)

Neat! What about the second series (with weekly smoothing)?
plotMovingAverage(currency, 7, plot_intervals=True, plot_anomalies=True)

Oh no! Here is the downside of our simple approach — it did not catch monthly seasonality in our data and marked almost all 30-day peaks as an anomaly. If you don’t want to have that many false alarms — it’s best to consider more complex models.

Weighted average is a simple modification of the moving average, inside of which observations have different weights summing up to one, usually more recent observations have greater weight.

Out: 98423.0

Exponential smoothing

And now let’s take a look at what happens if instead of weighting the last nn values of the time series we start weighting all available observations while exponentially decreasing weights as we move further back in historical data. There’s a formula of the simple exponential smoothing that will help us in that:

Here the model value is a weighted average between the current true value and the previous model values. The α weight is called a smoothing factor. It defines how quickly we will “forget” the last available true observation. The less α is the more influence previous model values have, and the smoother the series is.

Exponentiality is hiding in the recursivity of the function — we multiply each time (1−α) by the previous model value which, in its turn, also containes (1−α) and so forth until the very beginning.

Double exponential smoothing

Until now all we could get from our methods in the best case was just a single future point prediction (and also some nice smoothing), that’s cool but not enough, so let’s extend exponential smoothing so that we can predict two future points (of course, we also get some smoothing).

Series decomposition should help us — we obtain two components: intercept (also, level) ℓ and trend (also, slope) b. We’ve learnt to predict intercept (or expected series value) using previous methods, and now we will apply the same exponential smoothing to the trend, believing naively or perhaps not that the future direction of the time series changes depends on the previous weighted changes.

As a result we get a set of functions. The first one describes intercept, as before it depends on the current value of the series, and the second term is now split into previous values of the level and of the trend. The second function describes trend — it depends on the level changes at the current step and on the previous value of the trend. In this case β coefficient is a weight in the exponential smoothing. The final prediction is the sum of the model values of the intercept and trend.

Now we have to tune two parameters — α and β. The former is responsible for the series smoothing around trend, and the latter for the smoothing of the trend itself. The bigger the values, the more weight the latest observations will have and the less smoothed the model series will be. Combinations of the parameters may produce really weird results, especially if set manually. We’ll look into choosing parameters automatically in a bit, immediately after triple exponential smoothing.

Triple exponential smoothing a.k.a. Holt-Winters

Hooray! We’ve successfully reached our next variant of exponential smoothing, this time triple.

The idea of this method is that we add another, third component — seasonality. This means we should’t use the method if our time series do not have seasonality, which is not the case in our example. Seasonal component in the model will explain repeated variations around intercept and trend, and it will be described by the length of the season, in other words by the period after which variations repeat. For each observation in the season there’s a separate component, for example, if the length of the season is 7 (weekly seasonality), we will have 7 seasonal components, one for each day of the week.

Now we get a new system:

Intercept now depends on the current value of the series minus corresponding seasonal component, trend stays unchanged, and the seasonal component depends on the current value of the series minus intercept and on the previous value of the component. Please take into account that the component is smoothed through all the available seasons, for example, if we have a Monday component then it will only be averaged with other Mondays. You can read more on how averaging works and how initial approximation of the trend and seasonal components is done here. Now that we have seasonal component we can predict not one and not even two but arbitrary mm future steps which is very encouraging.

Below is the code for a triple exponential smoothing model, also known by the last names of its creators — Charles Holt and his student Peter Winters. Additionally Brutlag method was included into the model to build confidence intervals:

where T is the length of the season, d is the predicted deviation, and the other parameters were taken from the triple exponential smoothing. You can read more about the method and its applicability to anomalies detection in time series here.

Time series cross validation

Before we start building model let’s talk first about how to estimate model parameters automatically.

There’s nothing unusual here, as always we have to choose a loss function suitable for the task, that will tell us how close the model approximates data. Then using cross-validation we will evaluate our chosen loss function for given model parameters, calculate gradient, adjust model parameters and so forth, bravely descending to the global minimum of error.

The question is how to do cross-validation on time series, because, you know, time series do have time structure and one just can’t randomly mix values in a fold without preserving this structure, otherwise all time dependencies between observations will be lost. That’s why we will have to use a bit more tricky approach to optimization of the model parameters, I don’t know if there’s an official name to it but on CrossValidated, where one can find all the answers but the Answer to the Ultimate Question of Life, the Universe, and Everything, “cross-validation on a rolling basis” was proposed as a name.

The idea is rather simple — we train our model on a small segment of the time series, from the beginning until some t, make predictions for the next t+n steps and calculate an error. Then we expand our training sample until t+n value and make predictions from t+n until t+2∗n, and we continue moving our test segment of the time series until we hit the last available observation. As a result we have as many folds as many n will fit between the initial training sample and the last observation.

Now, knowing how to set cross-validation, we will find optimal parameters for the Holt-Winters model, recall that we have daily seasonality in ads, hence the slen=24 parameter

In the Holt-Winters model, as well as in the other models of exponential smoothing, there’s a constraint on how big smoothing parameters could be, each of them is in the range from 0 to 1, therefore to minimize loss function we have to choose an algorithm that supports constraints on model parameters, in our case — Truncated Newton conjugate gradient.

0.11652680227350454 0.002677697431105852 0.05820973606789237
CPU times: user 1.96 s, sys: 17.3 ms, total: 1.98 s
Wall time: 2 s

Chart rendering code

Judging by the chart, our model was able to successfully approximate the initial time series, catching daily seasonality, overall downwards trend and even some anomalies. If you take a look at the modeled deviation, you can clearly see that the model reacts quite sharply to the changes in the structure of the series but then quickly returns deviation to the normal values, “forgetting” the past. This feature of the model allows us to quickly build anomaly detection systems even for quite noisy series without spending too much time and money on preparing data and training the model.

We’ll apply the same algorithm for the second series which, as we know, has trend and 30-day seasonality

0.012841445048055122 0.04883371471892228 0.00943678056045777
CPU times: user 3.03 s, sys: 24.8 ms, total: 3.05 s
Wall time: 3.11 s

Looks quite adequate, model has caught both upwards trend and seasonal spikes and overall fits our values nicely

Econometric approach

Stationarity

Before we start modeling we should mention such an important property of time series as stationarity.

If the process is stationary that means it doesn’t change its statistical properties over time, namely mean and variance do not change over time (constancy of variance is also called homoscedasticity), also covariance function does not depend on the time (should only depend on the distance between observations). You can see this visually on the pictures from the post of Sean Abu:

• The red graph below is not stationary because the mean increases over time.

• We were unlucky with the variance, see the varying spread of values over time

• Finally, the covariance of the i th term and the (i + m) th term should not be a function of time. In the following graph, you will notice the spread becomes closer as the time increases. Hence, the covariance is not constant with time for the right chart.

So why stationarity is so important? Because it’s easy to make predictions on the stationary series as we assume that the future statistical properties will not be different from the currently observed. Most of the time series models in one way or the other model and predict those properties (mean or variance, for example), that’s why predictions would be wrong if the original series were not stationary. Unfortunately most of the time series we see outside of textbooks are non-stationary but we can (and should) change this.

So, to fight non-stationarity we have to know our enemy so to say. Let’s see how to detect it. To do that we will now take a look at the white noise and random walks and we will learn how to get from one to another for free, without registration and SMS.

White noise chart:

So the process generated by standard normal distribution is stationary and oscillates around 0 with with deviation of 1. Now based on this process we will generate a new one where each next value will depend on the previous one: x(t)=ρ*x(t−1)+e(t)

Chart rendering code

On the first chart you can see the same stationary white noise you’ve seen before. On the second one the value of ρρ increased to 0.6, as a result wider cycles appeared on the chart but overall it is still stationary. The third chart deviates even more from the 0 mean but still oscillates around it. Finally, the value of ρ equal to 1 gives us a random walk process — non-stationary time series.

This happens because after reaching the critical value the series x(t)=ρ*x(t−1)+e(t) does not return to its mean value. If we subtract x(t−1) from the left and the right side we will get x(t)−x(t−1)=(ρ−1)*x(t−1)+e(t), where the expression on the left is called the first difference. If ρ=1 then the first difference gives us stationary white noise e(t). This fact is the main idea of the Dickey-Fuller test for the stationarity of time series (presence of a unit root). If we can get stationary series from non-stationary using the first difference we call those series integrated of order 1. Null hypothesis of the test — time series is non-stationary, was rejected on the first three charts and was accepted on the last one. We’ve got to say that the first difference is not always enough to get stationary series as the process might be integrated of order d, d > 1 (and have multiple unit roots), in such cases the augmented Dickey-Fuller test is used that checks multiple lags at once.

We can fight non-stationarity using different approaches — various order differences, trend and seasonality removal, smoothing, also using transformations like Box-Cox or logarithmic.

Getting rid of non-stationarity and building SARIMA

Now let’s build an ARIMA model by walking through all the circles of hell stages of making series stationary.

Chart rendering code

Surprisingly, initial series are stationary, Dickey-Fuller test rejected null hypothesis that a unit root is present. Actually, it can be seen on the plot itself — we don’t have a visible trend, so mean is constant, variance is pretty much stable throughout the series. The only thing left is seasonality which we have to deal with before modelling. To do so let’s take “seasonal difference” which means a simple subtraction of series from itself with a lag that equals the seasonal period.

That’s better, visible seasonality is gone, however autocorrelation function still has too many significant lags. To remove them we’ll take first differences — subtraction of series from itself with lag 1

Perfect! Our series now look like something undescribable, oscillating around zero, Dickey-Fuller indicates that it’s stationary and the number of significant peaks in ACF has dropped. We can finally start modelling!

ARIMA-family Crash-Course

A few words about the model. Letter by letter we’ll build the full name — SARIMA(p,d,q)(P,D,Q,s), Seasonal Autoregression Moving Average model:

• AR(p) — autoregression model, i.e., regression of the time series onto itself. Basic assumption — current series values depend on its previous values with some lag (or several lags). The maximum lag in the model is referred to as p. To determine the initial p you need to have a look at PACF plot — find the biggest significant lag, after which most other lags are becoming not significant.
• MA(q) — moving average model. Without going into detail it models the error of the time series, again the assumption is — current error depends on the previous with some lag, which is referred to as q. Initial value can be found on ACF plot with the same logic.

Let’s have a small break and combine the first 4 letters:

AR(p) + MA(q) = ARMA(p,q)

What we have here is the Autoregressive–moving-average model! If the series is stationary, it can be approximated with those 4 letters. Shall we continue?

• I(d)— order of integration. It is simply the number of nonseasonal differences needed for making the series stationary. In our case it’s just 1, because we used first differences.

Adding this letter to four previous gives us ARIMA model which knows how to handle non-stationary data with the help of nonseasonal differences. Awesome, last letter left!

• S(s) — this letter is responsible for seasonality and equals the season period length of the series

After attaching the last letter we find out that instead of one additional parameter we get three in a row — (P,D,Q)

• P — order of autoregression for seasonal component of the model, again can be derived from PACF, but this time you need to look at the number of significant lags, which are the multiples of the season period length, for example, if the period equals 24 and looking at PACF we see 24-th and 48-th lags are significant, that means initial P should be 2.
• Q — same logic, but for the moving average model of the seasonal component, use ACF plot
• D — order of seasonal integration. Can be equal to 1 or 0, depending on whether seasonal differences were applied or not

Now, knowing how to set initial parameters, let’s have a look at the final plot once again and set the parameters:
tsplot(ads_diff[24+1:], lags=60)

• p is most probably 4, since it’s the last significant lag on PACF after which most others are becoming not significant.
• d just equals 1, because we had first differences
• q should be somewhere around 4 as well as seen on ACF
• P might be 2, since 24-th and 48-th lags are somewhat significant on PACF
• D again equals 1 — we performed seasonal differentiation
• Q is probably 1, 24-th lag on ACF is significant, while 48-th is not

Now we want to test various models and see which one is better

Out: 54

CPU times: user 41.3 s, sys: 2.76 s, total: 44.1 s
Wall time: 37.4 s

Let’s inspect the residuals of the model
tsplot(best_model.resid[24+1:], lags=60)

Well, it’s clear that the residuals are stationary, there are no apparent autocorrelations, let’s make predictions using our model

In the end we got quite adequate predictions, our model on average was wrong by 4.01%, which is very very good, but overall costs of preparing data, making series stationary and brute-force parameters selecting might not be worth this accuracy.

Linear (and not quite) models on time series

Small lyrical digression again. Often in my job I have to build models with the only principle guiding me known as fast, good, cheap. That means some of the models will never be “production ready” as they demand too much time for the data preparation (for example, SARIMA), or require frequent re-training on new data (again, SARIMA), or are difficult to tune (good example — SARIMA), so it’s very often much easier to select a couple of features from the existing time series and build a simple linear regression or, say, a random forest. Good and cheap.

Maybe this approach is not backed up by theory, breaks different assumptions (like, Gauss-Markov theorem, especially about the errors being uncorrelated), but it’s very useful in practice and quite frequently used in machine learning competitions.

Feature exctraction

Alright, model needs features and all we have is a 1-dimentional time series to work with. What features can we exctract?

Lags of time series, of course

Window statistics:

• Max/min value of series in a window
• Average/median value in a window
• Window variance
• etc.

Date and time features:

• Minute of an hour, hour of a day, day of the week, you get it
• Is this day a holiday? Maybe something special happened? Make it a boolean feature

Target encoding

Forecasts from other models (though we can lose the speed of prediction this way)

Let’s run through some of the methods and see what we can extract from our ads series

Lags of time series

Shifting the series n steps back we get a feature column where the current value of time series is aligned with its value at the time t−n. If we make a 1 lag shift and train a model on that feature, the model will be able to forecast 1 step ahead having observed current state of the series. Increasing the lag, say, up to 6 will allow the model to make predictions 6 steps ahead, however it will use data, observed 6 steps back. If something fundamentally changes the series during that unobserved period, the model will not catch the changes and will return forecasts with big error. So, during the initial lag selection one has to find a balance between the optimal prediction quality and the length of forecasting horizon.

Wonderful, we got ourselves a dataset here, why don’t we train a model?

Well, simple lags and linear regression gave us predictions that are not that far from SARIMA in quality. There are lot’s of unnecessary features, but we’ll do feature selection a bit later. Now let’s continue engineering!

We’ll add into our dataset hour, day of the week and boolean for the weekend. To do so we need to transform current dataframe index into datetime format and exctract hour and weekday out of it.

We can visualize the resulting features

Since now we have different scales of variables — thousands for lag features and tens for categorical, it’s reasonable to transform them into same scale to continue exploring feature importances and later — regularization.

Test error goes down a little bit and judging by the coefficients plot we can say that weekday and is_weekend are rather useful features

Target encoding

I’d like to add another variant of encoding categorical variables — by mean value. If it’s undesirable to explode dataset by using tons of dummy variables that can lead to the loss of information about the distance, and if they can’t be used as real values because of the conflicts like “0 hours < 23 hours”, then it’s possible to encode a variable with a little bit more interpretable values. Natural idea is to encode with the mean value of the target variable. In our example every day of the week and every hour of the day can be encoded by the corresponding average number of ads watched during that day or hour. It’s very important to make sure that the mean value is calculated over train set only (or over current cross-validation fold only), so that the model is not aware of the future.

Let’s have a look at hour averages

Finally, put all the transformations together in a single function

Here comes overfitting! Hour_average variable was so great on train dataset that the model decided to concentrate all its forces on it - as a result the quality of prediction dropped. This problem can be approached in a variety of ways, for example, we can calculate target encoding not for the whole train set, but for some window instead, that way encodings from the last observed window will probably describe current series state better. Or we can just drop it manually, since we're sure here it makes things only worse.

Regularization and feature selection

As we already know, not all features are equally healthy, some may lead to overfitting and should be removed. Besides manual inspecting we can apply regularization. Two most popular regression models with regularization are Ridge and Lasso regressions. They both add some more constrains to our loss function.

In case of Ridge regression — those constrains are the sum of squares of coefficients, multiplied by the regularization coefficient. I.e. the bigger coefficient feature has — the bigger our loss will be, hence we will try to optimize the model while keeping coefficients fairly low.

As a result of such regularization which has a proud name L2 we’ll have higher bias and lower variance, so the model will generalize better (at least that’s what we hope will happen).

Second model — Lasso regression, here we add to the loss function not squares but absolute values of the coefficients, as a result during the optimization process coefficients of unimportant features may become zeroes, so Lasso regression allows for automated feature selection. This regularization type is called L1.

First, make sure we have things to drop and data truly has highly correlated features

We can clearly see how coefficients are getting closer and closer to zero (thought never actually reach it) as their importance in the model drops

Lasso regression turned out to be more conservative and removed 23-rd lag from most important features (and also dropped 5 features completely) which only made the quality of prediction better.

Boosting

Why not try XGBoost now?

Here is the winner! The smallest error on the test set among all the models we’ve tried so far.

Yet this victory is decieving and it might not be the brightest idea to fit xgboost as soon as you get your hands over time series data. Generally tree-based models poorly handle trends in data, compared to linear models, so you have to detrend your series first or use some tricks to make the magic happen. Ideally — make the series stationary and then use XGBoost, for example, you can forecast trend separately with a linear model and then add predictions from xgboost to get final forecast.

Conclusion

We got acquainted with different time series analysis and prediction methods and approaches. Unfortunately, or maybe luckily, there’s no silver bullet to solve this kind of problems. Methods developed in the 60s of the last century (and some even in the beginning of the XIX century) are still popular along with the LSTM and RNN (not covered in this article). Partially this is related to the fact that the prediction task as any other data related task is creative in so many aspects and definitely requires research. In spite of the large number of formal quality metrics and approaches to parameters estimation, it’s often required to seek and try something different for each time series. Last but not least the balance between quality and cost is important. As a good example SARIMA model mentioned here not once or twice can produce spectacular results after due tuning but might require many hours of tambourine dancing time series manipulation, as in the same time simple linear regression model can be build in 10 minutes giving more or less comparable results.

My socials:

• Data Wondering blog [English]
• Data Wondering blog [Russian]

Useful resources

• Online textbook of the advanced statistical forecasting course of the Duke University — covers in details various smoothing techniques, linear and ARIMA models
• Comparison of ARIMA and Random Forest time series models for prediction of avian influenza H5N1 outbreaks — one of a few where random forest applicability to the tasks of time series forecasting is actively defended
• Time Series Analysis (TSA) in Python — Linear Models to GARCH ARIMA models family and their applicability to the task of modeling financial indicators (Brian Christopher)

Author: Dmitry Sergeyev. Translated and edited by Borys Zibrov, and Yuanyuan Pao.

Open Machine Learning Course. Topic 9. Part 1. Time series analysis in Python was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

This week, we’ll cover two reasons for Vowpal Wabbit’s exceptional training speed, namely, online learning and hashing trick, in both theory and practice. We will try it out with news, movie reviews, and StackOverflow questions.

Outline

1. Stochastic gradient descent and online learning

• 1.1. SGD
• 1.2. Online approach to learning

2. Categorical data processing: Label Encoding, One-Hot Encoding, Hashing trick

• 2.1. Label Encoding
• 2.2. One-Hot Encoding
• 2.3. Hashing trick

3. Vowpal Wabbit

• 3.1. News. Binary classification
• 3.2. News. Multiclass classification
• 3.3. IMDB reviews
• 3.4. Classifying gigabytes of StackOverflow questions

4. Assignment #8

5. Useful resources

The following content is better viewed and reproduced as a Jupyter-notebook or a Kaggle kernel.

1. Stochastic gradient descent and online learning

1.1. Stochastic gradient descent

Despite the fact that gradient descent is one of the first things learned in machine learning and optimization courses, it is one of its modifications, Stochastic Gradient Descent (SGD), that is hard to top.

Recall that the idea of gradient descent is to minimize some function by making small steps in the direction of the fastest decrease. This method was named due to the following fact from calculus: vector

of partial derivatives of the function f(x) = f(x1, …, xn) points to the direction of the fastest function growth. It means that, by moving in the opposite direction (antigradient), it is possible to decrease the function value with the fastest rate.

In addition to advertising beautiful landscapes, this picture depicts the idea of gradient descent. If you want to ride as fast as possible, you need to choose the path of steepest descent. Calculating antigradients can be seen as evaluating the slope at various spots.

Example

The paired regression problem can be solved with gradient descent. Let us predict one variable using another: height with weight. Assume that these variables are linearly dependent. We will use the SOCR dataset.

Here we have a vector x of dimension ℓ (weight of every person i.e. training sample) and y, a vector containing the height of every person in the dataset.

The task is the following: find weights w0 and w1 such that predicting height as yi = w0 + w_1 xi (where yi is i-th height value, xi is i-th weight value) minimizes the squared error (as well as mean squared error since 1/ℓ doesn’t make any difference):

We will use gradient descent, utilizing the partial derivatives of SE(w0,w1) over weights w0 and w1. An iterative training procedure is then defined by simple update formulas (we change model weights in small steps, proportional to a small constant η, towards the antigradient of the function SE(w0, w1)):

Computing the partial derivatives, we get the following:

This math works quite well as long as the amount of data is not large (we will not discuss issues with local minima, saddle points, choosing the learning rate, moments and other stuff – these topics are covered very thoroughly in the Numeric Computation chapter in “Deep Learning”). There is an issue with batch gradient descent — the gradient evaluation requires the summation of a number of values for every object from the training set. In other words, the algorithm requires a lot of iterations, and every iteration recomputes weights with formula which contains a sum over the whole training set. What happens when we have billions of training samples?

Hence the motivation for stochastic gradient descent! Simply put, we throw away the summation sign and update the weights only over single training samples (or a small number of them). In our case, we have the following:

With this approach, there is no guarantee that we will move in best possible direction at every iteration. Therefore, we may need many more iterations, but we get much faster weight updates.

Andrew Ng has a good illustration of this in his machine learning course. Let’s take a look.

These are the contour plots for some function, and we want to find the global minimum of this function. The red curve shows weight changes (in this picture, θ0 and θ1 correspond to our w0 and w1). According to the properties of a gradient, the direction of change at every point is orthogonal to contour plots. With stochastic gradient descent, weights are changing in a less predictable manner, and it even may seem that some steps are wrong by leading away from minima; however, both procedures converge to the same solution.

1.2. Online approach to learning

Stochastic gradient descent gives us practical guidance for training both classifiers and regressors with large amounts of data up to hundreds of GBs (depending on computational resources).

Considering the case of paired regression, we can store the training data set (X,y) in HDD without loading it into RAM (where it simply won’t fit), read objects one by one, and update the weights of our model:

After working through the whole training dataset, our loss function (for example, quadratic squared root error in regression or logistic loss in classification) will decrease, but it usually takes dozens of passes over the training set to make the loss small enough.

This approach to learning is called online learning, and this name emerged even before machine learning MOOC-s turned mainstream.

We did not discuss many specifics about SGD here. If you want dive into theory, I highly recommend “Convex Optimization” by Stephen Boyd. Now, we will introduce the Vowpal Wabbit library, which is good for training simple models with huge data sets thanks to stochastic optimization and another trick, feature hashing.

In scikit-learn, classifiers and regressors trained with SGD are named SGDClassifier and SGDRegressor in sklearn.linear_model. These are nice implementations of SGD, but we'll focus on VW since it is more performant than sklearn's SGD models in many aspects.

2. Categorical feature processing: Label Encoding, One-Hot Encoding, and Hashing trick

Many classification and regression algorithms operate in Euclidean or metric space, implying that data is represented with vectors of real numbers. However, in real data, we often have categorical features with discrete values such as yes/no or January/February/…/December. We will see how to process this kind of data, particularly with linear models, and how to deal with many categorial features even when they have many unique values.

Let’s explore the UCI bank marketing dataset where most of features are categorial. Actually, more features are present in this dataset, explore this Jupyter notebook.

We can see that most of features are not represented by numbers. This poses a problem because we cannot use most machine learning methods (at least those implemented in scikit-learn) out-of-the-box.

Let’s dive into the “education” feature.

df['education'].value_counts().plot.barh();

The most straightforward solution is to map each value of this feature into a unique number. For example, we can map university.degree to 0, basic.9y to 1, and so on. You can use sklearn.preprocessing.LabelEncoder to perform this mapping.

label_encoder = LabelEncoder()

The fit method of this class finds all unique values and builds the actual mapping between categories and numbers, and the transform method converts the categories into numbers. After fit is executed, label_encoder will have the classes_ attribute with all unique values of the feature. Let us count them to make sure the transformation was correct.

# Output:
 
{0: 'basic.4y', 1: 'basic.6y', 2: 'basic.9y', 3: 'high.school', 4: 'illiterate', 5: 'professional.course', 6: 'university.degree', 7: 'unknown'}

df['education'] = mapped_education df.head()

Let’s apply the transformation to other columns of type object.

The main issue with this approach is that we have now introduced some relative ordering where it might not exist.

For example, we implicitly introduced algebra over the values of the job feature where we can now substract the job of client #2 from the job of client #1 :

df.loc[1].job - df.loc[2].job
# -1.0

Does this operation make any sense? Not really. Let’s try to train logisitic regression with this feature transformation.

We can see that logistic regression never predicts class 1. In order to use linear models with categorial features, we will use a different approach: One-Hot Encoding.

2.2. One-Hot Encoding

Suppose that some feature can have one of 10 unique values. One-hot encoding creates 10 new features corresponding to these unique values, all of them except one are zeros.

This idea is implemented in the OneHotEncoder class from sklearn.preprocessing. By default OneHotEncoder transforms data into a sparse matrix to save memory space because most of the values are zeroes and because we do not want to take up more RAM. However, in this particular example, we do not encounter such problems, so we are going to use a "dense" matrix representation.

We have 53 columns that correspond to the number of unique values of categorical features in our data set. When transformed with One-Hot Encoding, this data can be used with linear models:

print(logistic_regression_accuracy_on(encoded_categorical_columns,  labels))

2.3. Hashing trick

Real data can be volatile, meaning we cannot guarantee that new values of categorial features will not occur. This issue hampers using a trained model on new data. Besides that, LabelEncoder requires preliminary analysis of the whole dataset and storage of constructed mappings in memory, which makes it difficult to work with large datasets.

There is a simple approach to vectorization of categorial data based on hashing and is known as, not-so-surprisingly, the hashing trick.

Hash functions can help us find unique codes for different feature values, for example:

# university.degree -> -6241459093488141593
# high.school -> 7728198035707179500
# illiterate -> -7360093633803373451

We will not use negative values or values of high magnitude, so we restrict the range of values for the hash function:

# university.degree -> 7
# high.school -> 0
# illiterate -> 24

Imagine that our data set contains a single (i.e. not married) student, who received a call on Monday. His feature vectors will be created similarly as in the case of One-Hot Encoding but in the space with fixed range for all features:

# job=student -> 20
# marital=single -> 23
# day_of_week=mon -> 9

We want to point out that we hash not only feature values but also pairs of feature name + feature value. It is important to do this so that we can distinguish the same values of different features.

Is it possible to have a collision when using hash codes? Sure, it is possible, but it is a rare case with large enough hashing spaces. Even if collision occurs, regression or classification metrics will not suffer much. In this case, hash collisions work as a form of regularization.

You may be saying “WTF?”; hashing seems counterintuitive. This is true, but these heuristics sometimes are, in fact, the only plausible approach to work with categorial data. Moreover, this technique has proven to just work. As you work more with data, you may see this for yourself.

3. Vowpal Wabbit

Vowpal Wabbit (VW) is one of the most widespread machine learning libraries used in industry. It is prominent for its training speed and support of many training modes, especially for online learning with big and high-dimentional data. This is one of the major merits of the library. Also, with the hashing trick implemented, Vowpal Wabbit is a perfect choice for working with text data.

Shell is the main interface for VW. Try this command:

!vw — help

Vowpal Wabbit reads data from files or from standard input stream (stdin) with the following format:

[Label] [Importance] [Tag]|Namespace Features |Namespace Features ... |Namespace Features

Namespace=String[:Value]

Features=(String[:Value] )*

here [] denotes non-mandatory elements, and (…)* means multiple inputs allowed.

• Label is a number. In the case of classification, it is usually 1 and -1; for regression, it is a real float value
• Importance is a number. It denotes the sample weight during training. Setting this helps when working with imbalanced data.
• Tag is a string without spaces. It is the “name” of the sample that VW saves upon prediction. In order to separate Tag from Importance, it is better to start Tag with the ‘ character.
• Namespace is for creating different feature spaces.
• Features are object features inside a given Namespace. Features have weight 1.0 by default, but it can be changed, for example feature:0.1.

The following string matches the VW format:

1 1.0 |Subject WHAT car is this |Organization University of Maryland:0.5 College Park

Let’s check the format by running VW with this training sample:

! echo '1 1.0 |Subject WHAT car is this |Organization University of Maryland:0.5 College Park' | vw

# Output:

Num weight bits = 18
learning rate = 0.5
initial_t = 0
power_t = 0.5
using no cache
Reading datafile = 
num sources = 1
average  since         example        example  current  current  current
loss     last          counter         weight    label  predict features
1.000000 1.000000            1            1.0   1.0000   0.0000       10

finished run
number of examples per pass = 1
passes used = 1
weighted example sum = 1.000000
weighted label sum = 1.000000
average loss = 1.000000
best constant = 1.000000
best constant's loss = 0.000000
total feature number = 10

VW is a wonderful tool for working with text data. We’ll illustrate it with the 20newsgroups dataset, which contains letters from 20 different newsletters.

3.1. News. Binary classification.

# load data with sklearn's function
newsgroups = fetch_20newsgroups(PATH_TO_ALL_DATA)
newsgroups['target_names']

The 20 topics are ['alt.atheism', 'comp.graphics', 'comp.os.ms-windows.misc', 'comp.sys.ibm.pc.hardware', 'comp.sys.mac.hardware', 'comp.windows.x', 'misc.forsale', 'rec.autos', 'rec.motorcycles', 'rec.sport.baseball', 'rec.sport.hockey', 'sci.crypt', 'sci.electronics', 'sci.med', 'sci.space', 'soc.religion.christian', 'talk.politics.guns', 'talk.politics.mideast', 'talk.politics.misc', 'talk.religion.misc']

Lets look at the first document in this collection:

# Output:

-----
rec.autos
-----
From: lerxst@wam.umd.edu (where's my thing)
Subject: WHAT car is this!?
Nntp-Posting-Host: rac3.wam.umd.edu
Organization: University of Maryland, College Park
Lines: 15

 I was wondering if anyone out there could enlighten me on this car I saw
the other day. It was a 2-door sports car, looked to be from the late 60s/
early 70s. It was called a Bricklin. The doors were really small. In addition,
the front bumper was separate from the rest of the body. This is 
all I know. If anyone can tellme a model name, engine specs, years
of production, where this car is made, history, or whatever info you
have on this funky looking car, please e-mail.

Thanks,
- IL
   ---- brought to you by your neighborhood Lerxst ----
----

Now we convert the data into something Vowpal Wabbit can understand. We will throw away words shorter than 3 symbols. Here, we will skip some important NLP stages such as stemming and lemmatization; however, we will later see that VW solves the problem even without these steps.

# Output:

'1 |text from lerxst wam umd edu where thing subject what car this nntp posting host rac3 wam umd edu organization university maryland college park lines was wondering anyone out there could enlighten this car saw the other day was door sports car looked from the late 60s early 70s was called bricklin the doors were really small addition the front bumper was separate from the rest the body this all know anyone can tellme model name engine specs years production where this car made history whatever info you have this funky looking car please mail thanks brought you your neighborhood lerxst\n'

We split the dataset into train and test and write these into separate files. We will consider a document as positive if it corresponds to rec.autos. Thus, we are constructing a model which distinguishes articles about cars from other topics:

Now, we pass the created training file to Vowpal Wabbit. We solve the classification problem with a hinge loss function (linear SVM). The trained model will be saved in the 20news_model.vw file:

!vw -d $PATH_TO_ALL_DATA/20news_train.vw \ --loss_function hinge -f $PATH_TO_ALL_DATA/20news_model.vw

# Output:

final_regressor = ../../data//20news_model.vw
Num weight bits = 18
learning rate = 0.5
initial_t = 0
power_t = 0.5
using no cache
Reading datafile = ../../data//20news_train.vw
num sources = 1
average  since         example        example  current  current  current
loss     last          counter         weight    label  predict features
1.000000 1.000000            1            1.0  -1.0000   0.0000      157
0.911276 0.822551            2            2.0  -1.0000  -0.1774      159
0.605793 0.300311            4            4.0  -1.0000  -0.3994       92
0.419594 0.233394            8            8.0  -1.0000  -0.8167      129
0.313998 0.208402           16           16.0  -1.0000  -0.6509      108
0.196014 0.078029           32           32.0  -1.0000  -1.0000      115
0.183158 0.170302           64           64.0  -1.0000  -0.7072      114
0.261046 0.338935          128          128.0   1.0000  -0.7900      110
0.262910 0.264774          256          256.0  -1.0000  -0.6425       44
0.216663 0.170415          512          512.0  -1.0000  -1.0000      160
0.176710 0.136757         1024         1024.0  -1.0000  -1.0000      194
0.134541 0.092371         2048         2048.0  -1.0000  -1.0000      438
0.104403 0.074266         4096         4096.0  -1.0000  -1.0000      644
0.081329 0.058255         8192         8192.0  -1.0000  -1.0000      174

finished run
number of examples per pass = 8485
passes used = 1
weighted example sum = 8485.000000
weighted label sum = -7555.000000
average loss = 0.079837
best constant = -1.000000
best constant's loss = 0.109605
total feature number = 2048932

VW prints a lot of interesting info while training (one can suppress it with the --quiet parameter). You can see documentation of the diagnostic output on GitHub. Note how average loss drops while training. For loss computation, VW uses samples it has never seen before, so this measure is usually accurate. Now, we apply our trained model to the test set, saving predictions into a file with the -p flag:

!vw -i $PATH_TO_ALL_DATA/20news_model.vw -t -d $PATH_TO_ALL_DATA/20news_test.vw \ -p $PATH_TO_ALL_DATA/20news_test_predictions.txt

Now we load our predictions, compute AUC, and plot the ROC curve:

The AUC value we get shows that we have achieved high classification quality.

3.2. News. Multiclass classification

We will use the same news dataset, but, this time, we will solve a multiclass classification problem. Vowpal Wabbit is a little picky – it wants labels starting from 1 till K, where K – is the number of classes in the classification task (20 in our case). So we will use LabelEncoder and add 1 afterwards (recall that LabelEncoder maps labels into range from 0 to K-1).

The data is the same, but we have changed the labels, train_labels_mult and test_labels_mult, into label vectors from 1 to 20.

We train Vowpal Wabbit in multiclass classification mode, passing the oaa parameter("one against all") with the number of classes. Also, let's see what parameters our model quality is dependent on (more info can be found in the official Vowpal Wabbit tutorial):

• learning rate (-l, 0.5 default) — rate of weight change on every step
• learning rate decay ( — power_t, 0.5 default) — it is proven in practice, that, if the learning rate drops with the number of steps in stochastic gradient descent, we approach the minimum loss better
• loss function ( — loss_function) — the entire training algorithm depends on it. See docs for loss functions
• Regularization (-l1) — note that VW calculates regularization for every object. That is why we usually set regularization values to about 1e-20.

Additionally, we can try automatic Vowpal Wabbit parameter tuning with Hyperopt.

!vw — oaa 20 $PATH_TO_ALL_DATA/20news_train_mult.vw -f $PATH_TO_ALL_DATA/20news_model_mult.vw \
 — loss_function=hinge

!vw -i $PATH_TO_ALL_DATA/20news_model_mult.vw -t -d $PATH_TO_ALL_DATA/20news_test_mult.vw \
-p $PATH_TO_ALL_DATA/20news_test_predictions_mult.txt

Here is how often the model misclassifies atheism with other topics:

# Output:

rec.autos 1
rec.sport.baseball 1
sci.med 1
soc.religion.christian 3
talk.religion.misc 5

3.3. IMDB movie reviews

In this part we will do binary classification of IMDB (International Movie DataBase) movie reviews. We will see how fast Vowpal Wabbit performs.

Using the load_files function from sklearn.datasets, we load the movie reviews here. If you want to reproduce the results, please download the archive, unzip it, and set the path to imdb_reviews (it already contains train and test subdirectories). Unpacking can take several minutes as there are 100k files. Both train and test sets hold 12.5k good and bad movie reviews. First, we split the texts and labels.

Take a look at examples of reviews and their corresponding labels.

text_train[0]

# Output:
b"Zero Day leads you to think, even re-think why two boys/young men would do what they did - commit mutual suicide via slaughtering their classmates. It captures what must be beyond a bizarre mode of being for two humans who have decided to withdraw from common civility in order to define their own/mutual world via coupled destruction.<br /><br />It is not a perfect movie but given what money/time the filmmaker and actors had - it is a remarkable product. In terms of explaining the motives and actions of the two young suicide/murderers it is better than 'Elephant' - in terms of being a film that gets under our 'rationalistic' skin it is a far, far better film than almost anything you are likely to see. <br /><br />Flawed but honest with a terrible honesty."

y_train[0] # good review

# Output:
1

text_train[1]

# Output:
b'Words can\'t describe how bad this movie is. I can\'t explain it by writing only. You have too see it for yourself to get at grip of how horrible a movie really can be. Not that I recommend you to do that. There are so many clich\xc3\xa9s, mistakes (and all other negative things you can imagine) here that will just make you cry. To start with the technical first, there are a LOT of mistakes regarding the airplane. I won\'t list them here, but just mention the coloring of the plane. They didn\'t even manage to show an airliner in the colors of a fictional airline, but instead used a 747 painted in the original Boeing livery. Very bad. The plot is stupid and has been done many times before, only much, much better. There are so many ridiculous moments here that i lost count of it really early. Also, I was on the bad guys\' side all the time in the movie, because the good guys were so stupid. "Executive Decision" should without a doubt be you\'re choice over this one, even the "Turbulence"-movies are better. In fact, every other movie in the world is better than this one.'

y_train[1] # bad review

# Output:
0

Now, we prepare training (movie_reviews_train.vw), validation (movie_reviews_valid.vw), and test (movie_reviews_test.vw) sets for Vowpal Wabbit. We will use 70% for training, 30% for the hold-out set.

Now, launch Vowpal Wabbit with the following arguments:

• -d, path to training set (corresponding .vw file)
• — loss_function — hinge (feel free to experiment here)
• -f — path to the output file (which can also be in the .vw format)

!vw -d $PATH_TO_ALL_DATA/movie_reviews_train.vw --loss_function hinge -f $PATH_TO_ALL_DATA/movie_reviews_model.vw --quiet

Next, make the hold-out prediction with the following VW arguments:

• -i –path to the trained model (.vw file)
• -t -d — path to hold-out set (.vw file)
• -p — path to a txt-file where the predictions will be stored

!vw -i $PATH_TO_ALL_DATA/movie_reviews_model.vw -t \ 
-d $PATH_TO_ALL_DATA/movie_reviews_valid.vw -p $PATH_TO_ALL_DATA/movie_valid_pred.txt --quiet

Read the predictions from the text file and estimate the accuracy and ROC AUC. Note that VW prints probability estimates of the +1 class. These estimates are distributed from -1 to 1, so we can convert these into binary answers, assuming that positive values belong to class 1.

Accuracy: 0.885
AUC: 0.942

Again, do the same for the test set.

!vw -i $PATH_TO_ALL_DATA/movie_reviews_model.vw -t \ -d $PATH_TO_ALL_DATA/movie_reviews_test.vw \ -p $PATH_TO_ALL_DATA/movie_test_pred.txt --quiet

Accuracy: 0.88
AUC: 0.94

Let’s try to achieve a higher accuracy by incorporating bigrams.

!vw -d $PATH_TO_ALL_DATA/movie_reviews_train.vw \ --loss_function hinge --ngram 2 -f $PATH_TO_ALL_DATA/movie_reviews_model2.vw --quiet

!vw -i$PATH_TO_ALL_DATA/movie_reviews_model2.vw -t -d $PATH_TO_ALL_DATA/movie_reviews_valid.vw \ -p $PATH_TO_ALL_DATA/movie_valid_pred2.txt --quiet

!vw -i $PATH_TO_ALL_DATA/movie_reviews_model2.vw -t -d $PATH_TO_ALL_DATA/movie_reviews_test.vw \ -p $PATH_TO_ALL_DATA/movie_test_pred2.txt --quiet

Adding bigrams really helped to improve our model!

# Validation 
Accuracy: 0.894
AUC: 0.954

# Test
Accuracy: 0.888
AUC: 0.952

3.4. Classifying gigabytes of StackOverflow questions

Now, let’s see Vowpal Wabbit work on large datasets. There is a 10GB dataset of StackOverflow questions here. The original dataset is comprised of 10 million questions; each question can have several tags. The data is quite clean, so don’t call it “Big Data”, even in a pub. :)

We chose only 10 tags: ‘javascript’, ‘java’, ‘python’, ‘ruby’, ‘php’, ‘c++’, ‘c#’, ‘go’, ‘scala’ and ‘swift’. Let’s solve the 10-class classification problem: we want to predict a tag corresponding to one of these 10 popular programming languages given only the text of the question. Your task in assignment is to

# change the path to data 
PATH_TO_STACKOVERFLOW_DATA = 'stackoverflow_10mln'

# Print the first 3 lines from a sample of the dataset.

!head -3 $PATH_TO_STACKOVERFLOW_DATA/stackoverflow_train.vw

# Output:

1 | i ve got some code in window scroll that checks if an element is visible then triggers another function however only the first section of code is firing both bits of code work in and of themselves if i swap their order whichever is on top fires correctly my code is as follows fn isonscreen function use strict var win window viewport top win scrolltop left win scrollleft bounds this offset viewport right viewport left + win width viewport bottom viewport top + win height bounds right bounds left + this outerwidth bounds bottom bounds top + this outerheight return viewport right lt bounds left viewport left gt bounds right viewport bottom lt bounds top viewport top gt bounds bottom window scroll function use strict var load_more_results ajax load_more_results isonscreen if load_more_results true loadmoreresults var load_more_staff ajax load_more_staff isonscreen if load_more_staff true loadmorestaff what am i doing wrong can you only fire one event from window scroll i assume not 
4 | redefining some constant in ruby ex foo bar generates the warning already initialized constant i m trying to write a sort of reallyconstants module where this code should have this behaviour reallyconstants define_constant foo bar gt sets the constant reallyconstants foo to bar reallyconstants foo gt bar reallyconstants foo foobar gt this should raise an exception that is constant redefinition should generate an exception is that possible 1 | in my form panel i added a checkbox setting stateful true stateid loginpanelremeberme then before sending form i save state calling this savestate on the panel all other componenets save their state and whe i reload the page they recall the previous state but checkbox alway start in unchecked state is there any way to force save value

After selecting our 10 tags, we have a 4.7G dataset that we divide into train and test.

!du -hs $PATH_TO_STACKOVERFLOW_DATA/stackoverflow_*.vw

# Output:

4,7G	
stackoverflow_10mln.vw 
1,6G	
stackoverflow_test.vw 
3,1G	
stackoverflow_train.vw

We will process the training set (3.1 GiB) with Vowpal Wabbit and the following arguments:

• -oaa 10 — for multiclass classification with 10 classes
• -d — path to data
• -f — path to output file of the trained model
• -b 28 — we will use 28 bits for hashing, resulting in a $2^{28}$-sized feature space
• fix random seed for reproducibility

%%time !vw --oaa 10 -d $PATH_TO_STACKOVERFLOW_DATA/stackoverflow_train.vw \ -f vw_model1_10mln.vw -b 28 --random_seed 17 --quiet

# Output:

CPU times: user 559 ms, sys: 171 ms, total: 730 ms Wall time: 38.5 s

%%time !vw -t -i vw_model1_10mln.vw -d $PATH_TO_STACKOVERFLOW_DATA/stackoverflow_test.vw \ -p vw_test_pred.csv --random_seed 17 --quiet

# Output:

CPU times: user 322 ms, sys: 97 ms, total: 420 ms Wall time: 22.8 s

# Output:

0.91728604842865913

The model has trained itself and made predictions in less than a minute (check it, these results are reported for a MacBook Pro, mid 2015, 2.2 GHz Intel Core i7, 16GB RAM). Its accuracy is almost 92%. All of this without a Hadoop cluster! :) Impressive, isn’t it?

4. Assignment #8

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

5. Useful resources

• Official VW documentation on Github
• “Numeric Computation” Chapter of the Deep Learning book
• “Convex Optimization” by Stephen Boyd
• “Command-line Tools can be 235x Faster than your Hadoop Cluster” post
• Benchmarking various ML algorithms on Criteo 1TB dataset on GitHub
• VW on FastML.com

Author: Yury Kashnitsky. Translated and edited by Serge Oreshkov, and Yuanyuan Pao.

Open Machine Learning Course. Topic 8. Vowpal Wabbit: Fast Learning with Gigabytes of Data was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

Welcome to the seventh part of our Open Machine Learning Course!

In this lesson, we will work with unsupervised learning methods such as Principal Component Analysis (PCA) and clustering. You will learn why and how we can reduce the dimensionality of the original data and what the main approaches are for grouping similar data points.

Article outline

1. Introduction
2. Principal Component Analysis

• Intuition, theories, and application issues
• Application examples

3. Cluster analysis

• K-means
• Affinity Propagation
• Spectral clustering
• Agglomerative clustering
• Accuracy metrics

4. Assignment #7

5. Useful links

1. Introduction

The main feature of unsupervised learning algorithms, when compared to classification and regression methods, is that input data are unlabeled (i.e. no labels or classes given) and that the algorithm learns the structure of the data without any assistance. This creates two main differences. First, it allows us to process large amounts of data because the data does not need to be manually labeled. Second, it is difficult to evaluate the quality of an unsupervised algorithm due to the absence of an explicit goodness metric as used in supervised learning.

One of the most common tasks in unsupervised learning is dimensionality reduction. On one hand, dimensionality reduction may help with data visualization (e.g. t-SNA method) while, on the other hand, it may help deal with the multicollinearity of your data and prepare the data for a supervised learning method (e.g. decision trees).

2. Principal Component Analysis (PCA)

Intuition, theories, and application issues

Principal Component Analysis is one of the easiest, most intuitive, and most frequently used methods for dimensionality reduction, projecting data onto its orthogonal feature subspace.

More generally speaking, all observations can be considered as an ellipsoid in a subspace of an initial feature space, and the new basis set in this subspace is aligned with the ellipsoid axes. This assumption lets us remove highly correlated features since basis set vectors are orthogonal. In the general case, the resulting ellipsoid dimensionality matches the initial space dimensionality, but the assumption that our data lies in a subspace with a smaller dimension allows us to cut off the “excessive” space with the new projection (subspace). We accomplish this in a ‘greedy’ fashion, sequentially selecting each of the ellipsoid axes by identifying where the dispersion is maximal.

“To deal with hyper-planes in a 14 dimensional space, visualize a 3D space and say ‘fourteen’ very loudly. Everyone does it.” — Geoffrey Hinton

Let’s take a look at the mathematical formulation of this process:

In order to decrease the dimensionality of our data from n to k with k ≤ n, we sort our list of axes in order of decreasing dispersion and take the top-k of them.

We begin by computing the dispersion and the covariance of the initial features. This is usually done with the covariance matrix. According to the covariance definition, the covariance of two features is computed as follows:

where µ is the expected value of the i-th feature. It is worth noting that the covariance is symmetric, and the covariance of a vector with itself is equal to its dispersion.

Therefore the covariance matrix is symmetric with the dispersion of the corresponding features on the diagonal. Non-diagonal values are the covariances of the corresponding pair of features. In terms of matrices where X is the matrix of observations, the covariance matrix is as follows:

Quick recap: matrices, as linear operators, have eigenvalues and eigenvectors. They are very convenient because they describe parts of our space that do not rotate and only stretch when we apply linear operators on them; eigenvectors remain in the same direction but are stretched by a corresponding eigenvalue. Formally, a matrix M with eigenvector w and eigenvalue λ satisfy this equation:

The covariance matrix for a sample X can be written as a product of a transposed matrix X and X itself. According to the Rayleigh quotient, the maximum variation of our sample lies along the eigenvector of this matrix and is consistent with the maximum eigenvalue. Therefore, the principal components we aim to retain from the data are just the eigenvectors corresponding to the top-k largest eigenvalues of the matrix.

The next steps are easier to digest. We multiply the matrix of our data X by these components to get the projection of our data onto the orthogonal basis of the chosen components. If the number of components was smaller than the initial space dimensionality, remember that we will lose some information upon applying this transformation.

Examples

Fisher’s iris dataset

Let’s start by uploading all of the essential modules and try out the iris example from the scikit-learn documentation.

Now let’s see how PCA will improve the results of a simple model that is not able to correctly fit all of the training data:

Accuracy: 0.88889

Let’s try this again, but, this time, let’s reduce the dimensionality to 2 dimensions:

Accuracy: 0.91111

The accuracy did not increase significantly in this case, but, with other datasets with a high number of dimensions, PCA can drastically improve the accuracy of decision trees and other ensemble methods.

Now let’s check out the percent of variance that can be explained by each of the selected components.

1 component: 92.46% of initial variance 0.362 x sepal length (cm) + -0.082 x sepal width (cm) + 0.857 x petal length (cm) + 0.359 x petal width (cm) 
2 component: 5.3% of initial variance 0.657 x sepal length (cm) + 0.730 x sepal width (cm) + -0.176 x petal length (cm) + -0.075 x petal width (cm)

Handwritten numbers dataset

Let’s look at the handwritten numbers dataset that we used before in the 3rd lesson.

Let’s start by visualizing our data. Fetch the first 10 numbers. The numbers are represented by 8 x 8 matrixes with the color intensity for each pixel. Every matrix is flattened into a vector of 64 numbers, so we get the feature version of the data.

Our data has 64 dimensions, but we are going to reduce it to only 2 and see that, even with just 2 dimensions, we can clearly see that digits separate into clusters.

Projecting 64-dimensional data to 2D

Indeed, with t-SNE, the picture looks better since PCA has a linear constraint while t-SNE does not. However, even with such a small dataset, the t-SNE algorithm takes significantly more time to complete than PCA.

CPU times: user 29.3 s, sys: 2.97 s, total: 32.3 s Wall time: 32.3 s

In practice, we would choose the number of principal components such that we can explain 90% of the initial data dispersion (via the explained_variance_ratio). Here, that means retaining 21 principal components; therefore, we reduce the dimensionality from 64 features to 21.

2. Clustering

The main idea behind clustering is pretty straightforward. Basically, we say to ourselves, “I have these points here, and I can see that they organize into groups. It would be nice to describe these things more concretely, and, when a new point comes in, assign it to the correct group.” This general idea encourages exploration and opens up a variety of algorithms for clustering.

The algorithms listed below do not cover all the clustering methods out there, but they are the most commonly used ones.

K-means

K-means algorithm is the most popular and yet simplest of all the clustering algorithms. Here is how it works:

1. Select the number of clusters k that you think is the optimal number.
2. Initialize k points as “centroids” randomly within the space of our data.
3. Attribute each observation to its closest centroid.
4. Update the centroids to the center of all the attributed set of observations.
5. Repeat steps 3 and 4 a fixed number of times or until all of the centroids are stable (i.e. no longer change in step 4).

This algorithm is easy to describe and visualize. Let’s take a look.

Here, we used Euclidean distance, but the algorithm will converge with any other metric. You can not only vary the number of steps or the convergence criteria but also the distance measure between the points and cluster centroids.

Another “feature” of this algorithm is its sensitivity to the initial positions of the cluster centroids. You can run the algorithm several times and then average all the centroid results.

Choosing the number of clusters for K-means

In contrast to the supervised learning tasks such as classification and regression, clustering requires more effort to choose the optimization criterion. Usually, when working with k-means, we optimize the sum of squared distances between the observations and their centroids.

where C — is a set of clusters with power K, µ is a centroid of a cluster.

This definition seems reasonable — we want our observations to be as close to their centroids as possible. But, there is a problem — the optimum is reached when the number of centroids is equal to the number of observations, so you would end up with every single observation as its own separate cluster.

In order to avoid that case, we should choose a number of clusters after which a function J(C) is decreasing less rapidly. More formally,

Let’s look at an example.

We see that J(C) decreases significantly until the number of clusters is 3 and then does not change as much anymore. This means that the optimal number of clusters is 3.

Issues

Inherently, K-means is NP-hard. For d dimensions, k clusters, and n observations, we will find a solution in O(n^(dk+1)) in time. There are some heuristics to deal with this; an example is MiniBatch K-means, which takes portions (batches) of data instead of fitting the whole dataset and then moves centroids by taking the average of the previous steps. Compare the implementation of K-means and MiniBatch K-means in the sckit-learn documentation.

The implemetation of the algorithm using scikit-learn has its benefits such as the possibility to state the number of initializations with the n_init function parameter, which enables us to identify more robust centroids. Moreover, these runs can be done in parallel to decrease the computation time.

Affinity Propagation

Affinity propagation is another example of a clustering algorithm. As opposed to K-means, this approach does not require us to set the number of clusters beforehand. The main idea here is that we would like to cluster our data based on the similarity of the observations (or how they “correspond” to each other).

Let’s define a similarity metric such that s(x, y) > s(x, z) if an observation x is more similar to observation y and less similar to observation z. A simple example of such a similarity metric is a negative square of distance s(x, y) = -||x — y||².

Now, let’s describe “correspondence” by making two zero matrices. One of them, r, determines how well the k-th observation is as a “role model” for the i-th observation with respect to all other possible “role models”. Another matrix, a determines how appropriate it would be for i-th observation to take the k-th observation as a “role model”. This may sound confusing, but it becomes more understandable with some hands-on practice.

The matrices are updated sequentially with the following rules:

Spectral clustering

Spectral clustering combines some of the approaches described above to create a stronger clustering method.

First of all, this algorithm requires us to define the similarity matrix for observations called the adjacency matrix. This can be done in a similar fashion as in the Affinity Propagation algorithm, so that matrix A hosts negative square of the distances between the corresponding points. This matrix describes a full graph with the observations as vertices and the estimated similarity value between a pair of observations as edge weights for that pair of vertices. For the metric defined above and two-dimensional observations, this is pretty intuitive — two observations are similar if the edge between them is shorter. We’d like to split up the graph into two subgraphs in such a way that each observation in each subgraph would be similar to another observation in that subgraph. Formally, this is a Normalized cuts problem; for more details, we recommend reading this paper.

Agglomerative clustering

The following algorithm is the simplest and easiest to understand among all the the clustering algorithms without a fixed number of clusters.

The algorithm is fairly simple:

1. We start by assigning each observation to its own cluster
2. Then sort the pairwise distances between the centers of clusters in descending order
3. Take the nearest two neigbor clusters and merge them together, and recompute the centers
4. Repeat steps 2 and 3 until all the data is merged into one cluster

The process of searching for the nearest cluster can be conducted with different methods of bounding the observations:

1. Single linkage — distance between centroids is the minimum distance between all the pairs of points where one is from the first cluster, and the second is from the other.
2. Complete linkage — distance between centroids is the maximum distance between all the pairs of points where one is from the first cluster, and the second is from the other.
3. Average linkage — distance between centroids is the adjusted sum of all of the distances between all the pairs of points where one is from the first cluster, and the second is from the other.
4. Centroid linkage — distance between centroids is the distance between their geometric centers of clusters.

The 3rd one is the most effective in computation time since it does not require recomputing the distances every time the clusters are merged.

The results can be visualized as a beautiful cluster tree (dendogram) to help recognize the moment the algorithm should be stopped to get optimal results. There are plenty of Python tools to build these dendograms for agglomerative clustering.

Let’s consider an example with the clusters we got from K-means:

Accuracy metrics

As opposed to classfication, it is difficult to assess the quality of results from clustering. Here, a metric cannot depend on the labels but only on the goodness of split. Secondly, we do not usually have true labels of the observations when we use clustering.

There are internal and external goodness metrics. External metrics use the information about the known true split while internal metrics do not use any external information and assess the goodness of clusters based only on the initial data. The optimal number of clusters is usually defined with respect to some internal metrics.

All the metrics described below are implemented in sklearn.metrics.

Adjusted Rand Index (ARI)

Here, we assume that the true labels of objects are known. This metric does not depend on the labels’ values but on the data cluster split. Let N be the number of observations in a sample. Let a to be the number of observation pairs with the same labels and located in the same cluster, and let b to be the number of observations with different labels and located in different clusters. The Rand Index can be calculated using the following formula: RI = 2(a + b)/n(n — 1) In other words, it evaluates a share of observations for which these splits (initial and clustering result) are consistent. The Rand Index (RI) evaluates the similarity of the two splits of the same sample. In order for this index to be close to zero for any clustering outcomes with any n and number of clusters, it is essential to scale it, hence the Adjusted Rand Index: ARI = (RI — E(RI))/(max(RI) — E(RI)).

This metric is symmetric and does not depend in the label permutation. Therefore, this index is a measure of distances between different sample splits. ARI takes on values in the [-1, 1] range. Negative values indicate the independence of splits, and positive values indicate that these splits are consistent (they match ARI = 1).

Adjusted Mutual Information (AMI)

This metric is similar to ARI. It is also symmetric and does not depend on the labels’ values and permutation. It is defined by the [entropy](https://en.wikipedia.org/wiki/Entropy_(information_theory) function and interprets a sample split as a discrete distribution (likelihood of assigning to a cluster is equal to the percent of objects in it). The MI index is defined as the mutual information for two distributions, corresponding to the sample split into clusters. Intuitively, the mutual information measures the share of information common for both clustering splits i.e. how information about one of them decreases the uncertainty of the other one.

Similarly to the ARI, the AMI is defined. This allows us to get rid of the MI index’s increase with the number of clusters. The AMI lies in the [0, 1] range. Values close to zero mean the splits are independent, and those close to 1 mean they are similar (with complete match at AMI = 1).

Homogeneity, completeness, V-measure

Formally, these metrics are also defined based on the entropy function and the conditional entropy function, interpreting the sample splits as discrete distributions:

where K is a clustering result and C is the initial split. Therefore, h evaluates whether each cluster is composed of same class objects, and c measures how well the same class objects fit the clusters. These metrics are not symmetric. Both lie in the [0, 1] range, and values closer to 1 indicate more accurate clustering results. These metrics’ values are not scaled as the ARI or AMI metrics are and thus depend on the number of clusters. A random clustering result will not have metrics’ values closer to zero when the number of clusters is big enough and the number of objects is small. In such a case, it would be more reasonable to use ARI. However, with a large number of observations (more than 100) and the number of clusters less than 10, this issue is less critical and can be ignored.

V-measure is a combination of h, and c and is their harmonic mean: v = (2hc)/(h + c). It is symmetric and measures how consistent two clustering results are.

Silhouette

In contrast to the metrics described above, this coefficient does not imply the knowledge about the true labels of the objects. It lets us estimate the quality of the clustering using only the initial, unlabeled sample and the clustering result. To start with, for each observation, the silhouette coefficient is computed. Let a be the mean of the distance between an object and other objects within one cluster and b be the mean distance from an object to an object from the nearest cluster (different from the one the object belongs to). Then the silhouette measure for this object is s = (b — a)/max(a, b).

The silhouette of a sample is a mean value of silhouette values from this sample. Therefore, the silhouette distance shows to which extent the distance between the objects of the same class differ from the mean distance between the objects from different clusters. This coefficient takes values in the [-1, 1] range. Values close to -1 correspond to bad clustering results while values closer to 1 correspond to dense, well-defined clusters. Therefore, the higher the silhouette value is, the better the results from clustering.

With the help of silhouette, we can identify the optimal number of clusters k (if we don’t know it already from the data) by taking the number of clusters that maximizes the silhouette coefficient.

To conclude, let’s take a look at how these metrics perform with the MNIST handwritten numbers dataset:

4. Assignment #7

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

5. Useful links

• Overview of clustering methods in the scikit-learn doc.
• Q&A for PCA with examples
• Notebook on k-means and the EM-algorithm

Author: Sergey Korolev, Software Engineer at Snap Inc. Translated and edited by Anna Golovchenko, Yury Kashnitskiy, and Yuanyuan Pao.

Open Machine Learning Course. Topic 7. Unsupervised Learning: PCA and Clustering was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

In this course, we have already seen several key machine learning algorithms. However, before moving on to the more fancy ones, we’d like to take a small detour and talk about data preparation. The well-known concept of “garbage in — garbage out” applies 100% to any task in machine learning. Any experienced professional can recall numerous times when a simple model trained on high-quality data was proven to be better than a complicated multi-model ensemble built on data that wasn’t clean.

To start, I wanted to review three similar but different tasks:

• feature extraction and feature engineering: transformation of raw data into features suitable for modeling;
• feature transformation: transformation of data to improve the accuracy of the algorithm;
• feature selection: removing unnecessary features.

This article will contain almost no math, but there will be a fair amount of code. (Jupyter version). Some examples will use the dataset from Renthop company, which is used in the Two Sigma Connect: Rental Listing Inquiries Kaggle competition. In this task, you need to predict the popularity of a new rental listing, i.e. classify the listing into three classes: ['low', 'medium' , 'high']. To evaluate the solutions, we will use the log loss metric (the smaller, the better). Those who do not have a Kaggle account, will have to register; you will also need to accept the rules of the competition in order to download the data.

import json
import pandas as pd

# Let's load the dataset from Renthop right away
with open('train.json', 'r') as raw_data:
    data = json.load(raw_data)
    df = pd.DataFrame(data)

Article outline

1. Feature Extraction
1.1. Texts
1.2. Images
1.3. Geospatial data
1.4. Date and time
1.5. Time series, web, etc.

2. Feature transformations
2.1. Normalization and changing distribution
2.2. Interactions
2.3. Filling in the missing values

3. Feature selection
3.1. Statistical approaches
3.2. Selection by modeling
3.3. Grid search

1. Feature Extraction

In practice, data rarely comes in the form of ready-to-use matrices. That's why every task begins with feature extraction. Sometimes, it can be enough to read the csv file and convert it into numpy.array, but this is a rare exception. Let's look at some of the popular types of data from which features can be extracted.

1.1. Texts

Text is a type of data that can come in different formats; there are so many text processing methods that cannot fit in a single article. Nevertheless, we will review the most popular ones.

Before working with text, one must tokenize it. Tokenization implies splitting the text into units (hence, tokens). Most simply, tokens are just the words. But splitting by word can lose some of the meaning -- "Santa Barbara" is one token, not two, but "rock'n'roll" should not be split into two tokens. There are ready-to-use tokenizers that take into account peculiarities of the language, but they make mistakes as well, especially when you work with specific sources of text (newspapers, slang, misspellings, typos).

After tokenization, you will normalize the data. For text, this is about stemming and/or lemmatization; these are similar processes used to process different forms of a word. One can read about the difference between them here.

So, now that we have turned the document into a sequence of words, we can represent it with vectors. The easiest approach is called Bag of Words: we create a vector with the length of the dictionary, compute the number of occurrences of each word in the text, and place that number of occurrences in the appropriate position in the vector. The process described looks simpler in code:

from functools import reduce 
import numpy as np

texts = [['i', 'have', 'a', 'cat'], 
        ['he', 'have', 'a', 'dog'], 
        ['he', 'and', 'i', 'have', 'a', 'cat', 'and', 'a', 'dog']]

dictionary = list(enumerate(set(list(reduce(lambda x, y: x + y, texts)))))

def vectorize(text): 
    vector = np.zeros(len(dictionary)) 
    for i, word in dictionary: 
        num = 0 
        for w in text: 
            if w == word: 
                num += 1 
        if num: 
            vector[i] = num 
    return vector

for t in texts: 
    print(vectorize(t))

Expected result:

[0. 1. 0. 1. 1. 0. 1.]
[0. 1. 1. 0. 1. 1. 0.]
[2. 1. 1. 1. 2. 1. 1.]

Here is an illustration of the process:

This is an extremely naive implementation. In practice, you need to consider stop words, the maximum length of the dictionary, more efficient data structures (usually text data is converted to a sparse vector), etc.

When using algorithms like Bag of Words, we lose the order of the words in the text, which means that the texts "i have no cows" and "no, i have cows" will appear identical after vectorization when, in fact, they have the opposite meaning. To avoid this problem, we can revisit our tokenization step and use N-grams (the sequence of N consecutive tokens) instead.

In : from sklearn.feature_extraction.text import CountVectorizer

In : vect = CountVectorizer(ngram_range=(1,1))

In : vect.fit_transform(['no i have cows', 'i have no cows']).toarray() 
Out: 
array([[1, 1, 1],
      [1, 1, 1]], dtype=int64)

In : vect.vocabulary_ 
Out: {'cows': 0, 'have': 1, 'no': 2}

In : vect = CountVectorizer(ngram_range=(1,2))

In : vect.fit_transform(['no i have cows', 'i have no cows']).toarray() 
Out: 
array([[1, 1, 1, 0, 1, 0, 1], 
      [1, 1, 0, 1, 1, 1, 0]], dtype=int64)

In : vect.vocabulary_ 
Out: {'cows': 0, 
      'have': 1,
      'have cows': 2,
      'have no': 3,
      'no': 4,
      'no cows': 5,
      'no have': 6}

Also note that one does not have to use only words. In some cases, it is possible to generate N-grams of characters. This approach would be able to account for similarity of related words or handle typos.

from scipy.spatial.distance import euclidean
from sklearn.feature_extraction.text import CountVectorizer

vect = CountVectorizer(ngram_range=(3,3), analyzer='char_wb')

n1, n2, n3, n4 = vect.fit_transform(['andersen', 'petersen', 'petrov', 'smith']).toarray()

euclidean(n1, n2), euclidean(n2, n3), euclidean(n3, n4)
# (2.8284271247461903, 3.1622776601683795, 3.3166247903554)

Adding onto the Bag of Words idea: words that are rarely found in the corpus (in all the documents of this dataset) but are present in this particular document might be more important. Then it makes sense to increase the weight of more domain-specific words to separate them out from common words. This approach is called TF-IDF (term frequency-inverse document frequency), which cannot be written in a few lines, so you should look into the details in references such as this wiki. The default option is as follows:

Analogs of Bag of Words can be found outside of text problems e.g. bag of sites in the Catch Me If You Can competition, bag of apps, bag of events, etc.

Using these algorithms, it is possible to obtain a working solution for a simple problem, which can serve as a baseline. However, for those who do not like the classics, there are new approaches. The most popular method in the new wave is Word2Vec, but there are a few alternatives as well (GloVe, Fasttext, etc.).

Word2Vec is a special case of the word embedding algorithms. Using Word2Vec and similar models, we can not only vectorize words in a high-dimensional space (typically a few hundred dimensions) but also compare their semantic similarity. This is a classic example of operations that can be performed on vectorized concepts: king - man + woman = queen.

It is worth noting that this model does not comprehend the meaning of the words but simply tries to position the vectors such that words used in common context are close to each other. If this is not taken into account, a lot of fun examples will come up.

Such models need to be trained on very large datasets in order for the vector coordinates to capture the semantics. A pretrained model for your own tasks can be downloaded here.

Similar methods are applied in other areas such as bioinformatics. An unexpected application is food2vec. You can probably think of a few other fresh ideas; the concept is universal enough.

1.2. Images

Working with images is easier and harder at the same time. It is easier because it is possible to just use one of the popular pretrained networks without much thinking but harder because, if you need to dig into the details, you may end up going really deep. Let's start from the beginning.

In a time when GPUs were weaker and the "renaissance of neural networks" had not happened yet, feature generation from images was its own complex field. One had to work at a low level, determining corners, borders of regions, color distributions statistics, and so on. Experienced specialists in computer vision could draw a lot of parallels between older approaches and neural networks; in particular, convolutional layers in today's networks are similar to Haar cascades. If you are interested in reading more, here are a couple of links to some interesting libraries: skimage and SimpleCV.

Often for problems associated with images, a convolutional neural network is used. You do not have to come up with the architecture and train a network from scratch. Instead, download a pretrained state-of-the-art network with the weights from public sources. Data scientists often do so-called fine-tuning to adapt these networks to their needs by "detaching" the last fully connected layers of the network, adding new layers chosen for a specific task, and then training the network on new data. If your task is to just vectorize the image (for example, to use some non-network classifier), you only need to remove the last layers and use the output from the previous layers:

from keras.applications.resnet50 import ResNet50, preprocess_input
from keras.preprocessing import image 
from scipy.misc import face 
import numpy as np

resnet_settings = {'include_top': False, 'weights': 'imagenet'}
resnet = ResNet50(**resnet_settings)

img = image.array_to_img(face())

# What a cute raccoon!

img = img.resize((224, 224))

# In real life, you may need to pay more attention to resizing

x = image.img_to_array(img) 
x = np.expand_dims(x, axis=0)
x = preprocess_input(x)

# Need an extra dimension because model is designed to work with an array of images - i.e. tensor shaped (batch_size, width, height, n_channels)

features = resnet.predict(x)

Here's a classifier trained on one dataset and adapted for a different one by "detaching" the last layer and adding a new one instead.

Nevertheless, we should not focus too much on neural network techniques. Features generated by hand are still very useful: for example, for predicting the popularity of a rental listing, we can assume that bright apartments attract more attention and create a feature such as "the average value of the pixel". You can find some inspiring examples in the documentation of relevant libraries.

If there is text on the image, you can read it without unraveling a complicated neural network. For example, check out pytesseract.

In : import pytesseract

In : from PIL import Image

In : import requests

In : from io import BytesIO

In : img = 'http://ohscurrent.org/wp-content/uploads/2015/09/domus-01-google.jpg'

##### Just a random picture from search

In : img = requests.get(img) 
...: img = Image.open(BytesIO(img.content)) 
...: text = pytesseract.image_to_string(img) 
...:

In : text 
Out: 'Google'

One must understand that pytesseract is not a solution for everything.

##### This time we take a picture from Renthop

In : img = requests.get('https://photos.renthop.com/2/8393298_6acaf11f030217d05f3a5604b9a2f70f.jpg') 
...: img = Image.open(BytesIO(img.content)) 
...: pytesseract.image_to_string(img) 
...:

Out: 'Cunveztible to 4}»'

Another case where neural networks cannot help is extracting features from meta-information. For images, EXIF stores many useful meta-information: manufacturer and camera model, resolution, use of the flash, geographic coordinates of shooting, software used to process image and more.

1.3. Geospatial data

Geographic data is not so often found in problems, but it is still useful to master the basic techniques for working with it, especially since there are quite a number of ready-to-use solutions in this field.

Geospatial data is often presented in the form of addresses or coordinates of (Latitude, Longitude). Depending on the task, you may need two mutually-inverse operations: geocoding (recovering a point from an address) and reverse geocoding (recovering an address from a point). Both operations are accessible in practice via external APIs from Google Maps or OpenStreetMap. Different geocoders have their own characteristics, and the quality varies from region to region. Fortunately, there are universal libraries like geopy that act as wrappers for these external services.

If you have a lot of data, you will quickly reach the limits of external API. Besides, it is not always the fastest to receive information via HTTP. Therefore, it is necessary to consider using a local version of OpenStreetMap.

If you have a small amount of data, enough time, and no desire to extract fancy features, you can use reverse_geocoder in lieu of OpenStreetMap:

In : import reverse_geocoder as revgc

In : revgc.search((df.latitude, df.longitude)) 
Loading formatted geocoded file... 
Out: [OrderedDict([('lat', '40.74482'), 
                   ('lon', '-73.94875'), 
                   ('name', 'Long Island City'), 
                   ('admin1', 'New York'), 
                   ('admin2', 'Queens County'), 
                   ('cc', 'US')])]

When working with geoсoding, we must not forget that addresses may contain typos, which makes the data cleaning step necessary. Coordinates contain fewer misprints, but its position can be incorrect due to GPS noise or bad accuracy in places like tunnels, downtown areas, etc. If the data source is a mobile device, the geolocation may not be determined by GPS but by WiFi networks in the area, which leads to holes in space and teleportation. While traveling along in Manhattan, there can suddenly be a WiFi location from Chicago.

WiFi location tracking is based on the combination of SSID and MAC-addresses, which may correspond to different points e.g. federal provider standardizes the firmware of routers up to MAC-address and places them in different cities. Even a company's move to another office with its routers can cause issues.

The point is usually located among infrastructure. Here, you can really unleash your imagination and invent features based on your life experience and domain knowledge: the proximity of a point to the subway, the number of stories in the building, the distance to the nearest store, the number of ATMs around, etc. For any task, you can easily come up with dozens of features and extract them from various external sources. For problems outside an urban environment, you may consider features from more specific sources e.g. the height above sea level.

If two or more points are interconnected, it may be worthwhile to extract features from the route between them. In that case, distances (great circle distance and road distance calculated by the routing graph), number of turns with the ratio of left to right turns, number of traffic lights, junctions, and bridges will be useful. In one of my own tasks, I generated a feature called "the complexity of the road", which computed the graph-calculated distance divided by the GCD.

1.4. Date and time

You would think that date and time are standardized because of their prevalence, but, nevertheless, some pitfalls remain.

Let's start with the day of the week, which are easy to turn into 7 dummy variables using one-hot encoding. In addition, we will also create a separate binary feature for the weekend called is_weekend.

df['dow'] = df['created'].apply(lambda x: x.date().weekday())
df['is_weekend'] = df['created'].apply(lambda x: 1 if x.date().weekday() in (5, 6) else 0)

Some tasks may require additional calendar features. For example, cash withdrawals can be linked to a pay day; the purchase of a metro card, to the beginning of the month. In general, when working with time series data, it is a good idea to have a calendar with public holidays, abnormal weather conditions, and other important events.

Q: What do Chinese New Year, the New York marathon, and the Trump inauguration have in common?

A: They all need to be put on the calendar of potential anomalies.

Dealing with hour (minute, day of the month ...) is not as simple as it seems. If you use the hour as a real variable, we slightly contradict the nature of data: 0<23 while 0:00:00 02.01> 01.01 23:00:00. For some problems, this can be critical. At the same time, if you encode them as categorical variables, you'll breed a large numbers of features and lose information about proximity -- the difference between 22 and 23 will be the same as the difference between 22 and 7.

There also exist some more esoteric approaches to such data like projecting the time onto a circle and using the two coordinates.

def make_harmonic_features(value, period=24):
    value *= 2 * np.pi / period 
    return np.cos(value), np.sin(value)

This transformation preserves the distance between points, which is important for algorithms that estimate distance (kNN, SVM, k-means ...)

In : from scipy.spatial import distance

In : euclidean(make_harmonic_features(23), make_harmonic_features(1)) 
Out: 0.5176380902050424

In : euclidean(make_harmonic_features(9), make_harmonic_features(11)) 
Out: 0.5176380902050414

In : euclidean(make_harmonic_features(9), make_harmonic_features(21)) 
Out: 2.0

However, the difference between such coding methods is down to the third decimal place in the metric.

1.5. Time series, web, etc.

Regarding time series — we will not go into too much detail here (mostly due to my personal lack of experience), but I will point you to a useful library that automatically generates features for time series.

If you are working with web data, then you usually have information about the user's User Agent. It is a wealth of information. First, one needs to extract the operating system from it. Secondly, make a feature is_mobile. Third, look at the browser.

In : ua = 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Ubuntu Chromium/56.0.2924.76 Chrome/ 
...: 56.0.2924.76 Safari/537.36'

In : import user_agents

In : ua = user_agents.parse(ua)

In : ua.is_bot 
Out: False

In : ua.is_mobile 
Out: False

In : ua.is_pc 
Out: True

In : ua.os.family 
Out: 'Ubuntu'

In : ua.os.version 
Out: ()

In : ua.browser.family 
Out: 'Chromium'

In : ua.os.version 
Out: ()

In : ua.browser.version 
Out: (56, 0, 2924)

As in other domains, you can come up with your own features based on intuition about the nature of the data. At the time of this writing, Chromium 56 was new, but, after some time, only users who haven't rebooted their browser for a long time will have this version. In this case, why not introduce a feature called "lag behind the latest version of the browser"?

In addition to the operating system and browser, you can look at the referrer (not always available), http_accept_language, and other meta information.

The next useful piece of information is the IP-address, from which you can extract the country and possibly the city, provider, and connection type (mobile/stationary). You need to understand that there is a variety of proxy and outdated databases, so this feature can contain noise. Network administration gurus may try to extract even fancier features like suggestions for using VPN. By the way, the data from the IP-address is well combined with http_accept_language: if the user is sitting at the Chilean proxies and browser locale is ru_RU, something is unclean and worth a look in the corresponding column in the table (is_traveler_or_proxy_user).

Any given area has so many specifics that it is too much for an individual to absorb completely. Therefore, I invite everyone to share their experiences and discuss feature extraction and generation in the comments section.

2. Feature transformations

2.1. Normalization and changing distribution

Monotonic feature transformation is critical for some algorithms and has no effect on others. This is one of the reasons for the increased popularity of decision trees and all its derivative algorithms (random forest, gradient boosting). Not everyone can or want to tinker with transformations, and these algorithms are robust to unusual distributions.

There are also purely engineering reasons: np.log is a way of dealing with large numbers that do not fit in np.float64. This is an exception rather than a rule; often it's driven by the desire to adapt the dataset to the requirements of the algorithm. Parametric methods usually require a minimum of symmetric and unimodal distribution of data, which is not always given in real data. There may be more stringent requirements; recall our earlier article about linear models.

However, data requirements are imposed not only by parametric methods; K nearest neighbors will predict complete nonsense if features are not normalized e.g. when one distribution is located in the vicinity of zero and does not go beyond (-1, 1) while the other’s range is on the order of hundreds of thousands.

A simple example: suppose that the task is to predict the cost of an apartment from two variables — the distance from city center and the number of rooms. The number of rooms rarely exceeds 5 whereas the distance from city center can easily be in the thousands of meters.

The simplest transformation is Standard Scaling (or Z-score normalization):

Note that Standard Scaling does not make the distribution normal in the strict sense.

In : from sklearn.preprocessing import StandardScaler

In : from scipy.stats import beta

In : from scipy.stats import shapiro

In : data = beta(1, 10).rvs(1000).reshape(-1, 1)

In : shapiro(data) 
Out: (0.8783774375915527, 3.0409122263582326e-27)

# Value of the statistic, p-value

In : shapiro(StandardScaler().fit_transform(data)) 
Out: (0.8783774375915527, 3.0409122263582326e-27)

# With such p-value we'd have to reject the null hypothesis of normality of the data

But, to some extent, it protects against outliers:

In : data = np.array([1, 1, 0, -1, 2, 1, 2, 3, -2, 4, 100]).reshape(-1, 1).astype(np.float64)

In : StandardScaler().fit_transform(data) 
Out: 
array([[-0.31922662], 
       [-0.31922662],
       [-0.35434155],
       [-0.38945648],
       [-0.28411169],
       [-0.31922662],
       [-0.28411169],
       [-0.24899676],
       [-0.42457141],
       [-0.21388184],
       [ 3.15715128]])

In : (data – data.mean()) / data.std() 
Out: 
array([[-0.31922662], 
       [-0.31922662],
       [-0.35434155],
       [-0.38945648],
       [-0.28411169],
       [-0.31922662],
       [-0.28411169],
       [-0.24899676],
       [-0.42457141],
       [-0.21388184],
       [ 3.15715128]])

Another fairly popular option is MinMax Scaling, which brings all the points within a predetermined interval (typically (0, 1)).

In : from sklearn.preprocessing import MinMaxScaler

In : MinMaxScaler().fit_transform(data) 
Out: 
array([[ 0.02941176],
       [ 0.02941176],
       [ 0.01960784],
       [ 0.00980392],
       [ 0.03921569],
       [ 0.02941176],
       [ 0.03921569],
       [ 0.04901961],
       [ 0. ],
       [ 0.05882353],
       [ 1. ]])

In : (data – data.min()) / (data.max() – data.min()) 
Out: 
array([[ 0.02941176],
       [ 0.02941176],
       [ 0.01960784],
       [ 0.00980392],
       [ 0.03921569],
       [ 0.02941176],
       [ 0.03921569],
       [ 0.04901961],
       [ 0. ],
       [ 0.05882353],
       [ 1. ]])

StandardScaling and MinMax Scaling have similar applications and are often more or less interchangeable. However, if the algorithm involves the calculation of distances between points or vectors, the default choice is StandardScaling. But MinMax Scaling is useful for visualization by bringing features within the interval (0, 255).

If we assume that some data is not normally distributed but is described by the log-normal distribution, it can easily be transformed to a normal distribution:

In : from scipy.stats import lognorm

In : data = lognorm(s=1).rvs(1000)

In : shapiro(data)
Out: (0.05714237689971924, 0.0)

In : shapiro(np.log(data))
Out: (0.9980740547180176, 0.3150389492511749)

The lognormal distribution is suitable for describing salaries, price of securities, urban population, number of comments on articles on the internet, etc. However, to apply this procedure, the underlying distribution does not necessarily have to be lognormal; you can try to apply this transformation to any distribution with a heavy right tail. Furthermore, one can try to use other similar transformations, formulating their own hypotheses on how to approximate the available distribution to a normal. Examples of such transformations are Box-Cox transformation (logarithm is a special case of the Box-Cox transformation) or Yeo-Johnson transformation (extends the range of applicability to negative numbers). In addition, you can also try adding a constant to the feature — np.log (x + const).

In the examples above, we have worked with synthetic data and strictly tested normality using the Shapiro-Wilk test. Let’s try to look at some real data and test for normality using a less formal method — Q-Q plot. For a normal distribution, it will look like a smooth diagonal line, and visual anomalies should be intuitively understandable.

Q-Q plot for lognormal distribution

Q-Q plot for the same distribution after taking the logarithm

# Let's draw plots!
In : import statsmodels.api as sm

Let's take the price feature from Renthop dataset and filter by hands the most extreme values for clarity

In : price = df.price[(df.price <= 20000) & (df.price > 500)]

In : price_log = np.log(price)

In : price_mm = MinMaxScaler().fit_transform(price.values.reshape(-1, 1).astype(np.float64)).flatten()

A lot of gestures so that sklearn didn't shower us with warnings

In : price_z = StandardScaler().fit_transform(price.values.reshape(-1, 1).astype(np.float64)).flatten()

In : sm.qqplot(price_log, loc=price_log.mean(), scale=price_log.std()).savefig('qq_price_log.png')

In : sm.qqplot(price_mm, loc=price_mm.mean(), scale=price_mm.std()).savefig('qq_price_mm.png')

In : sm.qqplot(price_z, loc=price_z.mean(), scale=price_z.std()).savefig('qq_price_z.png')

Q-Q plot of the initial feature

Q-Q plot after StandardScaler. Shape doesn’t change

Q-Q plot after MinMaxScaler. Shape doesn’t change

Q-Q plot after taking the logarithm. Things are getting better!

Let’s see whether transformations can somehow help the real model. There is no silver bullet here.

2.2. Interactions

If previous transformations seemed rather math-driven, this part is more about the nature of the data; it can be attributed to both feature transformations and feature creation.

Let’s come back again to the Two Sigma Connect: Rental Listing Inquiries problem. Among the features in this problem are the number of rooms and the price. Logic suggests that the cost per single room is more indicative than the total cost, so we can generate such a feature.

rooms = df["bedrooms"].apply(lambda x: max(x, .5))

# Avoid division by zero; .5 is chosen more or less arbitrarily

df["price_per_bedroom"] = df["price"] / rooms

You should limit yourself in this process. If there are a limited number of features, it is possible to generate all the possible interactions and then weed out the unnecessary ones using the techniques described in the next section. In addition, not all interactions between features must have a physical meaning; for example, polynomial features (see sklearn.preprocessing.PolynomialFeatures) are often used in linear models and are almost impossible to interpret.

2.3. Filling in the missing values

Not many algorithms can work with missing values, and the real world often provides data with gaps. Fortunately, this is one of the tasks for which one doesn’t need any creativity. Both key python libraries for data analysis provide easy-to-use solutions: pandas.DataFrame.fillna and sklearn.preprocessing.Imputer.

These solutions do not have any magic happening behind the scenes. Approaches to handling missing values are pretty straightforward:

• encode missing values with a separate blank value like "n/a" (for categorical variables);
• use the most probable value of the feature (mean or median for the numerical variables, the most common value for categorical variables);
• or, conversely, encode with some extreme value (good for decision-tree models since it allows the model to make a partition between the missing and non-missing values);
• for ordered data (e.g. time series), take the adjacent value — next or previous.

Easy-to-use library solutions sometimes suggest sticking to something like df = df.fillna(0) and not sweat the gaps. But this is not the best solution: data preparation takes more time than building models, so thoughtless gap-filling may hide a bug in processing and damage the model.

3. Feature selection

Why would it even be necessary to select features? To some, this idea may seem counterintuitive, but there are at least two important reasons to get rid of unimportant features. The first is clear to every engineer: the more data, the higher the computational complexity. As long as we work with toy datasets, the size of the data is not a problem, but, for real loaded production systems, hundreds of extra features will be quite tangible. The second reason is that some algorithms take noise (non-informative features) as a signal and overfit.

3.1. Statistical approaches

The most obvious candidate for removal is a feature whose value remains unchanged, i.e., it contains no information at all. If we build on this thought, it is reasonable to say that features with low variance are worse than those with high variance. So, one can consider cutting features with variance below a certain threshold.

In : from sklearn.feature_selection import VarianceThreshold

In : from sklearn.datasets import make_classification

In : x_data_generated, y_data_generated = make_classification()

In : x_data_generated.shape Out: (100, 20)

In : VarianceThreshold(.7).fit_transform(x_data_generated).shape Out: (100, 19)

In : VarianceThreshold(.8).fit_transform(x_data_generated).shape Out: (100, 18)

In : VarianceThreshold(.9).fit_transform(x_data_generated).shape Out: (100, 15)

There are other ways that are also based on classical statistics.

In : from sklearn.feature_selection import SelectKBest, f_classif

In : x_data_kbest = SelectKBest(f_classif, k=5).fit_transform(x_data_generated, y_data_generated)

In : x_data_varth = VarianceThreshold(.9).fit_transform(x_data_generated)

In : from sklearn.linear_model import LogisticRegression

In : from sklearn.model_selection import cross_val_score

In : cross_val_score(LogisticRegression(), x_data_generated, y_data_generated, scoring='neg_log_loss').mean() Out: -0.45367136377981693

In : cross_val_score(LogisticRegression(), x_data_kbest, y_data_generated, scoring='neg_log_loss').mean() Out: -0.35775228616521798

In : cross_val_score(LogisticRegression(), x_data_varth, y_data_generated, scoring='neg_log_loss').mean() Out: -0.44033042718359772

We can see that our selected features have improved the quality of the classifier. Of course, this example is purely artificial; however, it is worth using for real problems.

3.2. Selection by modeling

Another approach is to use some baseline model for feature evaluation because the model will clearly show the importance of the features. Two types of models are usually used: some “wooden” composition such as Random Forest or a linear model with Lasso regularization so that it is prone to nullify weights of weak features. The logic is intuitive: if features are clearly useless in a simple model, there is no need to drag them to a more complex one.

# Synthetic example

from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectFromModel
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline

x_data_generated, y_data_generated = make_classification()

pipe = make_pipeline(SelectFromModel(estimator=RandomForestClassifier()), LogisticRegression())

lr = LogisticRegression()
rf = RandomForestClassifier()

print(cross_val_score(lr, x_data_generated, y_data_generated, scoring='neg_log_loss').mean())
print(cross_val_score(rf, x_data_generated, y_data_generated, scoring='neg_log_loss').mean())
print(cross_val_score(pipe, x_data_generated, y_data_generated, scoring='neg_log_loss').mean())

# -0.455661299521484
# -0.6817083161800542
# -0.28572375767268404

We must not forget that this is not a silver bullet again — it can make the performance worse.

x_data, y_data = get_data() x_data = x_data.values

pipe1 = make_pipeline(StandardScaler(), SelectFromModel(estimator=RandomForestClassifier()), LogisticRegression())

pipe2 = make_pipeline(StandardScaler(), LogisticRegression())

rf = RandomForestClassifier()

print('LR + selection: ', cross_val_score(pipe1, x_data, y_data, scoring='neg_log_loss').mean())
print('LR: ', cross_val_score(pipe2, x_data, y_data, scoring='neg_log_loss').mean())
print('RF: ', cross_val_score(rf, x_data, y_data, scoring='neg_log_loss').mean())

# LR + selection: -0.714208124619
# LR: -0.669572736183

# It got worse!

# RF: -2.13486716798

3.3. Grid search

Finally, we get to the most reliable method, which is also the most computationally complex: trivial grid search. Train a model on a subset of features, store results, repeat for different subsets, and compare the quality of models to identify the best feature set. This approach is called Exhaustive Feature Selection.

Searching all combinations usually takes too long, so you can try to reduce the search space. Fix a small number N, iterate through all combinations of N features, choose the best combination, and then iterate through the combinations of (N + 1) features so that the previous best combination of features is fixed and only a single new feature is considered. It is possible to iterate until we hit a maximum number of characteristics or until the quality of the model ceases to increase significantly. This algorithm is called Sequential Feature Selection.

This algorithm can be reversed: start with the complete feature space and remove features one by one until it does not impair the quality of the model or until the desired number of features is reached.

In : selector = SequentialFeatureSelector(LogisticRegression(), scoring='neg_log_loss', verbose=2, k_features=3, forward=False, n_jobs=-1)

In : selector.fit(x_data_scaled, y_data)

In : selector.fit(x_data_scaled, y_data)

[2017-03-30 01:42:24] Features: 45/3 -- score: -0.682830838803 [2017-03-30 01:44:40] Features: 44/3 -- score: -0.682779463265 [2017-03-30 01:46:47] Features: 43/3 -- score: -0.682727480522 [2017-03-30 01:48:54] Features: 42/3 -- score: -0.682680521828 [2017-03-30 01:50:52] Features: 41/3 -- score: -0.68264297879
[2017-03-30 01:52:46] Features: 40/3 -- score: -0.682607753617 [2017-03-30 01:54:37] Features: 39/3 -- score: -0.682570678346 [2017-03-30 01:56:21] Features: 38/3 -- score: -0.682536314625 [2017-03-30 01:58:02] Features: 37/3 -- score: -0.682520258804 [2017-03-30 01:59:39] Features: 36/3 -- score: -0.68250862986
[2017-03-30 02:01:17] Features: 35/3 -- score: -0.682498213174
...
[2017-03-30 02:21:09] Features: 10/3 -- score: -0.68657335969
[2017-03-30 02:21:18] Features: 9/3 -- score: -0.688405548594
[2017-03-30 02:21:26] Features: 8/3 -- score: -0.690213724719
[2017-03-30 02:21:32] Features: 7/3 -- score: -0.692383588303
[2017-03-30 02:21:36] Features: 6/3 -- score: -0.695321584506
[2017-03-30 02:21:40] Features: 5/3 -- score: -0.698519960477
[2017-03-30 02:21:42] Features: 4/3 -- score: -0.704095390444
[2017-03-30 02:21:44] Features: 3/3 -- score: -0.713788301404

But improvement couldn’t last forever

Take a look how this approach was done in one simple yet elegant Kaggle kernel.

4. Assignment #6

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer. We don’t have a separate assignment on feature engineering (this is done in Kaggle Inclass competitions) so this demo assignment is on regression due.

Open Data Science wishes you good luck with the assignment as well as clean data and useful features in your real work!

Author: Arseny Kravchenko. Translated and edited by Christina Butsko, Yury Kashnitskiy, Anna Larionova, Anastasia Manokhina, Evgeny Sushko, Egor Polusmak, and Yuanyuan Pao.

Open Machine Learning Course. Topic 6. Feature Engineering and Feature Selection was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

In the previous articles, you saw different classification algorithms as well as techniques for how to properly validate and evaluate the quality of your models.

Now, suppose that you have chosen the best possible model for a particular problem and are struggling to further improve its accuracy. In this case, you would need to apply some more advanced machine learning techniques that are collectively referred to as ensembles.

An ensemble is a set of elements that collectively contribute to a whole. A familiar example is a musical ensemble, which blends the sounds of several musical instruments to create a beautiful harmony, or architectural ensembles, which are a set of buildings designed as a unit. In ensembles, the (whole) harmonious outcome is more important than the performance of any individual part.

Article outline

1. Ensembles
2. Bootstraping
3. Bagging
4. Out-of-bag error
5. Random Forest
6. Feature importance
7. Assignment #5
8. Useful resources

1. Ensembles

Condorcet’s jury theorem (1784) is about an ensemble in some sense. It states that, if each member of the jury makes an independent judgement and the probability of the correct decision by each juror is more than 0.5, then the probability of the correct decision by the whole jury increases with the total number of jurors and tends to one. On the other hand, if the probability of being right is less than 0.5 for each juror, then the probability of the correct decision by the whole jury decreases with the number of jurors and tends to zero.

One can find an analytic expression for this theorem in a Jupyter notebook accompanying this article.

Let’s look at another example of ensembles: an observation known as Wisdom of the crowd. In 1906, Francis Galton visited a country fair in Plymouth where he saw a contest being held for farmers. 800 participants tried to estimate the weight of a slaughtered bull. The real weight of the bull was 1198 pounds. Although none of the farmers could guess the exact weight of the animal, the average of their predictions was 1197 pounds.

A similar idea for error reduction was adopted in the field of Machine Learning.

2. Bootstraping

Bagging (also known as Bootstrap aggregation) is one of the first and most basic ensemble techniques. It was proposed by Leo Breiman in 1994. Bagging is based on the statistical method of bootstraping, which makes the evaluation of many statistics of complex models feasible.

The bootstrap method goes as follows. Let there be a sample X of size N. We can make a new sample from the original sample by drawing N elements from the latter randomly and uniformly, with replacement. In other words, we select a random element from the original sample of size N and do this N times. All elements are equally likely to be selected, thus each element is drawn with the equal probability 1/N.

Let’s say we are drawing balls from a bag one at a time. At each step, the selected ball is put back into the bag so that the next selection is made equiprobably i.e. from the same number of balls N. Note that, because we put the balls back, there may be duplicates in the new sample. Let’s call this new sample X1.

By repeating this procedure M times, we create M bootstrap samples X1, …, XM. In the end, we have a sufficient number of samples and can compute various statistics of the original distribution.

For our example, we’ll use the familiar telecom_churn dataset. Previously, when we discussed feature importance, we saw that one of the most important features in this dataset is the number of calls to customer service. Let's visualize the data and look at the distribution of this feature.

As you can see, loyal customers make fewer calls to customer service than those who eventually left. Now, it might be a good idea to estimate the average number of customer service calls in each group. Since our dataset is small, we would not get a good estimate by simply calculating the mean of the original sample. We will be better off applying the bootstrap method. Let’s generate 1000 new bootstrap samples from our original population and produce an interval estimate of the mean.

Service calls from loyal: mean interval [1.4077193 1.49473684] # Service calls from churn: mean interval [2.0621118 2.39761905]

In the end, we see that, with 95% probability, the average number of customer service calls from loyal customers lies between 1.4 and 1.49 while the churned clients called 2.06 through 2.40 times on average. Also, note that the interval for the loyal customers is narrower, which is reasonable since they make fewer calls (0, 1 or 2) in comparison with the churned clients who called until they became fed up and switched providers.

3. Bagging

Now that you’ve grasped the idea of bootstrapping, we can move on to bagging.

Suppose that we have a training set X. Using bootstrapping, we generate samples X1, …, XM. Now, for each bootstrap sample, we train its own classifier ai(x). The final classifier will average the outputs from all these individual classifiers. In the case of classification, this technique corresponds to voting:

The picture below illustrates this algorithm:

In a regression problem, by averaging individual answers, bagging reduces the mean squared error by a factor of M, the number of regressors. See the mathematical proof in our Jupyter notebook.

From our previous lesson, let’s recall the components that make up the total out-of-sample error:

Bagging reduces the variance of a classifier by decreasing the difference in error when we train the model on different datasets. In other words, bagging prevents overfitting. The efficacy of bagging comes from the fact that individual models are quite different due to the different training data and their errors cancel out during voting. Additionally, outliers are likely omitted in some of the training bootstap samples.

The scikit-learn library supports bagging with meta-estimators BaggingRegressor and BaggingClassifier. You can use most of the algorithms as a base.

Let’s examine how bagging works in practice and compare it with the decision tree. For this, we will use an example from sklearn’s documentation.

The error for the decision tree:

0.0255 = 0.0003 (bias²)+ 0.0152 (variance) + 0.0098 (σ²)

The error when using bagging:

0.0196 = 0.0004 (bias²) + 0.0092 (variance) + 0.0098 (σ²)

As you can see from the graph above, the variance in the error is much lower for bagging. Remember that we have already proved this theoretically.

Bagging is effective on small datasets. Dropping even a small part of training data leads to constructing substantially different base classifiers. If you have a large dataset, you would generate bootstrap samples of a much smaller size.

The example above is unlikely to be applicable to any real work. This is because we made a strong assumption that our individual errors are uncorrelated. More often than not, this is way too optimistic for real-world applications. When this assumption is false, the reduction in error will not be as significant. In the following lectures, we will discuss some more sophisticated ensemble methods, which enable more accurate predictions in real-world problems.

4. Out-of-bag error

Looking ahead, in case of Random Forest, there is no need to use cross-validation or hold-out samples in order to get an unbiased error estimation. Why? Because, in ensemble techniques, the error estimation takes place internally.

Random trees are constructed using different bootstrap samples of the original dataset. Approximately 37% of inputs are left out of a particular bootstrap sample and are not used in the construction of the K-th tree. The proof is easy and elegant, it’s done in our Jupyter notebook.

Let’s visualize how Out-of-Bag Error (or OOBE) estimation works:

The top part of the figure above represents our original dataset. We split it into the training (left) and test (right) sets. In the left image, we draw a grid that perfectly divides our dataset according to classes. Now, we use the same grid to estimate the share of the correct answers on our test set. We can see that our classifier gave incorrect answers in those 4 cases that have not been used during training (on the left). Hence, the accuracy of our classifier is 11/15*100% = 73.33%.

To sum up, each base algorithm is trained on ~ 63% of the original examples. It can be validated on the remaining ~37%. The Out-of-Bag estimate is nothing more than the mean estimate of the base algorithms on those ~37% of inputs that were left out of training.

5. Random Forest

Leo Breiman managed to apply bootstrapping not only in statistics but also in machine learning. He, along with Adel Cutler, extended and improved the Random Forest algorithm proposed by Tin Kam Ho. They combined the construction of uncorrelated trees using CART, bagging, and the random subspace method.

Decision trees are a good choice for the base classifier in bagging because they are quite sophisticated and can achieve zero classification error on any sample. The random subspace method reduces the correlation between the trees and thus prevents overfitting. With bagging, the base algorithms are trained on different random subsets of the original feature set.

The following algorithm constructs an ensemble of models using the random subspace method:

1. Let the number of instance be equal to n, and the number of feature dimensions be equal to d.
2. Choose M as the number of individual models in the ensemble.
3. For each model m, choose the number of features dm < d. As a rule, the same value of dm is used for all the models.
4. For each model m, create a training set by selecting dm features at random from the whole set of d features.
5. Train each model.
6. Apply the resulting ensemble model to a new input by combining the results from all the models in M. You can use either majority voting or aggregation of the posterior probabilities.

Outline of part 5

1. Algorithm
2. Comparison with Decision Trees and Bagging
3. Parameters
4. Extremely Randomized Trees
5. Transformation of a dataset into a high-dimensional representation
6. Pros and cons of Random Forests

5.1. Algorithm

The algorithm for constructing a random forest of N trees goes as follows:

The final classifier is defined by:

We use the majority voting for classification and the mean for regression.

For classification problems, it is advisable to set m to be equal the square root of d. For regression problems, we usually take m = d/3, where d is the number of features. It is recommended to build each tree until all of its leaves contain only 1 instance for classification and 5 instances for regression.

You can see Random Forest as bagging of decision trees with the modification of selecting a random subset of features at each split.

5.2. Comparison with Decision Trees and Bagging

See the code in this Jupyter notebook

As we can see from our graphs and the MSE values above, a Random Forest of 10 trees achieves a better result than a single decision tree and is comparable to bagging with 10 trees. The main difference between Random Forests and bagging is that, in a Random Forest, the best feature for a split is selected from a random subset of the available features while, in bagging, all features are considered for the next best split.

We can also look at the advantages of Random Forests and bagging in classification problems.

See the code in this Jupyter notebook

The figures above show that the decision boundary of the decision tree is quite jagged and has a lot of acute angles that suggest overfitting and a weak ability to generalize. We would have trouble making reliable predictions on new test data. In contrast, the bagging algorithm has a rather smooth boundary and has no obvious signs of overfitting.

Now, let’s investigate some parameters which can help us increase the model accuracy.

5.3. Parameters

The scikit-learn library implements random forests by providing two estimators: RandomForestClassifier and RandomForestRegressor.

Below are the parameters which we need to pay attention to when we are building a new model:

• n_estimators is the number of trees in the forest;
• criterion is the function used to measure the quality of a split;
• max_features is the number of features to consider when looking for the best split;
• min_samples_leaf is the minimum number of samples required to be at a leaf node;
• max_depth is the maximum depth of the tree.

Please refer to this Jupyter notebook to find an example of Random Forest tuning in a telecom churn problem. We get the following validation curves.

The most important fact about Random Forests is that it’s accuracy doesn’t decrease when we add trees, so the number of trees is not a compexity hyperparameter as opposed to max_depth and min_samples_leaf. This means that you can tune hyperparameters with, say, 10 trees, and then increase the number of trees till 500 and be sure that accuracy will only get better.

In this notebook you’ll also find more analysis of bias and variance for Random Forests and it’s analogies with k Nearset Neighbors.

5.4. Extremely Randomized Trees

Extremely Randomized Trees employ a greater degree of randomization at the cut-point choice when splitting a tree node. As in Random Forests, a random subset of features is used. But, instead of the search for optimal thresholds, their values are selected at random for each possible feature, and the best one among these randomly generated thresholds is used as the best rule to split the node. This usually trades off a slight reduction in the model variance with a small increase of bias.

In the scikit-learn library, there are 2 implementations of Extremely Randomized Trees: ExtraTreesClassifier and ExtraTreesRegressor.

This method should be used if you have greatly overfit with Random Forests or gradient boosting.

5.5. Transformation of a dataset into a high-dimensional representation

Random forests are mostly used in supervised learning, but there is a way to apply them in the unsupervised setting.

Using the scikit-learn method RandomTreesEmbedding, we can transform our dataset into a high-dimensional, sparse representation. We first build extremely randomized trees and then use the index of the leaf containing the example as a new feature.

For example, if the input appears in the first leaf, we assign 1 as the feature value; if not, we assign 0. This is a so-called binary coding. We can control the number of features and the sparseness of data by increasing or decreasing the number of trees and their depth. Because nearby data points are likely to fall into the same leaf, this transformation provides an implicit nonparametric estimate of their density.

5.6. Pros and cons of Random Forests

Pros:

• High prediction accuracy; will perform better than linear algorithms in most problems; the accuracy is comparable with that of boosting;
• Robust to outliers, thanks to random sampling;
• Insensitive to feature scaling as well as any other monotonic transformations due to the random subspace selection;
• Doesn’t require fine-grained parameter tuning, works quite well out-of-the-box. With tuning, it is possible to achieve a 0.5–3% gain in accuracy, depending on the problem setting and data;
• Efficient for datasets with a large number of features and classes;
• Handles both continuous and discrete variables equally well;
• Rarely overfits. In practice, an increase in the tree number almost always improves the composition. But, after reaching a certain number of trees, the learning curve is very close to the asymptote;
• There are developed methods to estimate feature importance;
• Works well with missing data and maintains good accuracy levels even when a large part of data is missing;
• Provides means to weight classes on the whole dataset as well as for each tree sample;
• Under the hood, calculates proximities between pairs of instamces that can subsequantly be used in clustering, outlier detection, or interesting data representations;
• The above functionality and properties may be extended to unlabeled data to enable unsupervised clustering, data visualization, and outlier detection;
• Easily parallelized and highly scalable.

Cons:

• In comparison with a single decision tree, Random Forest’s output is more difficult to interpret.
• There are no formal p-values for feature importance estimation.
• Performs worse than linear methods in the case of sparse data: text inputs, bag of words, etc.;
• Unlike linear regression, Random Forest is unable to extrapolate. But, this can be also regarded as an advantage because outliers do not cause extreme values in Random Forests;
• Prone to overfitting in some problems, especially, when dealing with noisy data;
• In the case of categorical variables with varying level numbers, Random Forests favor variables with a greater number of levels. The tree will fit more towards a feature with many levels because this gains greater accuracy;
• If a dataset contains groups of correlated features with similar importance for predicted classes, then the preference will be given to smaller groups;
• The resulting model is large and requires a lot of RAM.

6. Feature importance

Pretty often, you would like to make out the exact reasons why the algorithm outputted a particular answer. Or, not being able to understand the algorithm completely, then st lest we’d like to find out which input features contributed the most to the result. With Random Forest, you can obtain such information quite easily.

6.1. Essence of the method

From the picture below, it is intuitively clear that, in our credit scoring problem, Age is much more important than Income. This can be formally explained using the concept of information gain.

In the case of many decision trees (i.e. a Random Forest), the closer the mean position of a feature over all trees to the root, the more important it is for a given classification or regression problem. Gains in the splitting criterion, such as the Gini impurity, obtained at each optimal split in every tree, is a measure of importance that is directly associated with the splitting feature. The value of this score is distinct for each feature and accumulates over all the trees.

We go a little deeper into details in this notebook.

6.2. Practical example

The code can also be found in a notebook.

Let’s consider the results of a survey given to visitors of hostels listed on Booking.com and TripAdvisor.com. Our features here are the average ratings for different categories include service quality, room condition, value for money, etc. Our target variable is the hostel’s overall rating on the website.

The picture above shows that, more often than not, customers pay great attention to staff and the price-quality ratio. This couple of factors affects the resulting overall rating the most. But, the difference between these two features and other festures is not very large. We can therefore conclude that exclusion of any of these features will lead to a reduction in model accuracy. Based on our analysis, we csn recommend hostel owners to focus primarily on staff training and price-to-quality ratio.

7. Assignment #5

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

8. Useful resources

• Chapter 15 of the Elements of Statistical Learning by Jerome H. Friedman, Robert Tibshirani, and Trevor Hastie.
• More about practical applications of random forests and other algorithms can be found in the official documentation of scikit-learn.
• For more in-depth discussion of variance and decorrelation of random forests, see the original paper.

Authors: Yury Kashnitskiy, and Vitaliy Radchenko. Translated and edited by Christina Butsko, Anna Shirshova, Anastasia Manokhina, Egor Polusmak, and Yuanyuan Pao.

Open Machine Learning Course. Topic 5. Bagging and Random Forest was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

Welcome to the 4-th week of our course! Now we will present our most important topic — linear models. If you have your data prepared and want to start training models, then you will most probably first try either linear or logistic regression, depending on your task (regression or classification).

This week’s material covers both theory of linear models and practical aspects of their usage in real-world tasks. There’s going to be a lot of math in this topic, and we won’t even try to render all the formulas on Medium. Instead, we provide a Jupyter Notebook for each part of this article. In the assignment, you’ll beat two simple benchmarks in a Kaggle competition solving a problem of identifying a user based on her session of visited websites.

Article outline

Part1. Regression, nbviewer

• Ordinary Least Squares
• Maximum Likelihood Estimation
• Bias-Variance Decomposition
• Regularization of Linear Regression

Part 2. Linear Classification, nbviewer

• Linear Classifier
• Logistic Regression as a Linear Classifier
• Maximum Likelihood Estimation and Logistic Regression
• L2-Regularization of Logistic Loss

Part 3. An Illustrative Example of Logistic Regression Regularization, nbviewer

Part 4. Where Logistic Regression Is Good and Where It’s Not, nbviewer

• Analysis of IMDB movie reviews
• XOR-Problem

Part 5. Validation and Learning Curves, nbviewer

• How much model complexity is needed?
• How much data is needed?

Part 6. Kaggle Inclass competition “Catch Me If You Can”

• Data Downloading and Transformation
• Sparse Matrices
• Training the first model

Assignment 4

Useful resources

Part 6. Kaggle Inclass competition “Catch Me If You Can”

Competition’s page

We will be solving the intruder detection problem analyzing users’ behavior on the Internet. It is a complicated and interesting problem combining data analysis and behavioral psychology. As an illustration of one of such tasks, Yandex solves the mailbox intruder detection problem based on the user’s behavior patterns. In a nutshell, intruder’s behavior patterns may differ from those of the mailbox owner:

• the breaker may not delete emails right after they are read as the mailbox owner may do;
• the intruder may mark emails and even move the cursor differently;
• etc.

So the intruder could be detected and thrown out from the mailbox forcing the user to authenticate via SMS-code.

Similar approaches are being developed in Google Analytics and described in scientific research papers. You can find more on this topic by searching “Traversal Pattern Mining” and “Sequential Pattern Mining”.

In this competition we are going to solve a similar problem: our algorithm is supposed to analyze the sequence of websites consequently visited by a particular person and predict whether this person is a user named Alice or an intruder (somebody else). As a metric, we will use ROC AUC.

Data Downloading and Transformation

Register on Kaggle, if you have not done it before. Go to the competition page and download the data.

First, load the training and test sets. Then explore the data and perform a couple of simple exercises:

The training dataset contains the following features:

• site1 — id of the first visited website in the session;
• time1 — visiting time for the first website in the session;
• …
• site10 — id of the tenth visited website in the session;
• time10 — visiting time for the tenth website in the session;
• target — target variable, takes the value of 1 for Alice’s sessions, and 0 for the other users’ sessions.

User sessions are chosen in such a way that they are not longer than half an hour and/or contain more than ten websites; i.e. a session is considered as ended either if the user has visited ten websites or if the session has lasted over thirty minutes.

There are some empty values in the table, which means that some sessions contain less than ten websites. Replace empty values with 0 and change column types to integer. Also, load the website dictionary and see how it looks:

# Websites total: 48371

In order to train our first model, we need to prepare the data. First of all, exclude the target variable from the training set. Now both training and test sets have the same number of columns, and we can aggregate them into one dataframe. Thus, all transformations will be performed simultaneously on both the training and test datasets.

On the one hand, it leads to the fact that both of our datasets have one feature space (so you don’t have to worry that you may have forgotten to transform a feature in one of the datasets). On the other hand, the processing time will increase. In case of enormously large sets, it may turn out that it is impossible to transform both datasets simultaneously (and sometimes you have to split your transformations into several stages, separately for the train/test dataset). In our case, we are going to perform all the transformations for the united dataframe at once; and, before training the model or making predictions, we will just use the corresponding part of it.

For the sake of simplicity, we will use only the visited websites in the session (and we will not take into account the timestamp features). The point behind this data selection is: Alice has her favorite sites, and the more often you see these sites in the session, the higher the probability that this is an Alice’s session, and vice versa.

Let’s prepare the data. We will keep only the features site1, site2, … , site10 in the dataframe. Keep in mind that missing values were replaced with zeros. Here is how the first rows of the dataframe look like:

Sessions are the sequences of website indices, and such representation of data is inconvenient for linear methods. According to our hypothesis (Alice has favorite websites), we need to transform this dataframe so that each website has the corresponding feature (column) which value is equal to the number of visits on this website within the session. All of this can be done in two lines:

If you understand what just happened here, then you can skip the next section (perhaps, you can handle logistic regression too?). If not, then let us figure it out.

Sparse Matrices

Let’s estimate how much memory it would require to store our data in the example above. Our united dataframe contains 336 thousand samples of 48 thousand integer features in each. It’s easy to calculate the required amount of memory, roughly:

336K * 48K * 8 bytes = 16M * 8 bytes = 128 GB

Obviously, mere mortals don’t have such volumes of memory (strictly speaking, Python may allow you to create such a matrix, but it would not be easy to do anything with it). An interesting fact is that most of the elements of our matrix are zeros. If we counted non-zero elements, then it would make out about 1.8 million, i.е. slightly more than 10% of all matrix elements. Such a matrix, where most elements are zeros, is called sparse, and the ratio between the number of zero elements and the total number of elements is called the sparseness of the matrix.

To work with such matrices, you can use scipy.sparse library, check the documentation to understand what possible types of sparse matrices are, how to work with them and in which cases their usage is most effective. You can learn how they are arranged, for example, be reading the Wikipedia article. Note that a sparse matrix contains only non-zero elements. Finally, you can get the allocated memory size (significant memory savings are obvious):

1866898 elements * 8 bytes = 14935184 bytes
sparse_matrix_size = 14935184 bytes

Let’s explore how the matrix with the websites has been formed using a mini example. Suppose we have the following table with user sessions:

There are 3 sessions, and no more than 3 websites in each. Users visited four different sites in total (there are numbers from 1 to 4 in the table cells). And let us assume that:

1. vk.com
2. habrahabr.ru
3. yandex.ru
4. ods.ai

If the user has visited less than 3 websites during the session, the last few values will be zero. We want to convert the original dataframe in such a way that each session would have the corresponding row which shows the number of visits on each particular site; i.e. we want to transform the previous table into the following form:

To do this, use the constructor: csr_matrix ((data, indices, indptr)) and create a frequency table (see examples, code and comments on the links above to see how it works). Here, we set all the parameters explicitly for greater clarity:

matrix([[2, 1, 0, 0, 0],
          [0, 2, 0, 1, 0],
          [0, 0, 1, 1, 1]])

As you might have noticed, the number of the columns in the resulting matrix is not four (by the number of different websites), but five. A zero column has been added, which shows on how many sites the session was shorter (in our mini example we took sessions of length 3). This column is excessive and it should be removed from the dataframe.

Another benefit of using sparse matrices is that there are special implementations of both matrix operations and machine learning algorithms for them, which sometimes allows to significantly accelerate operations due to the data structure peculiarities. This applies to logistic regression as well. Now, everything is ready to build our first model.

3. Training the first model

Let’s build our first model, using logistic regression implementation from Sklearn with default parameters. We will use the first 90% of the data for training (the training data set is sorted by time), and the remaining 10% for validation. Let's write a simple function that returns the quality of the model, and then train our first classifier:

0.9198622553850315
CPU times: user 138 ms, sys: 77.1 ms, total: 216 ms
Wall time: 2.74 s

The first model demonstrated the quality of approximately 0.92 ROC AUC on the validation set. Let’s take it as the first baseline and a starting point. To make a prediction on the test set, we need to train the model again on the entire training dataset (until this moment, our model used only part of the data for training), which will increase its generalizing ability:

If you follow these steps and upload the answer to the competition page, then you should get the quality of ROC AUC = 0.91707 on the public leaderboard.

Assignment #4

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

Useful resources

• A nice and concise overview of linear models is given in the book “Deep Learning” (I. Goodfellow, Y. Bengio, and A. Courville).
• Linear models are covered practically in every ML book. We recommend “Pattern Recognition and Machine Learning” (C. Bishop) and “Machine Learning: A Probabilistic Perspective” (K. Murphy).
• If you prefer a thorough overview of linear model from a statistician’s viewpoint, then look at “The elements of statistical learning” (T. Hastie, R. Tibshirani, and J. Friedman).
• The book “Machine Learning in Action” (P. Harrington) will walk you through implementations of classic ML algorithms in pure Python.
• Scikit-learn library. These guys work hard on writing really clear documentation.
• Scipy 2017 scikit-learn tutorial by Alex Gramfort and Andreas Mueller.
• One more ML course with very good materials.
• Implementations of many ML algorithms. Search for linear regression and logistic regression.
• And many others, feel free to share in comments.

Author: Yury Kashnitskiy. Translated and edited by Yuriy Isakov, Christina Butsko, Nerses Bagiyan, Yulia Klimushina, Anastasia Manokhina, Evgeniy Mashkin, Egor Polusmak, and Yuanyuan Pao.

Open Machine Learning Course. Topic 4. Linear Classification and Regression was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

Hi! This is the third article in our series. Today we finally reach machine learning. This is going to be exciting!

Article outline

1. Introduction
2. Decision Tree

• How to Build a Decision Tree
• Tree-building Algorithm
• Other Quality Criteria for Splits in Classification Problems
• How a Decision Tree Works with Numerical Features
• Crucial Tree Parameters
• Class DecisionTreeClassifier in Scikit-learn
• Decision Tree in a Regression Problem

3. Nearest Neighbors Method

• Nearest Neighbors Method in Real Applications
• Class KNeighborsClassifier in Scikit-learn

4. Choosing Model Parameters and Cross-Validation

5. Application Examples and Complex Cases

• Decision trees and nearest neighbors method in a customer churn prediction task
• Complex Case for Decision Trees
• Decision Trees and k-NN in a Task of MNIST Handwritten Digits Recognition
• Complex Case for the Nearest Neighbors Method

6. Pros and Cons of Decision Trees and the Nearest Neighbors Method

7. Assignment #3

8. Useful resources

The following material is better viewed as a Jupyter notebook and can be reproduced locally with Jupyter if you clone the course repository.

1. Introduction

Before we dive into the material for this week’s article, let’s talk about the kind of problem that we are going to solve and its place in the exciting field of machine learning. T. Mitchell’s book Machine Learning (1997) gives a classic, general definition of machine learning as follows:

A computer program is said to learn from experience E with respect to some class of tasks T and performance measure P, if its performance at tasks in T, as measured by P, improves with experience E.

In the various problem settings T, P, and E can refer to completely different things. Some of the most popular tasks T in machine learning are the following:

• classification of an instance to one of the categories based on its features;
• regression — prediction of a numerical target feature based on other features of an instance;
• clustering — identifying partitions of instances based on the features of these instances so that the members within the groups are more similar to each other than those in the other groups;
• anomaly detection — search for instances that are “greatly dissimilar” to the rest of the sample or to some group of instances;
• and so many more.

A good overview is provided in the “Machine Learning basics” chapter of “Deep Learning” (by Ian Goodfellow, Yoshua Bengio, Aaron Courville, 2016).

Experience E refers to data (we can’t go anywhere without it). Machine learning algorithms can be divided into those that are trained in supervised or unsupervised manner. In unsupervised learning tasks, one has a set consisting of instances described by a set of features. In supervised learning problems, there’s also a target variable, which is what we would like to be able to predict, known for each instance in a training set.

Example

Classification and regression are supervised learning problems. For example, as a credit institution, we may want to predict loan defaults based on the data accumulated about our clients. Here, the experience E is the available training data: a set of instances (clients), a collection of features (such as age, salary, type of loan, past loan defaults, etc.) for each, and a target variable (whether they defaulted on the loan). This target variable is just a fact of loan default (1 or 0), so recall that this is a (binary) classification problem. If you were instead predicting by how much time the loan payment is overdue, this would become a regression problem.

Finally, the third term used in the definition of machine learning is a metric of the algorithm’s performance evaluation P. Such metrics differ for various problems and algorithms, and we’ll discuss them as we study new algorithms. For now, we’ll refer to a simple metric for classification algorithms, the proportion of correct answers — accuracy — on the test set.

Let’s take a look at two supervised learning problems: classification and regression.

2. Decision Tree

We begin our overview of classification and regression methods with one of the most popular ones — a decision tree. Decision trees are used in everyday life decisions, not just in machine learning. Flow diagrams are actually visual representations of decision trees. For example, Higher School of Economics publishes information diagrams to make the lives of its employees easier. Here is a snippet of instructions for publishing a paper on the Institution portal.

In terms of machine learning, one can see it as a simple classifier that determines the appropriate form of publication (book, article, chapter of the book, preprint, publication in the “Higher School of Economics and the Media”) based on the content (book, pamphlet, paper), type of journal, original publication type (scientific journal, proceedings), etc.

A decision tree is often a generalization of the experts’ experience, a means of sharing knowledge of a particular process. For example, before the introduction of scalable machine learning algorithms, the credit scoring task in the banking sector was solved by experts. The decision to grant a loan was made on the basis of some intuitively (or empirically) derived rules that could be represented as a decision tree.

In our next case, we solve a binary classification problem (approve/deny a loan) on the grounds of “Age”, “Home-ownership”, “Income” and “Education”.

The decision tree as a machine learning algorithm is essentially the same thing as the diagram shown above; we incorporate a stream of logical rules of the form “feature a value is less than x and feature b value is less than y … => Category 1” into a tree-like data structure. The advantage of this algorithm is that they are easily interpretable. For example, using the above scheme, the bank can explain to the client why they were denied for a loan: e.g the client does not own a house and her income is less than 5,000.

As we’ll see later, many other models, although more accurate, do not have this property and can be regarded as more of a “black box” approach, where it is harder to interpret how the input data was transformed into the output. Due to this “understandability” and similarity to human decision-making (you can easily explain your model to your boss), decision trees have gained immense popularity. C4.5, a representative of this group of classification methods, is even the first in the list of 10 best data mining algorithms (“Top 10 Algorithms in Data Mining”, Knowledge and Information Systems, 2008. PDF).

How to Build a Decision Tree

Earlier, we saw that the decision to grant a loan is made based on age, assets, income, and other variables. But what variable to look at first? Let’s discuss a simple example where all the variables are binary.

Recall the game of “20 Questions”, which is often referenced when introducing decision trees. You’ve probably played this game — one person thinks of a celebrity while the other tries to guess by asking only “Yes” or “No” questions. What question will the guesser ask first? Of course, they will ask the one that narrows down the number of the remaining options the most. Asking “Is it Angelina Jolie?” would, in the case of a negative response, leave all but one celebrity in the realm of possibility. In contrast, asking “Is the celebrity a woman?” would reduce the possibilities to roughly half. That is to say, the “gender” feature separates the celebrity dataset much better than other features like “Angelina Jolie”, “Spanish”, or “loves football.” This reasoning corresponds to the concept of information gain based on entropy.

Entropy

Shannon’s entropy is defined for a system with N possible states as follows:

where Pi is the probability of finding the system in the i-th state. This is a very important concept used in physics, information theory, and other areas. Entropy can be described as the degree of chaos in the system. The higher the entropy, the less ordered the system and vice versa. This will help us formalize “effective data splitting”, which we alluded to in the context of “20 Questions”.

Toy Example

To illustrate how entropy can help us identify good features for building a decision tree, let’s look at a toy example. We will predict the color of the ball based on its position.

There are 9 blue balls and 11 yellow balls. If we randomly pull out a ball, then it will be blue with probability p1 = 9/20 and yellow with probability p2 = 11/20, which gives us an entropy S0 = -9/20 log2(9/20) - 11/20 log2(11/20) ≈ 1. This value by itself may not tell us much, but let’s see how the value changes if we were to break the balls into two groups: with the position less than or equal to 12 and greater than 12.

The left group has 13 balls, 8 blue and 5 yellow. The entropy of this group is S1 = -5/13 log2(5/13) - 8/13 log2(8/13) ≈ 0.96. The right group has 7 balls, 1 blue and 6 yellow. The entropy of the right group is S2 = -1/7 log2(1/7) - 6/7 log2(6/7) ≈ 0.6. As you can see, entropy has decreased in both groups, more so in the right group. Since entropy is, in fact, the degree of chaos (or uncertainty) in the system, the reduction in entropy is called information gain. Formally, the information gain (IG) for a split based on the variable Q(in this example it’s a variable “x ≤ 12”) is defined as

where q is the number of groups after the split, Ni is number of objects from the sample in which variable Q is equal to the i-th value. In our example, our split yielded two groups (q = 2), one with 13 elements (N1 = 13), the other with 7 (N2 = 7). Therefore, we can compute the information gain as

It turns out that dividing the balls into two groups by splitting on “coordinate is less than or equal to 12” gave us a more ordered system. Let’s continue to divide them into groups until the balls in each group are all of the same color.

For the right group, we can easily see that we only need one extra partition using “coordinate less than or equal to 18”. But, for the left group, we need three more. Note that the entropy of a group where all of the balls are the same color is equal to 0 (log2(1) = 0).

We have successfully constructed a decision tree that predicts ball color based on its position. This decision tree may not work well if we add any balls because it has perfectly fit to the training set (initial 20 balls). If we wanted to do well in that case, a tree with fewer “questions” or splits would be more accurate, even if it does not perfectly fit the training set. We will discuss the problem of overfitting later.

Tree-building Algorithm

We can make sure that the tree built in the previous example is optimal: it took only 5 “questions” (conditioned on the variable x) to perfectly fit a decision tree to the training set. Under other split conditions, the resulting tree would be deeper, i.e. take more “questions” to reach an answer.

At the heart of the popular algorithms for decision tree construction, such as ID3 or C4.5, lies the principle of greedy maximization of information gain: at each step, the algorithm chooses the variable that gives the greatest information gain upon splitting. Then the procedure is repeated recursively until the entropy is zero (or some small value to account for overfitting). Different algorithms use different heuristics for “early stopping” or “cut-off” to avoid constructing an overfitted tree.

Other Quality Criteria for Splits in Classification Problems

We discussed how entropy allows us to formalize partitions in a tree. But this is only one heuristic; there exist others.

Maximizing this criterion can be interpreted as the maximization of the number of pairs of objects of the same class that are in the same subtree (not to be confused with the Gini index).

In practice, misclassification error is almost never used, and Gini uncertainty and information gain work similarly.

For binary classification, entropy and Gini uncertainty take the following form:

where (p+ is the probability of an object having a label +).

If we plot these two functions against the argument p+, we will see that the entropy plot is very close to the plot of Gini uncertainty, doubled. Therefore, in practice, these two criteria are almost identical.

Example

Let’s consider fitting a decision tree to some synthetic data. We will generate samples from two classes, both normal distributions but with different means.

Let’s plot the data. Informally, the classification problem in this case is to build some “good” boundary separating the two classes (the red dots from the yellow). Machine learning for this case boils down to choosing a good separating border. A straight line will be too simple while some complex curve snaking by each red dot will be too complex and will lead us to making mistakes on new samples. Intuitively, some smooth boundary, or, at least, a straight line or a hyperplane, would work well on new data.

Let’s try to separate these two classes by training an Sklearn decision tree. We will use max_depth parameter that limits the depth of the tree. Let's visualize the resulting separating boundary.

And how does the tree itself look? We see that the tree “cuts” the space into 8 rectangles, i.e. the tree has 8 leaves. Within each rectangle, the tree will make the prediction according to the majority label of the objects inside it.

How can we “read” such a tree?

In the beginning, there were 200 samples (instances), 100 of each class. The entropy of the initial state was maximal, S=1. Then, the first partition of the samples into 2 groups was made by comparing the value of x2 with 1.211 (find this part of the border in the picture above). With that, the entropy of both left and right groups decreased. The process continues up to depth 3. In this visualization, the more samples of the first class, the darker the orange color of the vertex; the more samples of the second class, the darker the blue. At the beginning, the number of samples from two classes is equal, so the root node of the tree is white.

How a Decision Tree Works with Numerical Features

Suppose we have a numeric feature “Age” that has a lot of unique values. A decision tree will look for the best (according to some criterion of information gain) split by checking binary attributes such as “Age <17”, “Age < 22.87”, and so on. But what if the age range is large? Or what if another quantitative variable, “salary”, can also be “cut” in many ways? There will be too many binary attributes to select from at each step during tree construction. To resolve this problem, heuristics are usually used to limit the number of thresholds to which we compare the quantitative variable.

Let’s consider an example. Suppose we have the following dataset:

We see that the tree used the following 5 values to evaluate by age: 43.5, 19, 22.5, 30, and 32 years. If you look closely, these are exactly the mean values between the ages at which the target class “switches” from 1 to 0 or 0 to 1. To illustrate further, 43.5 is the average of 38 and 49 years; a 38-year-old customer failed to return the loan whereas the 49-year-old did. The tree looks for the values at which the target class switches its value as a threshold for “cutting” a quantitative variable.

Given this information, why do you think it makes no sense here to consider a feature like “Age <17.5”?

Let’s consider a more complex example by adding the “Salary” variable (in the thousands of dollars per year).

If we sort by age, the target class ( “loan default”) switches (from 1 to 0 or vice versa) 5 times. And if we sort by salary, it switches 7 times. How will the tree choose features now? Let’s see.

We see that the tree partitioned by both salary and age. Moreover, the thresholds for feature comparisons are 43.5 and 22.5 years of age and 95k and 30.5k per year. Again, we see that 95 is the average between 88 and 102; the individual with a salary of 88k proved to be “bad” while the one with 102k was “good”. The same goes for 30.5k. That is, only a few values for comparisons by age and salary were searched. Why did the tree choose these features? Because they gave better partitioning (according to Gini uncertainty).

Conclusion: the simplest heuristics for handling numeric features in a decision tree is to sort its values in ascending order and check only those thresholds where the value of the target variable changes.

Furthermore, when there are a lot of numeric features in a dataset, each with a lot of unique values, only the top-N of the thresholds described above are selected, i.e. only use the top-N that give maximum gain. The process is to construct a tree of depth 1, compute the entropy (or Gini uncertainty), and select the best thresholds for comparison.

To illustrate, if we split by “Salary ≤ 34.5”, the left subgroup will have the entropy of 0 (all clients are “bad”), and the right one will have the entropy of 0.954 (3 “bad” and 5 “good”, you can check this yourself as it will be part of the assignment). The information gain is roughly 0.3. If we split by “Salary ≤ 95”, the left subgroup will have an entropy of 0.97 (6 “bad” and 4 “good”), and the right one will have an entropy of 0 (a group containing only one object). The information gain is about 0.11. If we calculate information gain for each partition in that manner, we can select the thresholds for comparison of each numeric feature before the construction of a large tree (using all features).

More examples of numeric feature discretization can be found in posts like this or this. One of the most prominent scientific papers on this subject is “On the handling of continuous-valued attributes in decision tree generation” (UM Fayyad. KB Irani, “Machine Learning”, 1992).

Crucial Tree Parameters

Technically, you can build a decision tree until each leaf has exactly one instance, but this is not common in practice when building a single tree because it will be overfitted, or too tuned to the training set, and will not predict labels for new data well. At the bottom of the tree, at some great depth, there will be partitions on less important features (e.g. whether a client came from Leeds or New York). We can exaggerate this story further and find that all four clients who came to the bank for a loan in green trousers did not return the loan. Even if that were true in training, we do not want our classification model to generate such specific rules.

There are two exceptions where the trees are built to the maximum depth:

• Random Forest (a group of trees) averages the responses from individual trees that are built to the maximum depth (we will talk later on why you should do this)
• Pruning trees. In this approach, the tree is first constructed to the maximum depth. Then, from the bottom up, some nodes of the tree are removed by comparing the quality of the tree with and without that partition (comparison is performed using cross-validation, more on this below).

The picture below is an example of a dividing border built in an overfitted tree.

The most common ways to deal with overfitting in decision trees are as follows:

• artificial limitation of the depth or a minimum number of samples in the leaves: the construction of a tree just stops at some point;
• pruning the tree.

Class DecisionTreeClassifier in Scikit-learn

The main parameters of the sklearn.tree.DecisionTreeClassifier class are:

• max_depth – the maximum depth of the tree;
• max_features – the maximum number of features with which to search for the best partition (this is necessary with a large number of features because it would be "expensive" to search for partitions for all features);
• min_samples_leaf – the minimum number of samples in a leaf. This parameter prevents creating trees where any leaf would have only a few members.

The parameters of the tree need to be set depending on input data, and it is usually done by means of cross-validation, more on this below.

Decision Tree in a Regression Problem

When predicting a numeric variable, the idea of a tree construction remains the same, but the quality criteria changes.

where n is the number of samples in a leaf, Yi is the value of the target variable. Simply put, by minimizing the variance around the mean, we look for features that divide the training set in such a way that the values of the target feature in each leaf are roughly equal.

Example

Let’s generate some data distributed by the function

with some noise. Then we will train a tree on it and show what predictions it makes.

We see that the decision tree approximates the data with a piecewise constant function.

3. Nearest Neighbors Method

The nearest neighbors method (k-Nearest Neighbors, or k-NN) is another very popular classification method that is also sometimes used in regression problems. This, like decision trees, is one of the most comprehensible approaches to classification. The underlying intuition is that you look like your neighbors. More formally, the method follows the compactness hypothesis: if the distance between the examples is measured well enough, then similar examples are much more likely to belong to the same class.

According to the nearest neighbors method, the green ball would be classified as “blue” rather than “red”.

For another example, if you do not know how to tag a Bluetooth-headset on an online listing, you can find 5 similar headsets, and, if 4 of them are tagged as “accessories” and only 1 as “Technology”, then you will also label it under “accessories”.

To classify each sample from the test set, one needs to perform the following operations in order:

1. Calculate the distance to each of the samples in the training set.
2. Select k samples from the training set with the minimal distance to them.
3. The class of the test sample will be the most frequent class among those k nearest neighbors.

The method adapts quite easily for the regression problem: on step 3, it returns not the class, but the number — a mean (or median) of the target variable among neighbors.

A notable feature of this approach is its laziness — calculations are only done during the prediction phase, when a test sample needs to be classified. No model is constructed from the training examples beforehand. In contrast, recall that for decision trees in the first half of this article the tree is constructed based on the training set, and the classification of test cases occurs relatively quickly by traversing through the tree.

Nearest neighbors is a well-studied approach. There exist many important theorems claiming that, on “endless” datasets, it is the optimal method of classification. The authors of the classic book “The Elements of Statistical Learning” consider k-NN to be a theoretically ideal algorithm which usage is only limited by computation power and the curse of dimensionality.

Nearest Neighbors Method in Real Applications

• k-NN can serve as a good starting point (baseline) in some cases;
• In Kaggle competitions, k-NN is often used for the construction of meta-features (i.e. k-NN predictions as input to other models) or for stacking/blending;
• The nearest neighbors method extends to other tasks like recommendation systems. The initial decision could be a recommendation of a product (or service) that is popular among the closest neighbors of the person for whom we want to make a recommendation;
• In practice, on large datasets, approximate methods of search are often used for nearest neighbors. There is a number of open source libraries that implement such algorithms; check out Spotify’s library Annoy.

The quality of classification/regression with k-NN depends on several parameters:

• The number of neighbors k.
• The distance measure between samples (common ones include Hamming, Euclidean, cosine, and Minkowski distances). Note that most of these metrics require data to be scaled. Simply speaking, we do not want the “salary” feature, which is on the order of thousands, to affect the distance more than “age”, which is generally less than 100.
• Weights of neighbors (each neighbor may contribute different weights; for example, the further the sample, the lower the weight).

Class KNeighborsClassifier in Scikit-learn

The main parameters of the class sklearn.neighbors.KNeighborsClassifier are:

• weights: uniform (all weights are equal), distance (the weight is inversely proportional to the distance from the test sample), or any other user-defined function;
• algorithm (optional): brute, ball_tree, KD_tree, or auto. In the first case, the nearest neighbors for each test case are computed by a grid search over the training set. In the second and third cases, the distances between the examples are stored in a tree to accelerate finding nearest neighbors. If you set this parameter to auto, the right way to find the neighbors will be automatically chosen based on the training set.
• leaf_size (optional): threshold for switching to grid search if the algorithm for finding neighbors is BallTree or KDTree;
• metric: minkowski, manhattan, euclidean, chebyshev, or other.

4. Choosing Model Parameters and Cross-Validation

The main task of learning algorithms is to be able to generalize to unseen data. Since we cannot immediately check the model performance on new, incoming data (because we do not know the true values of the target variable yet), it is necessary to sacrifice a small portion of the data to check the quality of the model on it.

This is often done in one of two ways:

• setting aside a part of the dataset (held-out/hold-out set). We reserve a fraction of the training set (typically from 20% to 40%), train the model on the remaining data (60–80% of the original set), and compute performance metrics for the model (e.g accuracy) on the hold-out set.
• cross-validation. The most frequent case here is k-fold cross-validation.

In k-fold cross-validation, the model is trained K times on different (K-1) subsets of the original dataset (in white) and checked on the remaining subset (each time a different one, shown above in orange). We obtain K model quality assessments that are usually averaged to give an overall average quality of classification/regression.

Cross-validation provides a better assessment of the model quality on new data compared to the hold-out set approach. However, cross-validation is computationally expensive when you have a lot of data.

Cross-validation is a very important technique in machine learning and can also be applied in statistics and econometrics. It helps with hyperparameter tuning, model comparison, feature evaluation, etc. More details can be found here (blog post by Sebastian Raschka) or in any classic textbook on machine (statistical) learning.

5. Application Examples and Complex Cases

Decision trees and nearest neighbors method in a customer churn prediction task

Let’s read data into a DataFrame and preprocess it. Store State in a separate Series object for now and remove it from the dataframe. We will train the first model without the State feature, and then we will see if it helps.

Let’s allocate 70% of the set for training (X_train, y_train) and 30% for the hold-out set (X_holdout, y_holdout). The hold-out set will not be involved in tuning the parameters of the models. We'll use it at the end, after tuning, to assess the quality of the resulting model. Let's train 2 models: a decision tree and k-NN. We do not know what parameters are good, so we will assume some random ones: a tree depth of 5 and the number of nearest neighbors equal 10.

Let’s assess prediction quality on our hold-out set with a simple metric — the proportion of correct answers (accuracy). The decision tree did better — percentage of correct answers is about 94% (decision tree) versus 88% (k-NN). Note that this performance is achieved by using random parameters.

Now, let’s identify the parameters for the tree using cross-validation. We’ll tune the maximum depth and the maximum number of features used at each split. Here is the essence of how the GridSearchCV works: for each unique pair of values of max_depth and max_features, compute model performance with 5-fold cross-validation, and then select the best combination of parameters.

Let’s list the best parameters and the corresponding mean accuracy from cross-validation.

Let’s draw the resulting tree. Due to the fact that it is not entirely a toy example (its maximum depth is 6), the picture is not that small, but you can “walk” over the tree if you locally open the corresponding picture downloaded from the course repo.

Now, let’s tune the number of neighbors k for k-NN:

Here, the tree proved to be better than the nearest neighbors algorithm: 94.2%/94.6% accuracy for cross-validation and hold-out respectively. Decision trees perform very well, and even random forest (let’s think of it for now as a bunch of trees that work better together) in this example cannot achieve better performance (95.1%/95.3%) despite being trained for much longer.

Complex Case for Decision Trees

To continue the discussion of the pros and cons of the methods in question, let’s consider a simple classification task, where a tree would perform well but does it in an “overly complicated” manner. Let’s create a set of points on a plane (2 features), each point will be one of two classes (+1 for red, or -1 for yellow). If you look at it as a classification problem, it seems very simple: the classes are separated by a line.

However, the border that the decision tree builds is too complicated; plus the tree itself is very deep. Also, imagine how badly the tree will generalize to the space beyond the 30 x 30 squares that frame the training set.

We got this overly complex construction, although the solution is just a straight line x1 = x2.

The method of one nearest neighbor does better than the tree but is still not as good as a linear classifier (our next topic).

Decision Trees and k-NN in a Task of MNIST Handwritten Digits Recognition

Now let’s have a look at how these 2 algorithms perform on a real-world task. We will use the sklearn built-in dataset on handwritten digits. This task is an example where k-NN works surprisingly well.

Pictures here are 8x8 matrices (intensity of white color for each pixel). Then each such matrix is ​​”unfolded” into a vector of length 64, and we obtain a feature description of an object.

Let’s draw some handwritten digits. We see that they are distinguishable.

Next, let’s do the same experiment as in the previous task, but, this time, let’s change the ranges for tunable parameters.

Let’s select 70% of the dataset for training (X_train, y_train) and 30% for holdout (X_holdout, y_holdout). The holdout set will not participate in model parameters tuning; we will use it at the end to check the quality of the resulting model.

Let’s train a decision tree and k-NN with our random parameters and make predictions on our holdout set. We can see that k-NN did much better, but note that this is with random parameters.

Now let’s tune our model parameters using cross-validation as before, but now we’ll take into account that we have more features than in the previous task: 64.

Let’s see the best parameters combination and the corresponding accuracy from cross-validation:

That has already passed 66% but not quite 97%. k-NN works better on this dataset. In the case of one nearest neighbor, we were able to reach 99% guesses on cross-validation.

Let’s train a random forest on the same dataset, it works better than k-NN on the majority of datasets. But we here have an exception.

You would be right to point out that we have not tuned any RandomForestClassifier parameters here. Even with tuning, the training accuracy doesn’t reach 98% as it did with one nearest neighbor.

The conclusion of this experiment (and general advice): first check simple models on your data: decision tree and nearest neighbors (next time we will also add logistic regression to this list). It might be the case that these methods already work well enough.

Complex Case for the Nearest Neighbors Method

Let’s consider another simple example. In the classification problem, one of the features will just be proportional to the vector of responses, but this won’t help for the nearest neighbors method.

As always, we will look at the accuracy for cross-validation and the hold-out set. Let’s construct curves reflecting the dependence of these quantities on the n_neighbors parameter in the method of nearest neighbors. These curves are called validation curves.

One can see that k-NN with the Euclidean distance does not work well on the problem, even when you vary the number of nearest neighbors over a wide range.

In contrast, the decision tree easily “detects” hidden dependencies in the data despite a restriction on the maximum depth.

In the second example, the tree solved the problem perfectly while k-NN experienced difficulties. However, this is more of a disadvantage of using Euclidian distance than of the method. It did not allow us to reveal that one feature was much better than the others.

6. Pros and Cons of Decision Trees and the Nearest Neighbors Method

Pros and cons of decision trees

Pros:

• Generation of clear human-understandable classification rules, e.g. “if age <25 and is interested in motorcycles, deny the loan”. This property is called interpretability of the model.
• Decision trees can be easily visualized, i.e. both the model itself (the tree) and prediction for a certain test object (a path in the tree) can “be interpreted”.
• Fast training and forecasting.
• Small number of model parameters.
• Supports both numerical and categorical features.

Cons:

• The trees are very sensitive to the noise in input data; the whole model could change if the training set is slightly modified (e.g. remove a feature, add some objects). This impairs the interpretability of the model.
• Separating border built by a decision tree has its limitations — it consists of hyperplanes perpendicular to one of the coordinate axes, which is inferior in quality to some other methods, in practice.
• We need to avoid overfitting by pruning, setting a minimum number of samples in each leaf, or defining a maximum depth for the tree. Note that overfitting is an issue for all machine learning methods.
• Instability. Small changes to the data can significantly change the decision tree. This problem is tackled with decision tree ensembles (discussed next time).
• The optimal decision tree search problem is NP-complete. Some heuristics are used in practice such as greedy search for a feature with maximum information gain, but it does not guarantee finding the globally optimal tree.
• Difficulties to support missing values in the data. Friedman estimated that it took about 50% of the code to support gaps in data in CART (an improved version of this algorithm is implemented in sklearn).
• The model can only interpolate but not extrapolate (the same is true for random forests and tree boosting). That is, a decision tree makes constant prediction for the objects that lie beyond the bounding box set by the training set in the feature space. In our example with the yellow and blue balls, it would mean that the model gives the same predictions for all balls with positions >19 or <0.

Pros and cons of the nearest neighbors method

Pros:

• Simple implementation.
• Well studied.
• Typically, the method is a good first solution not only for classification or regression, but also recommendations.
• It can be adapted to a certain problem by choosing the right metrics or kernel (in a nutshell, the kernel may set the similarity operation for complex objects such as graphs while keeping the k-NN approach the same). By the way, Alexander Dyakonov, a former top-1 kaggler, loves the simplest k-NN but with the tuned object similarity metric.
• Good interpretability. There are exceptions: if the number of neighbors is large, the interpretability deteriorates (“We did not give him a loan, because he is similar to the 350 clients, of which 70 are the bad, and that is 12% higher than the average for the dataset”).

Cons:

• Method considered fast in comparison with compositions of algorithms, but the number of neighbors used for classification is usually large (100–150) in real life, in which case the algorithm will not operate as fast as a decision tree.
• If a dataset has many variables, it is difficult to find the right weights and to determine which features are not important for classification/regression.
• Dependency on the selected distance metric between the objects. Selecting the Euclidean distance by default is often unfounded. You can find a good solution by grid searching over parameters, but this becomes very time consuming for large datasets.
• There are no theoretical ways to choose the number of neighbors — only grid search (though this is often true for all hyperparameters of all models). In the case of a small number of neighbors, the method is sensitive to outliers, that is, it is inclined to overfit.
• As a rule, it does not work well when there are a lot of features due to the “curse of dimensionality”. Professor Pedro Domingos, a well-known member in the ML community, talks about it here in his popular paper, “A Few Useful Things to Know about Machine Learning”; also “the curse of dimensionality” is described in the Deep Learning book in this chapter.

This is a lot of information, but, hopefully, this article will be a great reference for you for a long time :)

7. Assignment #3

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

8. Useful resources

• Decision trees and k Nearest Neighbors are covered practically in every ML book. We recommend “Pattern Recognition and Machine Learning” (C. Bishop) and “Machine Learning: A Probabilistic Perspective” (K. Murphy).
• The book “Machine Learning in Action” (P. Harrington) will walk you through implementations of classic ML algorithms in pure Python.
• Scikit-learn library. These guys work hard on writing really clear documentation.
• Scipy 2017 scikit-learn tutorial by Alex Gramfort and Andreas Mueller.
• One more ML course with very good materials.
• Implementations of many ML algorithms. Good to search for decision trees and k-NN.
• And many others, feel free to share in comments.

Author: Yury Kashnitskiy. Translated and edited by Christina Butsko, Gleb Filatov, and Yuanyuan Pao.

Open Machine Learning Course. Topic 3. Classification, Decision Trees and k Nearest Neighbors was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.

In the field of Machine Learning, data visualization is not just making fancy graphics for reports; it is used extensively in day-to-day work for all phases of a project.

To start with, visual exploration of data is the first thing one tends to do when dealing with a new task. We do preliminary checks and analysis using graphics and tables to summarize the data and leave out the less important details. It is much more convenient for us, humans, to grasp the main points this way than by reading many lines of raw data. It is amazing how much insight can be gained from seemingly simple charts created with available visualization tools.

Next, when we analyze the performance of a model or report results, we also often use charts and images. Sometimes, for interpreting a complex model, we need to project high-dimensional spaces onto more visually intelligible 2D or 3D figures.

All in all, visualization is a relatively fast way to learn something new about your data. Thus, it is vital to learn its most useful techniques and make them part of your everyday ML toolbox.

In this article, we are going to get hands-on experience with visual exploration of data using popular libraries such as pandas, matplotlib and seaborn.

Article outline

1. Dataset
2. Univariate visualization: quantitative and categorical distributions
3. Multivariate visualization: interactions between variables
4. Whole dataset: peeking into high-dimensional spaces
5. Assignment #2
6. Useful resources

The following material is better viewed as a Jupyter notebook and can be reproduced locally with Jupyter if you clone the course repository.

1. Dataset

Before we get to the data, let’s initialize our environment:

import numpy as np
import pandas as pd
pd.options.display.max_columns = 12

# Disable warnings in Anaconda
import warnings
warnings.simplefilter('ignore')

# We will display plots right inside Jupyter Notebook
%matplotlib inline
import matplotlib.pyplot as plt

# We will use the Seaborn library
import seaborn as sns
sns.set()

# Graphics in SVG format are more sharp and legible
%config InlineBackend.figure_format = 'svg'

# Increase the default plot size
from pylab import rcParams
rcParams['figure.figsize'] = 5, 4

In the first article, we looked at the data on customer churn for a telecom operator. We will load again that dataset into a DataFrame:

df = pd.read_csv('../../data/telecom_churn.csv')

To get acquainted with our data, let’s look at the first 5 entries using head():

df.head()

Here is the description of our features:

The last data column, Churn, is our target variable. It is binary: True indicates that that the company eventually lost this customer, and False indicates that the customer was retained. Later, we will build models that predict this variable based on the remaining features. This is why we call it a target.

2. Univariate visualization

Univariate analysis looks at one variable at a time. When we analyze a feature independently, we are usually mostly interested in the distribution of its values and ignore the other variables in the dataset.

Below, we will consider different statistical types of variables and the corresponding tools for their individual visual analysis.

2.1 Quantitative features

Quantitative features take on ordered numerical values. Those values can be discrete, like integers, or continuous, like real numbers, and usually express a count or a measurement.

Histograms and density plots

The easiest way to take a look at the distribution of a numerical variable is to plot its histogram using the DataFrame's method hist().

features = ['Total day minutes', 'Total intl calls']
df[features].hist(figsize=(12, 4));

A histogram groups values into bins of equal value range. The shape of the histogram may contain clues about the underlying distribution type: Gaussian, exponential etc. You can also spot any skewness in its shape when the distribution is nearly regular but has some anomalies. Knowing the distribution of the feature values becomes important when you use Machine Learning methods that assume a particular type of it, most often Gaussian.

In the above plot, we see that the variable Total day minutes is normally distributed, while Total intl calls is prominently skewed right (its tail is longer on the right).

There is also another, often clearer, way to grasp the distribution: density plots or, more formally, Kernel Density Plots. They can be considered a smoothed version of the histogram. Their main advantage over the latter is that they do not depend on the size of the bins. Let’s create density plots for the same two variables:

df[features].plot(kind='density', subplots=True, 
                  layout=(1, 2), sharex=False, figsize=(12, 4));

It is also possible to plot a distribution of observations with seaborn's distplot(). For example, let's look at the distribution of Total day minutes. By default, the plot displays both the histogram with the kernel density estimate (KDE) on top.

sns.distplot(df['Total intl calls']);

The height of the histogram bars here is normed and shows the density rather than the number of examples in each bin.

Box plot

Another useful type of visualization is a box plot. seaborn does a great job here:

_, ax = plt.subplots(figsize=(3, 4))

sns.boxplot(data=df['Total intl calls'], ax=ax);

Let’s see how to interpret a box plot. Its components are a box (obviously, this is why it is called a box plot), the so-called whiskers, and a number of individual points (outliers).

The box by itself illustrates the interquartile spread of the distribution; its length is determined by the 25th (Q1) and 75th (Q3) percentiles. The vertical line inside the box marks the median (50%) of the distribution.

The whiskers are the lines extending from the box. They represent the entire scatter of data points, specifically the points that fall within the interval (Q1−1.5⋅IQR, Q3+1.5⋅IQR), where IQR=Q3−Q1 is the interquartile range.

Outliers that fall out of the range bounded by the whiskers are plotted individually as black points along the central axis.

We can see that a large number of international calls is quite rare in our data.

Violin plot

The last type of distribution plots that we will consider is a violin plot.

Look at the figures below. On the left, we see the already familiar box plot. To the right, there is a violin plot with the kernel density estimate on both sides.

_, axes = plt.subplots(1, 2, sharey=True, figsize=(6, 4))

sns.boxplot(data=df['Total intl calls'], ax=axes[0]); sns.violinplot(data=df['Total intl calls'], ax=axes[1]);

The difference between the box and violin plots is that the former illustrates certain statistics concerning individual examples in a dataset while the violin plot concentrates more on the smoothed distribution as a whole.

In our case, the violin plot does not contribute any additional information about the data as everything is clear from the box plot alone.

describe()

In addition to graphical tools, in order to get the exact numerical statistics of the distribution, we can use the method describe() of a DataFrame:

df[features].describe()

Its output is mostly self-explanatory. 25%, 50% and 75% are the corresponding percentiles.

2.2 Categorical and binary features

Categorical features take on a fixed number of values. Each of these values assigns an observation to a corresponding group, known as a category, which reflects some qualitative property of this example. Binary variables are an important special case of categorical variables when the number of possible values is exactly 2. If the values of a categorical variable are ordered, it is called ordinal.

Frequency table

Let’s check the class balance in our dataset by looking at the distribution of the target variable: the churn rate. First, we will get a frequency table, which shows how frequent each value of the categorical variable is. For this, we will use the value_counts() method:

df['Churn'].value_counts()

False    2850
True      483
Name: Churn, dtype: int64

By default, the entries in the output are sorted from the most to the least frequently-occurring values.

In our case, the data is not balanced; that is, our two target classes, loyal and disloyal customers, are not represented equally in the dataset. Only a small part of the clients canceled their subscription to the telecom service. As we will see in the following articles, this fact may imply some restrictions on measuring the classification performance, and, in the future, we may want to additionaly penalize our model errors in predicting the minority “Churn” class.

Bar plot

The bar plot is a graphical representation of the frequency table. The easiest way to create it is to use the seaborn's function countplot(). There is another function in seaborn that is somewhat confusingly called barplot() and is mostly used for representation of some basic statistics of a numerical variable grouped by a categorical feature.

Let’s plot the distributions for two categorical variables:

_, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 4))

sns.countplot(x='Churn', data=df, ax=axes[0]);
sns.countplot(x='Customer service calls', data=df, ax=axes[1]);

While the histograms, discussed above, and bar plots may look similar, there are several differences between them:

1. Histograms are best suited for looking at the distribution of numerical variables while bar plots are used for categorical features.
2. The values on the X-axis in the histogram are numerical; a bar plot can have any type of values on the X-axis: numbers, strings, booleans.
3. The histogram’s X-axis is a Cartesian coordinate axis along which values cannot be changed; the ordering of the bars is not predefined. Still, it is useful to note that the bars are often sorted by height, that is, the frequency of the values. Also, when we consider ordinal variables (like Customer service calls in our data), the bars are usually ordered by variable value.

The left chart above vividly illustrates the imbalance in our target variable. The bar plot for Customer service calls on the right gives a hint that the majority of customers resolve their problems in maximum 2–3 calls. But, as we want to be able to predict the minority class, we may be more interested in how the fewer dissatisfied customers behave. It may well be that the tail of that bar plot contains most of our churn. These are just hypotheses for now, so let’s move on to some more interesting and powerful visual techniques.

3. Multivariate visualization

Multivariate plots allow us to see relationships between two and more different variables, all in one figure. Just as in the case of univariate plots, the specific type of visualization will depend on the types of the variables being analyzed.

3.1 Quantitative–Quantitative

We are going to start with the interaction between quantitative variables.

Correlation matrix

Let’s look at the correlations among the numerical variables in our dataset. This information is important to know as there are Machine Learning algorithms (for example, linear and logistic regression) that do not handle highly correlated input variables well.

First, we will use the method corr() on a DataFrame that calculates the correlation between each pair of features. Then, we pass the resulting correlation matrix to heatmap() from seaborn, which renders a color-coded matrix for the provided values:

# Drop non-numerical variables
numerical = list(set(df.columns) - 
                 set(['State', 'International plan', 
                      'Voice mail plan', 'Area code', 'Churn',
                      'Customer service calls']))

# Calculate and plot
corr_matrix = df[numerical].corr()
sns.heatmap(corr_matrix);

From the colored correlation matrix generated above, we can see that there are 4 variables such as Total day charge that have been calculated directly from the number of minutes spent on phone calls (Total day minutes). These are called dependent variables and can therefore be left out since they do not contribute any additional information. Let’s get rid of them:

numerical = list(set(numerical) - 
                 set(['Total day charge', 'Total eve charge', 
                      'Total night charge', 'Total intl charge']))

Scatter plot

The scatter plot displays values of two numerical variables as Cartesian coordinates in 2D space. Scatter plots in 3D are also possible.

Let’s try out the function scatter() from the matplotlib library:

plt.scatter(df['Total day minutes'], df['Total night minutes']);

We get an uninteresting picture of two normally distributed variables. Also, it seems that these features are uncorrelated because the ellpise-like shape is aligned with the axes.

There is a slightly fancier option to create a scatter plot with the seaborn library:

sns.jointplot(x='Total day minutes', y='Total night minutes',
              data=df, kind='scatter');

The function jointplot() plots two histograms that may be useful in some cases. Using the same function, we can also get a smoothed version of our bivariate distribution:

sns.jointplot('Total day minutes', 'Total night minutes',
              data=df, kind="kde", color="g");

This is basically a bivariate version of the Kernel Density Plot discussed earlier.

Scatterplot matrix

In some cases, we may want to plot a scatterplot matrix such as the one shown below. Its diagonal contains the distributions of the corresponding variables, and the scatter plots for each pair of variables fill the rest of the matrix.

# pairplot may become very slow with the SVG format
%config InlineBackend.figure_format = 'png' sns.pairplot(df[numerical]);

Sometimes, such visualization may help draw conclusions about data; but, in this case, everything is pretty clear with no surprises.

3.2 Quantitative–Categorical

In this section, we will make our simple quantitative plots a little more exciting. We will try to gain new insights for churn prediction from the interactions between the numerical and categorical features.

More specifically, let’s see how the input variables are related to the target variable Churn.

Previously, you learned about scatter plots. Additionally, their points can be color or size coded so that the values of a third categorical variable are also presented in the same figure. We can achieve this with the scatter() function seen above, but, let's try a new function called lmplot() and use the parameter hue to indicate our categorical feature of interest:

sns.lmplot('Total day minutes', 'Total night minutes', data=df, 
           hue='Churn', fit_reg=False);

It seems that our small proportion of disloyal customers lean towards the top-right corner; that is, such customers tend to spend more time on the phone during both day and night. But this is not absolutely clear, and we won’t make any definitive conclusions from this chart.

Now, let’s create box plots to visualize the distribution statistics of the numerical variables in two disjoint groups: the loyal customers (Churn=False) and those who left (Churn=True).

# Sometimes you can analyze an ordinal variable as numerical one
numerical.append('Customer service calls') 

fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(10, 7))

for idx, feat in enumerate(numerical): 
    ax = axes[int(idx / 4), idx % 4] 
    sns.boxplot(x='Churn', y=feat, data=df, ax=ax) 
    ax.set_xlabel('') 
    ax.set_ylabel(feat) 
fig.tight_layout();

From this chart, we can see that the greatest discrepancy in distribution between the two groups is for three variables: Total day minutes, Customer service calls, and Number vmail messages. Later in this course, we will learn how to determine feature importance in classification using Random Forest or Gradient Boosting; there, we will see that the first two features are indeed very important for churn prediction.

Let’s look at the distribution of day minutes spoken for the loyal and disloyal customers separately. We will create box and violin plots for Total day minutes grouped by the target variable.

_, axes = plt.subplots(1, 2, sharey=True, figsize=(10, 4)) 

sns.boxplot(x='Churn', y='Total day minutes',
            data=df, ax=axes[0]);
sns.violinplot(x='Churn', y='Total day minutes', 
               data=df, ax=axes[1]);

In this case, the violin plot does not contribute any additional information about our data as everything is clear from the box plot alone: disloyal customers tend to talk on the phone more.

An interesting observation: on average, customers that discontinue their contracts are more active users of communication services. Perhaps they are unhappy with the tariffs, so a possible measure to prevent churn could be a reduction in call rates. The company will need to undertake additional economic analysis to find out whether such measures would be beneficial.

When we want to analyze a quantitative variable in two categorical dimensions at once, there is a suitable function for this in the seaborn library called factorplot(). For example, let's visualize the interaction between Total day minutes and two categorical variables in the same plot:

sns.factorplot(x='Churn', y='Total day minutes',
               col='Customer service calls',
               data=df[df['Customer service calls'] < 8], 
               kind="box", col_wrap=4, size=3, aspect=.8);

From this, we could conclude that, starting with 4 calls, Total day minutes may no longer be the main factor for customer churn. Perhaps, in addition to our previous guess about the tariffs, there are customers that are dissatisfied with the service due to other problems, which might lead to fewer number of day minutes spent on calls.

3.3 Categorical–Categorical

As we saw earlier in this article, the variable Customer service calls has few unique values and, thus, can be considered either numerical or ordinal. We have already seen its distribution with a count plot. Now, we are interested in the relationship between this ordinal feature and the target variable Churn.

Let’s look at the distribution of the number of calls to the customer service, again using a count plot. This time, let’s also pass the parameter hue=Churn that adds a categorical dimension to the plot:

sns.countplot(x='Customer service calls', hue='Churn', data=df);

An observation: the churn rate increases significantly after 4 or more calls to the customer service.

Now, let’s look at the relationship between Churn and the binary features, International plan and Voice mail plan.

_, axes = plt.subplots(1, 2, sharey=True, figsize=(10, 4)) 

sns.countplot(x='International plan', hue='Churn',
              data=df, ax=axes[0]);
sns.countplot(x='Voice mail plan', hue='Churn',
              data=df, ax=axes[1]);

An observation: when International Plan is enabled, the churn rate is much higher; the usage of the international plan by the customer is a strong feature. We do not observe the same effect with Voice mail plan.

Contingency table

In addition to using graphical means for categorical analysis, there is a traditional tool from statistics: a contingency table, also called a cross tabulation. It represents multivariate frequency distribution of categorical variables in tabular form. In particular, it allows us to see the distribution of one variable conditional on the other by looking along a column or row.

Let’s try to see how Churn is related to the categorical variable State by creating a cross tabulation:

pd.crosstab(df['State'], df['Churn']).T

In the case of State, the number of distinct values is rather high: 51. We see that there are only a few data points available for each individual state — only 3 to 17 customers in each state abandoned the operator. Let’s ignore that for a second and calculate the churn rate for each state, sorting it from high to low:

df.groupby(['State'])['Churn'].
   agg([np.mean]).
   sort_values(by='mean', ascending=False).T

At first glance, it seems that the churn rate in New Jersey and California are above 25% and less than 6% for Hawaii and Alaska. However, these conclusions are based on too few examples, and our observation could be a mere property of our particular dataset. We can confirm this with the Matthews and Cramer correlation hypotheses, but this would be beyond the scope of this article.

4. Whole dataset

4.1 Naive approach

We have been looking at different facets of our dataset by guessing interesting features and selecting a small number of them at a time for visualization. We have only delt with two to three variables at once and were easily able to observe the structure and relationships in data. But, what if we want to display all the features and still be able to interpret the resulting visualization?

We could use hist() or create a scatterplot matrix with pairplot() for the whole dataset to look at all of our features simultaneously. But, when the number of features is high enough, this kind of visual analysis quickly becomes slow and inefficient. Besides, we would still be analyzing our variables in a pairwise fashion, not all at once.

4.2 Dimensionality reduction

Most real-world datasets have many features, sometimes, many thousands of them. Each of them can be considered as a dimension in the space of data points. Consequently, more often than not, we deal with high-dimensional datasets, where entire visualization is quite hard.

To look at a dataset as a whole, we need to decrease the number of dimensions used in visualization without losing much information about data. This task is called dimensionality reduction and is an example of an unsupervised learning problem because we need to derive new, low-dimensional features from the data itself, without any supervised input.

One of the well-known dimensionality reduction methods is Principal Component Analysis (PCA), which we will study later in this course. Its limitation is that it is a linear algorithm that implies certain restrictions on the data.

There are also many non-linear methods, collectively called Manifold Learning. One of the best-known of them is t-SNE.

4.3 t-SNE

Let’s create a t-SNE representation of the same churn data we have been using.

The name of the method looks complex and a bit intimidating: t-distributed Stohastic Neighbor Embedding. Its math is also impressive (we will not delve into it here, but, if you feel brave, here is the original article by Laurens van der Maaten and Geoffrey Hinton from JMLR). Its basic idea is simple: to find a projection for a high-dimensional feature space onto a plane (or a 3D hyperplane, but it is almost always 2D) such that those points that were far apart in the initial n-dimensional space will end up far apart on the plane. Those that were originally close would remain close to each other.

Essentially, neighbor embedding is a search for a new and less-dimensional data representation that preserves neighborship of examples.

Now, let’s do some practice. First, we need to import some additional classes:

from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler

We will leave out the State and Churn features and convert the values “Yes”/”No” of the binary features into numerical values using pandas.Series.map():

X = df.drop(['Churn', 'State'], axis=1)
X['International plan'] = X['International plan'].
                          map({'Yes': 1, 'No': 0})
X['Voice mail plan'] = X['Voice mail plan'].
                       map({'Yes': 1, 'No': 0})

We also need to normalize the data. For this, we will subtract the mean from each variable and divide it by its standard deviation. All of this can be done with StandardScaler.

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

Now, let’s build a t-SNE representation:

%%time tsne = TSNE(random_state=17)
tsne_repr = tsne.fit_transform(X_scaled)

CPU times: user 1min 32s, sys: 7.86 s, total: 1min 39s
Wall time: 1min 39s

and plot it:

plt.scatter(tsne_repr[:, 0], tsne_repr[:, 1]);

Let’s color this t-SNE representation according to the churn (green for loyal customers, and red for those who left).

plt.scatter(tsne_repr[:, 0], tsne_repr[:, 1], 
            c=df['Churn'].map({False: 'green', True: 'red'}));

We can see that customers who churned are concentrated in a few areas of the lower dimensional feature space.

To better understand the picture, we can also color it with the remaining binary features: International Plan and Voicemail. The green dots here indicate the objects that are positive for the corresponding binary feature.

_, axes = plt.subplots(1, 2, sharey=True, figsize=(12, 5))

for i, name in enumerate(['International plan', 'Voice mail plan']):
    axes[i].scatter(tsne_repr[:, 0], tsne_repr[:, 1],
                    c=df[name].map({'Yes': 'green', 'No': 'red'}))
    axes[i].set_title(name)

Now it is clear that, for example, many dissatisfied customers who canceled their subscription are crowded together in the most south-west cluster that represents the people with the international plan but no voice mail.

Finally, let’s note some disadvantages of t-SNE:

• High computational complexity. The implementation in scikit-learn is unlikely to be feasible in a real task. If you have a large number of samples, you should try Multicore-TSNE instead.
• The plot can change a great deal depending on the random seed, which complicates interpretation. Here is a good tutorial on t-SNE. In general, you shouldn’t make any far-reaching conclusions based on such graphs because it can equate to plain guessing. Of course, some findings in t-SNE pictures can inspire an idea and be confirmed through more thorough research down the line, but that does not happen very often.

Occasionally, using t-SNE, you can get a really good intuition for the data. The following is a good paper that shows an example of this for handwritten digits: Visualizing MNIST.

Sometimes t-SNE really helps to understand something, and sometimes you can just draw a Christmas tree toy :-)

5. Assignment #2

Full versions of assignments are announced each week in a new run of the course (October 1, 2018). Meanwhile, you can practice with a demo version: Kaggle Kernel, nbviewer.

6. Useful resources

• Here is the official documentation for the libraries we used: matplotlib, seaborn and pandas.
• The gallery of sample charts created with seaborn is a very good resource.
• Also, see the documentation on Manifold Learning in scikit-learn
• Efficient t-SNE implementation Multicore-TSNE.
• “How to Use t-SNE Effectively”, Distill.pub.

Authors: Egor Polusmak, and Yury Kashnitskiy. Translated and edited by Yuanyuan Pao, Christina Butsko, Anastasia Manokhina, Inga Kaydanova, and Artem Trunov.

Open Machine Learning Course. Topic 2. Visual Data Analysis with Python was originally published in Open Machine Learning Course on Medium, where people are continuing the conversation by highlighting and responding to this story.