PolyResids

Description

PolyResids returns four types of residuals in polynomial regression. These residuals are useful for assessing model fit and for identifying outliers.

Syntax

PolyResids(Y, X, [p = 1], [pmod = FALSE], [S]), where

• Y and X are the dependent and independent variables as row or column vectors
• p is the polynomial order, defaulting to 1 (best-fit line)
• pmod is a boolean that controls whether to return Bonferroni-corrected p-values for the studentized deleted residuals, providing a more objective method of labeling outliers. By default pmod = FALSE.
• S is a vector of standard deviations of the Y values. If provided then weighted least-squares is used.

The function returns 4-5 columns, depending on the value of pmod:

1. The raw residuals, $e_i = y_i - \hat{y}_i$
2. The standardized residuals, ${e_i}/{s(e_i)}$, appropriately weighted if S is provided in the function call.
3. The deleted (“leave-one-out”) residuals, $e_{-i}$, which are the residuals for a given value $y_i$ if that value were not used in the regression fitting.
4. The studentized residuals, $e_{-i}/s(e_{-i})$, appropriately weighted if S values are provided and WLS is used. These quantities follow a t-distribution if the underlying errors are Gaussian.
5. If pmod is TRUE then a fifth column exists that returns the Bonferroni-corrected p-values for the studentized residuals. If a value falls below 0.05 (or whatever significance level you are comfortable with) then you have some objective justification for omitting the point from the fit. At the very least you can delete this point and observe the effect on the fit and other calculations based on it.

Example

The example below shows how residuals can provide a visual clue to help choose the appropriate regression model. With a proper model the residuals should be randomly distributed around the fit. For the first and second orders below there is still some systematic structure to the residuals, implying some unmodeled dependence on the X values. But the residuals for the third and fourth order fits look randomized. It is best to choose the simplest model that adequately describes the functional dependence of Y on X, and in this case that would be a cubic (third order) fit.

Residuals can help identify outliers as well, though there don’t appear to be any in these data. The image below shows the entire output of PolyResids for the third-order fit, with the point with the smallest p-value highlighted. However its Bonferroni-corrected value is still nowhere near any reasonable significance level, meaning the justification for omitting this data point from modeling is weak.

Download: PolyResids.xlsx

Source Code

PolyResids = LAMBDA(
    Y, // row/column vector of dependent values
    X, // row/column vector of independent values
    [order], // optional polynomial order (default 1)
    [pmod], // optional flag to return Bonferroni-corrected p-values
    [S], // optional std devs of Y for weighted LS
    LET(
        order, IF(ISOMITTED(order), 1, order),
        pmod, IF(ISOMITTED(pmod), FALSE, pmod),
        p, order + 1, // number of regression params
        // filter out rows with non-numeric values
        NumFilter, IF(ISOMITTED(S), ISNUMBER(X) * ISNUMBER(Y), ISNUMBER(X) * ISNUMBER(Y) * ISNUMBER(S)),
        Yn, TOCOL(FILTER(Y, NumFilter)),
        Xc, FILTER(X, NumFilter),
        n, COUNT(Xc),
        // weights (all equal to 1 for ordinary LS)
        sdn, IF(ISOMITTED(S), 1, TOCOL(FILTER(S, NumFilter))),
        wt, 1 / sdn ^ 2,
        Xn, HSTACK(Ones(n), TOCOL(Xc) ^ SEQUENCE(, order)),
        // calculate the projection ("Hat") matrix, yhat = H*Y
        inner, MINVERSE(MMULT(TRANSPOSE(Xn * wt), Xn)),
        H, MMULT(MMULT(Xn, inner), TRANSPOSE(Xn * wt)),
        h_ii, MDiag(H),
        // calculate residuals and s_res
        resids, Yn - MMULT(H, Yn),
        sres, SQRT(SUM(wt * resids ^ 2) / (n - p)),
        // calculate the standardized residuals
        s_resids, SQRT((1 - h_ii) / wt) * sres,
        resids_std, resids / s_resids,
        // now the deleted and studentized residuals
        delresids, resids / (1 - h_ii),
        s_del, SQRT(((n - p) * sres ^ 2 - wt * resids ^ 2 / (1 - h_ii)) / (n - p - 1)),
        resids_stud, resids / (s_del * SQRT((1 - h_ii) / wt)),
        output, IF( //return Bonferroni-corrected modified p-values if requested
            pmod,
            HSTACK(resids, resids_std, delresids, resids_stud, T.DIST.2T(ABS(resids_stud), n-p-1) * n),
            HSTACK(resids, resids_std, delresids, resids_stud)
        ),
        IF(OR(sdn <= 0), "ERROR: std devs must be positive", output)
    )
);