[Gstat-info] recovering trend coefficients for GLS

[Edzer Pebesma]

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 
> estimates?
>   
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.

Best wishes,
--
Edzer

# how to get the BLUE trend coefficients out of a predict.gstat call?
# prepare data
library(gstat)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
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.)
newdat
out = as.data.frame(predict(g, newdat, BLUE = TRUE))[,3:4]
rownames(out) = c("Intercept", "slope")
colnames(out) = c("BLUE", "Variance")
out