6

I'm trying to interpolate the values of an attribute (called "disease_rate") that i've sampled at a set of spatial points and i'd like to add 'weights' to the resulting kriging estimates.

Basically, I have about 1500 points (a random sample of 50 is shown below). I have an attribute vector (the rate of a particular disease) and a vector of absolute population sizes ("population"). I'd like to interpolate about 10,000 new points for "disease_rate", and at the same time weight the estimates by absolute population size (as areas with higher absolute populations have absolutely more people with the disease, and should probably have greater influence on the kriging estimates).

I'm using the packages "automap" and "gstat", with the following code:

library(automap)
vgm <- autofitVariogram(disease_rate ~ 1, , input_data = dat)
plot(vgm)

kr.ok <- autoKrige(disease_rate ~ 1, input_data = dat, model = "Exp",
    verbose = TRUE, maxdist = Inf) 
plot(kr.ok)

I'm stuck on two points:

1) how to add the weights? Is there an argument for this within the autoKrige() or krige() functions? I don't see a way to do it.

2) when using autokrige() from the package "automap" a set of 5000 point locations for interpolated values is generated by default, by using I believe a convex hull. In "gstat", how do I get a similar sample of point locations? At the moment I'm doing this:

prediction_locations <- spsample(SpatialPolygonsDataFrame, n = 10000,
    type = "regular")
prediction_locations@coords <- jitter(prediction_locations@coords)

the object "SpatialPolygonsDataFrame" is too big to include here, but provides a reference area within which to sample the points. However, I always get the following error:

"chfactor.c", line 130: singular matrix in function LDLfactor()
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
  LDLfactor

Any help/suggestions for the above two problems is much appreciated.

# sample of data (50 points)
dat <- new("SpatialPointsDataFrame"
, data = structure(list(population = structure(c(13L, 37L, 9L, 25L, 5L, 
2L, 40L, 31L, 10L, 43L, 22L, 7L, 32L, 19L, 12L, 28L, 23L, 30L, 
44L, 6L, 39L, 1L, 4L, 35L, 15L, 16L, 20L, 21L, 27L, 24L, 18L, 
41L, 42L, 26L, 14L, 45L, 47L, 33L, 36L, 3L, 46L, 4L, 4L, 11L, 
29L, 38L, 48L, 34L, 8L, 17L), .Label = c("10323", "1070", "1142", 
"1218", "1226", "1228", "12858", "1583", "1619", "1645", "1737", 
"1832", "1944", "2003", "2125", "2140", "2248", "273", "2816", 
"2868", "2986", "3093", "322", "3338", "3707", "371", "3841", 
"4151", "4184", "4189", "4198", "434", "4528", "491", "546", 
"560", "65", "654", "7138", "735", "776", "82", "845", "879", 
"91027", "943", "948", "989"), class = "factor"), disease_rate = c(0, 
0, 154.4163064, 140.2751551, 228.3849918, 0, 0, 71.46260124, 
0, 0, 226.3174911, 107.5091726, 0, 159.8011364, 0, 120.4529029, 
0, 45.35688709, 227.5312856, 0, 77.59122357, 33.5322916, 0, 0, 
0, 93.45794393, 239.0914525, 0, 0, 0, 0, 0, 0, 0, 0, 189.7113215, 
0, 0, 0, 175.1313485, 0, 369.4581281, 0, 0, 103.5691523, 0, 0, 
0, 0, 111.2099644)), .Names = c("population", "disease_rate"), class = "data.frame",     row.names = c(351L, 
251L, 477L, 349L, 56L, 434L, 299L, 248L, 301L, 262L, 115L, 325L, 
402L, 238L, 73L, 376L, 454L, 470L, 6L, 359L, 473L, 233L, 448L, 
 121L, 35L, 253L, 324L, 369L, 492L, 463L, 400L, 365L, 245L, 157L, 
264L, 101L, 367L, 142L, 109L, 389L, 158L, 254L, 457L, 68L, 219L, 
395L, 450L, 380L, 474L, 199L))
, coords.nrs = c(4L, 3L)
, coords = structure(c(634873.420766408, 541790.625187116, 649805.450219927, 
536333.978570239, 581957.093009378, 607202.361088042, 617240.615621174, 
505627.220834258, 620543.999367496, 592414.016643991, 579995.994396474, 
636851.153664358, 593891.033349101, 594650.473384273, 607103.681876253, 
659509.454474378, 604222.534339966, 592766.490473358, 625515.430661865, 
545043.782766906, 663659.275576225, 626566.467170621, 687406.871345609, 
579464.084578232, 594577.343169983, 650081.936891944, 661853.556699763, 
620028.294881262, 622881.614147092, 714608.702621971, 620991.249432193, 
531188.877201674, 654835.610827768, 527241.805893339, 604428.976031172, 
688457.08368477, 613888.674693582, 616783.067305461, 500082.617682936, 
609515.575391672, 611609.283708876, 689412.437797476, 675561.054104031, 
601548.525491549, 656484.375080421, 677341.365311871, 616862.413097344, 
644065.67844776, 475677.631939326, 617193.359923705, 5144526.77927643, 
5328308.33322789, 5216468.95632244, 5278378.21530125, 5242407.33311755, 
5107204.72128753, 5104424.13808425, 5346313.09064684, 5291125.71228228, 
5260872.67506553, 5296414.86191041, 5153372.04174738, 5107141.20462901, 
5126046.0142604, 5188835.84829657, 5241495.91305329, 5120799.1233737, 
5077612.68066607, 5173630.04242822, 5351027.28834978, 5231260.93265272, 
5321245.1139323, 5239845.72123907, 5257533.19157413, 5093198.30043019, 
5176787.1924563, 5208692.02318177, 5301924.71404294, 5307513.36239273, 
5170246.64046221, 5129845.77701215, 5284862.55527842, 5147875.41139468, 
5295710.09388113, 5293345.39302848, 5189299.58825685, 5165031.11113162, 
5136185.92439977, 5334176.39641511, 5166895.32812857, 5279273.17144778, 
5236974.38111571, 5222457.79639944, 5279521.25697787, 5171514.42644342, 
5244570.31251435, 5113000.12179841, 5170833.89650194, 5346485.72226854, 
5068270.63886515), .Dim = c(50L, 2L), .Dimnames = list(NULL, 
c("long.proj", "lat.proj")))
, bbox = structure(c(475677.631939326, 5068270.63886515, 714608.702621971, 
5351027.28834978), .Dim = c(2L, 2L), .Dimnames = list(c("long.proj", 
"lat.proj"), c("min", "max")))
, proj4string = new("CRS"
, projargs = NA_character_
)
)

closed as too broad by Andre Silva, BERA, Jochen Schwarze, xunilk, Kersten Dec 19 '18 at 13:13

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    There seem to be some misconceptions here. First, why should population have any correspondence at all with disease rate? Second, the BLUP criteria for kriging determine the weights. If you were to force different weights, you wouldn't be doing kriging. Thus it makes no sense for kriging to accept weights as an argument. It sounds like you might need a slightly more sophisticated approach, such as a kriging-based GLM like that offered in GeoRGLM – whuber Mar 21 '12 at 16:44
2

Maybe it is too late but this error:

"chfactor.c", line 130: singular matrix in function LDLfactor()
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
LDLfactor

It is sometimes related with some points at nearly the same position.
Did you try zerodist(dat) to check that?

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