Partitioning the data

sparklyr gives us an easy way to partiton a Spark DataFrame. Here we’re creating a copy of df in Spark with compute(), before partitioning this copy into test and training data. We can then access the partitions with, e.g., df_part$test. We can have as many partitions as we want, so in practice we’d probably want to do k-fold cross validation (see below).

df_part <- df %>%
  compute("df_part") %>%
  sdf_partition(test = 0.2, train = 0.8, seed = 2017)

Fit the model

From sparklyr we can access a range of models available to use with Spark. They are all prefixed with ml_. To run a random forest we can use ml_random_forest(). One thing that is a bit confusing if you’re working with the MLlib documentation is that some of the parameter names are quite different for the sparklyr functions compared to what they’re called by MLlib. For example, the proportion of the predictors to use for making splits of a random forest is called featureSubsetStrategy in the MLlib documentation but col.sample.rate by sparklyr. Another example is that MLlib calls the number of trees for gradient boosted trees numInterations whereas sparklyr calls it num.trees.

fit_rf <- ml_random_forest(
  df_part$train, # the training partion
  response = "y",
  features = colnames(df_part$train)[1:16], # the names minus the outcome
  impurity = "entropy",
  max.bins = 32L, # default = 32L
  max.depth = 5L, # default = 5L
  num.trees = 100L,  
  learn.rate = 0.1, # default = 0.1
  col.sample.rate = 0.25, # aka featureSubsetStrategy
  type = "classification",
  seed = 2017
)

Evaluate the model

We can get out the feature importance for the model by providing our Spark connection and the model object.

ml_tree_feature_importance(sc = sc, model = fit_rf)

##             importance   feature
## 1    0.513011341044725  duration
## 2   0.0802082038509115  poutcome
## 3    0.059010844561907       age
## 4   0.0538740319352658     pdays
## 5   0.0500510804148416     month
## 6   0.0494220628070373   contact
## 7   0.0437808856900149  previous
## 8   0.0356504629964395       day
## 9   0.0294190550796405   balance
## 10  0.0207771897883123       job
## 11  0.0196039047639807   housing
## 12    0.01626807994962  campaign
## 13  0.0109601402422286   marital
## 14 0.00964582937505621      loan
## 15 0.00699162595376989 education
## 16 0.00132526154624899   default

We can make predictions from a model for some data with sdf_predict(). MLlib also has built in evaluation for models, which can be accessed through sparklyr. Below we are making predictions for the training data, then getting the F1 score. Note that the data we provide to sdf_predict() is in Spark not R, we’re just providing the R connection to the Spark data. Assuming we’re using Spark 2.X, we could alternatively set the metric to weightedPrecision, weightedRecall or accuracy. More metrics are available for Spark 1.6.

fit_rf %>%
  # predict from the model for the training data
  sdf_predict(df_part$train) %>%
  ml_classification_eval(label = "y",
                         predicted_lbl = "prediction",
                         metric = "f1") # default, F1 score 

## [1] 0.858982

Next we can do predictions for the test data

test_rf <- sdf_predict(fit_rf, df_part$test)

And calculate the F1 score again.

(test_rf_f1 <- test_rf %>%
    ml_classification_eval(label = "y",
                           predicted_lbl = "prediction",
                           metric = "f1")) # default, F1 score 

## [1] 0.8640745

The random forest model returns probababilities as well as predicted categories, unlike other models (e.g. gradient boosted trees). As our outcome is also binary we can use ml_binary_classification_eval() to get back the area under the ROC curve for our model.

(rf_test_auc <- test_rf %>%
    ml_binary_classification_eval(label = "y",
                                  score = "probability")) 

## [1] 0.8911765

This information is a little limited though. Data size permitting, we could collect() our DataFrame into R to get a tibble of the data used to make the predictions, with prediction and probability columns added.

rf_test_df <- collect(test_rf)

We can then use caret::confusionMatrix() to get a more exhaustive set of performance metrics. From the output it seems the RF is no better than just predicting the majority category for everyone (the No Information Rate).

confusionMatrix(data = rf_test_df$prediction, reference = rf_test_df$y)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 822  89
##          1   8  13
##                                           
##                Accuracy : 0.8959          
##                  95% CI : (0.8745, 0.9148)
##     No Information Rate : 0.8906          
##     P-Value [Acc > NIR] : 0.3222          
##                                           
##                   Kappa : 0.1808          
##  Mcnemar's Test P-Value : 4.557e-16       
##                                           
##             Sensitivity : 0.9904          
##             Specificity : 0.1275          
##          Pos Pred Value : 0.9023          
##          Neg Pred Value : 0.6190          
##              Prevalence : 0.8906          
##          Detection Rate : 0.8820          
##    Detection Prevalence : 0.9775          
##       Balanced Accuracy : 0.5589          
##                                           
##        'Positive' Class : 0               
## 

We could also look as those predicted probabilities.

rf_test_df <- rf_test_df %>%
  mutate(prob_no = unlist(map(probability, 1)),
         prob_yes = unlist(map(probability, 2)))

ggplot(rf_test_df, aes(x = prob_no, y = duration, fill = as.factor(y))) +
  geom_point(alpha = .5, color = "black", shape = 21)

K-fold cross validiation

Clearly a better way to evaluate a model is k-fold cross validation. This is pretty easy to do with sparklyr, and I’ve provided a function below to automate the process, at least for classification. The first thing we’ll need to do is create some partitions.

k = 5 # number of fold

# the proportion of the data to go in each fold
weights <- rep(1/k, times = k)
# our vector of weights needs to be named to work below
names(weights) <- paste0("fold", as.character(1:k))

# create a copy of df then partion it into 5 partitions 
# with 0.2 of the data in each one
df_cv <- df %>%
  compute("df_cv") %>%
  sdf_partition(weights = weights,
                seed = 2017)

Now, lets train a model on partitions 1 to 4, and test on 5, printing out the AUC at the end.

train_1 <- sdf_bind_rows(df_cv[1:4])

fit_rf_1 <- ml_random_forest(
  train_1,
  response = "y",
  features = colnames(df_cv$fold1)[1:16],
  impurity = "entropy",
  max.bins = 32L, # default = 32L
  max.depth = 6L, # default = 5L
  num.trees = 100L,  
  col.sample.rate = 0.25, # aka featureSubsetStrategy
  type = "classification",
  seed = 2017
)

# Get the AUC for the left out fold
sdf_predict(fit_rf_1, df_cv$fold5) %>%
  ml_binary_classification_eval(label = "y",
                                score = "probability")

## [1] 0.8876745

It would be a lot better to create a function to do cross validation where we have control over \(k\) and the model that is used.

ml_classification_cross_validation <- function(
  data,            # a tbl_spark
  response,        # the response column as a character
  features = NULL, # character vector or else all the columns except response
  model_fun,       # the modelling function to use (unquoted)
  k,               # number of folds
  ...) {           # additional arguments
  
  # first off create weights for partitioning the data
  weights <- rep(1 / k, times = k)
  # name the elements of weights
  names(weights) <- paste0("fold", as.character(1:k))
  
  # partition the data using weights
  data_cv <- sdf_partition(data, weights = weights)
  
  # get the indicies for our different training sets
  # e.g. a dataFrame for partitions 1, 2, 3, and 4
  K <- 1:k
  indices <- purrr::map(K, ~ K[-.x])
  
  # create our training sets by binding together partitions of data_cv
  # We get back a list where each element is a dataFrame of the partions
  # indexed by each element of indices combined together
  data_splits <- purrr::map(indices, ~ sdf_bind_rows(data_cv[.x]))
  
  # If a vector of feature names hasn't been specified
  if (is.null(features)) {
  # Get the column names of the data
  columns <- colnames(data_splits[[1]])
  # Create the feature names by using all the columns except the response
  features <- columns[columns != response]
  }
  
  # Map the specified model_fun over each of our training sets (the elements of
  # data splits)
  fits <- purrr::map(.x = data_splits,
                     .f = purrr::as_mapper(model_fun),
               response = response,
               features = features,
                   type = "classification",
                       ...)
  
  # Get some predictions for each model for the partiton that wasn't used
  # to train that model.
  preds <- purrr::map2(fits, K, ~ sdf_predict(.x, data_cv[[.y]]))
  
  # Checking whether a model that outputs probabilities was used.
  # If so the relevant evaluation function is used.
  # (Should also test for whether the response is binary but doesn't atm)
  if (any(stringr::str_detect(colnames(preds[[1]]), "probability"))) {
    
  evals <- purrr::map(preds,
                      ml_binary_classification_eval,
                      label = response,
                      score = "probability")
  } else {
    
  evals <- purrr::map(preds,
                      ml_classification_eval,
                      predicted_lbl = "prediction",
                      label = response)
  }
  # This is what will be returned
  list(fits = fits,
  predictions = preds,
  evals = evals)
  
}

Now this function can be used to do classification with whatever model we like. E.g. gradient boosted trees:

gbt_5_fold_cv <- 
  ml_classification_cross_validation(df, 
                                     response = "y",
                                     model_fun = ml_gradient_boosted_trees,
                                     k = 5,
                                     max.depth = 3L,
                                     num.trees = 10L,
                                     seed = 2017)

We can then get the average of our evaluation metric. Here it will be

mean(unlist(gbt_5_fold_cv$evals))

## [1] 0.8680992

Disconnect

Disconnect from our Spark cluster.

spark_disconnect(sc)

Further reading

Hopefully this post has been useful for providing an idea of what a sparklyr workflow could look like. In particular, the function to do cross validation should hopefully be useful for some people. Clearly, there’s loads of stuff I didn’t cover in this post. Most obviously, I haven’t gone into using dplyr verbs with a tbl_spark object. For examples of what this looks like you can see here. DataCamp also have a course on sparklyr, which I found very useful. The MLlib Main Guide is a good place to start for an overview of the features it offers, but you’ll want to look at the API Guide to get more information on the models.