Best Combined Spatial Prediction
Line 10:  Line 10:  
<geshi lang=R lines=0>  <geshi lang=R lines=0>  
−  +  > data(SAGA_pal)  
−  +  
−  +  
> data(meuse)  > data(meuse)  
> coordinates(meuse) < ~x+y  > coordinates(meuse) < ~x+y  
−  # load grids:  +  > ## load grids: 
> data(meuse.grid)  > data(meuse.grid)  
> meuse.grid$X < meuse.grid$x  > meuse.grid$X < meuse.grid$x  
Line 21:  Line 19:  
> meuse.grid$soil < as.factor(meuse.grid$soil)  > meuse.grid$soil < as.factor(meuse.grid$soil)  
> coordinates(meuse.grid) < ~x+y  > coordinates(meuse.grid) < ~x+y  
−  > gridded(meuse.grid) <  +  > gridded(meuse.grid) < T 
−  > fullgrid(meuse.grid) <  +  > fullgrid(meuse.grid) < T 
−  # get values from grids to points:  +  > # get values from grids to points: 
−  > meuse.ov <  +  > meuse.ov < over(meuse, meuse.grid) 
−  >  +  > meuse.ov < cbind(as.data.frame(meuse), meuse.ov) 
−  +  
−  +  
</geshi>  </geshi>  
Line 52:  Line 48:  
# GLM prediction (all possible predictors):  # GLM prediction (all possible predictors):  
−  > m.glm < glm(zinc~dist+soil+ffreq  +  > m.glm < glm(zinc~sqrt(dist)+soil+ffreq, meuse.ov, family=gaussian(link=log)) 
−  >  +  > summary(m.glm) 
> zinc.glm < as(meuse.grid["soil"], "SpatialPointsDataFrame")  > zinc.glm < as(meuse.grid["soil"], "SpatialPointsDataFrame")  
> zinc.glm$var1.pred < p.glm$fit  > zinc.glm$var1.pred < p.glm$fit  
−  > zinc.glm$var1.var < p.glm$se.fit  +  > zinc.glm$var1.var < p.glm$se.fit^2 + p.glm[["residual.scale"]]^2 
−  > zinc.glm$svar <  +  > zinc.glm$svar < zinc.glm$var1.var/(m.glm$null.deviance/m.glm$df.null) 
−  > gridded(zinc.glm) <  +  > gridded(zinc.glm) < TRUE 
−  > fullgrid(zinc.glm) <  +  > fullgrid(zinc.glm) < TRUE 
# regressionkriging:  # regressionkriging:  
−  > vr.fit < fit.variogram(variogram(log1p(zinc)~sqrt(dist), meuse), vgm(1, "Exp", 300, 1))  +  > vr.fit < fit.variogram(variogram(log1p(zinc)~sqrt(dist)+soil+ffreq, meuse), vgm(1, "Exp", 300, 1)) 
−  > plot(variogram(log1p(zinc)~sqrt(dist), meuse), vr.fit)  +  > plot(variogram(log1p(zinc)~sqrt(dist)+soil+ffreq, meuse), vr.fit) 
−  > zinc.rk < krige(log1p(zinc)~sqrt(dist), meuse, meuse.grid, vr.fit)  +  > zinc.rk < krige(log1p(zinc)~sqrt(dist)+soil+ffreq, meuse, meuse.grid, vr.fit) 
−  > zinc.rk$zinc.pred < expm1(zinc.rk$var1.pred)  +  > zinc.rk$zinc.pred < expm1(zinc.rk$var1.pred) 
</geshi>  </geshi>  
Line 71:  Line 67:  
<geshi lang=R lines=0>  <geshi lang=R lines=0>  
+  > zinc.rk$zinc.pred < expm1(zinc.rk$var1.pred)  
> meuse.grid$zinc.id < zinc.id$var1.pred  > meuse.grid$zinc.id < zinc.id$var1.pred  
> meuse.grid$zinc.tr < zinc.tr$zinc.pred  > meuse.grid$zinc.tr < zinc.tr$zinc.pred  
Line 77:  Line 74:  
> meuse.grid$zinc.glm < zinc.glm$var1.pred  > meuse.grid$zinc.glm < zinc.glm$var1.pred  
> meuse.grid$zinc.rk < zinc.rk$zinc.pred  > meuse.grid$zinc.rk < zinc.rk$zinc.pred  
−  > spplot(meuse.grid[c("zinc.  +  > spplot(meuse.grid[c("zinc.ok", "zinc.glm", "zinc.rk", "zinc.id", "zinc.tr", "zinc.lm")], col.regions=SAGA_pal[[1]]) 
</geshi>  </geshi>  
Line 90:  Line 87:  
<geshi lang=R lines=0>  <geshi lang=R lines=0>  
> pc.preds < prcomp(~zinc.id+zinc.tr+zinc.lm+zinc.ok+zinc.glm+zinc.rk, scale=F, meuse.grid)  > pc.preds < prcomp(~zinc.id+zinc.tr+zinc.lm+zinc.ok+zinc.glm+zinc.rk, scale=F, meuse.grid)  
−  > biplot(pc.preds, cex=0.9, xlabs=rep(".", length(pc.preds$x[,1])))  +  > biplot(pc.preds, cex=0.9, xlabs=rep(".", length(pc.preds$x[,1])), xlim=c(0.06,0.025), ylim=c(0.04,0.03)) 
</geshi>  </geshi>  
Revision as of 20:45, 28 May 2013
patial prediction (also known as spatial interpolation) is a process of estimating values of some target variable at a new location, using sampled values and auxiliary data (maps) that can be used to explain variation of the target variable. A technique (algorithm) that generates spatial predictions given a set of inputs can be called a spatial predictor. A common objective of somebody interested to generate quality maps using geostatistics is to compare a variety of predictors and then select the one that is: (a) the most accurate (crossvalidation), (b) least biased, (c) most robust (sensitivity analysis), (d) most reliable, and (e) fastest possible. There are many spatial predictors described in the literature. The Spatial Interpolation Comparison game, for example, lists over 35 techniques. Sometimes you will discover that there are not as many possibilities. In fact, there are only dozen truly independent spatial predictors (e.g. mechanical, splinesbased, krigingbased, regressionbased and machine learningbased predictors). Many predictors produce fairly similar outputs (even if their names seem to indicate completely different algorithms).
This article demonstrates how to make a quick comparison of the most popular spatial predictors (inverse distances, ordinary kriging, regressionbased predictors) implemented in the gstat package (read more on obtaining and installing the software used below). A simple approach is suggested to generate the Best Combined Spatial Prediction using predictions produced variety of methods. This assumes that the predictors are mutually independent and that they address different features of interest.

For this purpose we will use the "meuse" dataset provided together with the gstat package:
UNIQ2d4bd38a4dba325ageshi00000000QINU
Once we have prepared all maps, we can proceed with running predictions. We consider six standard techniques: inverse distances (ID), regression on coordinates (TR), simple linear regression (LM), ordinary kriging (OK), GLMbased regression, and regressionkriging (RK):
UNIQ2d4bd38a4dba325ageshi00000001QINU
In principle, all techniques can be automated. We can finally plot all predictions next to each other:
UNIQ2d4bd38a4dba325ageshi00000002QINU
To compare various predictions (maps), we can use the Principal Component Analysis (PCA):
UNIQ2d4bd38a4dba325ageshi00000003QINU
This shows that ID and OK, and LM and GLM are highly correlated (as expected); RK lies inbetween (as expected); all predictors are positively correlated (of course).
Different predictors can be also compared by estimating the time needed to generate predictions for all grid nodes:
UNIQ2d4bd38a4dba325ageshi00000004QINU
In this case the difference is fairly small. In reality (large datasets) some simple local predictor can be few thousands time faster than global predictors.
The Best Combined Spatial Prediction (BSCP) can be produced by combining the predictions of multiple prediction techniques. Assuming that a series of prediction techniques are mutually independent (i.e. they do not use the same model parameters; they treat different parts of spatial variation), the BCSP can be generated as a weighted average from multiple predictions:
where is the prediction error estimated by the model (prediction variance), and is the number of predictors. For example, in the case study from above (meuse), we can generate a combined prediction using OK and GLMregression:
UNIQ2d4bd38a4dba325ageshi00000008QINU
Note that the prediction will in some parts of the study are look more as OK, in others more as GLM, which actually depicts extrapolation areas of both methods (this map is probably very similar to predictions produced using regressionkriging). In general, the map in the middle looks more as the GLMregression map because this map is about 23 times more precise than the OK map (read more about visualization of uncertainty in maps). It is important to emphasize that, in order to combine various predictors, we do need to have an estimate of the prediction uncertainty, otherwise we are not able to assign the weights. In principle, linear combination of statistical techniques using the equation above should be avoided if a theoretical basis exists that incorporates such combination.
In the example above (GLM+OK), we assume that the the predictions/prediction errors are independent, and they are probably not. In addition, a statistical theory exists that supports a combination of regression and kriging (see e.g. Diggle and Ribeiro, 2007), so there is no need to run the predictions separately and derive an unrealistic measure of model error. The BCSP can be only interesting for situations where there are indeed several predictors possible, where no theory exists that reflects their combination, and/or where fitting of individual models is faster and less problematic than fitting of a hybrid model.
The combined prediction error of a BCSP can be estimated by summing up simulated values using the same formula from above. This might be computationally intensive as we need at least 20 simulations to get a reasonable estimate of the additive prediction variance.
<Rating> Rate this article: 1 Poor and in error 2 Far from ready 3 Needs improvements 4 Minor corrections 5 Very useful!! </Rating>