PolySres

Description

PolySres returns the standard deviation of the residuals and the standard deviation of the deleted (leave-one-out) residuals. These can be useful in choosing an appropriate calibration model.

The companion function Sres is a more generalized version of this function, applying to any linear regression model.

Syntax

PolySres(Y, X, [order = 1], [intcpt = TRUE], [S]) where

• Y and X are the dependent and independent variables, respectively, as row or column vectors
• order is the polynomial order, defaulting to order = 1 (best-fit line). It can also be a vector of orders, such as {1;2;3} (see below).
• intcpt is a boolean that governs whether to include the intercept in the model (defaults to TRUE)
• S is a vector of standard deviations of the Y values. If provided then weighted least-squares is used (see Sres).

The function returns two values as a row vector: the standard deviation of the residuals, $s_\text{res}$, and the standard deviation of the deleted residuals, $s_\text{del}$. A “deleted” residual is the residual of a data point when it is omitted from the fit.

If order is a vector of orders, the output has one row ($s_\text{res}$, $s_\text{del}$) per order, all fit to the same data points. This makes it easy to compare several orders in a single formula.

Example

It is well known that the $R^2$ value is guaranteed never to decrease when the order of a polynomial fit increases, hence it is not a useful statistic in helping to determine what order is appropriate to use when fitting. The standard deviation of the residuals is a better statistic to use since it penalizes you for increasing the order since the degrees of freedom (the denominator in calculating $s_\text{res}$) decreases. The standard deviation of the deleted residuals is an attempt to measure the uncertainty in predicting values which is often a major use case of a regression model, and it often penalizes one even more for overfitting.

Choosing an appropriate model is often best done with a mixture of statistics (like those returned by PolySres) and visualization. The example below illustrates this point by using both to conclude that a cubic model is the best model to fit a data set. (The statistics for all of the orders can also be returned by a single formula, such as =L.PolySres(Y, X, SEQUENCE(4)).)

The above case is unusually clear-cut, with both $s_\text{res}$ and $s_\text{del}$ declaring a clear winner by reaching a minimum. It is common that either (or both) values will decrease drastically and then plateau, bouncing around some minimum value. And sometimes they do not give the same answer, particularly for smaller sample sizes. That is where data visualization (eg, residual plots, as shown) can help, as well as a parsimonious preference for simpler models (ie ones with fewer parameters) when the evidence is mixed.

Download: PolySres.xlsx

Source Code

PolySres = LAMBDA(
    Y, // dependent values as a column/row vector
    X, // independent values as a column/row vector
    // X and Y must be same length but can contain non-numbers and empty cells
    [order], // polynomial order, or a vector of orders (defaults to 1)
    [intcpt], // optional flag, if TRUE (the default) an intercept is included in the model
    [S], // optional std devs of Y for weighted LS
    LET(
        order, IF(ISOMITTED(order), 1, order),
        intcpt, IF(ISOMITTED(intcpt), TRUE, intcpt),
        weighted, NOT(ISOMITTED(S)),
        // filter out rows with non-numeric values
        NumFilter, IF(weighted, ISNUMBER(X) * ISNUMBER(Y) * ISNUMBER(S), ISNUMBER(X) * ISNUMBER(Y)),
        Xn, TOCOL(FILTER(X, NumFilter)),
        Yn, TOCOL(FILTER(Y, NumFilter)),
        Sn, IF(weighted, TOCOL(FILTER(S, NumFilter)), 1),
        n, ROWS(Yn),
        // s_res and s_del for one order k, or an error message row
        one_order, LAMBDA(k,
            IF(
                NOT(ISNUMBER(k)),
                {"ERROR: order must be a whole number >= 1", ""},
                IF(
                    OR(k < 1, k <> INT(k)),
                    {"ERROR: order must be a whole number >= 1", ""},
                    IF(
                        n - k - intcpt - 1 < 1,
                        {"ERROR: too few points for this order", ""},
                        IF(weighted, Sres(Yn, Xn ^ SEQUENCE(, k), intcpt, Sn), Sres(Yn, Xn ^ SEQUENCE(, k), intcpt))
                    )
                )
            )
        ),
        // one row per order (REDUCE because MAP/BYROW can't return a row per element)
        rows_, REDUCE("", TOCOL(order), LAMBDA(acc, k, VSTACK(acc, one_order(k)))),
        IF(OR(Sn <= 0), "ERROR: std devs must be positive", DROP(rows_, 1))
    )
);