Spatial Modeling Using Statistical Learning Techniques

Linear Discriminant Analysis (LDA)

LDA is simple and fast, and often performs surprisingly well if the problem at hand is ‘linear enough’. As a start, let’s fit a model with all predictors and using all available data:

library(MASS)
fit <- lda(fo, data = maipo)

Predict the croptype with the fitted model and calculate the misclassification error rate (MER) on the training sample:

pred <- predict(fit, newdata = maipo)$class
mean(pred != maipo$croptype)

## [1] 0.0437

But remember that this result is over-optimistic because we are re-using the training sample for model evaluation. We will soon show you how to do better with cross-validation.

We can also take a look at the confusion matrix but again, this result is overly optimistic:

table(pred = pred, obs = maipo$croptype)

##        obs
## pred    crop1 crop2 crop3 crop4
##   crop1  1294     8     4    37
##   crop2    50  1054     4    44
##   crop3     0     0  1935     6
##   crop4    45   110    29  3093

RandomForest

Bagging, bundling and random forests build upon the CART technique by fitting many trees on bootstrap resamples of the original data set (Breiman 1996) (Breiman 2001) (Hothorn and Lausen 2005). They differ in that random forest also samples from the predictors, and bundling adds an ancillary classifier for improved classification. We will use the widely used random forest in its implementation in the ranger package.

library("ranger")

fit <- ranger(fo, data = maipo)
fit

## Ranger result
## 
## Call:
##  ranger(fo, data = maipo) 
## 
## Type:                             Classification 
## Number of trees:                  500 
## Sample size:                      7713 
## Number of independent variables:  64 
## Mtry:                             8 
## Target node size:                 1 
## Variable importance mode:         none 
## Splitrule:                        gini 
## OOB prediction error:             0.62 %

Let’s take a look at the MER achieved on the training sample:

pred <- predict(fit, data = maipo, type = "response")
mean(pred$predictions != maipo$croptype)

## [1] 0

table(pred = pred$predictions, obs = maipo$croptype)

##        obs
## pred    crop1 crop2 crop3 crop4
##   crop1  1389     0     0     0
##   crop2     0  1172     0     0
##   crop3     0     0  1972     0
##   crop4     0     0     0  3180

Isn’t this amazing? All grid cells were correctly classified! Even the OOB (out-of-bag) estimate of the error rate is < 1%. Too good to be true? We’ll see…

Linear Discriminant Analysis (LDA)

First, we need to create a wrapper predict method for LDA for sperrorest(). This is necessary in order to accommodate the majority filter, and also because class predictions from lda’s predict method are hidden in the $class component of the returned object.

lda_predfun <- function(object, newdata, fac = NULL) {

  library(nnet)
  majority <- function(x) {
    levels(x)[which.is.max(table(x))]
  }

  majority_filter <- function(x, fac) {
    for (lev in levels(fac)) {
      x[fac == lev] <- majority(x[fac == lev])
    }
    x
  }

  pred <- predict(object, newdata = newdata)$class
  if (!is.null(fac)) pred <- majority_filter(pred, newdata[, fac])
  return(pred)
}

To ensure that custom predict-functions will work with sperrorest(), we need to wrap all custom functions in one single function. Otherwise, sperrorest() might fail during execution.

Finally, we can run sperrorest() with a non-spatial sampling setting (partition_cv()). In this example we use a ‘3 repetitions - 5 folds’ setup. In reality, we recommend to use 100 repetitions to reduce the influence of random partitioning.

res_lda_nsp <- sperrorest(fo,
  data = maipo, coords = c("utmx", "utmy"),
  model_fun = lda,
  pred_fun = lda_predfun,
  pred_args = list(fac = "field"),
  smp_fun = partition_cv,
  smp_args = list(repetition = 1:3, nfold = 5),
  mode_rep = "sequential",
  progress = FALSE
)

Note that we’re running this in sequential mode since it’s not worth the trouble parallelizing the execution of three repetitions with a short runtime.

So what have we got:

summary(res_lda_nsp$error_rep)

##                    mean       sd   median      IQR
## train_error    3.38e-02 0.000567 3.38e-02 0.000567
## train_accuracy 9.66e-01 0.000567 9.66e-01 0.000567
## train_events   4.69e+03 0.000000 4.69e+03 0.000000
## train_count    3.09e+04 0.000000 3.09e+04 0.000000
## test_error     4.06e-02 0.001162 4.12e-02 0.001037
## test_accuracy  9.59e-01 0.001162 9.59e-01 0.001037
## test_events    1.17e+03 0.000000 1.17e+03 0.000000
## test_count     7.71e+03 0.000000 7.71e+03 0.000000

To run a spatial cross-validation at the field level, we can use partition_factor_cv() as the sampling function. Since we are using 5 folds, we get a coarse 80/20 split of our data. 80% will be used for training, 20% for testing our trained model.

To take a look where our training and tests sets will be partitioned on each fold, we can plot them. The red colored points represent the test set in each fold, the black colored points the training set. Note that because we plotted over 7000 points, overplotting occurs and since the red crosses are plotted after the black ones, it seems visually that way more than ~20% of red points exist than it is really the case.

resamp <- partition_factor_cv(maipo, nfold = 5, repetition = 1:1, fac = "field")
plot(resamp, maipo, coords = c("utmx", "utmy"))

Subsequently, we have to specify the location of the fields (fac = "field") in the prediction arguments (pred_args) and sampling arguments (smp_args) in sperrorest().

res_lda_sp <- sperrorest(fo,
  data = maipo, coords = c("utmx", "utmy"),
  model_fun = lda,
  pred_fun = lda_predfun,
  pred_args = list(fac = "field"),
  smp_fun = partition_factor_cv,
  smp_args = list(fac = "field", repetition = 1:3, nfold = 5),
  mode_rep = "sequential",
  benchmark = TRUE, progress = FALSE
)
res_lda_sp$benchmark$runtime_performance

## Time difference of 10.2 secs

summary(res_lda_sp$error_rep)

##                    mean       sd   median     IQR
## train_error    2.98e-02 0.000471 2.98e-02 0.00047
## train_accuracy 9.70e-01 0.000471 9.70e-01 0.00047
## train_events   4.69e+03 0.000000 4.69e+03 0.00000
## train_count    3.09e+04 0.000000 3.09e+04 0.00000
## test_error     6.17e-02 0.009400 6.16e-02 0.00940
## test_accuracy  9.38e-01 0.009400 9.38e-01 0.00940
## test_events    1.17e+03 0.000000 1.17e+03 0.00000
## test_count     7.71e+03 0.000000 7.71e+03 0.00000

RandomForest

In the case of Random Forest, the customized pred_fun looks as follows; it is only required because of the majority filter, without it, we could just omit the pred_fun and pred_args arguments below.

rf_predfun <- function(object, newdata, fac = NULL) {

  library(nnet)
  majority <- function(x) {
    levels(x)[which.is.max(table(x))]
  }

  majority_filter <- function(x, fac) {
    for (lev in levels(fac)) {
      x[fac == lev] <- majority(x[fac == lev])
    }
    x
  }

  pred <- predict(object, data = newdata)
  if (!is.null(fac)) pred <- majority_filter(pred$predictions, newdata[, fac])
  return(pred)
}

res_rf_sp <- sperrorest(fo,
  data = maipo, coords = c("utmx", "utmy"),
  model_fun = ranger,
  pred_fun = rf_predfun,
  pred_args = list(fac = "field"),
  smp_fun = partition_factor_cv,
  smp_args = list(
    fac = "field",
    repetition = 1:3, nfold = 5
  ),
  mode_rep = "sequential",
  benchmark = TRUE, progress = 2
)

summary(res_rf_sp$error_rep)

##                    mean      sd   median    IQR
## train_error    0.00e+00 0.00000 0.00e+00 0.0000
## train_accuracy 1.00e+00 0.00000 1.00e+00 0.0000
## train_events   4.69e+03 0.00000 4.69e+03 0.0000
## train_count    3.09e+04 0.00000 3.09e+04 0.0000
## test_error     9.56e-02 0.00494 9.70e-02 0.0048
## test_accuracy  9.04e-01 0.00494 9.03e-01 0.0048
## test_events    1.17e+03 0.00000 1.17e+03 0.0000
## test_count     7.71e+03 0.00000 7.71e+03 0.0000

summary(res_rf_sp$error_rep)["test_accuracy",]

##                mean      sd median    IQR
## test_accuracy 0.904 0.00494  0.903 0.0048

What a surprise! {ranger}‘s classification is not that good after all, if we acknowledge that in ’real life’ we wouldn’t be making predictions in situations where the class membership of other grid cells in the same field is known in the training stage. So spatial dependence does matter.