Op 26-10-2010 9:17, Edzer Pebesma schreef:
On 10/26/2010 09:02 AM, eric.ja...@ujf-grenoble.fr wrote:
Hello,
I would like to know if there is a way to get the trend coefficents when
doing kriging with external drift. I would like to compare with the
coefficients that I can have when doing temperature = f(altitude) for a
sample of stations.
Thanks for your answer(s).
Eric
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Eric, try
library(gstat)
demo(blue)
shows you one way how to do this. It is a bit tricky, because you're
using prediction equations (that disregard the residual, if BLUE=TRUE)
to get them. But then, in the end if your regression model is
beta0 * Intercept + beta1 * altitude,
setting Intercept to 0 and altitude to 1 gives you beta1, with its
estimation error.
Another, maybe more direct way could be find by using the mixed models
in package nlme.
You might want to take a look at how geoR package does this. geoR does
print by default the regression coefficients (which I think any RK/KED
system should always do) e.g.:
> zinc.rvgm2 <- likfit(zinc.geo, lambda=0,
trend=step.zinc$call$formula[c(1,3)],
+ messages=FALSE, ini=c(var(residuals(step.zinc)),500),
cov.model="exponential")
> zinc.rvgm2
likfit: estimated model parameters:
beta0 beta1 beta2 beta3 beta4 beta5
" 5.6919" " -0.4028" " -0.1203" " -0.0176" " 0.0090" " -0.3199"
tausq sigmasq phi
" 0.0526" " 0.1894" "499.9983"
Practical Range with cor=0.05 for asymptotic range: 1498
likfit: maximised log-likelihood = -975
See also p.143 in [http://spatial-analyst.net/book/GstatIntro]
T. Hengl
http://www.wewur.wur.nl/popups/vcard.aspx?id=HENGL001
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
