[GIS] Random Forest out of the bag and confusion matrix

classificationconfusion matrixrrandom forestremote sensing

I have implemented Random Forest classifier to classify remote sensing data in R. The original code comes from here:

How to perform Random Forest land classification?

Everything works fine, but what I need is to obtain confusion matrix and out of the bag error for fast classification accuracy assessment. I have created junk dataset for validating functionality of the whole script itself (only 2 classes). I added parallel computing to speed up computations.

I am wondering if my implementation is correct and If I will be able to obtain confusion matrix and out of the bag error when I will add much more features describing classes later. Here is part of my code:

sdata <- readOGR(dsn = vyber_shp, layer = "trenink")

# extraction from raster to be classified

rdata <- data.frame(extract(r, sdata))

# create random forest model in parallel

cl <- makeCluster(detectCores())

registerDoParallel(cl)

rf.mdl <- foreach(ntree=rep(125, 4), .combine = combine, .packages = 'randomForest') 
%dopar% {randomForest(x=rdata, y=sdata$Class, ntree=500, 
proximity=TRUE, importance=TRUE, confusion=TRUE, 
do.trace=TRUE, err.rate=TRUE) 
    }

stopCluster(cl)

# classify raster in parallel

beginCluster()

predikce <- clusterR(r, raster::predict, args=list(model=rf.mdl))

endCluster()

# write and save result of classification as raster image

klasifikace <- writeRaster(predikce, filename='Klasifikace_RandomForest', 
format='HFA', options='INTERLEAVE=BSQ', datatype='INT2S', overwrite=TRUE)

varImpPlot(rf.mdl)

rf.mdl$confusion

Best Answer

The combine function (which is used in foreach) does not store the relevant components into the final randomForest object. See ?randomForest::combine:

The confusion, err.rate, mse and rsq components (as well as the corresponding components in the test component, if exist) of the combined object will be NULL.

But the predict method returns OOB predictions if the newdata argument is omitted:

library("doSNOW")
library("foreach")

cl <- makeCluster(2, type="SOCK")
registerDoSNOW(cl)

set.seed(1)
x <- subset(iris, select=-Species)
y <- iris$Species

rf <- foreach(ntree=rep(250, 4), .combine=combine, .packages="randomForest") %dopar%
      randomForest(x=x, y=y, ntree=ntree, norm.votes=FALSE)

rf
# Call:
#  randomForest(x = x, y = y, ntree = ntree, norm.votes = FALSE) 
#                Type of random forest: classification
#                      Number of trees: 1000
# No. of variables tried at each split: 2

# `predict.randomForest` returns OOB predictions if `newdata` is not given
rf_pred <- predict(rf)

caret::confusionMatrix(rf_pred, y)
# Confusion Matrix and Statistics
#
#             Reference
# Prediction   setosa versicolor virginica
#   setosa         50          0         0
#   versicolor      0         47         3
#   virginica       0          3        47
#
# Overall Statistics
#
#                Accuracy : 0.96
#                  95% CI : (0.915, 0.985)
#     No Information Rate : 0.333
#     P-Value [Acc > NIR] : <2e-16
#
#                   Kappa : 0.94
#  Mcnemar's Test P-Value : NA
#
# Statistics by Class:
#
#                      Class: setosa Class: versicolor Class: virginica
# Sensitivity                  1.000             0.940            0.940
# Specificity                  1.000             0.970            0.970
# Pos Pred Value               1.000             0.940            0.940
# Neg Pred Value               1.000             0.970            0.970
# Prevalence                   0.333             0.333            0.333
# Detection Rate               0.333             0.313            0.313
# Detection Prevalence         0.333             0.333            0.333
# Balanced Accuracy            1.000             0.955            0.955

Related Solutions

[GIS] Computing unsupervised random forest classification in R

Random Forests in unlabeled (unsupervised) mode does not return explicit classes but, rather something analogous to scaled multivariate distances which is based on node proximities. Without the proximity matrix, you do not have a usable unlabeled model. And yes, for large problems, even using a sparse matrix, the very nature of the approach causes the proximity matrix to get huge. This may very well be the reason that you have not seen published approaches using Random Forests in unsupervised remote sensing.

Based on the proximities an approach that I have seen, to derive/test clusters, is to use a modified K-means on the proximity matrix. Alternately, you may be able to trick the imputation function, using the random forests option, in the yaImpute package to perform a matrix imputation which would return something analogs to a k nearest neighbor (kNN) that could then be assigned to clusters based on a similarity matrix.

It is nothing near as straightforward as what you are thinking and I would encourage you to research this approach before jumping in with both feet.

**** Edit 12/14/2018 A few versions ago I added an unsupervised random forests function to the rfUtilities package. I would not recommend it on large data such as rasters but it is a useful clustering method. Here is a simple example.

library(rfUtilities)
library(sp)

data(meuse)
  meuse <- na.omit(meuse)

n = 6  
clust.meuse <- rf.unsupervised(meuse, n=n, proximity = TRUE, 
                               silhouettes = TRUE)
( meuse$k <- clust.meuse$k )

mds <- stats:::cmdscale(clust.meuse$distances, eig=TRUE, k=n)
  colnames(mds$points) <- paste("Dim", 1:n)
  mds.col <- ifelse(clust.meuse$k == 1, rainbow(6)[1],
               ifelse(clust.meuse$k == 2, rainbow(6)[2],
                 ifelse(clust.meuse$k == 3, rainbow(6)[3],
                   ifelse(clust.meuse$k == 4, rainbow(6)[4],
                    ifelse(clust.meuse$k == 5, rainbow(6)[5],
                     ifelse(clust.meuse$k == 6, rainbow(6)[6], NA))))))
plot(mds$points[,1:2],col=mds.col, pch=20)                 
pairs(mds$points, col=mds.col, pch=20)

coordinates(meuse) <- ~x+y
plot(meuse, col=mds.col, pch=19)
  box()

[GIS] Random Forest land cover classification in ArcMap

I am the main developer of MGET.

The first step in your problem is to obtain values of the covariates that you will use to fit the model to your 90 GPS points. It sounds like you want to use the 8 bands as your covariates. You need to add 8 fields to your shapefile (one for each band) and populate them using a tool such as Extract Multi Values to Points from recent versions of ArcGIS or Interpolate Raster Values at Points from MGET (equivalent to what Arc provides but developed before the Arc tool existed).

After that, you need to fit a classification model to the GPS points, using the field containing the known cover type as the response variable and the 8 band fields as the covariates (a.k.a. predictor variables). After that you can obtain some performance statistics for your model and then predict it on rasters representing the covariates.

You can see a basic overview of MGET's modeling workflow for this here. The example is somewhat dated--not all of the tool parameters will look exactly like what you see there--but the basic workflow is the same: fit the model to a table of data, predict it against the table to get some performance statistics, and predict it on a stack of rasters. In MGET, the procedure is the same regardless of which modeling framework you use--MGET currently provides GLM, GAM, trees (a.k.a. CARTs), and random forest--so you can try different kinds of models with very similar workflows.

I'm sorry I don't have detailed instructions about this workflow written up. So far, we have not had funding to develop a complete manual. All MGET tools have documentation within ArcGIS, so be sure you click the Show Help >> button on the tool dialogs if you have not done so already.

Regarding Jeffrey Evans' speculation that MGET does not utilize the R raster package. That is correct. The code in MGET that performs raster predictions was developed before the raster package was released to CRAN (R's distribution system for R packages), thus it does not rely on that package. But it is not correct that MGET will crash due to memory limitations. MGET's raster prediction code was written specifically to handle the situation you're facing, by performing predictions in blocks similar to how the raster package does it. Prior to the raster package being developed, MGET was one of the only readily-available tools that could handle prediction of large rasters. MGET users have done this, for example, using large bathymetry rasters with 5m resolution.

All of that said, if you believe you will be performing a lot of modeling as your career progresses, I encourage you to learn how to do it in R directly, and about modeling and statistics more generally, independent of software. In a sense, MGET's modeling tools are a "gateway drug" to R. MGET's tools are just as robust as R--they utilize R to perform the actual model fitting and prediction--but they expose only a limited subset of what is possible in R itself. As you continue to do more modeling projects, eventually you may face a situation in which MGET is not enough and you need the full flexibility of R.

Best Answer

Related Solutions

[GIS] Computing unsupervised random forest classification in R

[GIS] Random Forest land cover classification in ArcMap

Related Question