Classification with SuperLearner

What I have come to understand in pulling this script together is that classification of discrete classes needs to be done step-wise. Here I mean that the training data needs to be converted to a binary question: Is this pixel part of ecotype_X?. Each discrete class is to evaluated individually. A composite map of all classes can be generated at the end using raster calculations … the easiest being raster::which.max().

Convert the question to a binary one

outcomes <- dat[,2]
table(outcomes)

## outcomes
## SBSwk1_01 SBSwk1_05 SBSwk1_06 SBSwk1_07 SBSwk1_08 SBSwk1_09 SBSwk1_11 SBSwk1_78 
##       118        31         9        18        37        54        24        65 
## SBSwk1_Wb sub-xeric     xeric 
##        30        13         9

Currently, I have outcomes that will be used to predict 11 discrete classes as present above. These need to be converted to 11 individual classifications. For an example let’s predict one of the most common classes SBSwk1_01.

outcomes_bin <- as.numeric(outcomes == "SBSwk1_01")
table(outcomes_bin)

## outcomes_bin
##   0   1 
## 290 118

Machine Learners

SuperLearner supports numerous learners use listWrappers() to view them. For this example I have chosen:

• SL.caret.rpart
• SL.ranger (a re-write of randomForest)
• SL.svm – Support vector machine
• SL.mean – Effectively no model – used for a baseline (should be the worst performer)

sl <- c("SL.caret.rpart", "SL.ranger", "SL.svm", "SL.mean")

Set hyper-parameters

The results show that ranger (aka randomForest) performs the best – distinctly outperforming the other models. However, how do other versions of `ranger compare when changing the hyperparameters?

below additional learners are created.

• Adjusting the number of trees within the random forest giving us four options to compare: using 250, 500 (default), 1000, and 2000trees.
• Adjusting the number of variables to split at each tree node is compared:

## Create Modifications to ranger's hyperparameters 

myRanger_250 <- create.Learner("SL.ranger", params = list(num.trees = 250),
                               name_prefix = "Ranger250_trees")

myRanger_500 <- create.Learner("SL.ranger", params = list(num.trees = 500),
                               name_prefix = "Ranger500_trees")

myRanger_1000 <- create.Learner("SL.ranger", params = list(num.trees = 1000),
                               name_prefix = "Ranger1000_trees")

myRanger_2000 <- create.Learner("SL.ranger", params = list(num.trees = 2000),
                               name_prefix = "Ranger2000_trees")



## Comparing var splits at nodes
# Let's try 3 multiplies of this default: 0.5, 1, and 2.
(mtry_seq = floor(sqrt(ncol(dat[,-c(1:2)])) * c(0.5, 1, 2)))

## [1]  4  8 16

myRanger_mtrys <- create.Learner("SL.ranger", tune = list(mtry = mtry_seq),
                                 name_prefix = "Ranger500_mtrys")

myRanger_1000_mtry4 <- create.Learner("SL.ranger", params = list(num.trees = 1000,
                                                                 mtry_seq = 4),
                               name_prefix = "Ranger1000_trees_mtry4")

Estimate of risk: CV.SuperLearner

The variations of SL.ranger along with SL.mean (for comparison) are tested.

Note the additional sending of these learners above to the multicore process

## Create vector of SuperLearners to use 
# sl2 <- c("my.ranger.250t", "SL.ranger", "my.ranger.1000t", "my.ranger.2000t", "SL.mean" )
sl2 <- c(myRanger_250$names, myRanger_500$names, myRanger_1000$names, 
         myRanger_2000$names, myRanger_mtrys$names, myRanger_1000_mtry4$names, 
         "SL.svm",  "SL.mean")

sl2

##  [1] "Ranger250_trees_1"        "Ranger500_trees_1"       
##  [3] "Ranger1000_trees_1"       "Ranger2000_trees_1"      
##  [5] "Ranger500_mtrys_1"        "Ranger500_mtrys_2"       
##  [7] "Ranger500_mtrys_3"        "Ranger1000_trees_mtry4_1"
##  [9] "SL.svm"                   "SL.mean"

modelName <- paste0("e:/workspace/tmp/SuperLearner/", "SiteSeries_rangers_01", ".rds")
## if statement -- don't run the model if it has already been done 
if(!file.exists(modelName)) {  ## if statement -- skip generation if it already exists 
  ## Create Function
  # sl_MultiCV <- function(y_train = y_train, x_train = x_train){
      library(snow)
      cluster = parallel::makeCluster(parallel::detectCores()-2) # use half cores
      parallel::clusterEvalQ(cluster, library(SuperLearner))
  
  ## IMPORTANT!! send additional learners to cluster -----------
      parallel::clusterExport(cluster, myRanger_250$names) 
      parallel::clusterExport(cluster, myRanger_500$names) 
      parallel::clusterExport(cluster, myRanger_1000$names)
      parallel::clusterExport(cluster, myRanger_2000$names)
      parallel::clusterExport(cluster, myRanger_mtrys$names)
      parallel::clusterExport(cluster, myRanger_1000_mtry4$names)
  
  ## Set seed -- common seed for all cores
      parallel::clusterSetRNGStream(cluster, 1) ## multicore type seed
      
  # with snow setup run the model
      system.time({
        cv_sl = CV.SuperLearner(Y = outcomes_bin, X = dat[, -c(1:2)], 
                                family = binomial(), ## use binomial for single classes
                                # For a real analysis we would use V = 10.
                                V = 3,
                                parallel = cluster,
                                SL.library = sl2
                                  )
      })
      parallel::stopCluster(cluster)
      # return(cv_sl) ## Turn on when converted to function.
  #} CLOSE FUNCTION
## Save the model object
saveRDS(cv_sl, modelName)
} else {
  cv_sl <- readRDS(modelName)
}
summary(cv_sl)

## 
## Call:  
## CV.SuperLearner(Y = outcomes_bin, X = dat[, -c(1:2)], V = 3, family = binomial(),  
##     SL.library = sl2, parallel = cluster) 
## 
## Risk is based on: Mean Squared Error
## 
## All risk estimates are based on V =  3 
## 
##                     Algorithm     Ave        se     Min     Max
##                 Super Learner 0.18604 0.0092195 0.16623 0.21043
##                   Discrete SL 0.18513 0.0094466 0.16167 0.21171
##         Ranger250_trees_1_All 0.18571 0.0098682 0.16320 0.21003
##         Ranger500_trees_1_All 0.18687 0.0096919 0.16427 0.21352
##        Ranger1000_trees_1_All 0.18570 0.0096900 0.16261 0.21156
##        Ranger2000_trees_1_All 0.18586 0.0097594 0.16315 0.21163
##         Ranger500_mtrys_1_All 0.18610 0.0094894 0.16457 0.21171
##         Ranger500_mtrys_2_All 0.18669 0.0097556 0.16331 0.21265
##         Ranger500_mtrys_3_All 0.18651 0.0098626 0.16381 0.21135
##  Ranger1000_trees_mtry4_1_All 0.18502 0.0096838 0.16167 0.21060
##                    SL.svm_All 0.19438 0.0098584 0.18525 0.21244
##                   SL.mean_All 0.20635 0.0095204 0.19599 0.22151

plot(cv_sl) + theme_minimal()

Tiling of rasters

As per the methods we learned in Victoria the data needs to be converted to a SpatialPointDataFrame. However, the current data set is much too large for this. Also, as emphasized in Victoria, this is where tiling of the raster becomes essential and to do that in parallel processing (while playing nice with the raster package)².

… long story short – I ended up building a wrapper function for SPAdes.tools::splitRaster which generates all the tiles needed using parallel processing – and is friendly with the raster package. I did not include that here – it will be included in the PEMgeneratr package that is in development.

Load a tile

As an example here one tile is loaded.

For full implementation this should be completed within a loop cycling through all the tiles – and done in parallel.

## Figure out which tile to load
# tiles <- "e:/tmpGIS/PEM_cvs/tiles.gpkg"
# tiles <- st_read(tiles)
# mapview::mapview(tiles, zcol = "index", alpha = 0.5) 

lt <- list.files("e:/tmp/cvs/tiles/aspect_025/", pattern = "tile109", full.names = TRUE)
lt <- lt[!grepl(".grd", lt)]

cv_tile <- stack(lt)

Convert to SpatialPixelDataFrame

system.time({
  spdf <- as(cv_tile, "SpatialPixelsDataFrame")
})

##    user  system elapsed 
##    3.82    0.15    4.04

plot(spdf[,"aspect_025"], main = "Aspect - 2.5m resolution")

Apply Prediction

Now all the data can be accessesed and sent through the model generated above.

tileOut <- "./SBSwk1.tif"
if(!file.exists(tileOut)){  ## if prediction is generated don't do it again.
new.data <- spdf@data ## just the data is used (as a data.fame)

## prediction data is generated 
pred <- predict(cv_sl_Prediction, new.data[,cv_sl_Prediction$varNames])

## prediction data is applied to raster 
spdf@data[,"SBSwk1_01"] <- pred ## puts the data back into the spatial dataframe -- as a new column 
r <- raster(spdf[,"SBSwk1_01"])
writeRaster(r, tileOut)
} else {
  r <- raster(tileOut)
}



my_pal <- colorRampPalette(brewer.pal(100, "Greens"))
my_pal <- my_pal((100))
plot(r, main = "Probability of 'SBSwk1_01'", col = my_pal)