Dear Pete, Edzer, If this is of any help, few years ago we have played with using auxiliary maps to interpolate categorical variables (multi-indicators?). The results are reported in:
Hengl T., Toomanian N., Reuter H.I., Malakouti M.J. 2007. Methods to interpolate soil categorical variables from profile observations: lessons from Iran. Geoderma, 140(4): 417-427. http://dx.doi.org/10.1016/j.geoderma.2007.04.022 Experiences were frustrating. Many (most of) classes come in isolated patches; hence getting something out from the variograms was difficult (close to pure nugget). In addition, if you go ahead and interpolate the 0/1 values (ignoring binomial properties), you will produce values <0 and >1 i.e. nonsense values. Then we tried with using memberships and things changed. First, memberships you can simply convert to logits, then interpolate and back-transform and then none of the values will be outside the 0-1 range. It was also much easier to fit the variograms. But we had a luxury that soils are fuzzy classes per definition and we had extra data to 'fuzzify' the A, B, C values to memberships (in fact we discovered that surveyors should have designated memberships from the beginning - this solves all problems of transition classes etc). Here is a similar type of approach (fit a GLM with a binomial model, interpolate the residuals using OK) with the meuse dataset: > data(meuse) > coordinates(meuse) <- ~x+y > data(meuse.grid) > coordinates(meuse.grid) <- ~x+y > gridded(meuse.grid) <- TRUE > fullgrid(meuse.grid) <- TRUE > meuse.ov <- overlay(meuse.grid, meuse) > meuse...@data <- cbind(meuse...@data, meuse[c("zinc", "lime")]...@data) # fit a GLM: > glm.lime <- glm(lime~dist+ffreq, meuse.ov, family=binomial(link="logit")) > step.lime <- step(glm.lime) # check if the predictions are within 0-1 range: > summary(round(step.lime$fitted.values, 2)) # convert to logits: > logits.lime <- log(step.lime$fitted.values/(1-step.lime$fitted.values)) # predict at all locations: > p.glm <- predict(glm.lime, newdata=meuse.grid, type="link", se.fit=T) # predictions as logits > str(p.glm) # attach spatial reference: > lime.glm <- as(meuse.grid, "SpatialPointsDataFrame") > lime.glm$lime <- p.glm$fit > gridded(lime.glm) <- TRUE > lime.glm <- as(lime.glm, "SpatialGridDataFrame") # fit a variogram for residuals: > lime.ivgm <- vgm(nugget=0, model="Exp", range=sqrt(diff(me...@bbox["x",])^2 + diff(me...@bbox["y",])^2)/4, psill=var(residuals(step.lime))) > lime.rvgm <- fit.variogram(variogram(residuals(step.lime)~1, meuse.ov), > model=lime.ivgm) # does not work - singular model; fix by hand: > lime.rvgm <- vgm(nugget=0.6, "Exp", psill=0.2, range=500) # interpolate residuals (logits): > lime.rk <- krige(residuals(step.lime)~1, meuse.ov, meuse.grid, lime.rvgm) # add predicted trend and residuals: > lime.rk$var1.rk <- lime.glm$lime + lime.rk$var1.pred # back-transform logits to the original response scale (0,1): > lime.rk$var1.rko <- exp(lime.rk$var1.rk)/(1+exp(lime.rk$var1.rk)) > spplot(lime.rk["var1.rko"], scales=list(draw=T), at=seq(0.05,1,0.05), col.regions=grey(rev(seq(0,0.95,0.05))), main="Liming requirements", sp.layout=list("sp.points", col="black", meuse)) > write.asciigrid(lime.rk["var1.rko"], "lime_rk.asc", na.value=-1) Another thing you should look at is the geoRglm package (it has a function "binom.krige" to work with binomial data), but as far as I remember - I could not get it running. Maybe you should dig into that and then contact the authors if you need more help. HTH, T. Hengl http://home.medewerker.uva.nl/t.hengl/ > -----Original Message----- > From: r-sig-geo-boun...@stat.math.ethz.ch > [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf > Of Edzer Pebesma > Sent: Sunday, November 29, 2009 9:10 PM > To: Pete Larson > Cc: r-sig-geo@stat.math.ethz.ch > Subject: Re: [R-sig-Geo] Indicator Kriging R/gstat > > Pete, for the universal kriging case I don't believe that indicator > theory allows what you want to do. See also > > https://stat.ethz.ch/pipermail/r-sig-geo/2009-May/005701.html > > If you disagree, could you please let me know why, or where you found > indications that this can be done? > > For ordinary kriging, usually values outside [0,1] are set back to their > nearest allowed value. The degree to which they will be outside [0,1] > also depends on the variogram used -- Gaussian variograms being > notorously suspect. > -- > Edzer > > Pete Larson wrote: > > Hello, > > > > I would like to do indicator kriging in R/gstat. I have dichotomous > > 0/1 data and have performed ordinary kriging and universal kriging, > > but get predctions that are far from 0 and 1. > > > > Am I doing something wrong? Here is the code I have been using: > > > > #### Ordinary Kriging with krige function ##### > > # ordinary kriging: > > x <- krige(z~1, ~x+y, model = v.fit, data = estand, newd = grd) > > > > ### universal block kriging: > > uk <- krige(z~x+y+DistHF+RiverDist+RiverDist2, ~x+y, model = v.fit, > > data = estand, newdata = > > grd) > > > > Any help would be appreciated. > > > > Pete > > > > _______________________________________________ > > R-sig-Geo mailing list > > R-sig-Geo@stat.math.ethz.ch > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > -- > Edzer Pebesma > Institute for Geoinformatics (ifgi), University of Münster > Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 > 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/ > http://www.springer.com/978-0-387-78170-9 e.pebe...@wwu.de > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo