# NOT RUN {
# generate response
y = rnorm(10)
# generate coordinates
x1 = runif(10); x2 = runif(10)
# data frame for observed data
data = data.frame(y, x1, x2)
coords = cbind(x1, x2)
d = as.matrix(dist(coords))
psill = 2 # partial sill
r = 4 # range parameter
evar = .1 # error variance
fvar = .1 # add finescale variance
# one can't generally distinguish between evar and fvar, but
# this is done for illustration purposes
# manually specify an exponential covariance model
v = psill * exp(-d/r) + (evar + fvar) * diag(10)
mod_man = cmod_man(v = v, evar = evar)
# coordinate names
cnames = c("x1", "x2")
# geolm for universal kriging
gearmod_uk = geolm(y ~ x1 + x2, data = data, mod = mod_man,
coordnames = cnames)
# newdata must have columns with prediction coordinates
# add 5 unsampled sites to sampled sites
newdata = data.frame(x1 = c(x1, runif(5)), x2 = c(x2, runif(5)))
newcoords = newdata[, cnames]
# create vop and vp using distances
dop = geodist(as.matrix(coords), as.matrix(newcoords))
dp = geodist(newcoords)
# manually create cross-covariance and covariance for
# prediction locations
vop = psill * exp(-dop/r) + fvar * (dop == 0)
vp = psill * exp(-dp/r) + fvar * diag(nrow(newcoords))
# prediction for universal kriging, with conditional simulation,
# using manual covariance matrices
pred_uk_man = predict(gearmod_uk, newdata, nsim = 2,
vop = vop, vp = vp, dmethod = "svd")
# do the same thing, but using cmod_std
# prediction for universal kriging, with conditional simulation
mod_std = cmod_std("exponential", psill = psill, r = r,
evar = evar, fvar = fvar)
gearmod_uk2 = geolm(y ~ x1 + x2, data = data, mod = mod_std,
coordnames = c("x1", "x2"))
pred_uk_std = predict(gearmod_uk2, newdata, nsim = 2, dmethod = "svd")
# compare results
all.equal(pred_uk_man$pred, pred_uk_std$pred)
all.equal(pred_uk_man$mspe, pred_uk_std$mspe)
# }
Run the code above in your browser using DataLab