Dear All, Finally I managed to develop this ad hoc procedure for the regression kriging cross validation; I tried to test the kriging part (the regression part has been already tested successfully) by comparison with the krige.cv function and I noticed that there are small discrepancies. I verified that they are due to the fact that in my case the variogram is recalculated for each step, while the function runs using the same initial variogram obtained on the complete set of data. What do you think about? Michele
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 #################### Cross validation REGRESSION KRIGING (leave on out) PP03.lm <- lm(reg, prec2) Nstazioni <- length(prec2$NOME_STAZ) for (i in 1:Nstazioni) { ### regression step for point “i” prec2.i <- prec2[-i,] md.i <- lm(reg, data=prec2.i) ypredlm.i <- predict(md.i,prec2[i,]) err[i] <- (ypredlm.i - prec2$PP03[i])^2 res[i] <- (ypredlm.i - prec2$PP03[i]) ### kriging step for residual point “i” PP03i.rsvar <- variogram(prec2.i$residuals~1, prec2.i) PP03i.ivgm <- vgm(nugget=0, model="Sph", range=sqrt(diff(prec2.i@bbox["x",])^2 + diff(prec2.i@bbox["y",])^2)/4, psill=var(prec2.i$residuals)) PP03i.rvgm <- fit.variogram(PP03i.rsvar, model=PP03i.ivgm) ypredok.i<- krige(md.i$residuals~1, prec2.i, prec2[i,], PP03i.rvgm) ypredok.i$var1.pred ### Regression + kriging predicted value for point “i” ypred.i <- ypredlm.i + ypredok.i$var1.pred err[i] <- (ypred.i - prec2$PP03[i])^2 ypred[i] <- ypred.i } err <-c(err) ypred <-c(ypred) ### mean squared error MSE <- sum(err)/(Nstazioni) MSE ### root mean square error RMSE.rk <- MSE^0.5 RMSE.rk Da: Moshood Agba Bakare [mailto:bak...@ualberta.ca] Inviato: giovedì 15 maggio 2014 17:39 A: Michele Fiori Cc: rubenfcasal; r-sig-geo@r-project.org Oggetto: Re: [R-sig-Geo] R: Regression Kriging cross validation Hi Michele, I have similar problem. I used ordinary kriging and inverse distance weighting method (IDW) to generate set of interpolated values from the same interpolation grid. I don't understand how cross validation can be done to come up with diagnostic statistics such as mse, rmse to use as basis for identifying the best interpolation method. I used krige.cv but I encountered error message. Please any advice on what to do please? ## Create grid for the interpolation through ordinary kriging and idw grid <- expand.grid(easting = seq(from = 299678, to = 301299, by=10), northing=seq(from = 5737278, to = 5738129, by=10)) ## convert the grid to SpatialPixel class to indicate gridded spatial data coordinates(grid)<-~easting + northing proj4string(grid)<-CRS("+proj=utm +zone=12 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0") gridded(grid)<- TRUE #### Ordinary kriging prok <- krige(id="yield",yield ~ 1, canmod.sp, newdata = grid, model=exp.mod,nmax=20,maxdist=33.0) ## Inverse Distance Weighting (IDW) Interpolation method yield.idw = idw(yield~1, canmod.sp, grid,nmax=20,maxdist=33.0,idp=1) Thanks Moshood On Thu, May 15, 2014 at 9:23 AM, Michele Fiori <mfi...@arpa.sardegna.it> wrote: 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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo