[R-sig-Geo] Indicator Kriging R/gstat

[Tomislav Hengl]

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.ov at data <- cbind(meuse.ov at 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(meuse at bbox["x",])^2 +
diff(meuse at 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-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
> Of Edzer Pebesma
> Sent: Sunday, November 29, 2009 9:10 PM
> To: Pete Larson
> Cc: r-sig-geo at 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 at 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.pebesma at wwu.de
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo