[R-sig-Geo] KED trend coefficients in Gstat with R

[Tomislav Hengl]

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