[R-sig-Geo] R: R: Regression Kriging cross validation
rubenfcasal
rubenfcasal at gmail.com
Fri May 16 14:32:22 CEST 2014
Hi Michele,
Kriging methods assume that the variogram is known, although it is not
in practice. With your procedure, the whole prediction method is
reproduced (adding the extra-variability due to the estimation of the
variogram), so a better approximation to the variability of the
prediction error should be (in principle) obtained.
Excuse the rush...
Ruben Fernandez-Casal
El 16/05/2014 9:50, Michele Fiori escribió:
> 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 at bbox["x",])^2 + diff(prec2.i at 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:bakare at ualberta.ca]
> Inviato: giovedì 15 maggio 2014 17:39
> A: Michele Fiori
> Cc: rubenfcasal; r-sig-geo at 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 <mfiori at 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-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
>
> _______________________________________________
> 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