R gstat krige() NA predictions increasing with larger nmax

Question

I've been playing around with a King County house dataset posted on Kaggle (link: https://www.kaggle.com/harlfoxem/housesalesprediction) and tried to predict house prices with Kriging. The entire coding is shown below:

library(gstat)
library(sp)
library(ggplot2)
setwd("~/Documentos/R/Kriging")
dados = read.csv('kc_house_data.csv')

dados$price_sqft = dados$price / dados$sqft_living

ggplot(dados, aes(x = long, y = lat, size = price_sqft)) + 
  geom_point(alpha = 0.2) + scale_size_continuous(range = c(0.5, 5))

set.seed(1337)

treino_ind <- sample(seq_len(nrow(dados)), size = 0.80*nrow(dados))
dados_treino = dados[treino_ind,]
dados_teste = dados[-treino_ind,]

coordinates(dados_treino)  = ~ long + lat
coordinates(dados_teste) = ~ long + lat

vrg = variogram(log(price_sqft)~1, (dados_treino))

vrg

vrg.fit = fit.variogram(vrg, model = vgm(1, "Sph", 1, 1))
vrg.fit

plot(vrg, vrg.fit)

vrg.krige = krige(log(price_sqft)~1, dados_treino, dados_teste, 
      model = vrg.fit, nmax = 1)

At first, I ran the code without the nmax = 1 bit, but saw some posts where the inclusion of some restriction on the amount of calculations was advised and finally added it. The thing with this is that if I increase the nmax parameter, kringe() function returns many NA. For instance, when nmax = 3, there was 216 NAs. With nmax = 10, there was 766 NAs. With nmax = 50, there was 2568 NAs (to give a little perspective, the entire test set is composed by about 4300 obs).

My question is: does that actually make sense or is it a bug?

Answer 1

I think this is because you have duplicated observation locations. If you jitter the coordinates slightly then the NAs are no longer returned.

This your data and it fails:

> vrg.krige = krige(log(price_sqft)~1, dados_treino, dados_teste, model=vrg.fit, nmax=3)
[using ordinary kriging]
There were 50 or more warnings (use warnings() to see the first 50)

Lets try making a tiny shift in all the training data:

> trainj = SpatialPointsDataFrame(jitter(coordinates(dados_treino),factor=0.2), dados_treino@data)

and all the test data...

> testj = SpatialPointsDataFrame(jitter(coordinates(dados_teste),factor=0.2), dados_teste@data)

Then kriging with nmax=3....

> vrg.krige = krige(log(price_sqft)~1, trainj, testj, model=vrg.fit, nmax=3)
[using ordinary kriging]
> 

gives no errors. I thought gstat would print warnings about colocated points, but I don't see any when running through your code.

zerodist is handy for finding coincident points:

> zerodist(dados_teste)
      [,1] [,2]
 [1,]  785  786
 [2,]  834  867
 [3,]  187 1496
 [4,]  905 1524

I don't understand why increasing nmax produces more NAs, since the first nearest neighbour of a point should always be a coincident point if it exists. It may be because as the nmax increases its hitting coincident points in the training data and picking up a pair of points rather than one.

So this is a partial answer which might help you on your way.