[R-sig-Geo] Regression Kriging cross validation

[rubenfcasal]

Hello Michele,

     Universal kriging is equivalent to Linear Regression (with the 
generalized-least-squaresestimator) + Simple Kriging of residuals (e.g. 
Cressie, 1993, section 3.4.5). The differences you observe are probably 
due to the use of ordinary least squares. If you use (leave-one-out) 
cross-validation with krige.cv (considering the UK model), the trend is 
also re-estimated at each prediction location. From my point of view, 
this would be the recommended way to proceed.
     As far as I know, there are no available implementations of the 
procedure you are suggesting.

     Best regards,
         Rubén.


El 29/04/2014 13:33, Michele Fiori escribió:
> Hi everyone,
> I am working on rainfall interpolation using regression kriging method and I
> need suggestions on how I can carry out a cross validation (leave-one-out)
> for this elaboration. At first I tried to apply directly Krige.cv, similarly
> to UK method (example for october: PP10uk.cv <- krige.cv(reg, prec2,
> PP10.vgm)), but unfortunately when I applied Universal Kriging on the same
> data, I realized that UK map was a little different from RK map.
> So my question is: How could I manage universal kriging in order to make it
> equivalent to regression kriging and use the above cross-validation, or is
> there
> another different method to apply cross validation (leave-one-out) on
> Regression Kriging interpolation?
> Below my code:
> Many thanks
>
> Michele Fiori
>
> ARPAS - Environmental Protection Agency of Sardinia
> MeteoClimatic Department - Meteorological Service
>
> Viale Porto Torres 119 - 07100 Sassari, Italy
> Tel + 39 079 258617
> Fax + 39 079 262681
> www.sardegnaambiente.it/arpas
>
> ####	Creating SpatialPixelDataFrame ("dem" - 250x250 m grid)
> 	....
> ####	Loading Precipitation data
>     	prec2 <- read.table("prec2.txt", sep="\t", header =TRUE)
>    	coordinates(prec2) <- c("x", "y")
>     	proj4string(prec2) <- CRS("+init=epsg:32632")
> ####	Linear regression Model
>    	mod.gen <- lm(PP10 ~ QUOTA_MARE + UTM_EST + UTM_NORD + DIST_MARE,
> prec2)
>    	step1 <- stepAIC(mod.gen, direction="both")
>    	reg <- formula(step1)
>     	PP10.lm <- lm(reg, prec2)
>    	summary(PP10.lm)
>     	prec2$residuals <- residuals(PP10.lm)
>      	dem$predlm <- predict(PP10.lm, dem)
> ####	Variogram of residuals
> 	PP10.vgm <- vgm(nugget=51.46, model="Sph", range=38038.89,
> psill=86.44)
> ####	Ordinary Kriging of residuals
> 	PP10.okr <- krige(PP10.lm$residuals ~ 1, prec2, dem, PP10.vgm,
> maxdist=Inf)
> 	dem$varokr <- PP10.okr$var1.pred
> ####	Regression Kriging (Linear Regression + Ordinary Kriging of
> residuals)
> 	dem$vark <- dem$predlm + dem$varokr
> ####	Universal kriging
> 	PP10.uk <- krige(reg, prec2, dem, PP10.vgm, maxdist=Inf)
> 	dem$varuk <- PP10.uk$var1.pred
> 	dem$difference <- dem$vark - dem$varuk
> 	spplot(dem[c("difference")], col.regions=terrain.colors(25),
> contour=FALSE, cuts = 15)
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>