Logistic regression

A shot with the foot from the penalty spot, not a penalty itself but a shot in open play with defenders nearby, becomes a goal about one time in five. From twenty-five metres, about one time in thirty. Where do these numbers come from? In the 2017/18 season the five biggest football leagues saw 40,461 shots in open play, and for each it is known where it was taken from and how it ended. But from any single spot on the pitch only a handful of shots were taken, and the share of goals “by spot” cannot be counted. Logistic regression reduces the whole pitch to four numbers: it adds up the shot’s features with weights and turns the sum into a probability. Below is a model trained on the autumn half of the season. Tap the pitch to “shoot” from there.

The thin lines join points with the same probability of a goal. They are not circles around the goal: from a tight angle goals come less often than from the same distance in the centre, and the model knows this from the feature “the angle at which the goal is seen”. Choose the model with distance alone, and the lines become arcs of circles centred on the middle of the goal: it can have no other shape. Press “learn from zero”: Newton’s method starts from zero weights, at which the probability is 50% everywhere, and reaches the answer in eight steps, and after the fifth the change can no longer be seen. The ring around the shot is a check on live data: at the penalty spot the model gives a shot with the foot 21.9%, and of the 513 foot shots of the whole season within two metres of it, 114 went in, 22.2%. The chart under the pitch checks the model on the spring, which it has not seen: shots it gave 5 to 10%, 7.4% on average, became goals 7.2% of the time in spring.

The idea

A shot has two outcomes, goal and no goal; let us write them as 1 and 0. The first thing that comes to mind is linear regression: draw a line through these zeros and ones by distance. The least squares line over the autumn shots: “probability” = 0.2805 − 0.00957·d, where d is the distance in metres. The trouble is visible at once: from 29.3 metres the “probability” is negative, and 11% of shots were taken from that far. Near the goal it is worse: at two metres the line promises 26%, while from three metres or closer 76% of shots went in. A probability is squeezed between 0 and 1 and changes with distance not at all along a line: first it falls steeply, then it slowly fades away.

Logistic regression unsqueezes this interval. Instead of the probability p it looks at the odds p / (1 − p), a probability of 0.2 corresponding to odds of 1 to 4, and at the odds through a logarithm. The logit ln(p / (1 − p)) runs over all numbers from −∞ to +∞, and it can be set equal to a weighted sum of features without fear of leaving the bounds. Back to a probability leads the sigmoid:

Logistic regression predicts not the outcome but the logarithm of its odds, as a weighted sum of features. The sigmoid turns the sum into a probability between 0 and 1. The weights are chosen so that the outcomes that actually happened are as probable as possible: the same as making the log loss smallest.

Our model has three features: the distance to the middle of the goal in metres, the angle at which the goal is seen (in radians), and a header (1 or 0). The trained model is:

logit = -1.1172 -0.1066·metres +1.5816·angle -0.9216·header
odds of a goal: x0.899 a metre farther, x0.398 for a header

The weights are read through odds. Each metre adds −0.1066 to the logit, that is, multiplies the odds by e^−0.1066 ≈ 0.899: ten per cent less per metre, anywhere on the pitch. A header multiplies the odds by 0.398: compared with a foot shot from the same spot, a header’s odds are two and a half times lower. This does not mean heading is harmful: headers come after crosses, in a challenge, at a ball in flight, and the feature “header” stands for all of this at once. The model describes shots as they happen, not what would have happened had the same player in the same situation shot differently.

Where does the model answer “goal”, and where “no goal”? Choose a threshold, say 20%: where p ≥ 0.2, the answer is “goal”. The condition p ≥ 0.2 is equivalent to z ≥ logit 0.2, and z is a linear function of the features, so the decision boundary is a hyperplane in feature space: in the axes “distance, angle” it is a line. On the pitch it is curved, because the features themselves are nonlinear functions of the shot’s position. That is how a model linear in its weights draws curved boundaries: it all depends on which columns it is given.

How to measure a miss

A model with weights w assigns each shot a probability of a goal p_i. How well does it agree with what actually happened? If shots are independent, the probability of seeing exactly this season under the model is the product of p_i over all goals and 1 − p_i over all misses. Good weights make this product large. Multiplying twenty thousand numbers below one is awkward, as the product quickly becomes smaller than any double, so we take the logarithm, flip the sign and divide by the number of shots:

where q_i is the probability the model gave to what happened: p_i for a goal and 1 − p_i for a miss. This is the log loss, also known as cross-entropy. The term −ln q_i is a measure of the model’s surprise at the outcome, introduced in the chapter on information in “Mathematics, the Queen of the Sciences”. Compare it with the squared loss (p − y)², familiar from linear regression, on a shot that became a goal:

The squared loss is bounded by one: it makes no difference to the model whether it called the goal that happened unlikely or impossible. The log loss grows without limit: confidence that was not warranted has to be paid for. But the main difference shows when the weights are sought by descent. The derivative of a shot’s log loss with respect to its logit z is simply p − y: how far the model missed. For the squared loss it is 2(p − y)·p(1 − p), and the factor p(1 − p) is almost zero exactly where the model is confident. A confident and wrong model hardly gets a signal that it should change.

Watch it live. Here the model is simpler, one curve by distance, logit = b + w·(distance − 20 m)/10 m, and it has only two weights, so all possible models fit on a plane. The colour is the relief of the loss, the crosses the bottom of each loss. From the starting point gradient descent with the same step goes over both losses; under the plane are the curves obtained on the way, over the real share of goals at each metre.

The start is deliberately bad: b = −6, w = 6, a model by which the farther the shot, the more dangerous it is: right at the goal there are no goals at all, and from 30 metres every other shot goes in. On the log loss the descent reaches the bottom in 39 steps. On the squared loss it first rolls quickly to the model “there are no goals”: a probability near zero at any distance. There its loss equals the share of goals, 0.1019, and the gradient almost vanishes, because p(1 − p) ≈ 0 everywhere. On this plateau the descent crawls for thousands of steps and reaches the bottom only on step 7896. Tap the plane in other places: on the log loss the descent arrives quickly from everywhere, on the squared loss it varies, depending on how far the start is from the plateau.

There is a second difference: the crosses stand in different places. These are two different answers to the question “which curve is better”. The squared loss draws a steeper curve (w = −1.92 against −1.48). The log loss answer is the maximum likelihood method, and it has a property that the next section is about.

The squared loss after a sigmoid cannot be convex, and nothing even needs computing for this: it is enclosed between 0 and 1 on the whole plane of weights, and a convex function bounded above on the whole plane must be constant. So it has plateaus and slopes, and on other data false pits too: in 1995 Peter Auer, Mark Herbster and Manfred Warmuth built examples where a single neuron with a sigmoid and squared loss has exponentially many local minima in the number of features. The least squares line from the start of the article and the descent from two starts on both losses are in one program:

Why its probabilities can be trusted

The log loss has a property that the line and many other models lack. At the minimum all derivatives are zero, and the derivative of the loss with respect to the weight w_j is Σ_i (p_i − y_i)·x_ij. So for a trained model

for every column j. Take the column of ones, for the intercept w₀: the sum of probabilities over all shots equals the number of goals. Take the column “header”: the sum of probabilities over headers equals the number of headed goals. The model cannot systematically overestimate or underestimate, either overall or in a group singled out by its feature. This is neither luck nor tuning but a consequence of the condition for the minimum:

autumn: expected goals 2094.0000, scored 2094

Such agreement of predicted probabilities with frequencies is called calibration. It is checked on data the model has not seen: the predictions are laid out in bins, and in each the mean prediction is compared with the share of goals. Here is the spring, 19,912 shots:

predicted   shots  goals  expected  scored
0.00-0.05    7103    190      2.8%    2.7%
0.05-0.10    5664    410      7.4%    7.2%
0.10-0.20    4822    725     13.8%   15.0%
0.20-0.30    1276    348     24.1%   27.3%
0.30-0.50     700    283     38.2%   40.4%
0.50-1.00     347    221     64.7%   63.7%

Shots to which the model gave 2.8% on average became goals 2.7% of the time, and shots with a prediction of 64.7% 63.7% of the time. In the middle the spring is a little luckier than predicted: the model expected 2081.8 goals, and 2177 were scored. This does not break the property, which holds on the data the model learned from; it is a reminder that the spring need not resemble the autumn: in spring 10.9% of shots went in against 10.2% in autumn. How to correct the model for a new season with one number is the subject of the article’s second problem.

Calibration is what makes logistic regression a workhorse where the probability itself matters, not only a “yes or no” answer. In search advertising the probability of a click on an ad, together with the advertiser’s bid, decides what to show and in what order; the engineers who described such a system in 2013 write plainly that without accurate, calibrated probabilities the auction does not work. In 2005 Alexandru Niculescu-Mizil and Rich Caruana compared the probabilities of ten learning methods on eight problems: logistic regression, neural networks and bagged trees turned out well calibrated from the start, while naive Bayes, because of its assumption that features are independent, pushed probabilities towards 0 and 1. And for methods whose scores are poorly calibrated there is a simple cure, Platt’s method: a logistic regression with one feature, the output of the original classifier.

The condition for the minimum promises agreement overall and for each column, but not in every bin. In the first scene choose the model by distance alone: shots it gives 30 to 50%, 34.8% on average, were scored 47.8% of the time in spring. A single logistic curve by distance cannot be both steep right at the goal and shallow far away, and near the goal it underestimates. Calibration in bins checks whether the model has enough features: with the angle and header, the bin furthest from the diagonal differs from reality by three percentage points, 24.1 against 27.3%.

Code

The minimum of the log loss cannot be written as a formula, but the function is convex and smooth, and Newton’s method finds it in a few steps. A Newton step is the minimum of a quadratic approximation of the loss: you need the gradient g = Σ(p_i − y_i)x_i and the matrix of second derivatives H = Σp_i(1 − p_i)x_ix_i^T, then the solution of the system Hδ = g; for four features it is a 4 × 4 system, solved by Gaussian elimination. Note that these are the normal equations of a weighted least squares problem with weights p_i(1 − p_i). So in statistics the method is called iteratively reweighted least squares, IRLS: each step is a linear regression with its own row weights.

Three details in the listing matter more than they seem. The sigmoid is computed by two formulas, for z ≥ 0 and z < 0, otherwise exp(-z) overflows at a large negative z. The loss is computed through softplus(z) = ln(1 + e^z) = max(z, 0) + ln(1 + e^−|z|), not through the logarithm of a computed probability: at z = 40 the probability 1 − σ(z) is already zero in double, while the honest loss is 40. The first problem is about this. And the third: there are 20,549 autumn shots, and 2664 distinct combinations of “spot, header”. Shots from one spot have the same features, so a row with n shots and g goals can be taken into account once with a weight: it contributes (n·p − g)·x to the gradient and n·p(1 − p)·xx^T to the matrix. Almost eight times less work, and the same answer.

The output starts with the course of training:

step 1: log loss 0.693147, largest change 2.0e+00
step 2: log loss 0.316471, largest change 4.4e-01
step 3: log loss 0.283604, largest change 6.4e-01
step 4: log loss 0.278656, largest change 3.3e-01
step 5: log loss 0.278264, largest change 3.4e-02
step 6: log loss 0.278264, largest change 3.5e-04
step 7: log loss 0.278264, largest change 3.8e-08
step 8: log loss 0.278264, largest change below 1e-12

The first line is the zero weights: a probability of 1/2 for every shot and a loss of ln 2 = 0.693. The last four steps show why Newton’s method is loved: the change goes 3.4·10⁻² → 3.5·10⁻⁴ → 3.8·10⁻⁸ → below 10⁻¹². The number of correct digits doubles with each step; this is quadratic convergence, the same as in Newton’s method for the roots of an equation: here it looks for a root of the gradient. At the end the listing compares models with different columns:

columns                     autumn  spring
distance                    0.2872  0.3030
angle                       0.2892  0.3053
distance + angle            0.2834  0.2995
distance + angle + header   0.2783  0.2922

The model “10.2% for any shot” gets a loss of 0.3454 on the spring. Distance alone brings it down to 0.3030, the angle adds a little more, the header noticeably more. On the autumn, where the model learned, the losses are lower than on the spring for all four models: the check has to be on data the model has not seen.

What it costs

A Newton step for n rows and d features costs O(n·d²): each row adds all d² pairwise products of its features to the matrix H. Another O(d³) goes on solving the system. The memory needed is O(d²) for the matrix, and the rows themselves can be read as a stream: only g and H build up over a pass. There are usually fewer than ten steps, because the convergence is quadratic. A gradient descent step is cheaper, O(n·d), as only the gradient is needed, but orders of magnitude more steps are required, and their number depends on how elongated the loss landscape is.

The slopes on logarithmic axes are about one for the descent step and about two for the Newton step: with eight features a Newton step costs several times more, with a hundred and twenty-eight almost a hundred times. For the xG model with four features it does not matter: eight Newton steps over 2664 rows take a fraction of a second even in Python. But in advertising and in texts there are millions of features, every word and every pair of “query, ad”, and a d × d matrix fits in no memory. There methods that do without it are used: L-BFGS remembers the last few steps and how the gradient changed over them, and approximates the curvature from them (this is scikit-learn’s default solver), while stochastic descent goes over one row or a batch of rows at a time. The advertising system from the 2013 paper learns exactly like this, online, on the stream of impressions, with a variant of descent called FTRL-Proximal, with a separate step size for each feature.

When there is no best answer

Try training the model by distance on a single match. One of the spring matches had 18 shots: two goals, from 4.7 and 6.3 metres, and all sixteen misses from 12.9 metres and farther. A curve that falls from one to zero somewhere between 6.3 and 12.9 metres explains this match perfectly, and the steeper its drop, the smaller the loss. The loss reaches zero at no finite weights; it only tends to it as the slope grows to infinity. Newton will step for ever, and the weights will keep growing.

This situation is called separation: there is a hyperplane with only goals on one side and only misses on the other. In 1984 Albert and Anderson showed that there is no finite minimum under quasi-complete separation either, when both goals and misses lie on the hyperplane itself, and that the minimum exists and is unique only when goals and misses overlap. With one feature it is seen from four numbers: the nearest and farthest goal, the nearest and farthest miss. Of the 1826 matches of the season, 330 have no best curve by distance: 184 matches had no goals, in 132 all goals were scored closer than all misses, in seven farther, and in seven more the farthest goal and the closest miss were taken from the same distance. The article’s third problem is about this.

Over a whole season four weights are in no such danger. But give the model many features, and separation comes back piece by piece. Cut the pitch into 4 × 4 metre cells and give each cell its own weight: logit = b + w_c. It is the same logistic regression, only the features are indicators “a shot from cell c”. In a cell from which nobody scored in autumn, the best weight is minus infinity: the model is sure that nobody scores from there. The very first spring goal from such a cell costs it an infinite loss.

The cure is an L2 penalty: λ/2·Σw_c² is added to the loss. Now each cell has a reason not to stray far from the common b, and the stronger the fewer shots it has: a cell with two shots has nothing to outweigh the penalty with, while a cell with a thousand has. The penalised loss is strictly convex and grows in all directions, so its minimum always exists, as long as the shots include both goals and misses.

Move λ. At λ = 0 the cells are just shares of goals, and red frames mark the “impossible” cells from which somebody scored in spring: with a hundred autumn matches there are 35 of them. A small penalty already saves the model, and λ around 0.5 gives the best loss on the spring. Penalise too hard, and all the cells merge into one constant. Switch the number of matches: the more data, the less the penalty matters and the lower the whole curve.

first 100 matches: 2263 shots in 107 cells, 58 cells without a goal; spring scores in 35 of them
  lambda 0: inf | 0.01: 0.3329 | 0.1: 0.3193 | 0.5: 0.3143 | 2: 0.3169 | 10: 0.3300 | 100: 0.3444 | 1000: 0.3468
first 909 matches: 20549 shots in 133 cells, 49 cells without a goal; spring scores in 10 of them
  lambda 0: inf | 0.01: 0.3029 | 0.1: 0.3021 | 0.5: 0.3017 | 2: 0.3019 | 10: 0.3057 | 100: 0.3273 | 1000: 0.3427

And the main lesson of this scene: even on the whole autumn and with the best λ, a hundred and thirty-three cells get 0.3017 on the spring, while the formula of distance and angle gets 0.2995 with three numbers. The formula manages with the first hundred matches too: trained on them, it gets 0.3004; the scene shows this as a dashed line. Knowledge of how a shot works turned out to be the best penalty possible: it keeps the model from fitting itself to chance. The penalty is chosen on held-out data, not on the data used for training: on training data a smaller penalty always looks better.

What the model does not know

• The boundary is a plane in the features. Logistic regression separates answers with a hyperplane, and anything more complex has to be built into the features. The model by distance alone draws only circles on the pitch; the angle, a feature devised by a person, gave it the shape of the real lines. If the contribution of one feature depends on another (a header from far away is more hopeless than from close by), a product column is needed. When there is nobody to devise features, models that build them themselves are used: trees and neural networks.
• Not everything is visible in the data. The model has no defenders, no goalkeeper and no idea how the ball arrived. Penalties are taken from the same spot from which, by our model, a foot shot in open play goes in 21.9% of the time, while in the same Wyscout data 393 of the 541 penalties in the five leagues were scored, 73%. The model averages everything it does not see over the situations that were in the training data, and for a different situation its answer is wrong.
• Association is not cause. The weight of the “header” says that headers go in less often, but not what would have happened had the player used the foot: those would have been different shots.
• More than two answers. If there are several outcomes, which of ten digits is written, softmax is used instead of the sigmoid: each class has its own weighted sum, and the probabilities are proportional to their exponentials. The loss is the same, −ln of the probability of the correct answer. This is how the output layer of the network that reads digits in the course chapter on machine learning works; logistic regression is such a network without hidden layers.

History

The curve itself is older than statistics. In 1838 Pierre François Verhulst published a short note on how a population grows: at first almost like an exponential, then ever more slowly, running up against a limit. In an 1845 paper he called this curve logistic and drew it next to the logarithmic one, as the exponential was then called. The work was forgotten, and in 1920 Raymond Pearl and Lowell Reed came to the same curve anew, also studying population growth.

The curve came into the statistics of probabilities from biological experiments: what share of test insects or animals dies at a given dose of poison. In 1934 Chester Bliss proposed the probit for this, a model with the normal distribution. In 1944 Joseph Berkson proposed taking the logistic function instead and, by analogy with the probit, called it the logit. The logit was considered simpler than the probit, but Berkson also argued against the maximum likelihood method, and the dispute with Ronald Fisher and his followers stretched over decades. Regression for a binary outcome with several features was described by David Cox in 1958, and in 1972 John Nelder and Robert Wedderburn included it in the general scheme of generalised linear models and described the solution method that R and statsmodels still use, iteratively reweighted least squares. The variant with several outcomes was linked to choice theory by Daniel McFadden, who received the Nobel Prize in economics in 2000 for his work on discrete choice.

Football shots were modelled with logistic regression as early as 1997: Richard Pollard and Charles Reep estimated the probability of a goal from the distance, angle and circumstances of the shot on 410 attacks ending in a shot, in 22 matches of the 1986 World Cup. Today such models are called expected goals, xG, and they are often shown in match statistics next to the score.

Problems

The first problem is about how to compute the log loss honestly when the model is very confident. The second is about calibration: how to correct a model trained on last season with one number. The third is about separation: when logistic regression has an answer at all.