ENH: Add confidence intervals / standard errors?

[nickeubank]

Looks like Gadfly.jl has implemented confidence intervals for linear models, but apparently they need this package to expose standard error estimates to add them for loess. Any chance that's feasible?

[tlnagy]

If you would like to take a stab at exposing the error estimates, I would be willing to take a look at your PR and merge it.

[nickeubank]

Unfortunately I don’t know how to calculate standard errors for loess... :( sorry!

[tpapp]

AFAIK the only way to do this for loess is with bootstrapping.

[juliohm]

As far as I remember, we can still have "deviation" bands as opposed to confidence bands in LOESS right? I've implemented it a long time ago in LocallyWeightedRegression.jl: https://nbviewer.jupyter.org/github/JuliaEarth/LocallyWeightedRegression.jl/blob/master/docs/Usage.ipynb

[Gkreindler]

If this is off-topic, I apologize in advance. I came here looking for confidence intervals for loess and ended up coding a bootstrap version.

It may be wrong (I always include the endpoints in every bootstrap run) but useful for quick data visualization.

Sharing in case it's useful for anyone else.

using Loess
using Random

function loess_bootstrap(xs, ys;
            nboot=200, step_x=0.1)

    nx = length(xs)

    # min and max of the x range
    idx_xmin = argmin(xs)
    idx_xmax = argmax(xs)

    # main model
    model = loess(xs, ys, span=0.5)
    us = range(extrema(xs)...; step = step_x)
    vs_main = Loess.predict(model, us)

    n_us = length(us)

    # Bootstrap
    vs_boot = zeros(nboot, n_us) # will hold all bootstrap predictions here
    for iboot=1:nboot
        print(".")

        # sample
        idx_boot = sample( (1:nx) |> Array |> vec, nx-2, replace=true)

        # always include x range endpoints (otherwise can't predict)
        # this point is probably sketchy
        push!(idx_boot, idx_xmin)
        push!(idx_boot, idx_xmax)

        # run model on bootstrap sample
        model = loess(xs[idx_boot], ys[idx_boot], span=0.5)

        vs_boot[iboot, :] .= Loess.predict(model, us)
    end

    # get CI percentiles for each value on the x axis
    vs_ci_low  = zeros(n_us)
    vs_ci_high = zeros(n_us)
    for icol=1:n_us
        vs_ci_low[icol]  = percentile(vs_boot[:, icol], 2.5)
        vs_ci_high[icol] = percentile(vs_boot[:, icol], 97.5)
    end

    # plot
    plot(us, vs_main, color=:black)
    plot!(us, vs_ci_low , color=:black, alpha=0.5, label="95% CI")
    plot!(us, vs_ci_high, color=:black, alpha=0.5, label="")

end