## ----setup, include=FALSE, cache=FALSE-----------------------------------
# set global chunk options
## this code is only required to set up knitr for reproducible documentation
## it has nothing to do with the exercise itself
## please ignore!
opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'), out.format='latex', use.highlight=TRUE)
opts_chunk$set(fig.path='kgraph/ck', fig.align='center', fig.show='hold', prompt=FALSE, comment=NA, fig.width=6, fig.height=6, out.width='.7\\linewidth', size='scriptsize')
# options to be read by formatR
options(replace.assign=TRUE,width=70)

## ----echo=FALSE, results='hide'------------------------------------------
rm(list=ls())

## ----ck1, child='CoKrigeR_1.Rnw'-----------------------------------------

## ----set-parent-1, echo=FALSE, cache=FALSE-------------------------------
set_parent('CoKrigeR.Rnw')

## ------------------------------------------------------------------------
library(sp)
library(gstat)
?meuse
data(meuse)
str(meuse)
?meuse.grid
data(meuse.grid)
str(meuse.grid)

## ----fig.width=10, fig.height=5------------------------------------------
require("lattice")
h1 <- histogram(~ lead , meuse, xlab="Pb", col="thistle3", nint=12)
h2 <- histogram(~ log10(lead), meuse, xlab="log10(Pb)", col="thistle3", nint=12)
print(h1, split = c(1,1,2,1), more=T)
print(h2, split = c(2,1,2,1), more=F)
rm(h1, h2)

## ------------------------------------------------------------------------
        # proportions higher than various thresholds
ph <- function(level) {
    round(100 * sum(meuse$lead > level)/length(meuse$lead), 1)
}
p <- NULL; lvls <- c(600, 300, 200, 100, 30)
for (l in lvls) p <- c(p, ph(l))
        # display a table of thresholds and proportions
(data.frame(cbind(level=lvls, percent.higher=p)))
rm(ph, l, lvls, p)

## ------------------------------------------------------------------------
h1 <- histogram(~ om, meuse, xlab="OM", col="lightblue", nint=12)
h2 <- histogram(~ zinc , meuse, xlab="Zn", col="red4", nint=12)
h3 <- histogram(~ log10(om), meuse, xlab="log10(OM)", col="lightblue", nint=12)
h4 <- histogram(~ log10(zinc) , meuse, xlab="log10(Zn)", col="red4", nint=12)
print(h1, split = c(1,1,2,2), more=T)
print(h3, split = c(2,1,2,2), more=T)
print(h2, split = c(1,2,2,2), more=T)
print(h4, split = c(2,2,2,2), more=F)
rm(h1, h2, h3, h4)

## ------------------------------------------------------------------------
meuse.pb <- meuse[ seq(1, length(meuse$lead), by=3),
                   c("x", "y", "lead", "om", "zinc")]
str(meuse.pb)

## ------------------------------------------------------------------------
rownames(meuse.pb)

## ------------------------------------------------------------------------
meuse.pb <- cbind(meuse.pb,
                  ltpb = log10(meuse.pb$lead),
                  ltom = log10(meuse.pb$om),
                  ltzn = log10(meuse.pb$zinc))
str(meuse.pb)

## ------------------------------------------------------------------------
meuse.extra <- meuse[setdiff(rownames(meuse), rownames(meuse.pb)),
                 c("x", "y", "lead")]
meuse.extra <- cbind(meuse.extra, ltpb = log10(meuse.extra$lead))
str(meuse.extra)

summary(log10(meuse$lead)); sd(log10(meuse$lead))
summary(meuse.pb$ltpb); sd(meuse.pb$ltpb)
summary(meuse.extra$ltpb); sd(meuse.extra$ltpb)

## ------------------------------------------------------------------------
class(meuse)
coordinates(meuse) <- ~ x + y  
    # alternate command format: coordinates(meuse) <- c("x", "y")
coordinates(meuse.pb) <- ~ x + y  
coordinates(meuse.extra) <- ~ x + y  
coordinates(meuse.grid) <- ~ x + y  
class(meuse)

## ------------------------------------------------------------------------
summary(meuse.pb)

## ------------------------------------------------------------------------
str(as.data.frame(meuse))

## ------------------------------------------------------------------------
xyplot(y ~ x, as.data.frame(meuse), asp="iso",
       panel = function(x, ...) {
         panel.points(coordinates(meuse),
                 cex=1.8*(log10(meuse$lead) - 1.3),
                 pch=1, col="blue");
         panel.points(coordinates(meuse.pb),
                 cex=1.8*(meuse.pb$ltpb - 1.3),
                 pch=20, col="red");
         panel.grid(h=-1, v=-1, col="darkgrey")
       })

## ----out.width='.5\\linewidth'-------------------------------------------
plot(v.ltpb.c <- variogram(ltpb  ~ 1, data=meuse.pb, cutoff=1800, cloud=T))

## ----out.width='.5\\linewidth'-------------------------------------------
plot(v.ltpb <- variogram(ltpb  ~ 1, data=meuse.pb, cutoff=1800, width=200), pl=T)

## ----out.width='.55\\linewidth'------------------------------------------
m.ltpb <- vgm(0.08,"Sph",800,0.03)
plot(v.ltpb, pl=T, model=m.ltpb)

## ----out.width='.5\\linewidth'-------------------------------------------
(m.ltpb.f <- fit.variogram(v.ltpb, m.ltpb))
plot(v.ltpb, pl=T, model=m.ltpb.f)
rm(v.ltpb.c, v.ltpb, m.ltpb)


## ----ck1a, child='CoKrigeR_1a.Rnw'---------------------------------------

## ----set-parent-1a, echo=FALSE, cache=FALSE------------------------------
set_parent('CoKrigeR.Rnw')

## ----krige-naive---------------------------------------------------------
        # interpolate
k.o <- krige(ltpb ~1, locations=meuse.pb, newdata=meuse.grid, model=m.ltpb.f)
        # summary statistics
summary(k.o)

## ----krige-naive-plot, out.width='\\linewidth', fig.height=6, fig.width=6, out.width='0.9\\linewidth', fig.align='center'----
source("ck_plotfns.R")
plot.kresults(k.o, "var1", meuse, meuse.pb, "lead", "log10(Pb), OK")

## ----krige-naive-extra-pts-----------------------------------------------
        # predict at the extra points
k <- krige(ltpb  ~ 1, meuse.pb, meuse.extra, m.ltpb.f)
        # compute and summarize evaluation errors
summary(k)
diff <- k$var1.pred - meuse.extra$ltpb
summary(diff)
sqrt(mean(sum(diff^2)))   # RMSE (precision)
mean(diff)           # mean error (bias)
median(meuse.extra$ltpb)         # median error

## ----bubble-extra-pts, out.width='\\linewidth', fig.height=6, fig.width=6----
diff.df <- as.data.frame(diff)
coordinates(diff.df) <- coordinates(meuse.extra)
bubble(diff.df, zcol="diff", pch=1,
       main="OK evaluation errors at undersampled points, log10(Pb)")

## ----xval, results='hide'------------------------------------------------
cv.o <- krige.cv(ltpb ~ 1, meuse.pb, model=m.ltpb.f, nfold=nrow(meuse.pb),
                 verbose=FALSE)

## ----xval-summary--------------------------------------------------------
summary(cv.o)
res <- as.data.frame(cv.o)$residual
sqrt(mean(res^2))
mean(res)
mean(res^2/as.data.frame(cv.o)$var1.var)

## ----xval-summary-plot, out.width='\\linewidth', fig.height=6, fig.width=12----
plot.valids(k, "var1", meuse.extra, "ltpb", cv.o, "OK")

## ----val-xval-errors-hist, out.width='0.5\\linewidth'--------------------
diff <- 10^(k$var1.pred) - 10^(meuse.extra$ltpb)
summary(diff)
histogram(diff, col="darkseagreen2", nint=12, type="count",
    main="Validation errors", xlab="Pb, ppm")
diff <- 10^(cv.o$var1.pred) - 10^(meuse.pb$ltpb)
summary(diff)
histogram(diff, col="lightblue2", nint=12, type="count",
    main="Cross-Validation errors", xlab="Pb, ppm")


## ----ck2, child='CoKrigeR_2.Rnw'-----------------------------------------

## ----set-parent-2, echo=FALSE, cache=FALSE-------------------------------
set_parent('CoKrigeR.Rnw')

## ----fig.height=6, fig.width=6, out.width='0.5\\linewidth'---------------
xyplot(ltpb ~ ltom,  data=meuse.pb@data, pch=20, cex=1.2,
       col="blue", ylab="log10(Pb)", xlab="log10(OM)")
with(meuse.pb@data, cor(ltom, ltpb))
sum(is.na(meuse.pb$om))
with(meuse.pb@data, cor(ltom, ltpb, use = "complete"))

## ----echo=FALSE, results='hide'------------------------------------------
tmp <- round(100*with(meuse.pb@data, cor(ltom, ltpb, use = "complete"))^2,2)

## ----fig.height=6, fig.width=6, out.width='0.5\\linewidth'---------------
        # all valid covariable observations, with coordinates
meuse.co <- subset(as.data.frame(meuse), !is.na(om), c(x, y, om))
        # add log10-transformed variables for convenience
meuse.co <- cbind(meuse.co, ltom = log10(meuse.co$om))
str(meuse.co)
        # convert to spatial object
coordinates(meuse.co) <- ~ x + y
        # experimental variogram
v.ltom <- variogram(ltom ~ 1, meuse.co, cutoff=1800)
plot(v.ltom, pl=T)
        # model by eye
m.ltom <- vgm(.035, "Sph", 800, .015)
        # fit
(m.ltom.f <- fit.variogram(v.ltom, m.ltom))
plot(v.ltom, pl=T, model=m.ltom.f)
        # compare variogram structure to target variable
m.ltom.f$range[2]; m.ltpb.f$range[2]
round(m.ltom.f$psill[1]/sum(m.ltom.f$psill),2)
round(m.ltpb.f$psill[1]/sum(m.ltpb.f$psill),2)

## ------------------------------------------------------------------------
(g <- gstat(NULL, id = "ltpb", form = ltpb ~ 1, data=meuse.pb))
(g <- gstat(g, id = "ltom", form = ltom ~ 1, data=meuse.co))

## ----fig.height=12, fig.width=12, out.width='\\linewidth'----------------
v.cross <- variogram(g)
str(v.cross)
plot(v.cross, pl=T)

## ------------------------------------------------------------------------
(g <- gstat(g, id = "ltpb", model = m.ltpb.f, fill.all=T))

## ----fig.height=12, fig.width=12, out.width='\\linewidth'----------------
(g <- fit.lmc(v.cross, g, fit.method=6, correct.diagonal=1.01))
plot(variogram(g), model=g$model)

## ------------------------------------------------------------------------
str(m.ltom.f)

## ------------------------------------------------------------------------
str(g, max.level = 1)  

## ------------------------------------------------------------------------
str(g$data, max.level = 1)  
str(g$model, max.level = 1)  

## ------------------------------------------------------------------------
str(g$data$ltpb, max.level = 1)  

## ------------------------------------------------------------------------
str(g$model$ltpb, max.level = 1)  

## ------------------------------------------------------------------------
g$model$ltom$psill - m.ltom.f$psill
sum(g$model$ltom$psill) - sum(m.ltom.f$psill)
sum(g$model$ltom$psill)

g$model$ltpb$psill - m.ltpb.f$psill
sum(g$model$ltpb$psill) - sum(m.ltpb.f$psill)
sum(g$model$ltpb$psill)


## ----ck2a, child='CoKrigeR_2a.Rnw'---------------------------------------

## ----set-parent-2a, echo=FALSE, cache=FALSE------------------------------
set_parent('CoKrigeR.Rnw')

## ------------------------------------------------------------------------
        # interpolate
k.c <- predict(g, meuse.grid)  
str(k.c)
        # summarize predictions and their errors
summary(k.c$ltpb.pred); summary(k.c$ltpb.var)

## ----out.width='\\linewidth', fig.height=6, fig.width=12, fig.align='center'----
plot.kresults(k.c, "ltpb", meuse, meuse.pb, "lead", "CK of Pb with OM covariable, ")

## ------------------------------------------------------------------------
        # predict at the extra points
k <- predict(g, meuse.extra)  
        # compute and summarize prediction errors
diff <- k$ltpb.pred - meuse.extra$ltpb
summary(diff)
sqrt(sum(diff^2)/length(diff))   # RMS error (precision)
sum(diff)/length(diff)               # mean error (bias)

## ----out.width='\\linewidth', fig.height=6, fig.width=12-----------------
diff <- as.data.frame(diff)
coordinates(diff) <- coordinates(meuse.extra)
bubble(diff, zcol="diff", pch=1,
       main="CK evaluation errors at undersampled points, log10(Pb)")

## ----results='hide'------------------------------------------------------
cv.c <- gstat.cv(g)

## ------------------------------------------------------------------------
str(cv.c)
summary(cv.c$residual)
sqrt(mean(cv.c$residual^2))
mean(cv.c$residual)
mean(cv.c$residual^2/cv.c$ltpb.var)

## ----out.width='\\linewidth', fig.height=6, fig.width=12-----------------
plot.valids(k, "ltpb", meuse.extra, "ltpb", cv.c, "CK, OM co-variable")


## ----ck3, child='CoKrigeR_3.Rnw'-----------------------------------------

## ----set-parent-3, echo=FALSE, cache=FALSE-------------------------------
set_parent('CoKrigeR.Rnw')

## ----out.width='\\linewidth', fig.height=12, fig.width=12, fig.align='center'----
plot.cf(k.o, k.c, "var1", "ltpb", meuse, meuse.pb, "pred",
            "log10(Pb), OK", "log10(Pb), CK with log10(OM)")

## ------------------------------------------------------------------------
summary(k.o$var1.pred); summary(k.c$ltpb.pred)
diff <- k.c$ltpb.pred - k.o$var1.pred
summary(diff); rm(diff)

## ----out.width='\\linewidth', fig.height=12, fig.width=12, fig.align='center'----
plot.cf(k.o, k.c, "var1", "ltpb", meuse, meuse.pb, "var",
            "log10(Pb), OK", "log10(Pb), CK with log10(OM)")

## ------------------------------------------------------------------------
summary(k.o$var1.var); summary(k.c$ltpb.var)
diff.var <- (k.c$ltpb.var - k.o$var1.var)
summary(diff.var)
round(100*sum(diff.var < 0)/length(diff.var))
rm(diff.var)

## ------------------------------------------------------------------------
kv.ok <- krige(ltpb ~ 1, meuse.pb, meuse.extra, m.ltpb.f)
kv.ck <- predict(g, meuse.extra, debug.level=0)  
kv.diff <- (kv.ck$ltpb.pred - kv.ok$var1.pred)
summary(kv.diff)
d <- SpatialPointsDataFrame(coordinates(kv.ok),
                            data=as.data.frame(kv.diff))
bubble(d, col = c("plum","seagreen"),
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
rm(kv.ok, kv.ck, kv.diff, d)  

## ------------------------------------------------------------------------
cv.diff <- data.frame(diff = cv.c$residual -
               cv.o$residual,
               better = (abs(cv.c$residual) < abs(cv.o$residual)))
summary(cv.diff$diff)
summary(cv.diff$better)
d <- SpatialPointsDataFrame(coordinates(cv.c),
                            data=as.data.frame(cv.diff$diff))
bubble(d,
        col=c("lightblue", "red4"),
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
rm(cv.diff)


## ----ck4, child='CoKrigeR_4.Rnw'-----------------------------------------

## ----set-parent-4, echo=FALSE, cache=FALSE-------------------------------
set_parent('CoKrigeR.Rnw')

## ----fig.height=6, fig.width=6, out.width='0.5\\linewidth'---------------
xyplot(ltpb ~ ltzn,  data=meuse.pb@data, pch=20, cex=1.2,
       col="red4", ylab="log10(Pb)", xlab="log10(Zn)")
with(meuse.pb@data,cor(ltzn, ltpb))

## ----fig.height=6, fig.width=6, out.width='0.5\\linewidth'---------------
        # all valid covariable observations, with coordinates
meuse.co <- as.data.frame(meuse)[, c("x", "y", "zinc")]
        # add log10-transformed variables for convenience
meuse.co <- cbind(meuse.co, ltzn = log10(meuse.co$zinc))
coordinates(meuse.co) <- ~ x + y
str(meuse.co)
        # experimental variogram
v.ltzn <- variogram(ltzn ~ 1, meuse.co, cutoff=1800)
plot(v.ltzn, pl=T)
        # model by eye
m.ltzn <- vgm(.11, "Sph", 1000, .02)
        # fit
(m.ltzn.f <- fit.variogram(v.ltzn, m.ltzn))
plot(v.ltzn, pl=T, model=m.ltzn.f)
        # compare variogram structure to target variable
m.ltzn.f$range[2]; m.ltpb.f$range[2]
round(m.ltzn.f$psill[1]/sum(m.ltzn.f$psill),2)
round(m.ltpb.f$psill[1]/sum(m.ltpb.f$psill),2)

## ----fig.height=6, fig.width=6, out.width='0.5\\linewidth'---------------
(g2 <- gstat(NULL, id = "ltpb", form = ltpb ~ 1, data = meuse.pb))
(g2 <- gstat(g2, id = "ltzn", form = ltzn ~ 1, data = meuse.co))

## ----fig.height=12, fig.width=12, out.width='\\linewidth'----------------
v.cross <- variogram(g2)
plot(v.cross, pl=T)

## ----fig.height=12, fig.width=12, out.width='\\linewidth'----------------
g2 <- gstat(g2, id = "ltpb", model = m.ltpb.f, fill.all=T)
(g2 <- fit.lmc(v.cross, g2, fit.method=6))
plot(variogram(g2), model=g2$model)

## ------------------------------------------------------------------------
        # zinc
g2$model$ltzn$psill - m.ltzn.f$psill
sum(g2$model$ltzn$psill) - sum(m.ltzn.f$psill)
sum(g2$model$ltzn$psill)
        # lead
g2$model$ltpb$psill - m.ltpb.f$psill
sum(g2$model$ltpb$psill) - sum(m.ltpb.f$psill)
sum(g2$model$ltpb$psill)

## ------------------------------------------------------------------------
g2$model$ltpb$psill - g$model$ltpb.$psill
sum(g2$model$ltpb$psill) - sum(g$model$ltpb$psill)

## ------------------------------------------------------------------------
k.c2 <- predict(g2, meuse.grid, debug.level=0)  
summary(k.c2$ltpb.pred); summary(k.c2$ltpb.var)

## ----out.width='\\linewidth', fig.height=6, fig.width=12, fig.align='center'----
plot.kresults(k.c2, "ltpb", meuse, meuse.pb, "lead", "CK with Zn covariable, ")

## ------------------------------------------------------------------------
        # predict at the extra points
k <- predict(g2, meuse.extra, debug.level=0)  
        # compute and summarize prediction errors
diff <- k$ltpb.pred - meuse.extra$ltpb
summary(diff)
sqrt(sum(diff^2)/length(diff))   # RMS error (precision)
sum(diff)/length(diff)               # mean error (bias)

## ----out.width='\\linewidth', fig.height=6, fig.width=12, fig.align='center'----
plot.valids(k, "ltpb", meuse.extra, "ltpb", cv.o, "CK, Zn co-variable")

## ------------------------------------------------------------------------
cv.c2 <- gstat.cv(g2, verbose=FALSE)
summary(cv.c2$residual)
sqrt(mean(cv.c2$residual^2))
mean(cv.c2$residual)
mean(cv.c2$residual^2/cv.c2$ltpb.var)

## ----out.width='\\linewidth', fig.height=6, fig.width=12, fig.align='center'----
plot.valids(k, "ltpb", meuse.extra, "ltpb", cv.o, "CK, Zn co-variable")

## ----out.width='\\linewidth', fig.height=12, fig.width=12, fig.align='center'----
plot.cf(k.o, k.c2, "var1", "ltpb", meuse, meuse.pb, "pred",
            "log10(Pb), OK", "log10(Pb), CK with log10(Zn)")

## ------------------------------------------------------------------------
summary(k.o$var1.pred); summary(k.c2$ltpb.pred)
diff <- k.c2$ltpb.pred - k.o$var1.pred
summary(diff); rm(diff)

## ------------------------------------------------------------------------
summary(k.o$var1.var); summary(k.c2$ltpb.var)
diff.var <- (k.c2$ltpb.var - k.o$var1.var)
summary(diff.var)

## ----out.width='\\linewidth', fig.height=12, fig.width=12, fig.align='center'----
plot.cf(k.o, k.c2, "var1", "ltpb", meuse, meuse.pb, "var",
            "log10(Pb), OK", "log10(Pb), CK with log10(Zn)")

## ------------------------------------------------------------------------
summary(k.o$var1.var); summary(k.c2$ltpb.var)
diff.var <- (k.c2$ltpb.var - k.o$var1.var)
summary(diff.var)

## ------------------------------------------------------------------------
kv.ok <- krige(ltpb ~ 1, meuse.pb, meuse.extra, m.ltpb.f)
kv.ck <- predict(g2, meuse.extra, debug.level=0)  
kv.diff <- data.frame(x = coordinates(kv.ck)[,"x"],
              y = coordinates(kv.ck)[,"y"],
              diff = kv.ck$ltpb.pred - kv.ok$var1.pred)
summary(kv.diff$diff)
coordinates(kv.diff) <- ~ x +y
bubble(kv.diff, z="diff",
        main="CK - OK predictions, validation points, co-variable Zn", 
        col = c("plum","seagreen"),
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
rm(kv.ok, kv.ck, kv.diff)  

## ------------------------------------------------------------------------
cv.diff <- data.frame(x = coordinates(cv.c2)[,"x"],
              y = coordinates(cv.c2)[,"y"],
              diff = cv.c2$residual -
               cv.o$residual, better = (abs(cv.c2$residual) < abs(cv.o$residual)));
str(cv.diff);
summary(cv.diff$better);
coordinates(cv.diff) <- ~ x +y
bubble(cv.diff, z="diff",
        main="CK - OK cross-validation residuals, co-variable Zn", 
        col = c("lightblue", "red4"),
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
rm(cv.diff)


## ----ck5, child='CoKrigeR_5.Rnw'-----------------------------------------

## ----set-parent-5, echo=FALSE, cache=FALSE-------------------------------
set_parent('CoKrigeR.Rnw')

## ----compare-all-preds---------------------------------------------------
        # common colour scale
range <- range(k.o$var1.pred, k.c$ltpb.pred, k.c2$ltpb.pred)
breaks <- seq(range[1], range[2], length=32)
p1 <- levelplot(var1.pred ~ x+y, as.data.frame(k.o), aspect="iso", 
    col.regions=bpy.colors(64), at=breaks,
    main="log10(Pb), OK predictions",
         panel = function(x, ...) {
            panel.levelplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")
})
p2 <- levelplot(ltpb.pred ~ x+y, as.data.frame(k.c), aspect="iso", 
    col.regions=bpy.colors(64), at=breaks,
    main="log10(Pb), CK predictions, OM covariable",
         panel = function(x, ...) {
            panel.levelplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")
})
p3 <- levelplot(ltpb.pred ~ x+y, as.data.frame(k.c2), aspect="iso", 
    col.regions=bpy.colors(64), at=breaks,
    main="log10(Pb), CK predictions, Zn covariable",
         panel = function(x, ...) {
            panel.levelplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")
})

## ----plot-all-preds,  out.width='\\linewidth', fig.height=6, fig.width=18, fig.align='center'----
print(p1, split=c(1,1,3,1), more=T)
print(p2, split=c(2,1,3,1), more=T)
print(p3, split=c(3,1,3,1), more=F)
rm(range, breaks)

## ----compare-all-pred-errors---------------------------------------------
        # common colour scale
range <- range(k.o$var1.var, k.c$ltpb.var, k.c2$ltpb.var)
breaks <- seq(range[1], range[2], length=32)
p1 <- levelplot(var1.var ~ x+y, as.data.frame(k.o), aspect="iso", 
	col.regions=cm.colors(64), at=breaks, colorkey=T,
	main="log10(Pb), OK errors",
         panel = function(x, ...) {
            panel.levelplot(x, ...);
            panel.points(coordinates(meuse), col="green", pch=20, cex=.6);
            panel.points(coordinates(meuse.pb), col="red", pch=20, cex=.8);
            panel.grid(h=-1, v=-1, col="darkgrey")
})
p2 <- levelplot(ltpb.var ~ x+y, as.data.frame(k.c), aspect="iso", 
	col.regions=cm.colors(64), at=breaks, colorkey=T,
	main="log10(Pb), CK errors, OM covariable",
         panel = function(x, ...) {
            panel.levelplot(x, ...);
            panel.points(coordinates(meuse), col="green", pch=20, cex=.6);
            panel.points(coordinates(meuse.pb), col="red", pch=20, cex=.8);
            panel.grid(h=-1, v=-1, col="darkgrey")
})
p3 <- levelplot(ltpb.var ~ x+y, as.data.frame(k.c2), aspect="iso", 
	col.regions=cm.colors(64), at=breaks, colorkey=T,
	main="log10(Pb), CK errors, Zn covariable",
         panel = function(x, ...) {
            panel.levelplot(x, ...);
            panel.points(coordinates(meuse), col="green", pch=20, cex=.6);
            panel.points(coordinates(meuse.pb), col="red", pch=20, cex=.8);
            panel.grid(h=-1, v=-1, col="darkgrey")
})
rm(range, breaks)

## ----plot-all-pred-errors, out.width='\\linewidth', fig.height=6, fig.width=18, fig.align='center'----
print(p1, split=c(1,1,3,1), more=T)
print(p2, split=c(2,1,3,1), more=T)
print(p3, split=c(3,1,3,1), more=F)

## ----compute-all-errors--------------------------------------------------
     # compute all the evaluations and x-validations, and differences
k.ok <- krige(ltpb  ~ 1, meuse.pb, meuse.extra, m.ltpb.f)
diff.ok <- k.ok$var1.pred - meuse.extra$ltpb
k.ck <- predict(g, meuse.extra)  
diff.ck <- k.ck$ltpb.pred - meuse.extra$ltpb
k.ck.zn <- predict(g2, meuse.extra)  
diff.ck.zn <- k.ck.zn$ltpb.pred - meuse.extra$ltpb
cv.o <- krige.cv(ltpb ~ 1, meuse.pb, model=m.ltpb.f, nfold=nrow(meuse.pb),
                 verbose=FALSE)
cv.c <- gstat.cv(g, verbose=FALSE)
cv.c2 <- gstat.cv(g2, verbose=FALSE)
        # value that corresponds to cex=3
extreme <- max(abs(range(diff.ok, diff.ck, diff.ck.zn,
        cv.o$residual, cv.c$residual, cv.c2$residual)))
        # max. cex in each case proportional to its largest


## ----compare-all-errors--------------------------------------------------
      # display all plots
d <- SpatialPointsDataFrame(coordinates(k.ok), data=as.data.frame(diff.ok))
b.ok <- bubble(d,
         main="OK evaluation errors, log10(Pb)",
         maxsize = 3 * (max(abs(range(diff.ok))))/extreme,
         panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
d <- SpatialPointsDataFrame(coordinates(k.ck), data=as.data.frame(diff.ck))
b.ck <- bubble(d,
          main="CK evaluation errors, log10(Pb), co-variable log10(OM)",
         maxsize = 3 * (max(abs(range(diff.ck))))/extreme,
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
d <- SpatialPointsDataFrame(coordinates(k.ck.zn), data=as.data.frame(diff.ck.zn))
b.ck.zn <- bubble(d,
          main="CK evaluation errors, log10(Pb), co-variable log10(Zn)",
         maxsize = 3 * (max(abs(range(diff.ck.zn))))/extreme,
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
d <- SpatialPointsDataFrame(coordinates(cv.o), data=as.data.frame(cv.o$residual))
b.x.ok <- bubble(d,
         main="OK cross-validation errors, log10(Pb)",
         maxsize = 3 * (max(abs(range(diff.ok))))/extreme, col=c(4,5),
         panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
d <- SpatialPointsDataFrame(coordinates(cv.c), data=as.data.frame(cv.c$residual))
b.x.ck <- bubble(d,
          main="CK cross-validation errors, log10(Pb), co-variable log10(OM)",
         maxsize = 3 * (max(abs(range(diff.ck))))/extreme, col=c(4,5),
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
d <- SpatialPointsDataFrame(coordinates(cv.c2), data=as.data.frame(cv.c2$residual))
b.x.ck.zn <- bubble(d,
          main="CK cross-validation errors, log10(Pb), co-variable log10(Zn)",
         maxsize = 3 * (max(abs(range(diff.ck.zn))))/extreme, col=c(4,5),
          panel = function(x, ...) {
            panel.xyplot(x, ...);
            panel.grid(h=-1, v=-1, col="darkgrey")})
rm(k.ok, k.ck, k.ck.zn)
rm(diff.ok, diff.ck, diff.ck.zn, extreme)

## ----plot-all-errors, fig.width=12, fig.height=18, out.width='\\linewidth'----
print(b.ok, split=c(1,1,2,3), more=T)
print(b.ck, split=c(1,2,2,3), more=T)
print(b.ck.zn, split=c(1,3,2,3), more=T)
print(b.x.ok, split=c(2,1,2,3), more=T)
print(b.x.ck, split=c(2,2,2,3), more=T)
print(b.x.ck.zn, split=c(2,3,2,3), more=F)
rm(b.ok, b.ck, b.ck.zn, b.x.ok, b.x.ck, b.x.ck.zn)