Regression kriging

[juliaDa]

Hi,

I am trying to perform regression kriging in ArcMap. Is there way to do this in ArcMap? Will appreciate any input!!!

Thanks in advance

[SteveLynch]

It is currently under development.

[SteveLynch]

What about Cokriging?

[EricKrause]

It depends on what kind of regression kriging you want to do.  If the explanatory variables are polynomials of the (x,y) coordinates, we call that Universal Kriging, and it is available in the Geostatistical Wizard.

If your covariates are not (x,y) polynomials, you can use the Ordinary Least Squares tool to create a regression equation, and then you can perform kriging on the residuals.  To get predictions, add the interpolated residuals to the prediction from the OLS model.

Ideally, the kriging parameters and the coefficients from OLS will be calculated simultaneously (rather than doing OLS then doing kriging), but there is currently no way to do this in ArcGIS.  Though that will likely change in the not-so-distant future.

[RoniLuo]

Hello Eric,

I am trying to carry out a regression kriging analysis for my groundwater project, but I am not sure if my approaches are correct and struggling to calculate RMSE and ME. Can you please give me some comments?

I have groundwater level(n=124) as dependent variable, elevation, temperature, precipitation and NDVI as independent variables.

1.run OLS in Geoda, the Moran`s I test showed that residuals are spatially autocorrelated

2.populate predicted values and residuals back to the groundwater layer

3.use ordinary kriging to interpolate the predicted values and residuals separately

4.convert the output interpolated surfaces to raster

5.use "Plus" tool in Raster to add the interpolated predicted values and residual raster together

thanks so much in advance,

Roni

[EricKrause]

There is an issue in your third step.  I'm assuming that your dependent variable and independent (explanatory) variables are coming from fields in a point feature class.

In this setup, there is no "correct" way to produce a prediction surface.  The explanatory variables must be measured in every location where you want to make a prediction.  I understand that this is often impractical, and what is often done in this case is to interpolate the explanatory variables.  It is not a good idea to interpolate the predicted values because you know that they depend on the explanatory variables.  However, when you interpolate the explanatory variables, all subsequent analysis will assume that the interpolated values are measured values when, in fact, they are predictions and have associated error.

In your case, you need to calculate the OLS equation from the input points, then interpolate the independent variables and the residuals, and plug the results into the OLS equation to get the predicted values.

Regarding cross validation, only the residuals (error terms) are needed.  When you are interpolating the residuals, note the cross validation statistics.  These are what you should report as the cross validation statistics for the regression kriging model.

However, as I alluded to in my first post, we are releasing "EBK Regression Prediction" in ArcGIS Pro 1.2, which does a regression kriging interpolation using Empirical Bayesian Kriging.  It will very shortly be available as a Beta through Esri's customer care website.  I am not sure exactly where to go to download it, but I can find out if you are interested.

[RoniLuo]

Hi Eric,

Thanks for your reply. Yes, I have point data. Based on your suggestions, I did the followings:

1. Run OLS and the results showed elevation and precipitation are significant, and residuals are spatially autocorrelated
2. Interpolate significant independent variables (elevation and precipitation) and residuals separately
3. Convert the interpolated surfaces to raster
4. Use “Plus” tool to add the two rasters together

I am afraid I am still not quite clear what you meant by “plug the results into the OLS equation to get the predicted values”. Is it the last step that I will get my prediction map? Can you please elaborate more? My guess was to add the interpolated independent variable and residual surfaces together, which I did, but predicted values are way too high and have negative values.

It would be very nice to have access to EBK Regression Prediction. Please help me find it if it becomes available.By the way, do I have to pay this new product?

Much appreciate your help!

Roni

[EricKrause]

This is what I mean by plugging the the rasters into the OLS equation:

When you run OLS, you will get coefficients for the explanatory variables along with an intercept value.  For example, let's say that the precipitation variable gets a coefficient of 2, elevation gets a coefficient of 3, and the intercept term is -4.

You would interpolate elevation, precipitation, and residuals to rasters.

You would then apply the following formula to get your prediction surface:

Prediction = -4 + 2 * (precipitation) + 3 * (elevation) + (residual)

The Raster Calculator tool is the easiest way to apply this formula to all raster cells.

For example, in a location where the precipitation is 500, elevation is 1000, and the residual is 50, the prediction would be:

Prediction = -4 + 2 * 500 + 3 * 1000 + 50 = 4046

So, your prediction at that location would be 4046.

I will look into how you can download ArcGIS Pro 1.2 beta.  It will be available as part of the ArcGIS for Desktop 10.4 release, so if you are on a maintenance contract that has this product available, you will be able to download the beta.

[RoniLuo]

Hello Eric,

Your example is very helpful. I was able to follow and calculate. The maximum value in my RK prediction map is 40, and the max in my cokriging prediction map is 109. Is the difference mainly due to the regression equation?

Also, I realize that my dependent variable, groundwater level, and independent variable, elevation are positively skewed. Should I transform them before I build OLS? However, field calculator only has "log()" function but not "ln()". For example, ln(GW) = intercept + b*ln(elevation) + residuals

If so, I will need to back transform the results after I fit the residuals back into the trend. Is that correct?

thanks again!

[DanPatterson_Retired]

python parser in the field calculator

math.log(some number, base)  if base is omitted, it is assumed to be base e

math.log10(some number)    base 10 log, usually more accurate than math.log(some number,10) according to help