Notes on errorbar enhancements

[mwaskom]

Here are some notes on planned changes to error bar specification.

Currently, what the error bars show is controlled through the ci parameter. This can be either a number, setting the width of a bootstrap confidence interval, or the string "sd", indicating that the error var covers +/- the standard deviation of the data around the estimate value.

Some problems with this have been routinely noted:

• There's no option for parametric confidence intervals/standard error
• There's no option for showing a measure of data spread other than +/- 1 sigma
• ci="sd" does not really make conceptual sense as a parametrizing (it is the result of a short-sighted API decision)

In effect, you can think of the options as having a 2D taxonomy defined by whether the error bars show a measure of estimate certainty or data spread and whether the computation is parametric or nonparametric. Currently, we occupy two cells in this matrix:

	Estimate certainty	Data spread
Parametric		ci="sd"
Nonparametric	ci=95, ci=68, etc.

I would like to fill out the matrix. But we have a few challenges:

• As mentioned, the current API is not great, and overloading the meaning of ci further is a nonstarter
• There is no centralized location in the code where this parameter is interpreted and used

Plans for the new API will involve a new parameter, probably called errorbar but possibly error, errbar or some other shorthand, that accepts a tuple of the form (kind, level). The first element determines what the error bars show, and the second parametrizes them. One proposal is to fill out the space like this:

	Estimate certainty	Data spread
Parametric	("se", scale)	("sd", scale)
Nonparametric	("ci", size)	("pi", size)

IMO, there is a lot of sense to this. You have four options for kind, each named using a bigram initialism. There are two kinds of level parameters:

• scale: multiplicatively scales a parametric error metric (e.g. ("sd", 3) gives you a 3-sigma error bar, ("se", 1.96) gives you a ~95% parametric confidence interval
• size sets the size of a nonparametric interval with percentiles (of the boostrap distribution for ci and the input data for pi) of (1 - size) / 2, 1 - (1 - size) / 2

There are also some potential drawbacks

• "pi" (i.e., "percentile interval") doesn't seem to be a commonly used term for a nonparametric measure of data spread. Actually I'm not sure there really is a term in the stats literature for such an interval, even though it's a very reasonable thing to plot (e.g. Support plotting quantiles of the data distribution with lineplot #1501)
• If you really want parametric 95% confidence intervals, this parametrization leaves you limited to a Z interval (and requires you to understand how to construct one from a standard error)

API decisions aside, the right implementation is going to take some thinking. Currently each module does its own errorbar computations. Most errorbars appear in the context of an aggregation-with-estimator operation. This can likely be abstracted. The other place they show up is in the regression module, where error bars are shown around the regression line. This needs to be handled differently, but statsmodels now has the get_prediction method which will do a lot of the work for us. We'll need a general enough implementation such that we can handle special cases (like logistic regression, where the SE/SD scaling should happen in logit space).

Here are some assorted open questions

• Should we accept simple strings (e.g. errorbar="sd") with a default level value used internally?
• This simple 4 option system is still fairly limiting; it may disappoint those who would like to be able to use a generic function to get error bars (e.g. Added option for ci to be a callable. #2332). What might that API look like?
• Is it a sensible API option for sd to correspond to the prediction interval in a regression model?
• Should standard error correspond to the estimator and raise if the one used doesn't have a defined standard error? In other words, what would we do with estimator="median", errorbar="se"? And if the estimator is a callable, should we use its name to associate with the correct standard error function?
• It would be nice to have seaborn support multiple error bars from a sequence of level parameters, e.g. 1-2-3 sigma bands or 68-95-99 CIs (e.g. [Feature request] Support lineplot with error bands of list of scaled standard deviations #1492). I like this kind of plot, but each plotting function will have to define its own logic for showing multiple error bars (e.g. layered alpha for error bands in lineplot, lines of diminishing width in pointplot). But still, if it's going to happen, we should at least plan for it here.
• What about additional arguments for bootstrapping (i.e. n_boot, seed?) It would be nice to reduce the number of parameters in the main function signatures, but I would like to keep the argument for errorbar a simple tuple and not a more complex object that could take optional parameters. I think...
• What about loess? (regplot: support confidence intervals with lowess model #552). Bootstrapping is still very slow, but statsmodels seems to still not have analytic confidence bands.
• What's the right order of operations for working on this? It should probably not be (fully) implemented until the categorical/regression modules can be refactored to use the core objects (where this should be handled).

[shoyer]

Hmm, it's puzzling indeed that nobody talks about "percentile intervals."

Perhaps it's more common to be interested in assymetric data quantiles rather than confidence bounds? In practice, people definitely do plot this sort of thing, I've just seen it explicitly called "the 5th and 95th percentiles" rather than "the 90th percentile intervals".

[shoyer]

To support generic error estimates, I would suggest supporting a pair of column names giving lower/upper bounds. This would be nice for cases where you want more control over error estimates for some reason, e.g., when Seaborn is used to plot statistics from a large dataset (so bootstrapping becomes very slow/impractical). This has come up for me a few times.

[mwaskom]

I would suggest supporting a pair of column names giving lower/upper bounds.

Do you mean in the case where you're skipping aggregation altogether? That's something to think about but a bit different from the question of whether to allow users to supply their own function for computing an error bar interval from a vector of data that is getting aggregated.

[mwaskom]

Regarding cases where the user has their own error estimates already, the question for me is about whether it's better for each function that can do aggregation to accept xmin/xmax/ymin/ymax (adding four parameters to many functions) or if there should be a separate intervalplot (or similar) function.

The downside is the latter will be more annoying to use in a faceted context.

[mwaskom]

Oh I guess the point is that if you make it easy to add error bars at specific values, users can just do their bespoke computation externally to seaborn, obviating the need for accepting a generic function.

[mwaskom]

A rough proof of principle for the plots that do aggregation with an estimator would use an object that looks like this:

def validate_errorbar_arg(arg):
    DEFAULT_LEVELS = {
        "ci": 95,
        "pi": 95,
        "sd": 1,
        "se": 1,
    }
    if callable(arg):
        return arg, None
    elif isinstance(arg, str):
        method = arg
        level = DEFAULT_LEVELS.get(arg, None)
    else:
        method, level = arg
        # TODO error handling on noniterable or iterable of wrong length

    sns.utils._check_argument("errorbar", DEFAULT_LEVELS, method)
    # TODO check that level is a number

    return method, level

class EstimateAggregator:

    def __init__(self, estimator, errorbar, level=None, **bootstrap_kws):
        
        self.estimator = estimator
        self.errorbar = errorbar
        self.level = level
        self.bootstrap_kws = bootstrap_kws

    def __call__(self, data):

        estimate = data.agg(self.estimator)

        # TODO handle case where length of data is <= 1
        
        if callable(self.errorbar):
            err_min, err_max = self.errorbar(data)
        elif self.errorbar == "se":
            # TODO error handling when estimator != "mean"?
            se = data.sem() * self.level
            err_min, err_max = estimate - se, estimate + se
        elif self.errorbar == "sd":
            sd = data.std() * self.level
            err_min, err_max = estimate - sd, estimate + sd
        elif self.errorbar == "pi":
            alpha = 100 - self.level
            pcts = alpha / 2, 100 - alpha / 2
            err_min, err_max = np.percentile(data, pcts)
        elif self.errorbar == "ci":
            units = data.get("units", None)
            boots = sns.algorithms.bootstrap(data, func=self.estimator, units=units, **self.bootstrap_kws)
            err_min, err_max = sns.utils.ci(boots, self.level)

        return pd.Series({"_": estimate, "_min": err_min, "_max": err_max})

The regression

Status of open questions:

• Leaning towards accepting a string argument specifying the method with a default level assigned internally
• Leaning towards accepting a generic callable that takes a vector and returns an interval
• Leaning against doing anything to facilitate multiple errorbars in a single function call
• Leaning towards keeping bootstrap arguments separate