[R] surf.ls

[Roger Bivand]

On Thu, 9 Dec 2004, m p wrote:

> My application of surf.ls is giving wrong results.
> CAn somebody point to the problem?

Your problem is of your own making: do read the book that this software 
supports, namely MASS, see ?surf.ls. I have already told you what it says 
about trend surface orders, now for some fairly standard trend surface 
numerical knowledge - the coefficients relate to coordinates that have 
been transformed to the -1,1 range, not to the input data range, because, 
when using say UTM coordinates typically for Northings in m north or south 
of the equator, y is in millions, y^2 is already very big indeed, and 
higher powers lead to meaningless fits. Try using lm() with UTM 
coordinates and you'll see what I mean. Prof. Ripley's code does this 
properly, and can be relied on to higher orders irrespective of the 
scaling of the coordinates. If you are in doubt, please read the source 
code. 

Beyond that, please do use the appropriate predict() function, which takes 
account of the change of range. Rolling your own can be good sometimes, 
but not unless you take account of all the underlying methods.

> Here is part of the code:
> 
> var is three columnd array filled with x,y,z values
> x,y are equally spaced
> zetaframe <- data.frame(x=var[,3],y=var[,2],z=var[,1])
> 
> #var is three columnd array filled with x,y,z values
> #part of zetaframe looks like
> 2211 10.000000  7.995732 -1.4149990
> 2212 10.000000  8.495732 -1.5649990
> 2213 10.000000  8.995732 -1.7249990
> 2214 10.000000  9.495732 -1.8849990
> 2215 10.000000  9.995732 -2.0549990
> 2216 10.000000 10.495730 -2.2249980
> 2217 10.000000 10.995730 -2.3949980
> 2218 10.000000 11.495730 -2.5749980
> 2219 10.000000 11.995730 -2.7549980
> 2220 10.000000 12.495730 -2.9349980
> 2221 10.000000 12.995730 -3.1249980
> 2222 10.000000 13.495730 -3.3149970
> 
> zeta.kr <- surf.ls(2,
> zetaframe$x,zetaframe$y,zetaframe$z)
> 
> #zeta.kr$beta =  
> #-2.0909907  0.9969686 -0.5092430 -2.0151540 
> 0.5047758 #-0.2155143
> #To plot  I fill an array with output from
> #zeta.kr for same x,y grid points
> 
> k <- 0
> for (i in 1:nx) {
>         for (j in 1:ny) {
>                 k <- k+1
>                 xi <- var[k,3]
>                 yi <- var[k,2]
>                 zetasurf[i,j] <-
>                 zeta.kr$beta[1]*xi^2 +
>                 zeta.kr$beta[2]*xi*yi +
>                 zeta.kr$beta[3]*yi^2 +
>                 zeta.kr$beta[4] +
>                 zeta.kr$beta[5]*xi +
>                 zeta.kr$beta[6]*yi
>         }
> }
> 
> 
> #When  I plot it with 
> 
> filled.contour(zetasurf,nlevels=11,
> color.palette=rainbow,
> xlab=xlabstring,ylab=ylabstring,cex.lab=1.2,cex.axis=1.2,
> xaxs = "i", yaxs = "i", las = 1,
> plot.axes={ axis(1); axis(2); points(10,10) })
> 
> # gives ridiculous results
> 
> Thanks,
> MArk
> 
> 
> 
> 
> --- Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 
> > On Thu, 9 Dec 2004, m p wrote:
> > 
> > > Hello,
> > > I am looking into description of surf.ls(spatial)
> > > and see under value $beta - the coefficients.
> > > When I use polynomial of degree 2 to fit surface
> > > I expect to get 4 coefficients:
> > > 
> > > z = a_1 x^2 + a_2 xy + a_3 y^2 + a_4
> > > 
> > > What do beta really stand for and why do I get
> > > $beta vector of length 6?
> > 
> > No, z = a_1 x^2 + a_2 xy + a_3 y^2 + a_4 + a_5 x +
> > a_6 y,
> > 
> > if you like, order 2 is linear + quadratic, see p.
> > 420 in MASS (4th 
> > edition), eq. 15.1:
> > 
> > f((x,y)) = \sum_{r+s \leq p} a_{rs} x^r y^s, 
> > 
> > with P = (p+1)(p+2)/2 coefficients for order p; for
> > p=2, P=6.
> > 
> > > 
> > > Thakns,
> > > Mark
> > > 
> > > ______________________________________________
> > > R-help at stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide!
> > http://www.R-project.org/posting-guide.html
> > > 
> > 
> > -- 
> > Roger Bivand
> > Economic Geography Section, Department of Economics,
> > Norwegian School of
> > Economics and Business Administration, Breiviksveien
> > 40, N-5045 Bergen,
> > Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
> > e-mail: Roger.Bivand at nhh.no
> > 
> > 
> >
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no