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

[Sandeep Pulla]

​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]]