[GIS] Backtransformation of kriging predictions and variances

gstatkrigingr

I'm not sure how to back-transform log-normal kriged results. This example using the "meuse" data shows how to make variogram and use it to get kriging predictions (and variances) using the popular 'gstat' package of R. The last few lines show backtransformation from log-space to original concentrations just using the 'exp()' function. I read some articles confirming biased outcome of this approach. I was wondering if anyone exactly knows how to convert the kriged predictions and variances to unbiased original mean concentrations and standard deviations using R. 

library(gstat)
library(sp)

data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)

#fit variogram
v = variogram(log(zinc)~1, data=meuse, meuse.grid)
m <- vgm(0.59,"Sph", 897, 0.05)
plot(v$dist, v$gamma, ylim=c(0,0.75), xlim=c(0,1600), xaxs="i", yaxs="i")
lines(variogramLine(m, maxdist=1600), lwd=2, lty=2, col=4) 

# do the prediction by kriging
gridded(meuse.grid) = ~x+y
a <- krige(log(zinc)~1, meuse, meuse.grid, m)

#back transformation to original concentration space?
#predictions
exp(a$var1.pred)
#standard deviations
exp(sqrt(a$var1.var))

Best Answer

When the prediction error of a model based on log(Y) is Gaussian like here, you normally should use the correction of Laurent (1963) to estimate back to the original scale.
In your case:
MeanY = exp(a$var1.pred + 0.5 * (a$var1.var))
and:
SdY = exp(2*a$var1.pred + a$var1.var)*(exp(a$var1.var)-1)

You can verify on yourself with simulation of a Gaussian distribution:

# Gaussian distribution
logY <- rnorm(100000, 0, 0.5)
mean(logY)
exp(mean(logY))
# Corrected Mu
exp(mean(logY) + 0.5 * var(logY))
# Corrected Var
exp(2*mean(logY) + var(logY))*(exp(var(logY))-1) 

# Verification with log-Gaussian distribution
Y <- exp(logY)
# Correct Mu
mean(Y)
# Correct Var
var(Y)
Related Question