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
