[STATSGRASS] gstat - trend over 3th polinominal

[Roger Bivand]

On Sat, 11 Nov 2006, Edzer J. Pebesma wrote:

> Jarek, we don't know what your "newdata" object is. It is a name of an 
> argument, but passed as data; I would expect something like
> 
> predict(g, newdata = meuse.grid)
> 
> I didn't know of that form in formulas, so thanks! Pretty much shorter 
> than writing the 4th order polynomial out.
> 
> It might work "in principle" but you may expect increasing inaccuracies 
> with higher order trends due to the large difference between say x^4 and 
> x. Rescaling the coordinates prior to the analysis (e.g. dividing by 
> 1000, going from m to km) might help. poly() doesn't help here, as it 
> will rescale differently for meuse and meuse.grid, afaik.

I have tried with:

unroll_trend <- function(np=1, x="x", y="y") {
    ip <- 0
    res <- character()
    if (np > 0) {
        for (i in (0:np)) {
            for (j in (0:(np - i))) {
                if ((i == 0) & (j == 0)) 
                    next
                else if (i == 0) {
                    ip <- ip + 1
                    res[ip] <- paste("I(", x, "^", j, ")", sep = "")
                } else if (j == 0) {
                    ip <- ip + 1
                    res[ip] <- paste("I(", y, "^", i, ")", sep = "")
                } else {
                    ip <- ip + 1
                    res[ip] <- paste("I(", x, "^", j, "*", 
                        y, "^", i, ")", sep = "")
                }
            }
        }
    }
    res
}

and:

> coefficients(lm(as.formula(paste("zinc ~ ", paste(unroll_trend(1), 
collapse="+"))), meuse))
  (Intercept)        I(x^1)        I(y^1) 
-2.244418e+04 -4.129057e-01  2.932104e-01 
> coefficients(lm(as.formula(paste("zinc ~ ", paste(unroll_trend(2), 
collapse="+"))), meuse))
  (Intercept)        I(x^1)        I(x^2)        I(y^1)  I(x^1 * y^1) 
 7.999712e+06  1.363638e+02  3.618219e-04 -1.223970e+02 -8.059965e-04 
       I(y^2) 
 4.040368e-04 
> coefficients(lm(as.formula(paste("zinc ~ ", paste(unroll_trend(3), 
collapse="+"))), meuse))
  (Intercept)        I(x^1)        I(x^2)        I(x^3)        I(y^1) 
 7.999712e+06  1.363638e+02  3.618219e-04            NA -1.223970e+02 
 I(x^1 * y^1)  I(x^2 * y^1)        I(y^2)  I(x^1 * y^2)        I(y^3) 
-8.059965e-04            NA  4.040368e-04            NA            NA 

but:

> lm(zinc ~ (x+y)^2, meuse)

Call:
lm(formula = zinc ~ (x + y)^2, data = meuse)

Coefficients:
(Intercept)            x            y          x:y  
 -2.994e+06    1.608e+01    9.254e+00   -4.974e-05  

Roger


> 
> I also wonder what, besides curiosity, the value of a 4-th or higher 
> order trend surface might be.
> --
> Edzer
> 
> Jarek Jasiewicz wrote:
> > Hi
> > I found that result of formula:
> >
> > gs=gstat(id="x", formula=zinc~(x+y)^2, data=meuse)
> > pgs=predict(gs, newdata, meuse.grid)
> >
> > is equvalent to:
> >
> >
> > gs=gstat(id="x", formula=zinc~x+y, data=meuse, degree=2)
> > pgs=predict(gs, newdata, meuse.grid)
> >
> >
> > so I tried:
> > gs=gstat(id="x", formula=zinc~(x+y)^4, data=meuse)
> > pgs=predict(gs, newdata, meuse.grid)
> >
> > and it worked!
> >
> > so it is possible to calculate trend over 3th polynominal (now its not 
> > a metter for what :))) in that way?
> >
> > regards
> > Jarek
> >
> >
> >
> > _______________________________________________
> > statsgrass mailing list
> > statsgrass at grass.itc.it
> > http://grass.itc.it/mailman/listinfo/statsgrass
> 
> _______________________________________________
> statsgrass mailing list
> statsgrass at grass.itc.it
> http://grass.itc.it/mailman/listinfo/statsgrass
> 

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