[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