Smoothing and spline bases

Smoothing and spline bases#

Scatterplot smoothing and natural cubic spline bases.

lowess() is Cleveland’s (1979) robust locally weighted regression, and matches R’s stats::lowess. The n points are sorted by x. The fit at x_i uses the q = max(2, min(n, floor(span * n + 1e-7))) points nearest to x_i, which lie within h_i, the distance from x_i to its q-th nearest neighbour (counting x_i itself). Points at distance r get the tricube weight (1 - (r / h_i)^3)^3, which is set to 1 for r <= 0.001 h_i and to 0 for r > 0.999 h_i; when h_i is zero, every point tied with x_i gets weight 1. The fitted value is that of a weighted least squares line, or of a weighted mean when the weighted standard deviation of x in the neighbourhood is at most 0.001 times the range of x. A neighbourhood whose weights are all zero keeps the observed value. Each robustness iteration refits with every weight multiplied by Tukey’s bisquare (1 - (e_k / s)^2)^2 of the residuals e_k of the previous fit, with s six times the median absolute residual; the bisquare is set to 1 for |e_k| <= 0.001 s and to 0 for |e_k| > 0.999 s. The iterations stop early when s is below 1e-7 times the mean absolute residual. To save time, points within delta of the last fitted point are not fitted but interpolated linearly, and points tied with a fitted point share its value.

choose_lowess_span() is limma’s rule for the span of a LOWESS trend through n points (limma’s chooseLowessSpan): 1 for n up to small_n points, and min_span + (1 - min_span) (small_n / n)^power for larger n. With limma’s settings (small_n = 50, min_span = 0.3, power = 1/3), the part of the span above 0.3 falls in proportion to n^(-1/3), as the width of a histogram’s bins does in the rules of Scott (1979) and Freedman and Diaconis (1981): the span is 1 at 50 points, 0.65 at 400 and 0.395 at 20,000.

natural_spline_basis() returns a basis of the natural cubic splines with given interior and boundary knots: cubic between the knots, twice continuously differentiable at simple knots, and linear beyond the boundary knots. Its column space matches that of R’s splines::ns. The basis is built from the cubic B-splines on the knot sequence with each boundary knot repeated four times (de Boor 1978), restricted to the subspace with zero second derivative at both boundary knots.

References

Cleveland WS (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association 74, 829-836. doi:10.1080/01621459.1979.10481038

Cleveland WS (1981). LOWESS: a program for smoothing scatterplots by robust locally weighted regression. The American Statistician 35, 54. doi:10.2307/2683591

de Boor C (1978). A Practical Guide to Splines. Springer, New York. doi:10.1007/978-1-4612-6333-3

Freedman D, Diaconis P (1981). On the histogram as a density estimator: L2 theory. Zeitschrift fuer Wahrscheinlichkeitstheorie und verwandte Gebiete 57, 453-476. doi:10.1007/BF01025868

Hastie TJ (1992). Generalized additive models. Chapter 7 of Statistical Models in S, eds Chambers JM and Hastie TJ. Wadsworth & Brooks/Cole, Pacific Grove.

Scott DW (1979). On optimal and data-based histograms. Biometrika 66, 605-610. doi:10.1093/biomet/66.3.605

cytepool.stats.smoothing.choose_lowess_span(n, *, small_n=50, min_span=0.3, power=1 / 3)#

The span of a LOWESS trend through n points, by limma’s rule.

limma’s chooseLowessSpan(n, small.n, min.span, power): the span is 1 for n up to small_n, and

min_span + (1 - min_span) (small_n / n)^power

for larger n, which falls toward min_span as n grows and stays above it. With the defaults, limma’s, the span is 1 up to 50 points, 0.65 at 400 and 0.395 at 20,000. limma tunes the rule for mean-variance trends, which are usually monotonic and so suit large spans; voom() chooses its span with it from the number of features, as limma’s voom does.

limma computes its formula, capped at 1, for any arguments. This function refuses arguments outside the ranges below, where limma gives NaN for a negative n, and spans below 1 for fewer than small_n points for a negative power or a min_span above 1. At 0 points with a min_span of 1 or a small_n of 0, limma’s formula gives NaN, and this function gives 1, as limma’s documentation does.

Parameters:

  • n – Number of points to smooth, a non-negative integer.

  • small_n – Number of points up to which the span is 1, a non-negative integer.

  • min_span – Span that the rule approaches as n grows, in (0, 1].

  • power – Rate at which the span falls with n, in [0, 1]; 0 gives a span of 1 for every n.

Returns:

The span, the fraction of the points in each local fit of lowess(): in [min_span, 1].

Raises:

  • TypeError – If n or small_n is not an integer, or min_span or power is not a real number.

  • ValueError – If n or small_n is negative, min_span is not in (0, 1], or power is not in [0, 1].

Explained in: Precision weights

cytepool.stats.smoothing.lowess(x, y, *, span=2.0 / 3.0, iterations=3, delta=None)#

Smooth a scatterplot by robust locally weighted linear regression.

Cleveland’s LOWESS, with the conventions of R’s lowess(x, y, f=span, iter=iterations, delta=delta); the module docstring gives the method. Unlike R, which returns the fitted values sorted by x, the values are returned in the order of the input points. Every point tied in x gets the same fitted value, so the result does not depend on the order of the input points (up to rounding when x has ties).

Parameters:

  • x – (n,) finite values of the predictor.

  • y – (n,) finite values of the response.

  • span – Fraction of the points in each local fit. A span of 1 or more uses all points, as in R; a span below 2 / n still uses two.

  • iterations – Number of robustness iterations after the first fit.

  • delta – Distance in x within which a point is interpolated from the neighbouring fits instead of being fitted. Defaults to 1% of the range of x, as in R; 0 fits every distinct x.

Returns:

(n,) fitted values, in the order of x. A single point is returned unchanged.

Raises:

  • TypeError – If an input is not real-valued or is a numpy masked array, or iterations is not an integer.

  • ValueError – If x and y are not 1-D of the same non-zero length or are not finite, span is not positive and finite, iterations is negative, or delta is negative or not finite.

Explained in: Precision weights

cytepool.stats.smoothing.natural_spline_basis(x, df=None, *, knots=None, intercept=False, boundary_knots=None)#

Basis of natural cubic splines evaluated at x.

The columns span the same space as R’s splines::ns(x, df=df, knots=knots, intercept=intercept, Boundary.knots=boundary_knots), so a least squares fit on either gives the same fitted values; the columns themselves differ from R’s. Without an intercept the splines are those that vanish at the lower boundary knot, as in R.

With df, the df - 1 - intercept interior knots are placed at quantiles of x (R’s default type 7) with probabilities k / (K + 1), k = 1, …, K, over the values of x within the boundary knots. Ties in x can make knots coincide, which lowers the smoothness there as for repeated knots. A knot that falls on a boundary knot is moved inside, as R’s ns does, by one eighth of its distance to the nearest other knot (first at the lower boundary, then at the upper one).

NaN in x gives a row of NaN, as R’s NA does; knots and boundary knots are computed from the other values. Where R’s ns warns and increases a df below 1 + intercept, or accepts a knot on the lower boundary knot (returning a basis with a degenerate column), this function raises an error.

Parameters:

  • x – (n,) values at which to evaluate the basis; NaN marks a missing value.

  • df – Number of columns. Give df or knots; given both, they must agree.

  • knots – Interior knots (one or a sequence), strictly inside the boundary knots; repeated knots are allowed. Defaults to none, or to quantiles when df is given.

  • intercept – Whether the basis spans constant functions.

  • boundary_knots – The two knots beyond which the splines are linear. Defaults to the range of the observed values of x; given in either order.

Returns:

(n, df) matrix, with df = len(knots) + 1 + intercept.

Raises:

  • TypeError – If an input is not real-valued or is a numpy masked array, df is not an integer or intercept is not a bool.

  • ValueError – If x is not 1-D or has infinite values, df is too small, knots are not finite or not strictly inside the boundary knots, the boundary knots are not two distinct finite values, or there are too few values of x to place the knots.

Explained in: Empirical Bayes moderation and TREAT