How do I implement 'weighted' kriging in R?

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