# IDW in R with cross validation

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.

• doesn't krige.cv already give you a data.frame? If you want to save all the data.frames from your loops, use a `list`. Then you will have a list with length(power)*length(neigh) data.frames. Very simple. Jun 17, 2015 at 15:29

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()`)

• Thank you very much @aaryno. I run it and it came up with the following warnings: "In if (d\$nmax == Inf) nmax = as.integer(-1) else nmax = as.integer(d\$nmax) : the condition has length > 1 and only the first element will be used" Also, when I do result[[power[i]_neigh[j]] it always returns NULL. Jun 17, 2015 at 16:04
• `for` loops were not formulated correctly. I changed the code to use the example `meuse` so I could run it through. Jun 17, 2015 at 16:48
• Thank you very much @aaryno. This is exactly what I wanted the script to do! Jun 17, 2015 at 17:59
• I fail to see how the power and neighbors are passed to the krige.cv call. If I don't specify a model in krige.cv it defaults to IDW? Mar 14, 2017 at 13:32
• Just found a way to set idp in krige.cv call (gis.stackexchange.com/questions/190558/…) Mar 14, 2017 at 14:21