swh:1:snp:33a53053e50f7abe7d281cc0c803be827debf4a3
Tip revision: c3530d114e0e58b870ff87ca29e22a54a7d5d5a7 authored by Edzer Pebesma on 14 March 2015, 18:40:41 UTC
version 1.0-22
version 1.0-22
Tip revision: c3530d1
blue.R
# how to get the BLUE trend coefficients out of a predict.gstat call?
# prepare data
library(sp)
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