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

[Edzer Pebesma]

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140821/fe6680fa/attachment.bin>