github linkedin email
Process of Data science - Numerical analysis
May 21, 2021
18 minutes read

Numerical analysis for error reduction

In the previous post, we picked a model class for a hospital readmission problem. The task was to predict which discharged patients are at high risk of returning to the hospital within thirty days, and we walked through why a particular model class fit the shape of the data and the tradeoff between bias and variance that comes with any modeling choice. Picking a model class is only half the job. Once you have chosen, say, a logistic regression or a gradient boosted tree for the readmission problem, there is still a lot of numerical work left to do before that model is trustworthy enough to sit in front of a clinician.

This post covers that numerical work. It is not about picking a model family, it is about tuning the model you already picked so that its errors are as small as possible, and just as important, so that the errors it does make are the right kind of errors for the problem at hand. We will look at gradient based optimization as the mechanical process that reduces error, why the loss function you optimize matters as much as the optimizer itself, regularization as a numerical technique that trades training error for generalization, and calibration, which asks whether the probabilities coming out of the model actually mean what they claim to mean.

  1. Hypothesis
  2. Measurement variables
  3. Latent or unobservable factors
  4. Experimental design (0 to 1)
    1. Controlling other factors to observe primary effect.
  5. Collection and analysis of data for pattern discovery
    1. Hypothesis driven Exploration
  6. Modeling of patterns for prediction
    1. Numerical Analysis for error reduction (discussed here)
    2. Qualitative modeling
  7. Generalizing or scaling the experiment (1 to n)
  8. Establishing a baseline
  9. Monitoring through controls and baselines
  10. Ethics and governance

What numerical analysis means here

Once a model class is fixed, that model is defined by a set of parameters. A logistic regression has a coefficient for every measurement variable and an intercept. A gradient boosted tree has the structure and leaf values of every tree in the ensemble. Numerical analysis, in the sense used in this post, is the set of techniques for choosing values of those parameters that make the model as accurate as possible, and for shaping what "as accurate as possible" even means for the problem at hand.

This is a different activity from the one covered in the next post in this series, qualitative modeling, which is about the structural choices a data scientist makes: which variables to include, how to encode categorical variables, whether an interaction term belongs in the model, whether a nonlinear transform of a variable is warranted. Qualitative modeling shapes the hypothesis space the model can express. Numerical analysis searches inside whatever hypothesis space the model class defines and finds a good point in it. Both matter, and in practice a data scientist moves back and forth between them, but they are worth separating conceptually because the techniques and the failure modes are different.

Four ideas do most of the work in numerical analysis for error reduction: how you search for good parameters (gradient based optimization), what you are actually minimizing while you search (the loss function), how you keep the search from overfitting to the training data (regularization), and how you verify that the resulting probabilities mean what they say (calibration). We will take each in turn and then tie them together in the hospital readmission example.

Gradient based optimization, intuitively

Most modern model classes, from logistic regression to deep neural networks, are fit using some variant of gradient based optimization. The intuition does not require calculus to grasp. Picture the loss, a single number that measures how wrong the model currently is on the training data, as the elevation on a landscape, where the position on the landscape is given by the current values of the model's parameters. Training a model is walking downhill on that landscape, one small step at a time, until you reach a low point.

At any given position, the gradient tells you the direction of steepest ascent. Taking a step in the opposite direction, scaled by a small learning rate, is gradient descent. Formally, if $\theta$ is the vector of model parameters and $L(\theta)$ is the loss computed on training data, each step of gradient descent moves the parameters according to

$$ \theta_{t+1} = \theta_t - \eta \nabla L(\theta_t) $$

where $\eta$ is the learning rate and $\nabla L(\theta_t)$ is the gradient of the loss with respect to the parameters, evaluated at the current parameter values. Repeat this update many times and, under reasonable conditions, the parameters converge to a point where the gradient is close to zero and the loss cannot be reduced much further by a small local move.

A few practical points tend to matter more than the formula itself.

The learning rate is a tradeoff

A learning rate that is too large overshoots the low point and can make the loss oscillate or diverge. A learning rate that is too small converges reliably but slowly, and can get stuck making imperceptible progress on a nearly flat stretch of the landscape long before it should stop. Most practical training procedures decay the learning rate over time, taking large exploratory steps early and small refining steps later.

The landscape is rarely a single bowl

For a plain logistic regression the loss landscape is convex, meaning it has a single low point and gradient descent is guaranteed to find it. For gradient boosted trees and neural networks the landscape can have many local low points, flat plateaus and narrow ravines. In practice this matters less than it sounds, because in high dimensional parameter spaces most local low points found by gradient based methods tend to have similar loss to each other, and techniques like momentum (carrying some of the previous step's direction forward) and stochastic sampling of training batches help the search avoid getting trapped in small dents in the landscape.

Stochastic rather than full batch

Computing the exact gradient on the entire training set at every step is expensive and, for the hospital readmission problem with a modest number of records, may not even be necessary. Stochastic gradient descent computes the gradient on a small random subset of the data at each step. This is noisier, but the noise itself has a mild regularizing effect, which is a preview of the next idea.

The important takeaway is that gradient based optimization is a mechanical, general purpose search procedure. It will faithfully find a low point of whatever loss function you hand it. It has no opinion about whether that loss function reflects what actually matters for the problem. That responsibility falls entirely on the person choosing the loss.

Why the loss function choice matters

A loss function turns a prediction and a true outcome into a single number representing how bad that prediction was. The most common loss for a binary classification problem like readmission or no readmission is log loss, also called cross entropy, defined for a single patient as

$$ \ell(y, \hat{p}) = -\Big[ y \log(\hat{p}) + (1 - y) \log(1 - \hat{p}) \Big] $$

where $y$ is $1$ if the patient was actually readmitted and $0$ otherwise, and $\hat{p}$ is the model's predicted probability of readmission. This loss treats a false negative, predicting a low probability for a patient who is readmitted, and a false positive, predicting a high probability for a patient who is not readmitted, as symmetric mistakes of equal severity once you account for the predicted probability involved. That symmetry is a modeling choice baked into the formula, not a law of nature, and it is worth questioning every time.

For hospital readmission risk, the symmetry is wrong. Missing a patient who will be readmitted, a false negative, means that patient leaves the hospital without the extra follow up call, the medication reconciliation visit or the home health referral that might have kept them out of the hospital. The cost is a possible readmission: more suffering for the patient and a real financial cost to the hospital, which under many reimbursement models is penalized for preventable readmissions. A false alarm, flagging a patient as high risk who was never going to be readmitted, costs a follow up call or a visit that turns out to be unnecessary. That is wasted staff time, not a missed opportunity to prevent harm. The two kinds of error are not remotely equivalent, and a loss function that treats them as equivalent will produce a model whose decision threshold is calibrated to the wrong tradeoff.

The numerical fix is to weight the two terms of the loss asymmetrically. A common approach is a weighted log loss

$$ \ell_w(y, \hat{p}) = -\Big[ w_1 \, y \log(\hat{p}) + w_0 \, (1 - y) \log(1 - \hat{p}) \Big] $$

with $w_1 > w_0$ so that missing a true readmission is penalized more heavily during training than raising a false alarm. Choosing $w_1$ and $w_0$ is itself a modeling decision, ideally informed by the actual costs involved: the expected cost of a missed readmission (patient harm plus any financial penalty) versus the expected cost of an unnecessary intervention (staff time for a follow up call). Getting these weights from real cost estimates, rather than picking a round number, is the difference between a loss function that reflects the business problem and one that reflects a default in a software library.

This idea generalizes well beyond weighted log loss. Any loss function can be shaped to reflect the true cost structure of a problem. For a model predicting a continuous quantity like length of stay, an asymmetric loss can penalize underprediction more than overprediction if running out of beds is worse than having spare capacity. For a demand forecasting model, the reverse might be true. The lesson is general: the loss function is not an implementation detail to accept by default, it is where the real world cost of being wrong gets encoded into the numbers that gradient descent will faithfully minimize.

Regularization: trading training error for generalization

The previous post introduced the bias and variance tradeoff: a model with too little flexibility underfits and has high bias, a model with too much flexibility overfits and has high variance, memorizing quirks of the training data that do not generalize. Regularization is the numerical technique for pulling a model back from the high variance end of that tradeoff without changing the model class itself.

Mechanically, regularization adds a penalty term to the loss function that grows with the size or complexity of the model's parameters. For a logistic regression, two common choices are L2 regularization (also called ridge), which penalizes the sum of squared coefficients, and L1 regularization (also called lasso), which penalizes the sum of absolute coefficients and tends to push some coefficients exactly to zero, effectively performing a kind of variable selection. The regularized loss looks like

$$ L_{reg}(\theta) = \frac{1}{n}\sum_{i=1}^{n} \ell_w(y_i, \hat{p}_i) \; + \; \lambda \sum_{j} \theta_j^2 $$

for L2 regularization, where $\lambda$ controls how strongly complexity is penalized. When $\lambda$ is zero, the model is free to fit the training data as closely as gradient descent allows, which for a small dataset with many measurement variables usually means memorizing noise. As $\lambda$ increases, the model is pushed toward smaller, smoother coefficients, which almost always increases the loss measured on the training data itself.

This last point is worth sitting with because it seems backwards at first. Regularization deliberately makes the model worse on the data it was trained on. The justification is that training error is not what we actually care about, generalization error, the error on new patients the model has never seen, is what matters. A model with no regularization can have very low training error and much higher error on new data, because some of what it learned was noise specific to the training sample rather than signal that generalizes. A modest amount of regularization raises training error slightly while lowering the gap between training error and new data error by more than that, producing a net improvement on the data that matters.

Gradient boosted trees and neural networks have their own analogues of this idea: limiting tree depth, shrinking the contribution of each tree in an ensemble, early stopping (halting gradient descent before it fully converges on the training set), and dropout in neural networks all serve the same purpose as L1 or L2 penalties, restraining the model's capacity to memorize the training sample. The right amount of regularization is almost never known in advance. It is chosen numerically, typically by trying a range of values for $\lambda$ and picking the one that minimizes error on a validation set the model was not trained on, rather than on the training set itself.

Calibration: do predicted probabilities mean what they say

A model can discriminate well, meaning it ranks high risk patients above low risk patients most of the time, while still producing probability values that are wrong in an absolute sense. Calibration asks a narrower and, for many applications, more important question than discrimination: among all the patients the model assigns a predicted probability of thirty percent, roughly thirty percent should actually go on to be readmitted. If the true rate among those patients is closer to fifteen percent or fifty percent, the model is miscalibrated, even if it still correctly ranks patients relative to each other.

Calibration matters enormously for a clinical decision tool because clinicians and hospital administrators do not just use the ranking, they use the number. A probability is often fed into a downstream decision rule: patients above a certain predicted risk get a home health referral, a follow up call within forty eight hours, or a slot in a transitional care program. If a model systematically reports thirty percent when the true rate is fifty percent, every one of those downstream decisions is miscalibrated in a way that is invisible unless someone specifically checks for it. Two models with identical discrimination, measured for instance by an area under the ROC curve, can have very different calibration, and only one of them is safe to hand a raw probability threshold to.

Miscalibration commonly appears after regularization or after training on class imbalanced data, both of which are common for a readmission problem where most patients are not readmitted. Regularization, by pulling coefficients toward zero, tends to pull predicted probabilities toward the middle of the range, a phenomenon sometimes described as the model being underconfident. Training on a rebalanced dataset, where readmitted patients were oversampled to help the model learn the minority class, shifts predicted probabilities away from the true base rate in the opposite direction.

Checking calibration is done empirically, using data the model was not trained on. A common tool is the reliability diagram: bucket predictions into ranges (zero to ten percent, ten to twenty percent, and so on), compute the average predicted probability and the observed readmission rate within each bucket, and plot one against the other. A perfectly calibrated model produces points that fall on the diagonal line where predicted equals observed. Systematic deviation from the diagonal, for instance if every bucket's observed rate sits above the predicted rate, indicates the model is underconfident and its probabilities need adjustment. A single summary number for calibration and discrimination together is the Brier score, the mean squared difference between predicted probability and the actual binary outcome, though a reliability diagram is more informative because it shows where along the probability range the miscalibration occurs.

When a model is discriminating well but miscalibrated, the usual fix is not to retrain the whole model but to apply a small correction on top of its output, fit on a held out set. Platt scaling fits a simple logistic function to remap the model's raw scores into calibrated probabilities. Isotonic regression fits a more flexible, monotonic remapping when there is enough held out data to support it. Both techniques leave the model's ranking of patients unchanged and only adjust the numeric probability values so that they match observed frequencies.

Worked example: tuning the hospital readmission model

Picking up the hospital readmission model from the previous post, assume the model class chosen was a regularized logistic regression over a modest set of measurement variables: prior admission count, discharge diagnosis category, length of stay, age, and a handful of lab values at discharge. The dataset is one hospital's worth of discharges over a few years, small by machine learning standards, perhaps a few thousand records with only a few hundred readmissions.

Step one: choosing an asymmetric loss

The hospital's clinical and finance teams estimate, in a targeted interview much like the ones described earlier in this series, that a missed readmission costs roughly four times what an unnecessary follow up call costs, once penalty exposure and patient harm are weighed against staff time. Rather than training on plain log loss, the weighted version from earlier is used with $w_1 = 4$ and $w_0 = 1$, so a false negative on a true readmission contributes four times as much to the loss gradient as a false positive of the same predicted probability magnitude. This single change shifts the model's effective decision threshold: it will tolerate more false alarms in exchange for catching more of the true readmissions, which matches what the hospital actually wants from the tool.

Step two: adding regularization for a small dataset

With only a few hundred positive examples and several dozen candidate measurement variables after encoding diagnosis categories, an unregularized logistic regression trained by gradient descent will happily drive some coefficients to large values chasing patterns that are specific to this hospital's small sample: perhaps a particular combination of diagnosis code and discharge day of week that appears correlated with readmission purely by chance in a few hundred cases. Adding L2 regularization, with the penalty strength $\lambda$ chosen by trying a range of values and selecting the one that minimizes weighted log loss on a validation split held out from training, shrinks the coefficients toward zero and away from those chance patterns. Training loss goes up slightly compared to the unregularized fit, exactly as expected, but validation loss goes down, which is the signal that the model will generalize better to next month's discharges than the unregularized version would have.

Step three: checking calibration before deployment

Before the model goes anywhere near a clinician, its predictions on a held out set, patients not used for training or for choosing $\lambda$, are bucketed into deciles of predicted probability and compared against observed readmission rates in each bucket. Suppose the reliability diagram shows that patients predicted at around sixty percent risk are actually readmitted only about forty percent of the time, an overconfident model in the high risk range, which is precisely the range the hospital cares most about since it drives the most expensive interventions. A Platt scaling correction is fit on that held out set, remapping raw model scores into calibrated probabilities without changing which patients rank above which other patients. After recalibration the reliability diagram sits much closer to the diagonal across the full range of predicted probabilities, and the resulting number, "this patient has roughly a forty percent chance of readmission," is now a number a clinician can trust in something close to its literal sense, rather than a rough proxy for relative risk.

None of these three steps required changing the model class. The logistic regression from the previous post is still a logistic regression. What changed is entirely numerical: which loss the gradient descent procedure was minimizing, how strongly the fitted coefficients were penalized for complexity, and a small correction applied to the raw output. Together they take a model that discriminates reasonably well and turn it into one whose errors are weighted the way the hospital actually experiences them, whose fit is less likely to be an artifact of a small sample, and whose probabilities can be read at face value.

Common pitfalls

A few mistakes recur often enough to call out directly. Optimizing the default loss function a library ships with, usually plain unweighted log loss or mean squared error, without asking whether it reflects the real asymmetry of the problem's errors. Tuning regularization strength by looking at training error instead of a held out validation or cross validation error, which defeats the entire purpose of regularization. Treating a high area under the ROC curve as proof that a model is ready for deployment without ever checking a reliability diagram, since discrimination and calibration are different properties and a model can have one without the other. And retraining an entire model from scratch to fix a calibration problem that a simple Platt scaling or isotonic regression correction on held out data would have fixed in a few lines of code.

All four ideas in this post, the optimizer, the loss function, the regularization penalty and the calibration check, are numerical knobs sitting on top of whatever model class was chosen. The next post in this series turns to the other half of modeling: the qualitative and structural choices about which variables to include and how to represent them, choices that shape the hypothesis space this post's numerical techniques search inside of.

References

[1] Platt, J. 1999. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in Large Margin Classifiers. 10, 3 (1999), 61 to 74.

[2] Niculescu-Mizil, A. and Caruana, R. 2005. Predicting good probabilities with supervised learning. Proceedings of the 22nd International Conference on Machine Learning. (2005), 625 to 632.

[3] Hastie, T., Tibshirani, R. and Friedman, J. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer.


Back to posts


comments powered by Disqus