[R-sig-Geo] Residual variogram from known trend

[Paul Hiemstra]

Hi,

You could try to (ab)use the standard lm function. First calculate the
regression using lm, than replace the fitted coefficients by the values
you want them to have and than use predict to estimate the values at
your measurement locations. Subtracting that value from the measurements
gives you the residuals. Than you can use variogram to create a
variogram. In code:

library(gstat)
loadMeuse()

meuse$logzinc = log(meuse$zinc)
lmfit = lm(log(zinc)~dist, meuse)
lmfit
lmfit$coefficients = c(10,-5)
lmfit
meuse$logzinctrend2 = predict(lmfit, meuse['dist'])
head(meuse at data[c('logzinc','logzinctrend')])
meuse$residuals = meuse$logzinc - meuse$logzinctrend
resvariogram = variogram(residuals~1, meuse)
plot(resvariogram)

cheers,
Paul

On 04/04/2011 02:02 PM, Tom Gottfried wrote:
> Dear list,
>
> is there a way in gstat to calculate a variogram from residuals of a
> known trend (that can not be estimated from the data) without passing
> the residuals to variogram() or gstat()? I tried the arguments
> trend.beta respectively beta. The first one gives an error, the second
> gives the same result as when I do not give it:
>
> library(gstat)
> data(meuse)
> # trying with variogram(..., trend.beta=...)
> var1 <- variogram(zinc~sqrt(dist), ~x+y, meuse)
> var2 <- variogram(zinc~sqrt(dist), ~x+y, meuse, trend.beta=c(1000,
> -1000))
>
> # trying with gstat(..., beta=...)
> coordinates(meuse) <- ~x+y
> meuse1 <- gstat(NULL, "zinc", zinc~sqrt(dist), meuse)
> var1 <- variogram(meuse1)
> meuse2 <- gstat(NULL, "zinc", zinc~sqrt(dist), meuse, beta=c(1000,
> -1000))
> var2 <- variogram(meuse2)
> identical(var1, var2) # both variograms are identical
>
> I see it's possible to calculate the residuals with
>
> residuals(lm(zinc~sqrt(dist), meuse))
>
> and pass them directly. But I wonder whether a more convenient direct
> passing of the trend coefficients is possible.
>
> Thanks!
> Tom
>


-- 
Paul Hiemstra, MSc
Global Climate Division
Royal Netherlands Meteorological Institute (KNMI)
Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39
P.O. Box 201 | 3730 AE | De Bilt
tel: +31 30 2206 494

http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770