[R-sig-Geo] [Gstat-info] LOCAL universal kriging with GLOBAL reg. coeff.estimation

[Edzer Pebesma]

Dave Depew wrote:
> A late follow up question to this thread....
>
> What if the local neighbor hood was restricted to the range of 
> autocorrelation? My impression was that values beyond that have little 
> weight in the interpolation anyways.
Often, yes. Three exceptions are (i) when the nugget forms a large part 
of the sill, or (ii) when you use a Gaussian variogram, or (iii) when 
you use a trend function.
> I realize that this is of course not statistically optimal, but for 
> those with large data sets (and lacking the computing power) would 
> this not give a reasonable approximation?
If the nugget is small, and you have say 500+ observations within the 
range, restricting it to a much smaller number may still be near-optimal.
--
Edzer
>
>
>
>
>
>
>
>
>
>
> >Dear Radim, Edzer,
>
> I was thinking about the same problem few years ago (I assume that you 
> work with auxiliary maps and
> not only coordinates).
>
> I think (have a feeling) that local and global Universal kriging 
> should be treated as two things
> (especially if you put a very narrow search radius). This is because, 
> in the case of KED algorithm,
> both regression and kriging part of predictions are solved at the same 
> time. Hence, if we limit the
> search window, but keep a constant variogram model, we could obtain 
> very different predictions then
> if we use a global (regression-kriging) model. Only if the variogram 
> of residuals if absolutely
> stationary (constant), then we can limit the search window to fit the 
> KED weights. In fact, I am
> confident that software packages should not allow UK/KED predictions 
> with a limited search radius,
> because this conflicts the model assumptions (a bit deeper discussion 
> about this problem can be
> found in my lecture notes section 2.2 "Localized versus local models").
>
> In the case of regression-kriging, this is less problematic because it 
> is computationally cheap to
> fit the regression part, and then you are allowed to limit the search 
> radius to interpolate the
> residuals (we still cheat but not so much). But then the problem is 
> that we do not have an estimate
> of the RK kriging variance (you could estimate the OLS prediction 
> variance and then OK-residuals
> prediction variance and then add them together - both can be done in 
> gstat - but then you completely
> ignore correlation of the two terms; although this is often 
> insignificant).
>
> If gstat already does something like this (global GLS simulations of 
> trend with local kriging), then
> this would be a legitime solution. But I still can not run much 
> geostatistics with big datasets
> (>>1000 points, >>10^6 grids) in R+gstat (on standard PCs), not to 
> mention geostatistical
> simulations (that are computationally much more demanding than 
> predictions)... if you would have to
> do it 100 times, this would really take a lot of processing.
>
> One alternative is to run UK/KED in SAGA software (I have tested it - 
> it runs about 5-10 times
> faster than R+gstat and has no memory limit problems), which is 
> possible through R also, e.g.:
>
> rsaga.get.usage("geostatistics_kriging", 3)
> rsaga.geoprocessor(lib="geostatistics_kriging", module=3,
> param=list(GRIDS="DEM.sgrd;SLOPE.sgrd;PLANC.sgrd;TWI.sgrd;SINS.sgrd;SMU1.sgrd;SMU3.sgrd;SMU4.sgrd;SM 
>
> U9.sgrd", GRID="SOLUM_rk.sgrd", SHAPES="baranja.shp", FIELD=0, 
> MODEL=1, NUGGET=0, SILL=200,
> RANGE=500, INTERPOL=0))
>
> see also http://geomorphometry.org/view_scripts.asp?id=24
>
> But the last version of SAGA had some bug in this module, so I 
> remember that I did not get correct
> results (I did not test the most recent version of SAGA).
>
> I personally think that we should simply think about implementing 
> local RK (regressions/variograms
> in moving window; maps of regression coefficients, variogram 
> parameters, R-square etc.), and this
> would then solve this problem (+give us much more insight into the data).
>
>
> Tom Hengl
> http://spatial-analyst.net
>
>
> -----Original Message-----
> From: gstat-info-bounces at geo.uu.nl [mailto:gstat-info-bounces at 
> geo.uu.nl] On Behalf Of Edzer Pebesma
> Sent: donderdag 21 augustus 2008 13:07
> To: Vasat, Radim
> Cc: Gstat-info at geo.uu.nl; sig-geo
> Subject: Re: [Gstat-info] LOCAL universal kriging with GLOBAL reg. 
> coeff.estimation
>
> Good question, I'll include r-sig-geo as well.
>
> actually the prediction equations you end up with are kind of funny, and
> I've never seen them written out. Two different covariance matrices
> being inversed.
>
> Package gstat can do the two-step approach: global BLUE, then simple
> kriging of residual, but that ignores the correlation of both terms,
> when added.
>
> I'm quite sure that if you do Gaussian simulation however with gstat,
> the trend is fitted globally (and simulated from the corresponding GLS
> mean & covariance), whereas the kriging is done locally, so that should
> result in realisations that have the variance you want. You just have to
> compute the average and variance, to get estimates of the kriging mean
> and variance you're looking for. Also, it might be computationally
> demanding, both in time and space (memory).
> -- 
> Edzer
>
> Vasat, Radim wrote:
> > Dear all,
> >
> > My question is how to performe LOCAL universal kriging with
> > estimation of regression coefficients from the whole area
> > (GLOBAL estimation) at once.
> >
> > For large data (I have more than 8000 records) is computationaly too
> > demanding to handle with all the data for UK prediction. Hence,
> > I prefer to do it LOCALLY (namx=100), but the reg. coeff.
> > would be estimated only with limited data in this case (not all
> > data are included).
> >
> > To estimate the reg. coeff. first (with BLUE=TRUE) and then perform
> > the simple kriging with beta equals to the reg. coeff. and nmax=100
> > might be the
> > solution. But the reg. coeff. estimation error is not included into the
> > final kriging variance in this case.
> >
> > Anybody knows how to join these two requirements (LOCAL UK and GLOBAL
> > estinmation of reg. coeff.) into one prediction???
> >
> > I will appreciate any coments!
> > Radim
> >
> > _______________________________________________
> > Gstat-info mailing list
> > Gstat-info at geo.uu.nl
> > http://mailman.geo.uu.nl/mailman/listinfo/gstat-info
> > 

-- 
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.springer.com/978-0-387-78170-9