Kriging using multiple spatial input variables and one spatial output (response) variable

Question

I have a data table with the following columns in a CSV file.

X    Y     Zinc      L    F     C
---  ---   -----     --   ---  ----

Of course, I can import the CSV file as a data frame

a <- read.csv('givendata.csv', header = T)

And convert this data into spatial data by

coordinates(a) = ~X + Y 

Now, I wish to run kriging by using L, F and C as spatial input variables and Zinc as the output variable. How do I achieve this?

Edit: I am new to Kriging.

I am adding some more explanation to the question, based on a comment that I do not have an idea about linear modelling.

If I do the following

A1 <- read.csv('givendata.csv', header = T)
mymodel <- lm(Zinc ~ L + F + C, data = A1) 

it gives the output, with all the regression coefficients and the constant. But, these observations have spatial inputs (see the meuse dataset), which is why it is referred to as a spatial variable (this is GIS.SE!).

Now suppose I were to take the following steps:

x.range <- (range(A1@coords[,1]))
y.range <- (range(A1@coords[,2]))

A1.grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=30),  
                    y=seq(from=y.range[1], to=y.range[2], by=30))

coordinates(A1.grd) <- ~x+y
gridded(A1.grd) <- TRUE

library(automap) 
kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd) 

I get the following error:

Error in eval(expr, envir, enclos) : object 'L' not found

which is funny because the lm was working and the Kriging does not work!

Answer 1

As I stated in the comments, you need values of your explanatory variables at the prediction locations in order to make predictions. Here is a fully reproducible example:

library(automap)

Make 100 random data points since we don't have your data file:

A1 = data.frame(x=runif(100),y=runif(100), Zinc=runif(100), L=runif(100), F=runif(100), C=runif(100))
coordinates(A1)=~x+y

Predict over this grid (note I've used len=20 instead of by=30 from your code because my numbers are different):

x.range <- (range(A1@coords[,1]))
y.range <- (range(A1@coords[,2]))

A1.grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], len=20),  
                    y=seq(from=y.range[1], to=y.range[2], len=20))
coordinates(A1.grd) <- ~x+y
gridded(A1.grd) <- TRUE

Now let's reproduce what you did:

kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd) 

That reproduces your error. Now I can see what your problem is, something impossible to tell from your initial post. Thank you for editing.

The problem is the model wants to predict Zinc at the grid locations, but you need values of L, F, and C at those locations, in order to do predictions. Let's just make some up to illustrate. I have a 20x20 grid so that's 400 points:

A1.grd$L=runif(400)
A1.grd$F=runif(400)
A1.grd$C=runif(400)

Moment of truth:

kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd)

And that works.

Now we have no idea what your L, F, and C represent. Normally they will be one of two kinds of things:

1. Spatially continuous quantities that you can get the value of at any (x,y) location, such as altitude, or distance to nearest road. Compute them over your (x,y) grid and do kriging as above.

2. Non-continuous quantities only valid at the sample level. For example if these are the amounts of Zinc in human blood samples you might have the age of the human sampled. How do you then do kriging predictions at a point? You don't have sample-specific measurements there because you haven't any samples there.

In this case you can literally make up some values. To continue my example you might choose the average age of people in your study, do the kriging with that, and then show a kriged map of "Zinc amount for 27 year olds". With more variables you may end up with something like "Zinc amount for unvaccinated 27 year old men" if your variables are vaccination status, age, and sex.

You could also do separate kriged maps for extreme values of your covariates, and show them on two plots, something like "Range of Zinc values for ages 12 (left figure) to 52 (right figure)".

Or you might find that your covariates aren't significant at all and just drop them.

Answer 2

There are two widely used packages to achieve Kriging interpolation and variogram model fitting in R, these are gstat and geoR. Depending on the package you use, you can build a native gstat object or a geoR object. Also, you can manage to do all your work using the known SpatialPointsDataFrame object.

Some general steps prior to achieve Kriging interpolation would be:

1. Representing Variation: describe the variation of your data by meanings of Frequency distribution (Histograms and Box-plots), Centre (mean, median and mode) and Dispersion (range, interquartile range, standard deviation, variance, coefficient of variation, skewness and kurtosis) analysis.

For Frequency distribution analysis you can use histogram and boxplot functions or if you want more attractive plots you can use ggplot2 library (geom_histogram and geom_boxplot functions). For Centre analysis you can use mean, median and this Mode function: Mode <- function(x) {ux = unique(x); ux[which.max(tabulate(match(x, ux)))]}. For Dispersion analysis you can use range, IQR, sd and var functions, calculate coefficient of variation as (sd/mean)*100 and skewness and kurtosis functions from moments library. Also, try to find outliers at this points and evaluate if it's better to remove them or not.

2. Transformation: attempt to transform the measured values to a new scale on which the distribution is more nearly normal to overcome difficulties arising from departures from normality.

To test normality of your data you can perform a Shapiro-Wilk test with shapiro.test function. If your data is not normal you can try simple transformations (log and sqrt) or see Box-Cox transformations using forecast library (BoxCox.lambda and BoxCox functions). Also, you can perform a "visual" analysis of normality using Q-Q plots, you can use ggQQ <- ggData + geom_qq(aes(sample = var)) + geom_abline(intercept = mean(var), slope = sd(var), col = "red") from ggplot2 library.

3. Trend analysis: discover trends between the response variable and co-variables using scatter plots with smoothed lines and/or correlations and/or correlograms.

To make a scatter plot you just can simply plot(x = var1, y = var2) and try a linear model abline(lm(var2 ~ var1, data = data), col = "red") or using ggplot2 library ggTrend_var1var2 <- ggData + geom_point(aes(y = var1, x = var2)) + geom_smooth(aes(y = var1, x = var2)). For correlation you can use Pearson correlation test: cor.test(x = var2, y = var1, method = "pearson"). To perform correlograms you can use the function corrgram from library corrgram. After analysing possible trends by means of correlations and plots, you can evaluate some trend models e.g.: multiple linear regression model with vars 2, 3 and 4 (model1 <- lm(formula = var1 ~ var2 + var3 + var4, data = data)) or a simple linear regression model with only var 2 (model2 <- lm(formula = var1 ~ var2, data = data)) and after that choose the best model with Akaike criterion AIC(model1, model2).

4. Find spatial outliers: try to find which points values are atypically different from their neighborhood point values, you can use Moran Index to quantify it and make a "posting" plot to explore your data. At this point you need to evaluate to remove or not (depending on your data) those spatial outliers points.

You can use moran.plot from spdep library to make a Moran Scatterplot and get influenced values. To perform a "posting" plot you can use the function spplot(obj = SpatialPointsDataFrame, zcol = "var1", colorkey = TRUE) from sp library or the great interactive map mapview(SpatialPointsDataFrame, zcol = "var1") from mapview library.

5. Empiric variogram: characterizing the spatial processes of your data.

Here you can use the function variog from geoR library or variogram function from gstat library. In geoR and gstat the trend can be included as a formula: trend = trend.spatial(var1 ~ var2 + var3) option in geoR::variog function and object = var1 ~ var2 + var3 option in 'gstat::variogram'. At this point you should consider find possible anisotropy process in your data. To check that, you have to perform directional variograms (typically 0, 45, 90 and 135 degrees). See direction, tolerance and unit.angle options in geoR::variog function or alpha option in gstat::variogram function.

6. Adjust theoretical variogram: adjust a theoretical variogram model to the empirical variogram. The choose of the model will depend on the best fit of a series of theoretic variogram models tried or one choose based on your previously knowledge of how the spatial structure of your data is best described.

See gstat::vgm(), geoR::likfit and geoR::variofit functions.

7. Kriging interpolation: there are different types of Kriging interpolations, you should consider which one is best for you (or your data!). At this point, you need your variogram model (output from vgm, likfit or variofit) and know which points positions in space you are interested to predict (it can be a grid, or some irregular point data).

One example of Ordinary Kriging using gstat could be gstat::krige(formula = var1 ~ var2 + var3, locations = data, newdata = pointsPositionsToPredict, model = variogramModel)

• Note: geoR and gstat packages are not compatible, these libraries work with different objects. For example, you can't use the variogram object model from geoR in a gstat::krige function. Try use one package first, but you can use both packages building manually the objects using the parameters derived from one package function to the other.

• For further information you can check this book Geostatistics for Environmental Scientists, 2nd Edition (Webster & Oliver).

• A code work example IntroductionToGeostatistics (only spanish version)