Skip to content

cross-validation

8 messages · Paul Hiemstra, Marta Rufino, Edzer Pebesma

Original
Marta Rufino
# 
Dear list members,

Is it possible to do cross-validation on multivariate kriging?
i.e. apply krige.cv to multivariate kriging in gstat?

thank you very much,
Best wishes,
Marta
Edzer Pebesma
# 
Yes, Martha; function gstat.cv works on multivariable objects; do read 
it's documentation.

Best regards,
--
Edzer
Marta Rufino wrote:
Marta Rufino
# 
Hello,

yes, I know it is suppose to do it, but I could not find how, because it 
gives me an error... for example:

require(gstat); require(lattice)
data(meuse)
coordinates(meuse) = ~x + y
data(meuse.grid)
gridded(meuse.grid) = ~x + y

meuse.g <- gstat(id = "zn", formula = log(zinc) ~ 1, data = meuse)
meuse.g <- gstat(meuse.g, "cu", log(copper) ~ 1, meuse)

meuse.g <- gstat(meuse.g, model = vgm(1, "Sph", 900, 1), fill.all = T)
x <- variogram(meuse.g, cutoff = 1000)
meuse.fit = fit.lmc(x, meuse.g)
plot(x, model = meuse.fit)
z <- predict(meuse.fit, newdata = meuse.grid)
spplot(z) #map
gstat.cv(meuse.g) #does not work...
gstat.cv(meuse.g, remove.all=T) #either
gstat.cv(meuse.g, all.residuals=T) #either
gstat.cv(object=meuse.g, formula = log(zinc) ~ 1, data = meuse, model = 
vgm(1, "Sph", 900, 1), nmax=40, verbose=F) #either :-(

#                             
# Intrinsic Correlation found. Good.
# [using ordinary cokriging]

# "chfactor.c", line 130: singular matrix in function LDLfactor()
# Error in predict.gstat(object, newdata = data[sel, ], ...) :
#         LDLfactor

Maybe an example on the help file would be nice (eheheh).. I
What am I missing?


Thank you very much in advance,
Marta
Paul Hiemstra
# 
Hi,

You should check if you have duplicate observations, duplicate 
observations lead to a singular matrix. Use the function zerodist() to 
check where the observations are and remove.duplicates() to remove them.

cheers,
Paul

Marta Rufino schreef:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
That was not the problem, the problem was that you used meuse.g instead 
of meuse.fit to pass on to gstat.cv. For meuse.g, you have perfect 
correlation between Cu and Zn, so that collocated observations (meaning 
a Zn and a Cu observation at each obs location) act as a duplicate in 
univarite kriging.

Try:

out = gstat.cv(object=meuse.fit, nmax=40, verbose=F, nfold=5)
--
Edzer

        
Paul Hiemstra wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Marta Rufino wrote:

            

      
      
      
    

        
Hi Marta,

I now have in the example section of gstat.cv help:

# multivariable; thanks to M. Rufino:
meuse.g <- gstat(id = "zn", formula = log(zinc) ~ 1, data = meuse)
meuse.g <- gstat(meuse.g, "cu", log(copper) ~ 1, meuse)
meuse.g <- gstat(meuse.g, model = vgm(1, "Sph", 900, 1), fill.all = T)
x <- variogram(meuse.g, cutoff = 1000)
meuse.fit = fit.lmc(x, meuse.g)
out = gstat.cv(meuse.fit, nmax = 40, nfold = 5)
summary(out)
# mean error, ideally 0:
mean(out$residual)
# MSPE, ideally small
mean(out$residual^2)
# Mean square normalized error, ideally close to 1
mean(out$zscore^2)
# correlation observed and predicted, ideally 1
cor(out$observed, out$observed - out$residual)
# correlation predicted and residual, ideally 0
cor(out$observed - out$residual, out$residual)

I'm a bit hesitant to write a dedicated function for such simple 
calculations. In your function below, I object to the use of "expected", 
where you mean "desired". It's definitely not expected in the 
statistical sense.

The difficult questions about the cross validation results in 
geostatistics are usually "can we attribute the difference of a MSNE of 
1.1 from the idealized value of 1 attribute to sampling error, or is it 
an indication of a misfitting model?" In the latter case: "should we 
worry?" (only if kriging errors are of importance, many people ignore them).

Also, if the data are nicely spread, say the shortest distance is 100m 
(stay with the example above), then CV is never going to tell anything 
about the variogram for distances up to 100m.

Leave-on-out vs. n-fold? I believe Hastie, Tibshirani and Friedman 
promote the 10-fold. The idea was that leave-one-out may be too 
optimistic and not reveal model error, as refitted models are almost 
identical for moderately sized or large data sets, in each fold. 
Leave-one-out will have smaller RMSE than say 5-fold, but this is not an 
indication of a better model nor of a substantially better procedure, IMO.

Best wishes,
--
Edzer

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Marta Rufino
    

    

      
      #
    

  


  

      
Excellent!

It was exactly this what I meant.
I also agree that it is not needed a function to calculate such simple 
stuff, but certainly would be desirable (and effective) to have it in 
the help file from krige.cv (to clarify the less statistically minded 
like myself :-))
I totally agree with 'desirable', instead of expected ;-)

About the 10-fold issue (thank you for the reference. It was exactly 
what I have heard about), maybe it would nice if the default was to do 
10-fold, because it is also substantially quicker to estimate. hihihi

In relation to the use of CV, I found an article that brought my 
attention to the subject again, where they compared the the mean 
estimated by kriging and sample mean through cross-validation, as a way 
to show how kriging is improving the estimates. (Sales MH, Souza JCM, 
Kyriakidis PC, Roberts DA, Vidal E (2007) Improving spatial distribution 
estimation of forest biomass with geostatistics: A case study for 
Rondonia, Brazil. Ecol Model 205:221-230). Thus this is one of the 
issues (for many users, I believe...) how to quantify/prove that kriging 
is actually doing better than other methods?
Does any one has any further alternatives?


Best wishes and thank you very much once again,
Marta

        
Edzer Pebesma wrote: