Linear Model

First I am training a Logit LASSO model as my benchmark model, a type of linear regression that uses shrinkage to reduce model complexity.

I am going to use different lambdas, which defines how much weight we want to add to penalization (the larger lambda the higher the penalty). Cross validation is used to check the best lambda.

The data is pre-processed before training, which includes centering and scaling the data. I am going to set metric to “ROC” to choose models based on AUC.

# tuning parameters
train_control <- trainControl(method = "cv", 
                         number = 5,
                         classProbs = TRUE,
                         summaryFunction = twoClassSummary)

lasso_tg <- expand.grid("alpha" = 1, 
                        "lambda" = 10^seq(-1, -4, length = 10))

set.seed(1102)
logit_lasso_model <- train(
    formula(paste0("is_popular_f", variables)),
    data = temp_train,
    method = "glmnet",
    preProcess = c("center", "scale"),
    family = "binomial",
    trControl = train_control,
    tuneGrid = lasso_tg,
    metric = "ROC"
  )

# print(logit_lasso_model)

tuned_logit_lasso_model <- logit_lasso_model$finalModel
best_lambda <- logit_lasso_model$bestTune$lambda
lasso_coeffs <- as.matrix(coef(tuned_logit_lasso_model, best_lambda))

The LASSO algorithm picked a regression with 49 predictor variables/interactions and lambda = 0.0004641589. I will use the trained model to predict the probability of an article being shared or not on the holdout dataset. The default threshold for prediction is 50%. If we were to increase the threshold for predicting something to be positive, we will have fewer cases that we label as positive and therefor decrease the false positives and true positives.

# Probability prediction
holdout_prediction_probs <- predict.train(logit_lasso_model,
                                          newdata = temp_holdout, 
                                          type = "prob")

The ROC curve summarizes how a binary classifier performs “overall”, taking into account all possible thresholds. It shows the trade-off between true positive rate (a.k.a sensitivity, # true positives / # all positives) and the false positive rate (a.k.a 1 - specificity, # false positive / # negatives).

# ROC curve
rocr_prediction <- prediction(holdout_prediction_probs$yes,
                              temp_holdout[["is_popular_f"]])

plot(performance(rocr_prediction, "tpr", "fpr"), colorize=TRUE) 

AUC is the “area under the (ROC) curve”. This is a number between 0 and 1. Higher AUC generally means better classification.

# AUC
holdout_AUC[["logit_lasso"]] <- performance(rocr_prediction, 
                                            measure = "auc")@y.values[[1]]
# print(holdout_AUC[["logit_lasso"]])

The bigger the area under the ROC curve is, the better is our prediction. In this case we have an AUC of 0.7111.

# PREDICT on test_data
test_prediction_probs <- predict.train(logit_lasso_model, 
                                       newdata = test_data, 
                                       type = "prob")

# WRITE prediction on test set
filename <- paste0("predictions/",round(holdout_AUC[["logit_lasso"]]*1000),"_lasso.csv")
df <- data.frame(article_id = test_data$article_id,
                   score = test_prediction_probs[,2])
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

Random Forest

Next, I am going to build a Probability Random Forest to predict y.

Random Forest combines multiple decision trees, each fit to a random sample of the original data. It builds a predictive model without intervention to select variables or functional form. By combining multiple trees it reduces variance with minimal increase in bias.

I am going to use the default option of growing 500 trees, with several different tuning parameters, such as number of variables randomly chosen for any split in the tree or the minimum number of observations in the terminal nodes. To verify the optimal parameters I use 5-fold cross-validation. Here data scaling was not applied during training.

# tuning parameters
train_control <- trainControl(method = "cv", 
                         number = 5,
                         classProbs = TRUE,
                         summaryFunction = twoClassSummary)

tune_grid <- expand.grid(
  .mtry = c(5, 7, 9),
  .splitrule = "gini",
  .min.node.size = c(10, 15)
)
# suggested mtry: square root of the number of predictor variables


# train model "ranger"
set.seed(654)
rf_model_p <- train(
  formula(paste0("is_popular_f", variables)),
  method = "ranger",
  data = temp_train,
  tuneGrid = tune_grid,
  trControl = train_control,
  metric = "ROC",
  verbose = FALSE
)

# rf_model_p

The performance of the random forest model changed only marginally with different tuning paramaters. The best results were produced with the following parameters: 5 features to consider for each split and minimum 10 observations in each terminal node.

# Probability prediction
holdout_prediction_probs <- predict.train(rf_model_p, 
                                          newdata = temp_holdout, 
                                          type = "prob")

# AUC
rocr_prediction <- prediction(holdout_prediction_probs$yes,
                              temp_holdout[["is_popular_f"]])
holdout_AUC[["random_forest"]] <- performance(rocr_prediction,
                                              measure = "auc")@y.values[[1]]
# print(holdout_AUC[["random_forest"]])

The random forest model produced an AUC of 0.7213, which is an improvement to the previous model (Logit Lasso AUC: 0.7111).

# PREDICT on test_data
test_prediction_probs <- predict.train(rf_model_p, 
                                       newdata = test_data, type = "prob")

# WRITE prediction on test set
filename <- paste0("predictions/",round(holdout_AUC[["random_forest"]]*1000),"_rf.csv")
df <- data.frame(article_id = test_data$article_id,
                 test_prediction_probs[,2])
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

Gradient Boosting

Gradient boosting machines are also an ensemble of trees, however, the method of building the trees is different. The idea is to gradually (step-by-step) improve the trees by using the residuals of the previous tree. For not overfitting the data, a shrinkage parameter is constantly added, as well as a learning rate.

I am going to use different hyperparameters, which include the number of trees (the number of gradient boosting iteration), maximum nodes per tree, the learning rate and the minimum number of observations in trees’ terminal nodes.

tune_grid <- expand.grid(n.trees = c(100, 500, 1000), 
                        interaction.depth = c(2, 3, 4, 6), 
                        shrinkage = c(0.005, 0.01, 0.1),
                        n.minobsinnode = c(5,10))

set.seed(1234)
gbm_model <- train(formula(paste0("is_popular_f", variables)),
                   method = "gbm",
                   data = temp_train,
                   trControl = train_control,
                   tuneGrid = tune_grid,
                   preProcess = c("center", "scale", "pca"),
                   metric = "ROC",
                   verbose = FALSE
                   )

# gbm_model

The final values used for the model were n.trees = 1000, interaction.depth = 6, shrinkage = 0.01 and n.minobsinnode = 10.

# Probability prediction
holdout_prediction_probs <- predict.train(gbm_model, 
                                          newdata = temp_holdout, 
                                          type = "prob")

# AUC
rocr_prediction <- prediction(holdout_prediction_probs$yes,
                              temp_holdout[["is_popular_f"]])

holdout_AUC[["gradient_boosting"]] <- performance(rocr_prediction, measure = "auc")@y.values[[1]]
# print(holdout_AUC[["gradient_boosting"]])

The gradient boosting model produced an AUC of 0.7067, which shows no improvement compared to the previous models (Logit Lasso AUC: 0.7111, Random Forest AUC: 0.7213).

I also tried to remove pca variable reduction, which improved the AUC to 0.7338. This might be an indicator that variables are not very much correlated, which I had already assumed from the correlation table in the DataExplorer report.

# PREDICT on test_data
test_prediction_probs <- predict.train(gbm_model, newdata = test_data, type = "prob")

# WRITE prediction on test set
filename <- paste0("predictions/",round(holdout_AUC[["gradient_boosting"]]*1000),"_gbm.csv")
df <- data.frame(article_id = test_data$article_id,
                   score = test_prediction_probs[,2])
names(df) <- c("article_id", "score")
write.csv(df,file=filename, row.names = FALSE)

Neural Network

I am going to train my neural network model with the caret package and tune the hyper-parameters using cross validation. The hyperparameters of my neural net model are size and decay. Size is the number of units in hidden layer (nnet fit a single hidden layer neural network) and decay is the regularization parameter to avoid over-fitting. Data is preprocessed, meaning centered, scaled and dimensions are reduced using PCA.

tune_grid <- expand.grid(size = c(3, 5, 7, 10, 15),
                         decay = c(1, 1.5, 2, 2.5, 5, 7, 10))

set.seed(857)
nnet_model <- train(formula(paste0("is_popular_f", variables)),
                   method = "nnet",
                   data = temp_train,
                   trControl = train_control,
                   tuneGrid = tune_grid,
                   preProcess = c("center", "scale", "pca"),
                   metric = "ROC",
                   trace = FALSE) # avoid extensive iteration output

nnet_model

The final values used for the model were size = 7 and decay = 5.

# Probability prediction
holdout_prediction_probs <- predict.train(nnet_model,
                                          newdata = temp_holdout, 
                                          type = "prob")

# AUC
rocr_prediction <- prediction(holdout_prediction_probs$yes,
                              temp_holdout[["is_popular_f"]])

holdout_AUC[["nnet_model"]] <- performance(rocr_prediction, 
                                           measure = "auc")@y.values[[1]]
# print(holdout_AUC[["nnet_model"]])

The neural network model’s AUC of 0.7137 is similar to the lasso model (Logit Lasso AUC: 0.7111, Random Forest AUC: 0.7213, Gradient Boosting AUC: 0.7067).

Again I tried to train the same model without using PCA as the variables, but it only marginally improved the result with an AUC of 0.7177.

# PREDICT on test_data
test_prediction_probs <- predict.train(nnet_model, newdata = test_data, type = "prob")

# WRITE prediction on test set
filename <- paste0("predictions/",round(holdout_AUC[["nnet_model"]]*1000),"_nnet.csv")
df <- data.frame(article_id = test_data$article_id,
                   score = test_prediction_probs[,2])
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

For comparison I am going to train a neural network model using H2O. H2O is an open source, in-memory, distributed, fast, and scalable machine learning and predictive analytics platform.

First I am going to connect to a local H2O instance.

# Try to connect to a local H2O instance that is already running.
# If not found, start a local H2O instance from R with 8 gigabytes of memory.
h2o.init(max_mem_size = '8g') 

##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         22 hours 16 minutes 
##     H2O cluster timezone:       Europe/Vienna 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.28.0.2 
##     H2O cluster version age:    2 months and 7 days  
##     H2O cluster name:           H2O_started_from_R_lisahlmsch_htr830 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   5.82 GB 
##     H2O cluster total cores:    4 
##     H2O cluster allowed cores:  4 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     H2O API Extensions:         Amazon S3, XGBoost, Algos, AutoML, Core V3, TargetEncoder, Core V4 
##     R Version:                  R version 3.6.1 (2019-07-05)

For this model I have chosen 3 hidden layers of 32 neurons each. This is a rather small network but it will run faster and hopefully still produces good results, as in the previous example also a small number of neurons was chosen by the algorithm. I have chosen 10000 epochs, but expect the misclassification rate to converge earlier, which will stop the training due to the early-stopping rule I will implement.

h2o_nnet_model <- h2o.deeplearning(y = "is_popular_f",
                             x = names(train_data[,-c(59,61)]),
                             training_frame = as.h2o(temp_train),
                             validation_frame = as.h2o(temp_holdout),
                             activation = "Rectifier",
                             hidden=c(32,32,32), 
                             epochs = 10000,
                             stopping_rounds=2,
                             stopping_metric="misclassification",
                             stopping_tolerance=0.01,
                             reproducible = TRUE,
                             nfolds = 5,
                             keep_cross_validation_predictions = TRUE,
                             seed = 123)

holdout_AUC[["h2o_nnet_model"]] <- h2o.auc(h2o_nnet_model, valid = TRUE)
# print(holdout_AUC[["h2o_nnet_model"]])

The h2o neural network model was not able to produce a better result than the neural net model trained with caret: 0.7031 (h2o) versus 0.7137 (caret).

test_prediction_probs <- h2o.predict(h2o_nnet_model, 
                                     newdata = as.h2o(test_data))

filename <- paste0("predictions/",round(holdout_AUC[["h2o_nnet_model"]]*1000),"_h2o_nnet.csv")
df <- data.frame("article_id" = test_data$article_id,
                 "score" = as.data.table(test_prediction_probs[,3]))
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

Stacked model

Finally I am going to build a stacked model using h2o.

I am going to combine four models of different families using cross validation for hyperparameter tuning. These are my base models. With these models I am going to approximate the target variable. With the output of each model I am going to use some kind of combining mechanism to combine the results: I am going to use the preditions of my base models as features to train a new model - my metalearner. With this combination I am going make my final prediction. The metalearner is supposed to highlight each base learner where it performs better and where it performs worse. In short: Stacking is a model ensembling technique used to combine information from multiple predictive models to produce a new one.

Model 1 - Elastic net:
In this model I am again using grid search with cross-validation to tune the hyperparamter alpha and then extract the best models in terms of cross validated AUC.

hyper_params <- list(alpha = c(0, .25, .5, .75, 1))

glm_grid <- h2o.grid(y = "is_popular_f",
                     x = names(train_data[,-c(59,61)]),
                     training_frame = as.h2o(temp_train),
                     validation_frame = as.h2o(temp_holdout),
                     algorithm = "glm",
                     lambda_search = TRUE,
                     nfolds = 5,
                     seed = 123,
                     family = "binomial",
                     hyper_params = hyper_params,
                     keep_cross_validation_predictions = TRUE,
                     grid_id = "1glm")

# Get the grid results, sorted by AUC
glm_gridperf <- h2o.getGrid("1glm",
                              sort_by = "auc",
                              decreasing = TRUE) 

# print(glm_gridperf)

# Grab the model_id for the top GLM model, chosen by best AUC
best_glm_model_id <- glm_gridperf@model_ids[[1]]
best_glm <- h2o.getModel(best_glm_model_id)

# get validation AUC
holdout_AUC[["h2o_glm"]] <- h2o.auc(best_glm, valid = TRUE)
# holdout_AUC[["h2o_glm"]]

The best model in terms of AUC was an elastic net model (alpha = 0.25) and produced an AUC of 0.7108 on the validation data set.

test_prediction_probs <- h2o.predict(best_glm, newdata = as.h2o(test_data))

filename <- paste0("predictions/",round(holdout_AUC[["h2o_glm"]]*1000),"_h2o_glm.csv")
df <- data.frame("article_id" = test_data$article_id,
                 "score" = as.data.table(test_prediction_probs[,3]))
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

Model 2 - Random Forest: As in the previous random forest example I will be growing 500 trees. The default in h2o is 50 but I have incresed it because I will let the early stopping criteria decide when the random forest is sufficiently accurate. The algorithm will stop fitting new trees when the 2-tree average is within 0.001 (default) of the prior two 2-tree averages. As part of hyperparameter tuning I am going to specify different number of columns to randomly select at each level.

rf_params <- list(ntrees = c(500),
                  mtries = c(3, 5, 7, 9))

rf_grid <- h2o.grid(y = "is_popular_f",
                    x = names(train_data[,-c(59,61)]),
                    training_frame = as.h2o(temp_train),
                    validation_frame = as.h2o(temp_holdout), 
                    algorithm = "randomForest", 
                    nfolds = 5,
                    seed = 123,
                    keep_cross_validation_predictions = TRUE,
                    hyper_params = rf_params,
                    stopping_rounds = 2,
                    grid_id = "2rf")


# Get the grid results, sorted by AUC
rf_gridperf <- h2o.getGrid("2rf",
                              sort_by = "auc",
                              decreasing = TRUE)
# print(rf_gridperf)

# Grab the model_id for the top RF model, chosen by validation AUC
best_rf_model_id <- rf_gridperf@model_ids[[1]]
best_rf <- h2o.getModel(best_rf_model_id)

# get validation AUC
holdout_AUC[["h2o_rf"]] <- h2o.auc(best_rf, valid = TRUE)
# holdout_AUC[["h2o_rf"]]

The final model produced an AUC of 0.7278 on the validation set. It used five randomly selected columns at each level (mtries = 5) and 259 trees.

test_prediction_probs <- h2o.predict(best_rf, newdata = as.h2o(test_data))

filename <- paste0("predictions/",round(holdout_AUC[["h2o_rf"]]*1000),"_h2o_rf.csv")
df <- data.frame("article_id" = test_data$article_id,
                 "score" = as.data.table(test_prediction_probs[,3]))
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

Model 3 - Gradient Boosting:

For this model I am going to adjust some of the default parameters and use grid search to find the optimal hyperparameters. I will increase the number of trees from the default of 50 to 300. I will try different learning rates. Increasing the learning rate means the contribution of each tree will be stronger, so the model will move further away from the overall mean. I will also try different values for depth which adjusts the “weakness” of each learner. Adding depth makes each tree fit the data closer.

gbm_params <- list(learn_rate = c(0.01, 0.05, 0.1, 0.2),
                    max_depth = c(2, 3, 5, 10),
                    sample_rate = c(0.5), # row sample rate
                    col_sample_rate = c(0.5, 1.0))

# Train and validate a cartesian grid of GBMs
gbm_grid <- h2o.grid(y = "is_popular_f",
                     x = names(train_data[,-c(59,61)]),
                     training_frame = as.h2o(temp_train),
                     validation_frame = as.h2o(temp_holdout), 
                     algorithm = "gbm", 
                     nfolds = 5,
                     seed = 123,
                     ntrees = 300,
                     keep_cross_validation_predictions = TRUE,
                     hyper_params = gbm_params,
                     stopping_rounds = 2,
                     grid_id = "3gbm")

# Get the grid results, sorted by AUC
gbm_gridperf <- h2o.getGrid("3gbm",
                              sort_by = "auc",
                              decreasing = TRUE)
# print(gbm_gridperf)

# Grab the model_id for the top GBM model, chosen by validation AUC
best_gbm_model_id <- gbm_gridperf@model_ids[[1]]
best_gbm <- h2o.getModel(best_gbm_model_id)

holdout_AUC[["h2o_gbm"]] <- h2o.auc(best_gbm, valid = TRUE)
# holdout_AUC[["h2o_gbm"]]

The best GBM model produced an AUC of 0.7307 on the validation set. It used a learning rate of 0.01, max_depth = 10, col_sample_rate of 0.5 and 295 number of trees.

test_prediction_probs <- h2o.predict(best_gbm, newdata = as.h2o(test_data))

filename <- paste0("predictions/",round(holdout_AUC[["h2o_gbm"]]*1000),"_h2o_gbm.csv")
df <- data.frame("article_id" = test_data$article_id,
                 "score" = as.data.table(test_prediction_probs[,3]))
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

Model 4 - Deep learning:
I am going to include the previously trained neural net in my stacked model.

Validation performances:
Here is a summary of performances (evaluated on the validation set) for each of the different models and its optimal hyperparameters.

kable(data.frame(holdout_AUC)) %>%
  kable_styling()

Stacked ensemble model:
Recall, stacking is a combination of predictive models. The idea is that combined models should produce something more stable and have larger predictive performance (decrease the variance).

To produce a stacked ensemble model I am going to define several independent models and use the same folds for training each model. Then I do an out of fold prediction (produce prediction for these left out folds) and record the value we predicted as new features. With another predictive model I use these new features to estimate the outcome based on these baseline predictive scores. I am stacking the out of fold predictions and produce one single score from those features.

First I start with the baseline meta-learner - a glm model with non-negative weights (coefficients are constraints to be non-negative).

ensemble_model_glm <- h2o.stackedEnsemble(
                     y = "is_popular_f",
                     x = names(train_data[,-c(59,61)]),
                     training_frame = as.h2o(temp_train),
                     validation_frame = as.h2o(temp_holdout), 
                     base_models = list(best_glm,
                                        best_rf,
                                        best_gbm,
                                        h2o_nnet_model))

Then I use “gbm” as meta-learner.

ensemble_model_gbm <- h2o.stackedEnsemble(
                     y = "is_popular_f",
                     x = names(train_data[,-c(59,61)]),
                     training_frame = as.h2o(temp_train),
                     validation_frame = as.h2o(temp_holdout), 
                     base_models = list(best_glm,
                                        best_rf,
                                        best_gbm,
                                        h2o_nnet_model),
                     metalearner_algorithm = "gbm")

Validation Performances:

validation_performances <- list()
validation_performances[["ensemble_glm"]] = h2o.auc(h2o.performance(ensemble_model_glm, valid = TRUE))
validation_performances[["ensemble_gbm"]] = h2o.auc(h2o.performance(ensemble_model_gbm, valid = TRUE))

kable(data.frame(validation_performances)) %>% 
  kable_styling(full_width = FALSE)
  # column_spec(1:2,bold = T, color = "white", background = "darkgrey")

The first ensemble, using glm as meta-learner performed better than elastic net, random forest, gbm and deep learning.

# PREDICT on test_data
test_prediction_probs <- h2o.predict(ensemble_model_glm, newdata = as.h2o(test_data))

# WRITE prediction on test set
filename <- paste0("predictions/",round(validation_performances[["ensemble_glm"]]*1000),"_stacking.csv")
df <- data.frame("article_id" = test_data$article_id,
                 "score" = as.data.table(test_prediction_probs[,3]))
names(df) <- c("article_id", "score")
write.csv(df,file=filename,row.names=FALSE)

For my final model I therefor choose a stacked ensemble model with glm as meta learner and train it on the entire training dataset to predict my final results: