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:

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


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.fit = fit.variogram(vrg, model = vgm(1, "Sph", 1, 1))

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?


1 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.

  • Well, that was a great "partial answer" actually. Worked just fine for me. Thank you so much. Oct 4, 2019 at 3:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.