Spark SQL transforms

These are feature transforms (aka mutating or filtering the columns/rows) using Spark SQL. This allows you to create new columns and modify existing columns while still employing the dplyr syntax. Here let’s modify 4 columns:

1. Family_Size - create number of siblings and parents
2. Pclass - format passenger class as character not numeric
3. Embarked - remove a small number of missing records
4. Age - impute missing age with average age

We use sdf_register() at the end of the operation to store the table in the Spark cluster.

titanic2_tbl <- titanic_tbl %>% 
  mutate(Family_Size = SibSp + Parch + 1L) %>% 
  mutate(Pclass = as.character(Pclass)) %>%
  filter(!is.na(Embarked)) %>%
  mutate(Age = if_else(is.na(Age), mean(Age), Age)) %>%
  sdf_register("titanic2")

Spark ML transforms

Spark also includes several functions to transform features. We can access several of them directly through sparklyr. For instance, to transform Family_Sizes into bins, use ft_bucketizer(). Because this function comes from Spark, it is used within sdf_mutate(), not mutate().

titanic_final_tbl <- titanic2_tbl %>%
  mutate(Family_Size = as.numeric(Family_size)) %>%
  sdf_mutate(
    Family_Sizes = ft_bucketizer(Family_Size, splits = c(1,2,5,12))
    ) %>%
  mutate(Family_Sizes = as.character(as.integer(Family_Sizes))) %>%
  sdf_register("titanic_final")

ft_bucketizer() is equivalent to cut() in R.

Train-validation split

Randomly partition the data into training/test sets.

# Partition the data
partition <- titanic_final_tbl %>% 
  mutate(Survived = as.numeric(Survived),
         SibSp = as.numeric(SibSp),
         Parch = as.numeric(Parch)) %>%
  select(Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked, Family_Sizes) %>%
  sdf_partition(train = 0.75, test = 0.25, seed = 1234)

# Create table references
train_tbl <- partition$train
test_tbl <- partition$test

Logistic regression

# Model survival as a function of several predictors
ml_formula <- formula(Survived ~ Pclass + Sex + Age + SibSp +
                        Parch + Fare + Embarked + Family_Sizes)

# Train a logistic regression model
(ml_log <- ml_logistic_regression(train_tbl, ml_formula))

## Call: Survived ~ Pclass_2 + Pclass_3 + Sex_male + Age + SibSp + Parch + Fare + Embarked_Q + Embarked_S + Family_Sizes_1 + Family_Sizes_2
## 
## Coefficients:
##    (Intercept)       Pclass_2       Pclass_3       Sex_male            Age 
##    3.696483969   -1.102340720   -2.053772367   -2.786065613   -0.035430315 
##          SibSp          Parch           Fare     Embarked_Q     Embarked_S 
##    0.158164571    0.457385475    0.002038998    0.214551594   -0.159123177 
## Family_Sizes_1 Family_Sizes_2 
##   -0.306388170   -3.749909900

Other machine learning algorithms

Run the same formula using the other machine learning algorithms. Notice that training times vary greatly between methods.

# Decision Tree
ml_dt <- ml_decision_tree(train_tbl, ml_formula)

# Random Forest
ml_rf <- ml_random_forest(train_tbl, ml_formula)

# Gradient Boosted Tree
ml_gbt <- ml_gradient_boosted_trees(train_tbl, ml_formula)

# Naive Bayes
ml_nb <- ml_naive_bayes(train_tbl, ml_formula)

# Neural Network
ml_nn <- ml_multilayer_perceptron(train_tbl, ml_formula, layers = c(11, 15, 2))

Validation data

# Bundle the models into a single list object
ml_models <- list(
  "Logistic" = ml_log,
  "Decision Tree" = ml_dt,
  "Random Forest" = ml_rf,
  "Gradient Boosted Trees" = ml_gbt,
  "Naive Bayes" = ml_nb,
  "Neural Net" = ml_nn
)

# Create a function for scoring
score_test_data <- function(model, data = test_tbl){
  pred <- sdf_predict(model, data)
  select(pred, Survived, prediction)
}

# Score all the models
ml_score <- map(ml_models, score_test_data)

Model lift

Lift compares how well the model predicts survival compared to random guessing. The function below calculates the model lift for each scored decile in the test data.

# Lift function
calculate_lift <- function(scored_data) {
  scored_data %>%
    mutate(bin = ntile(desc(prediction), 10)) %>% 
    group_by(bin) %>% 
    summarize(count = sum(Survived)) %>% 
    mutate(prop = count / sum(count)) %>% 
    arrange(bin) %>% 
    mutate(prop = cumsum(prop)) %>% 
    select(-count) %>% 
    collect() %>% 
    as.data.frame()
}

# Initialize results
ml_gains <- data_frame(
  bin = 1:10,
  prop = seq(0, 1, len = 10),
  model = "Base"
)

# Calculate lift
for(i in names(ml_score)){
  ml_gains <- ml_score[[i]] %>%
    calculate_lift %>%
    mutate(model = i) %>%
    bind_rows(ml_gains, .)
}

# Plot results
ggplot(ml_gains, aes(x = bin, y = prop, color = model)) +
  geom_point() +
  geom_line() +
  scale_color_brewer(type = "qual") +
  labs(title = "Lift Chart for Predicting Survival",
       subtitle = "Test Data Set",
       x = NULL,
       y = NULL)

The lift chart suggests the tree-based models (random forest, gradient boosted trees, and decision tree) provide the best prediction.

Accuracy and AUC

Receiver operating characteristic (ROC) curves are graphical plots that illustrate the performance of a binary classifier. They visualize the relationship between the true positive rate (TPR) against the false positive rate (FPR).

From Receiver operating characteristic

From Receiver operating characteristic

The ideal model perfectly classifies all positive outcomes as true and all negative outcomes as false (i.e. TPR = 1 and FPR = 0). The line on the second graph is made by calculating predicted outcomes at different cutpoint thresholds (i.e. \(.1, .2, .5, .8\)) and connecting the dots. The diagonal line indicates expected true/false positive rates if you guessed at random. The area under the curve (AUC) summarizes how good the model is across these threshold points simultaneously. An area of 1 indicates that for any threshold value, the model always makes perfect preditions. This will almost never occur in real life. Good AUC values are between \(.6\) and \(.8\). While we cannot draw the ROC graph using Spark, we can extract the AUC values based on the predictions.

# Function for calculating accuracy
calc_accuracy <- function(data, cutpoint = 0.5){
  data %>% 
    mutate(prediction = if_else(prediction > cutpoint, 1.0, 0.0)) %>%
    ml_classification_eval("prediction", "Survived", "accuracy")
}

# Calculate AUC and accuracy
perf_metrics <- data_frame(
  model = names(ml_score),
  AUC = 100 * map_dbl(ml_score, ml_binary_classification_eval, "Survived", "prediction"),
  Accuracy = 100 * map_dbl(ml_score, calc_accuracy)
  )
perf_metrics

## # A tibble: 6 x 3
##                    model      AUC Accuracy
##                    <chr>    <dbl>    <dbl>
## 1               Logistic 81.00928 82.43902
## 2          Decision Tree 87.78392 79.51220
## 3          Random Forest 86.90500 82.43902
## 4 Gradient Boosted Trees 85.89769 79.51220
## 5            Naive Bayes 66.47738 69.26829
## 6             Neural Net 77.32076 78.04878

# Plot results
gather(perf_metrics, metric, value, AUC, Accuracy) %>%
  ggplot(aes(reorder(model, value), value, fill = metric)) + 
  geom_bar(stat = "identity", position = "dodge") + 
  coord_flip() +
  labs(title = "Performance metrics",
       x = NULL,
       y = "Percent")

Overall it appears the tree-based models performed the best - they had the highest accuracy rates and AUC values.

Feature importance

It is also interesting to compare the features that were identified by each model as being important predictors for survival. The tree models implement feature importance metrics (a la randomForest::varImpPlot(). Sex, fare, and age are some of the most important features.

# Initialize results
feature_importance <- data_frame()

# Calculate feature importance
for(i in c("Decision Tree", "Random Forest", "Gradient Boosted Trees")){
  feature_importance <- ml_tree_feature_importance(sc, ml_models[[i]]) %>%
    mutate(Model = i) %>%
    mutate(importance = as.numeric(levels(importance))[importance]) %>%
    mutate(feature = as.character(feature)) %>%
    rbind(feature_importance, .)
}

# Plot results
feature_importance %>%
  ggplot(aes(reorder(feature, importance), importance, fill = Model)) + 
  facet_wrap(~Model) +
  geom_bar(stat = "identity") + 
  coord_flip() +
  labs(title = "Feature importance",
       x = NULL) +
  theme(legend.position = "none")

H2O and logistic regression

As a quick demonstration, let’s estimate a logistic regression model with 10-fold CV using H2O. First we need to load some additional packages:

library(rsparkling)
library(h2o)

We will reuse the previously modified Titanic table titanic_final_tbl. However to use it with H2O functions, we need to convert it to an H2O data frame:

titanic_h2o <- titanic_final_tbl %>% 
  mutate(Survived = as.numeric(Survived),
         SibSp = as.numeric(SibSp),
         Parch = as.numeric(Parch)) %>%
  select(Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked, Family_Sizes) %>%
  as_h2o_frame(sc, ., strict_version_check = FALSE)

Next we can estimate the logistic regression model using h2o.glm().

• This function does not use a formula to pass in the indepenent and dependent variables; instead they are passed as character vector arguments to x and y
• family = "binomial" - ensure we run logistic regression, not linear regression for continuous dependent variables
• training_frame - data frame containing the training set (here we use the entire data frame because we also use cross-validation)
• lambda_search = TRUE - argument for the optimizer function to calculate the parameter values
• nfolds = 10 - estimate the model using 10-fold cross-validation

glm_model <- h2o.glm(x = c("Pclass", "Sex", "Age", "SibSp", "Parch",
                           "Fare", "Embarked", "Family_Sizes"), 
                     y = "Survived",
                     family = "binomial",
                     training_frame = titanic_h2o,
                     lambda_search = TRUE,
                     nfolds = 10)

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=                                                                |   1%
  |                                                                       
  |========================                                         |  37%
  |                                                                       
  |=================================================================| 100%

glm_model

## Model Details:
## ==============
## 
## H2OBinomialModel: glm
## Model ID:  GLM_model_R_1503421695031_1 
## GLM Model: summary
##     family  link                                regularization
## 1 binomial logit Elastic Net (alpha = 0.5, lambda = 0.003129 )
##                                                                    lambda_search
## 1 nlambda = 100, lambda.max = 0.248, lambda.min = 0.003129, lambda.1se = 0.03515
##   number_of_predictors_total number_of_active_predictors
## 1                          4                           4
##   number_of_iterations                                   training_frame
## 1                   68 frame_rdd_13566_a40ade0b999e58d83ef2c902891449bb
## 
## Coefficients: glm coefficients
##       names coefficients standardized_coefficients
## 1 Intercept    -0.236119                 -0.454621
## 2       Age    -0.022417                 -0.290716
## 3     SibSp    -0.289895                 -0.319958
## 4     Parch     0.101851                  0.082169
## 5      Fare     0.017416                  0.865539
## 
## H2OBinomialMetrics: glm
## ** Reported on training data. **
## 
## MSE:  0.2084469
## RMSE:  0.4565598
## LogLoss:  0.6122621
## Mean Per-Class Error:  0.3641755
## AUC:  0.7074815
## Gini:  0.414963
## R^2:  0.117433
## Residual Deviance:  1088.602
## AIC:  1098.602
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##          0   1    Error      Rate
## 0      288 261 0.475410  =261/549
## 1       86 254 0.252941   =86/340
## Totals 374 515 0.390326  =347/889
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.322308 0.594152 268
## 2                       max f2  0.172284 0.759268 395
## 3                 max f0point5  0.400116 0.645032 173
## 4                 max accuracy  0.400116 0.724409 173
## 5                max precision  0.999640 1.000000   0
## 6                   max recall  0.172284 1.000000 395
## 7              max specificity  0.999640 1.000000   0
## 8             max absolute_mcc  0.399184 0.393883 174
## 9   max min_per_class_accuracy  0.346658 0.632353 236
## 10 max mean_per_class_accuracy  0.399184 0.677215 174
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## 
## H2OBinomialMetrics: glm
## ** Reported on cross-validation data. **
## ** 10-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
## 
## MSE:  0.2108965
## RMSE:  0.4592347
## LogLoss:  0.6189053
## Mean Per-Class Error:  0.4220026
## AUC:  0.6925587
## Gini:  0.3851173
## R^2:  0.1070613
## Residual Deviance:  1100.414
## AIC:  1110.414
## 
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
##          0   1    Error      Rate
## 0      147 402 0.732240  =402/549
## 1       38 302 0.111765   =38/340
## Totals 185 704 0.494938  =440/889
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.297723 0.578544 304
## 2                       max f2  0.152770 0.757238 395
## 3                 max f0point5  0.413114 0.625000 151
## 4                 max accuracy  0.413114 0.713161 151
## 5                max precision  0.999707 1.000000   0
## 6                   max recall  0.152770 1.000000 395
## 7              max specificity  0.999707 1.000000   0
## 8             max absolute_mcc  0.399255 0.368296 167
## 9   max min_per_class_accuracy  0.346688 0.623529 228
## 10 max mean_per_class_accuracy  0.392200 0.670556 173
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## Cross-Validation Metrics Summary: 
##                 mean          sd cv_1_valid cv_2_valid cv_3_valid
## accuracy   0.6500479 0.041881736  0.6363636  0.6091954 0.71428573
## auc       0.71170306 0.032053903 0.74575895  0.6638825 0.70077866
## err        0.3499521 0.041881736 0.36363637  0.3908046  0.2857143
## err_count       31.0   3.7549968       36.0       34.0       26.0
## f0point5   0.5726497 0.045476418 0.53505534  0.5106383 0.56213015
##           cv_4_valid cv_5_valid cv_6_valid cv_7_valid cv_8_valid
## accuracy  0.55172414  0.6835443   0.682243  0.5833333  0.7582418
## auc       0.70297295  0.7615283  0.7119403  0.6064706 0.74852073
## err       0.44827586  0.3164557   0.317757 0.41666666 0.24175824
## err_count       39.0       25.0       34.0       35.0       22.0
## f0point5    0.542522 0.62211984 0.58189654  0.5252101 0.72192514
##           cv_9_valid cv_10_valid
## accuracy  0.61904764      0.6625
## auc         0.719496   0.7556818
## err        0.3809524      0.3375
## err_count       32.0        27.0
## f0point5         0.5       0.625
## 
## ---
##                         mean          sd cv_1_valid cv_2_valid  cv_3_valid
## precision          0.5381394 0.055250388  0.4915254 0.47058824  0.54285717
## r2                 0.0975009 0.031127375  0.1564119 0.10831892 0.086596675
## recall             0.8112482 0.086222775 0.82857144  0.7741935   0.6551724
## residual_deviance  109.25886   5.4304786   114.6542  104.17874   106.71968
## rmse              0.45955113 0.008883104  0.4390905 0.45223105  0.44533285
## specificity        0.5446867  0.11514253    0.53125 0.51785713   0.7419355
##                   cv_4_valid cv_5_valid cv_6_valid  cv_7_valid  cv_8_valid
## precision          0.4868421  0.5869565     0.5625  0.49019608   0.7297297
## r2                 0.1278114 0.09742886 0.13443607 0.027979773 0.121132664
## recall                   1.0  0.8181818      0.675   0.7352941   0.6923077
## residual_deviance  108.03258  101.55076  128.32715  111.066345   112.77616
## rmse               0.4617127 0.46854305 0.45012507  0.48393032  0.46393192
## specificity             0.22  0.5869565  0.6865672        0.48   0.8076923
##                    cv_9_valid cv_10_valid
## precision          0.44642857  0.57377046
## r2                0.008515219  0.10637745
## recall             0.96153843   0.9722222
## residual_deviance   103.23352  102.049416
## rmse               0.46032485   0.4702888
## specificity        0.46551725   0.4090909

We get lots of information back about the model. Many of these statistics can be extracted and stored as tidy data frames or used to create visualizations.