[Gstat-info] recovering trend coefficients for GLS
edzer.pebesma at uni-muenster.de
Sun Apr 27 20:46:51 CEST 2008
Ashton Shortridge wrote:
> Hello list,
> It's straightforward to generate GLS trend surfaces with gstat; an example is
> given in the man page for predict.gstat in R. This is useful for calculating
> residual variograms for universal kriging. The result is a set of predictions
> and variances, just like the output when kriging.
> However, I'm also interested in the trend model coefficients (betas), and
> conducting some diagnostics on that trend. I played with the debug.level
> settings but couldn't determine much from those. How can I obtain the beta
By manipulating the predictor values. If you have a linear trend in a
single external variable, pass the values c(0,1) as predictor values to
get the slope (or second coefficient); the variance then is the variance
corresponding to that coefficient.
The example below should work. The hack is ugly because it need to make
the intercept explicit instead of automatic, which is the usual case in
R, leading to the -1 in the formula's. If it's automatic, there is no
way of setting the corresponding predictor to 0 (if there is, tell me!).
This example is now in demo(blue) in the R package gstat, but might not
be in the CRAN version yet; in that case it'll be in the next one.
Note the ugliness of the hack, and note that it needs careful attention
to do it for more complex models. I might work it into a more generic
function one day, I can see the use; this is asked with some frequency.
# how to get the BLUE trend coefficients out of a predict.gstat call?
# prepare data
coordinates(meuse) = ~x+y
gridded(meuse.grid) = ~x+y
# create a manual, non-automatic intercept
meuse$Int = rep(1, length(meuse$zinc))
meuse.grid$Int = rep(1, length(meuse.grid$dist))
# create a gstat object without an automatic intercept:
g = gstat(formula = log(zinc)~-1+Int+sqrt(dist), data=meuse, model =
vgm(1, "Exp", 300))
newdat = meuse.grid[1:2,c("Int", "dist")]
gridded(newdat) = FALSE
newdat$dist = c(0,1)
newdat$Int = c(1,0)
# (in the more general case of n predictors, make sure that the matrix with
# predictors, or design matrix, equals the identity matrix.)
out = as.data.frame(predict(g, newdat, BLUE = TRUE))[,3:4]
rownames(out) = c("Intercept", "slope")
colnames(out) = c("BLUE", "Variance")
More information about the Gstat-info