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-boun...@r-project.org [mailto:r-sig-geo-boun...@r-project.org] Per conto di rubenfcasal Inviato: martedì 29 aprile 2014 19:49 A: r-sig-geo@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@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
