gstat trend beta ignored?
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>