Improved Spread-Location Visualization

A.I. McLeod

Abstract

An improved spread-location plot for checking for lack of constant variance is derived.  The improved plot allows for p-th power transformations on the absolute residuals as well as a boxplot to check for symmetry in the plot.  This plot is useful in visually checking the adequacy of many types of statistical models.  In practice, a common inadequacy of many statistical models is that the variability of the residuals increase or decrease monotonically with the fit.  This type of statistical model inadequacy is called montone spread and can usually be corrected by a suitable power transformation.  The use of a Mathematica package, slplot.m,  which implements this improved plot is illustrated. This is an electronic supplement to McLeod (1999).

Keywords:  Model criticism; Monotone spread; Power transformation; Statistical visualization; Symbolic skewness computations in Mathematica

Plotting residuals vs fitted values is a widely used model diagnostic.  Most books on regression illustrate the use of this plot to detect wedge shape patterns which may indicate an increase in spread of the residuals with fit.

The Spread-Location (SL) Plot introduced by Cleveland (1979) is an improved version of this plot.  Rather than plot the residual, Cleveland (1979) suggested the absolute residual vs the fitted value be plotted and a robust loess fit be used to assess if there is a monotonic change in spread.  Notice that in the original version a so-called wedge shape pattern  could in some cases just be due to more data at higher fitted values.  Cleveland (1993) suggested that rather than simply plot the absolute residual it would be better to plot the square root of the absolute residual since in the Gaussian case, the distribution of the absolute residual is quite skewed while that of the square-root absolute residual is not.  However we have found that for many non-Gaussian distributions the square-root of the absolute residual is still appreciably skewed.  Therefore we recommend choosing a suitable transformation to remove skewness in the SL plot.

First we present some skewness computations which demonstrate the possible need for an improved visualization.  Then we conclude by illustrating the use of Mathematica package which implements the improved SL plot and demonstrates its usefulness with real data.

Skewness Derivations

Given a random variable X with the probability density, f(x), defined on (-∞, ∞), we compute the skewness coefficient of  \(TraditionalForm\`| X  |\^p\). We look at plots of the skewness coefficient vs p for some specific distributions.  These distributions include the normal distribution as well as various thick-tailed distributions such as (Subscript[t, ν]) and contaminated normal.  Additionally we look at the centered exponential as an example of a non-symmetric error distribution.

First we derive some general functions which we can use to compute the skewness coefficient of  \( | X\( | \^p\)\).  Symbolic or exact computation can be done using Integrate but this sometimes does not provide an answer.  So we switch to NIntegrate.

General Methods

Given  X and its probabity density function f(x), defined on (-∞, ∞), the probabilty density function of | X |  is given by h(x)  =   f(x) + f(-x) with support on (0, ∞) .

h[x_]:=f[x]+f[-x]

More generally for p ≠ 0, the density of  ((| X)^p) |   is (Subscript[f, p])(x) =  q h((x^q)) (x^(q - 1)) where q = 1/p .

f[x_, p_ /; p ≠ 0]  :=   
      (1/p) h[x^(1/p)] x^(1/p - 1)

The case p = 0, corresponds to the natural log transformation, log( | X |  ) and so in this case the density function can be written (Subscript[f, p])(x) = (e^x) f((e^x)).

f[x_, p_ /; p == 0]  :=    (E^x)  h[ (E^x) ]

The following Mathematica functions use symbolic integration to obtain the skewness coefficient for X with pdf  h(x) and given specified range of support.

\(Ex[h_, \ support_] := Integrate[x*h[x], support]; 

    Var[\ h_, \ support_] := 
      Integrate[\((x - Ex[\ h, \ support])\)^2\ *\ h[x], \ support]; 

    g1[h_, \ support_] := \ 
      Integrate[\((x - Ex[\ h, \ support])\)^3*h[x], \ support]/
        \((Var[\ h, \ support])\)^\((3/2)\); \)

Sometimes an analytical solution is not available so we use numerical integatration functions,

\(NEx[h_, \ support_] := NIntegrate[x*h[x], support]; 

    NVar[h_, \ support_] := 
\ \ 
      Module[{mean = NEx[h, \ support]}, 
        NIntegrate[\((x - mean)\)^2\ *\ h[x], \ support]]; 

    Ng1[h_, \ support_] := 
\ \ 
      Module[{mean = NEx[h, \ support]}, \ 
\ \ \ \ 
        NIntegrate[\((x - mean)\)^3*h[x], \ support]/
          \((NVar[\ h, \ support])\)^\((3/2)\)]; \)

Standard Normal Distribution

f[x_]  :=  1/Sqrt[2 Pi]*Exp[-x^2/2]

g1N0 = {0, Ng1[f[#, 0]&, {x, -∞, ∞}]};

g1N0

In many situations, the warning message about underflow can be ignored.  The underflow message means that a very small number was produced during the computation and that it was set to zero.  Since the end numerical result of this last integration was about -1.53514, it is almost certainly fully accurate.  Also the fact that the skewness when p = 1/2 is nearly 0 verifies the result of Cleveland (1993).  Also if numerical error were important, it is unlikely that the plot below of g1N would look so reasonable.  It some cases we can verify using the analytical functions that the result is correct.  For example, use Integrate instead of NIntegrate produces other warning messages but the result is the same since Mathematica using sensible default actions which produce a valid answer in many cases.  Further experimentation could be done using other quadrature methods available in the Standard Library NumericalMath`GaussianQuadrature`.

Our viewpoint here is scientific rather than trying derive an absolutely rigorously mathematical justification.

g1[f[#, 0]&, {x, -∞, ∞}] // N

g1N = 
      Table[{p = 1 - i/10, g1[f[#, p]&, {x, 0, ∞}]}, {i, 0, 9}];

g1N0

g1N = Join[g1N // N, { g1N0}]

ListPlot[g1N,  PlotJoined -> True,  Frame -> True];

Student t on 5 df

Needs["Statistics`ContinuousDistributions`"];

f[x_]  := PDF[StudentTDistribution[5],  x]

g1t0 = {0, Ng1[f[#, 0]&, {x, -∞, ∞}]};

g1t = 
      Table[{p = 1 - i/10, g1[f[#, p]&, {x, 0, ∞}]}, {i, 0, 9}];

g1T = Join[g1t // N, { g1t0}];

Needs["Graphics`MultipleListPlot`"]

MultipleListPlot[g1T, g1N,  PlotJoined -> True,  Frame -> True,  
      SymbolShape -> None, FrameLabel -> {"p", "g"}];

Student t on 4 df

Needs["Statistics`ContinuousDistributions`"];

f[x_]  := PDF[StudentTDistribution[4],  x]

g1t04 = {0, Ng1[f[#, 0]&, {x, -∞, ∞}]};

g1t4 = 
      Table[{p = 1 - i/10, Ng1[f[#, p]&, {x, 0, ∞}]}, {i, 0, 9}];

g1T4 = Join[g1t4 // N, { g1t04}];

ListPlot[g1T4,  PlotJoined -> True,  Frame -> True];

MultipleListPlot[g1L, g1N,  PlotJoined -> True,  Frame -> True,  
      SymbolShape -> None, FrameLabel -> {"p", "g"}];

MultipleListPlot[Drop[Drop[g1T4, 4], -2], 
      Drop[Drop[g1N, 4], -2],  PlotJoined -> True,  Frame -> True,  
      SymbolShape -> None, FrameLabel -> {"p", "g"} ];

Laplace Distribution

LaplacePDF[x_,  β_]  :=  
      (1/(2 β))  E^(-Abs[x]/β ); 
    f[x_]  =  LaplacePDF[x,  1/Sqrt[2]];

g1L0 = {0, Ng1[f[#, 0]&, {x, -∞, ∞}]};

g1L = 
      Table[{p = 1 - i/10, g1[f[#, p]&, {x, 0, ∞}]}, {i, 0, 9}];

g1L = Join[g1L // N, { g1L0}];

ListPlot[g1L,  PlotJoined -> True,  Frame -> True]

MultipleListPlot[g1L, g1N,  PlotJoined -> True,  Frame -> True,  
      SymbolShape -> None, FrameLabel -> {"p", "g"}];

Contaminated Normal Distribution

\(TukeyPDF[x_, \ p_, \ c_]\  := 
      Sqrt[1 - p + p*c^2]*
        \((\((1 - p)\)*\ 
              PDF[NormalDistribution[0, 1], \ x*Sqrt[1 - p + p*c^2]] + 
            p*PDF[NormalDistribution[0, \ c], \ x*Sqrt[1 - p + p*c^2]])\); 

    f[x_]\  := \ TukeyPDF[x, \ 0.05, \ 3]; \)

g1mix0 = {0, Ng1[f[#, 0]&, {x, -∞, ∞}]};

g1mix = 
      Table[{p = 1 - i/10, Ng1[f[#, p]&, {x, 0, ∞}]}, {i, 0, 9}];

g1mixtab = Join[g1mix // N, { g1mix0}];

MultipleListPlot[g1mixtab, g1N,  PlotJoined -> True,  
      Frame -> True,  SymbolShape -> None, 
      FrameLabel -> {"p", "g"}];

Centered Exponential Distribution

Needs["Statistics`ContinuousDistributions`"];

f[x_]  :=  PDF[GammaDistribution[1,  1], x + 1]

Integrate[x f[x],{x,-1,∞}]

g1e0 = {0, Ng1[f[#, 0]&, {x, -∞, ∞}]};

g1e = 
      Table[{p = 1 - i/10, Ng1[f[#, p]&, {x, 0, ∞}]}, {i, 0, 9}];

g1E = Join[g1e // N, { g1e0}];

MultipleListPlot[g1T, g1N,  PlotJoined -> True,  Frame -> True,  
      SymbolShape -> None, FrameLabel -> {"p", "g"}];

Illustrative Examples

SetDirectory["d:\\p\\math\\1999\\writalk"];

<< "slplot.m"

?SLPlot

Options[SLPlot]

Example 1: An adequate model

\(Needs["Statistics`ContinuousDistributions`"]; 

    Needs["Statistics`LinearRegression`"]; \)

\(SeedRandom[3]; 
X = Range[50]; 

    error = RandomArray[NormalDistribution[0, 1], 50]; 

    Y = 25 + 0.5*X\  + \ error; \)

{res,fits} = {StudentizedResiduals,PredictedResponse}/.Regress[Transpose[{X,Y}], {1, x^2}, x,RegressionReport->{StudentizedResiduals,PredictedResponse} ];

Note that one could also use StandizedResiduals or FitResiduals.  It usually won't make much difference in the plot.

SLPlot[fits, res];

The SL-plot looks fine as it should.  We can redo the SL-plot with a boxplot added to check for symmetry.

SLPlot[fits, res, BoxPlot -> True];

Example 2: Montone Spread

In this example, we construct a mis-specified model which has monotone spread in the residuals.

\(Y = 1 + 0.1*X\  + \ error/4; 
Y = Exp[Y]; 
    
{res, fits}\  = 
      \ {StudentizedResiduals, PredictedResponse} /. 
        Regress[Transpose[{X, Y}], \ {1, \ x^2}, \ x, 
          RegressionReport -> {StudentizedResiduals, PredictedResponse}\ ]; 
    \)

SLPlot[fits, res];

This diagnostic plot indicates that a power transformation is needed. Since the spread is increasing, we should try going down the ladder of power transformation.  If we use a log transformation, we get an adequate model.

\({res, fits}\  = 
      \ {StudentizedResiduals, PredictedResponse} /. 
        Regress[Transpose[{X, Log[Y]}], \ {1, \ x^2}, \ x, 
          RegressionReport -> {StudentizedResiduals, PredictedResponse}\ ]; 

    SLPlot[fits, res]; \)

Example 3: Another Mis-specified model

This example shows the importance of visualization.  There is no monotone spread problem, but the SL-plot clearly indicates an inadequate model since there is a systematic component present.  For an adequate model, the SL-plot should not have any non-random looking features.

\(Y = 25 + 0.5*X\  + \ error/4; 
Y = Sqrt[Y]; 
    
{res, fits}\  = 
      \ {StudentizedResiduals, PredictedResponse} /. 
        Regress[Transpose[{X, Y}], \ {1, \ x^2}, \ x, 
          RegressionReport -> {StudentizedResiduals, PredictedResponse}\ ]; 
    \)

SLPlot[fits, res];

Example 4: Ethanol Data

We illustrate with the famous ethanol data series which was given first in the paper of Brinkman (1981) and is available electronically from the statistics archive Statlib as part of a collection of datasets from the book of Cleveland (1993).  The Statlib archive file is available from http://www.stat.cmu.edu/S/visualizing.tables (Cleveland, 1997). For convenience, the files ethanol.dat and ethanol.txt used below are also available from my homepage (http://www.stats.uwo.ca/mcleod).  The files ethanol.dat includes the original data as well as the incorrect model used for illustration by Cleveland (1993).  Cleveland (1993, 1997) showed how mis-specifying the data transformation as a log when in fact no transformation was needed induced monotonic spread.  This example shows how Cleveland's visualization (Cleveland, 1993, Figure 4.33 and McLeod, 1999, Figure 2) can be improved on by noting the obvious problem with skewness.

ethanol = 
      ReadList["ethanol.dat", {Real, Real, Real, Real, Real}];

The variables in this data are described in the file ethanol.txt,

ReadList["ethanol.txt", {String}] // TableForm

{c, e, nox, fit, res} = Transpose[ethanol];

Now take a log transformation and fit a 5-th order polynomial with all interactions.  We regress log NOx on polynomials in C, denoted by x, and E, denoted by y.

SLPlot[fit, res, BoxPlot -> True];

SLPlot[fit, res, 0, 1,  BoxPlot -> True];

References

Brinkman, N.D. (1981), Ethanol Fuel -- A single-cylinder engine study of efficiency and exhaust emissions, Society of Automobile Engineers, 80, 1410-1424.

Cleveland, W.S. (1979),  Robust locally weighted regression and smoothing scatterplots,  Journal of the American Statistical Association,  74, 829-836.

Cleveland, W.S. (1993),  Visualizing Data. Summit: Hobart Press.

Cleveland, W.S. (1997),  Visualizing Scripts,  http://www.stat.cmu.edu/S/visualizing.scripts

McLeod, A.I. (1999),  Improved spread-location visualization,  Journal of  Computational and Graphical Statistics 8, 135-141.