Robust, Quantile, and Count Regression

Choosing a regression family begins with the decision, not with the fact that the target is numeric. One column of numbers can carry a mean question, a median question, a tail question and a rate question. The loss is what decides which of them gets answered.

A server that usually answers in 10ms can still be the reason a service is slow. Say its 99th-percentile latency is one second. On its own it makes only one user request in 100 slow. Then it is asked to work in parallel, and in 2013 Jeffrey Dean and Luiz André Barroso did that arithmetic: “If a user request must collect responses from 100 such servers in parallel, then 63% of user requests will take more than one second (marked “x” in the figure).” Push the slow responses down to one in 10,000 and a 2,000-server service still sees almost one in five requests exceed one second.

Their Table 1, measured on a real Google fan-out service, shows one percentile answering three different questions. The 99th-percentile latency is 10ms for one random leaf request, 70ms for 95% of all leaf requests, and 140ms for 100% of all leaf requests (50th percentiles 1ms/12ms/40ms; 95th percentiles 5ms/32ms/87ms). A model fitted to mean latency would have reported none of that.

The fitted function is defined partly by the objective, even when the features stay unchanged. Squared loss targets the conditional mean and responds strongly to large residuals. Absolute loss targets a conditional median with a linear residual penalty, nondifferentiable at zero but optimizable. Pinball loss targets a chosen conditional quantile through asymmetric error weights, and produces a separate function for every quantile. Poisson deviance targets a nonnegative mean under a count-oriented likelihood, commonly through a log link with an exposure offset and an assumed mean–variance relationship.

Pinball loss is not a classroom device. The M5 “Uncertainty” competition made it the operational scoring rule at scale. Entrants had to forecast 42,840 hierarchical unit-sales time series from Walmart, and not as point predictions. The organisers' 2022 report on the competition states the requirement: “To do so, the competition required the prediction of nine different quantiles (0.005, 0.025, 0.165, 0.250, 0.500, 0.750, 0.835, 0.975, and 0.995), that can sufficiently describe the complete distributions of future sales.” Nine quantiles per series means nine fitted functions per series, scored with a scaled pinball loss so that series of different magnitudes could be compared. The choice of loss was the choice of what a forecast even is.

Huber loss behaves quadratically for small residuals and linearly beyond a threshold. Peter J. Huber introduced the M-estimator behind it in 1964. The threshold is not left abstract in practice either. R's MASS::rlm defaults psi.huber's tuning constant to k = 1.345, statsmodels' HuberT defaults t = 1.345, and scikit-learn's HuberRegressor defaults epsilon = 1.35. The three most-used implementations put the quadratic-to-linear crossover at roughly 1.35 standard-scaled residuals. scikit-learn writes the split out — squared loss where |(y - Xw - c) / sigma| < epsilon, absolute loss where it exceeds epsilon — on the reasoning that “The Huber loss function has the advantage of not being heavily influenced by the outliers while not completely ignoring their effect.”

RANSAC takes the other route. It repeatedly fits candidate models from subsets and searches for a consensus set.

Neither method decides whether an extreme row is an error, a rare valid case, or the population that matters most, so robustness choices must match the operational objective. A default of 1.345 is a decision about which rows stop counting fully. It was made by someone who never saw your data.

Several quantiles can summarize asymmetry and spread without assuming one constant error variance.

The family has a single origin and a dated one. Roger Koenker and Gilbert Bassett set it out in Econometrica in January 1978, opening with the generalisation: “A simple minimization problem yielding the ordinary sample quantiles in the location model is shown to generalize naturally to the linear model generating a new class of statistics we term ‘regression quantiles.’” Least absolute deviations traces the conditional median. It appears there as “an important special case” of that class, not as the class itself. Their claim for the family is comparative, and it is the sentence to hold onto whenever squared error is reached for by default. These estimators “have comparable efficiency to least squares for Gaussian linear models while substantially out-performing the least-squares estimator over a wide class of non-Gaussian error distributions.”

Fitted together, the pieces are readable. The median at q=0.50 says half of comparable conditional outcomes are expected below the predicted value. A lower quantile at q=0.10 marks a low conditional boundary, an upper quantile at q=0.90 a high one. The gap between two predicted quantiles reflects conditional dispersion, not a guaranteed coverage interval.

A count cannot be negative. Its variance often changes with its mean. Poisson regression commonly models the log of the conditional mean as a linear predictor, so exponentiated coefficients act multiplicatively, and time or units at risk enter as an exposure term.

That exposure term is not a conceptual flourish. scikit-learn's worked example on Poisson regression and non-normal loss fits claim frequency on freMTPL2freq, the French motor third-party liability data uploaded to OpenML on 2018-11-05 under CC0: 678,013 rows, 12 features, from Dutang and Charpentier's CASdatasets, which describes it as risk features for 677,991 policies observed mostly over one year around 2011–2013. The example models y = ClaimNb / Exposure and carries Exposure as sample_weight. On test data it reports a mean Poisson deviance of 0.594 for PoissonRegressor and 0.575 for a HistGradientBoostingRegressor trained with Poisson loss. The squared-error alternative fitted to the same rows makes the wrong distributional assumption, and the example says plainly what that costs: “The `Ridge` regression model can predict very low expected frequencies that do not match the data. It can therefore severely under-estimate the risk for some policyholders.”

Real counts can be overdispersed, zero-inflated, truncated, or dependent, and the size of that failure is measurable rather than rhetorical. In 2008 Achim Zeileis and two co-authors fitted competing count models to US National Medical Expenditure Survey data on physician office visits for 4,406 individuals aged 66 and over. “The observed number of zeros is 683 in 4406 observations.” The maximum-likelihood Poisson GLM expects 47. A negative binomial gives 608, a negative-binomial hurdle model 683 by construction, and a zero-inflated negative binomial 709. A fourteen-fold shortfall on the single quantity the model is supposed to describe is not a subtlety you catch in a summary table of coefficients.

Which count model you pick also decides which rows the fit listens to. Jay M. Ver Hoef and Peter L. Boveng put it plainly in Ecology in November 2007: “The variance of a quasi-Poisson model is a linear function of the mean while the variance of a negative binomial model is a quadratic function of the mean. These variance relationships affect the weights in the iteratively weighted least-squares algorithm of fitting models to data. Because the variance is a function of the mean, large and small counts get weighted differently in quasi-Poisson and negative binomial regression.” Their worked example is harbour seal counts from aerial surveys. They report “a data set that showed a dramatic difference on estimating abundance of harbor seals when using quasi-Poisson vs. negative binomial regression”. Two models with equal numbers of parameters, one data set, and an answer that moved because of the weights.

Use the target distribution and the decision asymmetry together.

State the estimand first: a mean, a median, a quantile, a count rate, or another functional. Describe the impossible outputs — nonnegativity, bounds, discreteness, exposure — before choosing a link. Examine tails and variance for skew, outliers, zero mass and scale-dependent spread; 683 observed zeros against 47 expected is exactly what that check exists to surface. Match a loss and link so that residuals and output constraints both enter training, rather than only one of them. Validate in decision units — asymmetric costs, quantile coverage, slice behaviour — instead of in whatever metric the library prints first. Then stress the rare regimes, and look hard at whether robustness has quietly hidden the very cases the system was built to handle.

A median model should be evaluated with median-oriented loss and residual summaries. A 90th-quantile model needs calibration of the fraction of outcomes below its predictions. A count model needs deviance, dispersion, zero behavior, and exposure checks: the scikit-learn insurance example reports a test-set mean Poisson deviance of 0.594 and 0.575 rather than an RMSE, and the four count models above are compared on the zero count itself — 683 observed against 47, 608, 683 and 709 expected.

Comparing every regression family only by RMSE can reward the wrong estimand. M5 ranked entrants on the Weighted Scaled Pinball Loss. The pinball loss at each of the nine quantiles is scaled by the mean absolute error of a one-step-ahead in-sample naive forecast, so series of different magnitudes are comparable. The scaled losses are then weighted by each series' dollar sales and aggregated across the twelve levels of the hierarchy. That is the machinery it takes to score a distributional claim honestly across 42,840 series. None of it reduces to a single squared-error number.