[R] polynomial fitting

[Prof Brian Ripley]

On 29 Apr 2003, Suchandra Thapa wrote:

> I'm trying to find a way to fit a polynomial of degree n in x  and y to
> a set of x, y, and z data that I have and obtain the coefficients for
> the terms of the fitted polynomial.  However, when I try to use the
> surf.ls function I'm getting odd results.  
> 
> > x <- seq(0, 10, length=50)
> > y <- x
> > f <- function (x, y) {x^2 + y}  
> > library(spatial)
> > test <- data.frame(x=x, y=y, z=f(x, y))
> > test.kr <- surf.ls(2, test)
> > test.kr$beta
> [1] 0 0 0 0 0 0
> 
> When I try the example from the help I get:
> 
> >      library(MASS)
> >      data(topo, package="MASS")
> >      topo.kr <- surf.ls(2, topo)
> > topo.kr$beta
> [1] 801.217617 -11.018887  68.229148 -73.992968   3.343573   8.342717
> 
> Why is my test data causing problems?  

They lie on a line in the x-y plane, so there is collinearity in the fit
and you cannot fit a 2-dimensional trend surface uniquely.

> Also it seems that the beta
> attribute from the object returned the surf.ls correspond with the terms
> of the fitted polynomial.  If this is correct, in what order are the
> coefficients for the fitted polynomial given?  

It is not fully correct, as there is scaling of the data going on before 
fitting the polynomial.

> Finally, the R source
> code for the spatial library indicate that both surf.ls and surf.gls are
> limited to polynomials of degree 6 or below.  Is there another function
> that will work with higher order polynomials.  

You can use the R function poly() to fit orthogonal polynomials using lm
(provided you have genuinely two-dimensional data).

> lm(z ~ poly(x,y, degree= 2), data=topo)
Coefficients:
              (Intercept)  poly(x, y, degree = 2)1.0  
                   827.07                     -20.54  
poly(x, y, degree = 2)2.0  poly(x, y, degree = 2)0.1  
                   167.27                    -337.40  
poly(x, y, degree = 2)1.1  poly(x, y, degree = 2)0.2  
                    67.26                      20.13  


Reminder: spatial is part of the VR bundle, whose DESCRIPTION file says it
is software to support a book.  It is not intended to be comprehensive,
just sufficient for the techniques covered in that book.  I have never
seen a real application of trend surfaces of degree more than 6.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595