Bivariate Loess

A.I. McLeod, University of Western Ontario
Copyright ©, A.I. McLeod, 1999.

This notebook illustrates the use of the bilo.m package for visualization and modelling using bivariate loess models.

Gaussian Case

Given bivariate observations ((Subscript[x, i]), (Subscript[y, i])),  i = 1, 2,   ...  ,  n, the basic model that can be fit may be written,

(Subscript[y, i]) = g((Subscript[x, i])) + (Subscript[ϵ, i]),  i = 1,   ...  , 
       n.

where (Subscript[e, i]) ~ NID( 0,  (σ^2)) and g(x) is a local linear polynomial of degree λ ≥ 0.  The local linear polynomial may be written,

\(TraditionalForm
    \`g(x) = \(β\_0\%\((x)\) + \(β\_1\%\((x)\)\) x\  + \  ... 
          \)\  + β\_λ\%\((x)\)\ x\^λ\)

The parameters (Subscript[β, 0]^((x))), (Subscript[β, 1]^((x))),   ...  , 
      (Subscript[β, λ]^((x))) are estimated by weighted least squares for each value of x.  The weight function weights the data, ((Subscript[x, i]), (Subscript[y, i])),  so that data values near to x have greater weight than those farther away from x.  Following Cleveland (1996) we use the tricube weight function,

T(z) =({{\(\((\(1 - \) | z\( | \^3\))\)\^3\),  | z | ,  ≤ 1}, {0,  | z | ,  > 1}})

to define the local neighborhood weights for the data at the point x,

\(TraditionalForm
      \`\(w\_i\)(x) = 
        T(\(\(Δ\_i\)(x)\)/\(Δ(x, α)\))
          \),

where \(TraditionalForm
      \`\(Δ\_i\)(x)\  = \ \( | x - x\_i | \)\) and Δ(x, α) controls the amount smoothing.  Larger values of  Δ(x, α) result in more smoothing.  As \(TraditionalForm
      \`Δ(x, α) -> ∞\), (Subscript[w, i])(x) -> 1 for each i = 1, 2,   ...  ,  n and the local linear model reduces to the standard parametric polynomial regression.  For 0 < α ≤ 1,  Δ (x, α) is the distance to the q-th nearest neighbor where q = [α n], where [●] is the integer part.  Hence, \(TraditionalForm
      \`Δ \((x, α)\) = 
        \(Δ\_\((q)\)\)(x)\), where (Subscript[Δ, (q)]) (x) denotes the q-th largest value of (Subscript[Δ, i]) (x),  i = 1,   ...  , n
      .  For α > 1, \(TraditionalForm
      \`Δ(x, α) = 
        α\ \(\(Δ\_\((n)\)\)(x)\)\).  It follows that as α -> ∞, the local linear model reduces to a parametric polynomial regression of degree λ.

Approximate statistical inference methods are available for Gaussian loess models (Cleveland and Devlin, 1988).

<< "c:\\p\\math\\1998\\jsm\\bilo.m"

?LoessFit

curve = 
      ReadList["c:\\p\\math\\1998\\jsm\\curve.dat", {Real, Real}];

ganglion = 
      ReadList["d:\\p\\math\\1998\\jsm\\ganglion.dat", {Real, Real}];

LoessPlot[ganglion, 0.5, 1];

LoessPlot[curve,  0.3,  2,  50, PlotStyle -> {PointSize[0.01]}];

Symmetric Errors Case

In the robust case it is only assumed that the errors are symmetric.  Hence for data ((Subscript[x, i]), (Subscript[y, i])),  i = 1, 2,   ...  ,  n, the basic model that can be fit may be written,

(Subscript[y, i]) = g((Subscript[x, i])) + (Subscript[ϵ, i]),  i = 1,   ...  , 
       n.

where (Subscript[e, i]) ~ IID( 0,  (σ^2)) with a symmetric distribution and g(x) is again a local linear polynomial of degree λ ≥ 0.

The model is fit to each of the n observed data values using a Iteratively Reweighted Least Squares (IRWLS) algorithm.  Each iteration fits the model using the current weights which are the product of the robustness weights times the local neighborhood weights.   So for fitting at the point (Subscript[x, j]) the weights are given by (Subscript[w, i])((Subscript[x, j])) (Subscript[r, j]),  j = 1,   ...  ,  n
      .

An additional set of weights are used to provide robustness.  Given the residuals from the last fit, \(TraditionalForm\`\(ϵ\&^\)\_i, \ i = 1, \  ... \ , \ n\), the robustness weights are then defined by

\(TraditionalForm\`r\_i = B(\(ϵ\&^\)\_i\/\(6\ s\))\)

where \(TraditionalForm
      \`s = median {\( | \(ϵ\&^\)\_1 | \), \  ... \ , 
            \( | \ \(ϵ\&^\)\_n | \)}\) and B(z) is the bisquare function,

B(z) =({{(((1 - (z^2)))^2), | z | ,  ≤ 1}, {0, | z | ,  > 1}})

Robust loess fits are often used for preliminary investigations of a dataset to guard against outliers and for the purpose of detecting outliers.  In various types of residual analysis plots with many types of other models such as in linear models it is useful to use a robust loess curve.  For example in a residual dependency plot the residuals are plotted against some covariate.  A trend in the loess curve in this plot would indicate that the model could be improved.

rubber = 
      ReadList["c:\\p\\math\\1998\\jsm\\rubber.dat", {Real, Real, Real}];

!! "c:\\p\\math\\1998\\jsm\\rubber.dat"

!! "c:\\p\\math\\1998\\jsm\\rubber.txt"

{hardness, tensile, abrasion} = Transpose[rubber];

Fit[rubber,  {1, h, t}, {h, t}]

fits = % /. {h -> hardness,  t -> tensile};

res = abrasion - fits;

?RobustLoessPlot

RobustLoessPlot[Transpose[{hardness,  res}], 1, 1, 
      FrameLabel -> {"hardness", "residual"}];

Salinity Data Interpolation

salinity = 
      ReadList["c:\p\datasets\gemaichen\salinity.dat", {Real,  Real}];

References

Chambers, J.M. and Hastie, T. (1991).  Statistical Models in S.  Wadsworth & Brooks/Cole.

Cleveland, W.S. (1995).  Visualizing Data.  New Jersey: Summit.

Cleveland, W.S. and Devlin, S.J. (1988). Locally-weighted regression: an approach to regression analysis by local fitting.  Journal of the American Statistical Association 83, 596-610.

Loader, C. (1998, forthcoming).   Local Regression.  New York: Springer-Verlag.