[R-sig-Geo] R: Regression Kriging cross validation

Michele Fiori mfiori at arpa.sardegna.it
Thu May 15 17:23:29 CEST 2014


Thank you for your kind reply
therefore as I have used the Osl method for regression, my result will never
match the universal kriging; However, in order to validate my method, I'm
trying to implement in the script a calculation loop witch runs n times (the
number of stations) regression + kriging without one station at a time.
Thank you again
Michele


-----Messaggio originale-----
Da: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org]
Per conto di rubenfcasal
Inviato: martedì 29 aprile 2014 19:49
A: r-sig-geo at r-project.org
Oggetto: Re: [R-sig-Geo] Regression Kriging cross validation

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
>

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



More information about the R-sig-Geo mailing list