IDW in R with cross validation

Question

I want to perform IDW interpolation using R, and more specifically, using the idw command from the gstat package. Also, I want to use the optimal number of neighbors (point locations) and distance power, which will be determined by what combination of "idp" and "nmax" produces the lowest RMSE in leave-one-out-cross validation. For this reason, after each run of IDW with different "nmax" and "idp" parameters, I perform the leave-one-out cross validation. In order to do this, I have written a double for loop in R that tests each "idp" parameter for each number of "nmax" neighbors, but my script can't produce the result I want.

power = seq(from = 1, to = 1.6, by = 0.1)
neigh = seq(from = 4, to = 5, by = 1)
# January
for(i in power) {
  for(j in 1:length(neigh)
    jan_idw = idw(Jan ~ 1, st, grd, nmax = power, idp = neigh)
    list_jan_idw = c(list_jan_idw,jan_idw )
    krige.cv(Jan ~ 1, st, nfold = 52)
  }
} 

For each run, I want the results to be stored in a data frame or matrix.

Answer 1

You could use a list to save all the krige.cv data.frames. I changed your example to the meuse data to make it reproducible. I'm assuming your data is similarly formatted.

require(gstat)
require(sp)
data(meuse)
data(meuse.grid)
coordinates(meuse) = ~x+y
coordinates(meuse.grid) = ~x+y 
gridded(meuse.grid) = TRUE

power = seq(from = 1, to = 1.6, by = 0.1)
neigh = seq(from = 4, to = 5, by = 1)
# January
results=list()
results.cv=list()
for (i in power) {
  for (j in neigh) {
    results[[paste0(i,"_",j)]] = idw (zinc ~ 1, meuse, meuse.grid, nmax = i, idp = j)
    results.cv[[paste0(i,"_",j)]] = krige.cv (zinc ~ 1, meuse, nfold = 52)
  }
} 

to get access to a specific list element, for example for power=1, neigh=4:

results[["1_4"]] # idw model

To see how the results are stored in the lists:

> str(results)
List of 14
 $ 1_4  :Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
  .. ..@ data       :'data.frame':  3103 obs. of  2 variables:
  .. .. ..$ var1.pred: num [1:3103] 634 713 654 604 857 ...
  .. .. ..$ var1.var : num [1:3103] NA NA NA NA NA NA NA NA NA NA ...
  .. ..@ coords.nrs : int [1:2] 1 2
  .. ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
  .. .. .. ..@ cellcentre.offset: Named num [1:2] 178460 329620
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. .. ..@ cellsize         : Named num [1:2] 40 40
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. .. .. ..@ cells.dim        : Named int [1:2] 78 104
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
  .. ..@ grid.index : int [1:3103] 69 146 147 148 223 224 225 226 300 301 ...
   ...

    >  str(results.cv)
List of 14
 $ 1_4  :Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  .. ..@ data       :'data.frame':  155 obs. of  6 variables:
  .. .. ..$ var1.pred: num [1:155] 796 738 606 481 369 ...
  .. .. ..$ var1.var : num [1:155] NA NA NA NA NA NA NA NA NA NA ...
  .. .. ..$ observed : num [1:155] 1022 1141 640 257 269 ...
  .. .. ..$ residual : num [1:155] 225.5 402.8 33.9 -223.6 -99.7 ...
  .. .. ..$ zscore   : num [1:155] NA NA NA NA NA NA NA NA NA NA ...
  .. .. ..$ fold     : int [1:155] 1 36 21 40 48 36 39 18 18 32 ...
  .. ..@ coords.nrs : num(0) 
  .. ..@ coords     : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
  .. .. .. ..$ : chr [1:2] "x" "y"
  .. ..@ bbox       : num [1:2, 1:2] 178605 329714 181390 333611
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. ..$ : chr [1:2] "min" "max"
  .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. ..@ projargs: chr NA

You can be more creative about how you name your list indices. You can also make a list of lists (e.g., results[[i]][[j]], in which case you would have to initialize each results list as a new list outside the inner loop, as in results[[i]]=list())