library("sperrorest")

data("maipo", package="sperrorest")

predictors<-colnames(maipo)[5:ncol(maipo)]# Construct a formula:fo<-as.formula(paste("croptype ~",paste(predictors, collapse="+")))

Linear Discriminant Analysis (LDA)

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

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

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

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

## [1] 0.0437

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

We can also take a look at the confusion matrix but again, thisresult 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

library("rpart")

fit<-rpart(fo, data=maipo)## optional: view the classiciation tree# par(xpd = TRUE)# plot(fit)# text(fit, use.n = TRUE)

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

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

RandomForest

Bagging, bundling and random forests build upon the CART technique byfitting many trees on bootstrap resamples of the original data set(Breiman 1996)(Breiman2001)(Hothorn and Lausen 2005). They differin that random forest also samples from the predictors, and bundlingadds an ancillary classifier for improved classification. We will usethe widely used random forest in its implementation in theranger 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.53 %

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! Eventhe OOB (out-of-bag) estimate of the error rate is < 1%. Too good tobe true? We’ll see…

Linear Discriminant Analysis (LDA)

First, we need to create a wrapper predict method for LDA forsperrorest(). This is necessary in order to accommodate themajority filter, and also because class predictions fromlda’s predict method are hidden in the$classcomponent 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(levinlevels(fac)){x[fac==lev]<-majority(x[fac==lev])}x}pred<-predict(object, newdata=newdata)$classif(!is.null(fac))pred<-majority_filter(pred,newdata[,fac])return(pred)}

To ensure that custom predict-functions will work withsperrorest(), we need to wrap all custom functions in onesingle function. Otherwise,sperrorest() might fail duringexecution.

Finally, we can runsperrorest() with a non-spatialsampling setting (partition_cv()). In this example we use a‘3 repetitions - 5 folds’ setup. In reality, we recommend to use 100repetitions 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 worththe trouble parallelizing the execution of three repetitions with ashort runtime.

So what have we got:

summary(res_lda_nsp$error_rep)

##                    mean       sd   median      IQR## train_error    3.35e-02 0.000357 3.33e-02 0.000324## train_accuracy 9.67e-01 0.000357 9.67e-01 0.000324## 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     3.93e-02 0.003323 3.97e-02 0.003306## test_accuracy  9.61e-01 0.003323 9.60e-01 0.003306## 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 usepartition_factor_cv() as the sampling function. Since weare using 5 folds, we get a coarse 80/20 split of our data. 80% will beused for training, 20% for testing our trained model.

To take a look where our training and tests sets will be partitionedon each fold, we can plot them. The red colored points represent thetest set in each fold, the black colored points the training set. Notethat because we plotted over 7000 points, overplotting occurs and sincethe red crosses are plotted after the black ones, it seems visually thatway 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)insperrorest().

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 6.45 secs

summary(res_lda_sp$error_rep)

##                    mean      sd   median     IQR## train_error    3.10e-02 0.00185 3.16e-02 0.00177## train_accuracy 9.69e-01 0.00185 9.68e-01 0.00177## train_events   4.69e+03 0.00000 4.69e+03 0.00000## train_count    3.09e+04 0.00000 3.09e+04 0.00000## test_error     6.23e-02 0.00724 6.47e-02 0.00694## test_accuracy  9.38e-01 0.00724 9.35e-01 0.00694## test_events    1.17e+03 0.00000 1.17e+03 0.00000## test_count     7.71e+03 0.00000 7.71e+03 0.00000

RandomForest

In the case of Random Forest, the customizedpred_funlooks as follows; it is only required because of the majority filter,without it, we could just omit thepred_fun andpred_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(levinlevels(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.01e-02 0.00943 9.09e-02 0.0094## test_accuracy  9.10e-01 0.00943 9.09e-01 0.0094## 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.91 0.00943  0.909 0.0094

What a surprise! {ranger}‘s classification is not that good afterall, if we acknowledge that in ’real life’ we wouldn’t be makingpredictions in situations where the class membership of other grid cellsin the same field is known in the training stage. So spatial dependencedoes matter.