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
KED trend coefficients in Gstat with R
3 messages · eric.jabot at ujf-grenoble.fr, Edzer Pebesma, Tomislav Hengl
On 10/26/2010 09:02 AM, eric.jabot at 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 at 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.
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.52north.org/geostatistics e.pebesma at wwu.de
Op 26-10-2010 9:17, Edzer Pebesma schreef:
On 10/26/2010 09:02 AM, eric.jabot at 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 at 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