[R-sig-Geo] gstat trend beta ignored?

[Sandeep Pulla]

Thanks for the clarifications, Edzer.

Actually I was passing the parameters ('beta' and 'gls'; the latter
inside 'set') to create gstat objects that are then passed on to
variogram() and predict(). Because variogram() respects 'set', I
assumed that would also be the case with predict() (the idea here was
to retrieve the trend predictions and beta that were used internally
by gstat prior to variogram estimation). If it is working as designed,
perhaps a note can be added to the documentation for predict()
indicating that GLS is always used? Alternately, and especially if
there is no other way to retrieve OLS trend predictions, perhaps
predict() could switch between using OLS/GLS depending on the 'gls'
flag?

Regarding documentation, I'm using gstat 1.0-19 for R. For 'set', the
user is referred to the gstat standalone user's manual (I'm using
v2.5.1 April 2014). I had interpreted the following to mean that
estimation by GLS or OLS is overridden when trend beta are specified:
   "set gls=1; use generalised least squares residuals instead of the
default ordinary least squares (OLS) residuals for sample variograms
or covariograms [0, use OLS or WLS residuals]"
   "trend.beta - vector with trend coefficients, in case they are
known. By default, trend coefficients are estimated from the data"
A minor change to the documentation would eliminate any doubt,
although it *could* be argued that always letting any specified beta
override trend estimation (regardless of GLS/OLS) would be less
surprising to a user.

Thanks,
Sandeep

On Thu, Aug 21, 2014 at 1:44 PM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
>
> Sandeep, before you ask whether certain bugs are known it is often more
> constructive to start the discussion about whether there is agreement
> that your findings point to bugs.
>
> As of your point 2, this is according to documentation: beta is not in
> the parameter list of predict.gstat, hence is absorbed by ..., and ...
> according to the documentation is ignored. If you expect, or wished for
> other behaviour, that would be a feature request.
>
> As of point 1, in the call to variogram() specifying gls=1 asks for
> generalized least squares _estimation_ of the trend coefficients,
> specifying beta means that you know, hence don't want to estimate these
> coefficients. This indicates some misunderstanding of what you could do,
> or what you expect the software to do. The right behaviour would be to
> generate an error message.
>
> If you could point me to the relevant sections in the documentation that
> made you believe that something useful should happen ("GLS + beta
> specified") then I would be happy to learn about that.
>
>
> On 08/21/2014 08:12 AM, Sandeep Pulla wrote:
> > There appear to be couple of problems in specifying known trend
> > coefficients in gstat. Specifically,
> > 1. The 'beta' parameter seems to be ignored when the 'gls' parameter is 1.
> > 2. The 'beta' parameter is ignored by predict.gstat() (regardless of the
> > 'gls' parameter).
> >
> > Here's a full example:
> >
> > library(gstat)
> > library(sp) # for coordinates()
> > library(geoR) # for grf()
> > n = 20 ^ 2 # number of observations
> > xmin = ymin = 0; xmax = ymax = 2 # rectangular domain
> > trend.beta = c(pi, exp(1), -exp(1)) # linear trend-surface coefficients
> > trend.surf = cbind(1, as.matrix(expand.grid(x = seq(xmin, xmax, len =
> > sqrt(n)),
> >   y = seq(xmin, ymax, len = sqrt(n))))) %*% trend.beta
> > set.seed(100)
> > field = grf(n, "reg", xlims = c(xmin, xmax), ylims = c(ymin, ymax),
> > cov.model = "mat",
> >   cov.pars = c(1, .1), kappa = 1, nugget = .5, mean = trend.surf)
> > dat = data.frame(var = field$data, field$coords, intercept = 1)
> > coordinates(dat) = ~ x + y
> > vgm.mod = vgm(1, "Mat", .5, 1)
> > # Trend surface predictions to retrieve beta (see demo(blue))
> > pred.grid = data.frame(intercept = c(1, 0, 0), x = c(0, 1, 0), y = c(0, 0,
> > 1))
> > coordinates(pred.grid) = ~ x + y
> > # OK: Trend beta are estimated using OLS
> > plot(plot(variogram(gstat(NULL, "var", var ~ x + y, data = dat,
> >   model = vgm.mod, set = list(gls = 0))), main = "OLS"))
> > # OK: Predefined trend beta are used
> > plot(plot(variogram(gstat(NULL, "var", var ~ x + y, data = dat,
> >   model = vgm.mod, set = list(gls = 0), beta = c(20, 10, -10))),
> >   main = "OLS + predefined"))
> > # OK: Trend beta are estimated using GLS
> > plot(plot(variogram(gstat(NULL, "var", var ~ x + y, data = dat,
> >   model = vgm.mod, set = list(gls = 1))), main = "GLS"))
> > # BUG: Predefined trend beta are ignored
> > plot(plot(variogram(gstat(NULL, "var", var ~ x + y, data = dat,
> >   model = vgm.mod, set = list(gls = 1), beta = c(20, 10, -10))),
> >   main = "GLS + predefined"))
> > # BUG: Predefined trend beta are ignored
> > print(predict(gstat(NULL, "var", var ~ x + y + intercept - 1, data = dat,
> >     model = vgm.mod, set = list(gls = 1), beta = c(20, 10, -10)),
> >   pred.grid, BLUE = T)@data$var.pred)
> > # BUG: 'gls' option is ignored (uses GLS instead of OLS for trend
> > estimation)
> > print(predict(gstat(NULL, "var", var ~ x + y + intercept - 1, data = dat,
> >     model = vgm.mod, set = list(gls = 0)),
> >   pred.grid, BLUE = T)@data$var.pred)
> >
> > Are these known bugs?
> >
> > Thanks,
> > Sandeep
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
> 83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>