[R-sig-Geo] pooled confidence interval on the total resource estimated by kriging

[Burton Shank]

I'm trying to calculate the confidence interval on the total resource
in a geographic area (i.e. get an estimate of the total quantity of
zinc in the region +/- a confidence interval).While this would seem
like a common product from kriging, this seems more complex than
expected from searching references.  Naturally, it is possible to
perform cross-validation to get  error estimates at sampled location
(i.e. using krige.cv in gstat).  However, I don't think this is useful
with highly aggregated sample sets with open spaces between sampled
areas.
For example, one can approximate sampling the meuse data along
buffered transects and examine the kriging predictions.

library(gstat); library(lattice);

data(meuse);
meuse=meuse[, c("x", "y", "zinc", "dist")]    ## trim dataset to useful columns
  meuse$logZinc=log(meuse$zinc)               ## log transform zinc

Transects=seq(from=330000, by=700, length=5); ## set positions for
sampling transects
Sampled=meuse[which(x=                        ## extract sample data
along transects
  (meuse$y>Transects[1]-100 & meuse$y< (Transects[1]+100) ) |
  (meuse$y>Transects[2]-100 & meuse$y< (Transects[2]+100) ) |
  (meuse$y>Transects[3]-100 & meuse$y< (Transects[3]+100) ) |
  (meuse$y>Transects[4]-100 & meuse$y< (Transects[4]+100) ) |
  (meuse$y>Transects[5]-100 & meuse$y< (Transects[5]+100) )
    ), ]

coordinates(meuse)=c("x", "y");
coordinates(Sampled)=c("x", "y");
plot(meuse)
points(Sampled, pch=20, cex=1.5)

W=W+10
Var=variogram(logZinc~sqrt(dist), Sampled, cutoff=1300, width=90);
plot(Var)## calculate variogram
Var.Fit=fit.variogram(Var, vgm(model="Sph", nugget=0.1, psill=0.3,
range=800)); ## fit variogram model
plot(Var, Var.Fit)

data(meuse.grid)  ## prepare prediction data set
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")

Pred=krige(logZinc~sqrt(dist), Sampled, meuse.grid, Var.Fit)#,
block=c(40,40)) ## predict by kriging
meuse.grid$Zinc.logPred=Pred$var1.pred; ## log-transformed predictions
meuse.grid$Zinc.Var=Pred$var1.var;  ## log-transformed variances
meuse.grid$Zinc.Pred.BT=exp(meuse.grid$Zinc.logPred+meuse.grid$Zinc.Var/2)
 ## bias-adjust and back-transform

trellis.par.set(sp.theme()) ## plot
spplot(meuse.grid["Zinc.logPred"], main="Kriging Prediction")
spplot(meuse.grid["Zinc.Var"], main="Kriging variance")

mean(meuse$zinc)
mean(meuse.grid$Zinc.Pred.BT)

The kriging variance shows a clear geographic pattern with respect to
the "transects". Thus, cross-validation is not useful for calculating
a confidence interval for the estimate for the entire geographic
region. Similarly, it is possible to get an estimate of the variance
of the mean from the likfit function in geoR, however, again, this can
not be applied to unsampled regions.
How would one calculate the pooled variance of the mean?

Any suggestions are appreciated.

thanks,
Burton