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

[Edzer Pebesma]

What is exactly wrong with the kriging variance of the block kriging
estimate, with as block the whole study area (besides that it gives the
mean, instead of the total resource)?

Of course you can never cross validate this in real world cases, but I'd
say you could with data simulated over a fine grid.

On 04/16/2012 04:49 PM, Burton Shank wrote:
> 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
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de