[R-sig-Geo] R: R: Regression Kriging cross validation
Michele Fiori
mfiori at arpa.sardegna.it
Fri May 16 09:50:29 CEST 2014
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