# Classification using Generalized Linear Models, Gradient Boosting Machines, Random Forests and Deep Learning in H2O

###### This tutorial demonstrates classification modeling in H2O using generalized linear models (GLM), gradient boosting machines (GBM), random forests and Deep Learning. It requires an installation of the h2o R package and its dependencies.

### Load the h2o R package and start an local H2O cluster

###### We will begin this tutorial by starting a local H2O cluster using the default heap size and as much compute as the operating system will allow.

    library(h2o)
    h2oServer <- h2o.init()

### Load the training and testing data into the H2O key-value store

###### This tutorial uses a 0.1% sample of the Person-Level 2013 Public Use Microdata Sample (PUMS) from United States Census Bureau with 75% of that sample being designated to the training data set and 25% to the test data set. This data set is intended to be an update to the [UCI Adult Data Set](https://archive.ics.uci.edu/ml/datasets/Adult).

    # we can skip this portion if continuing from Basics.Rmd!!!
    rootdir <- "/Users/spencer/meetups/introR"
    datadir <- file.path(rootdir, "/data")
    pumsdir <- file.path(datadir, "h2o-training", "pums2013")
    trainfile <- "adult_2013_train.csv.gz"
    testfile  <- "adult_2013_test.csv.gz"
    adult_2013_train <- h2o.importFile(h2oServer,
                                       path = file.path(pumsdir, trainfile),
                                       destination_frame = "adult_2013_train", sep = ",")
    adult_2013_test <- h2o.importFile(h2oServer,
                                      path = file.path(pumsdir, testfile),
                                      destination_frame = "adult_2013_test", sep = ",")
    dim(adult_2013_train)
    dim(adult_2013_test)
    
    adult_2013_train[,"TOP2_WAGP"] <- as.factor(adult_2013_train[,"TOP2_WAGP"])
    adult_2013_test[, "TOP2_WAGP"] <- as.factor(adult_2013_test[,"TOP2_WAGP"])
    

###### For the purposes of validation, we will create a single column data set containing only the target variable `TOP2_WAGP` from the test data set.

    actual_top2_wagp <- h2o.assign(adult_2013_test[, "TOP2_WAGP"],
                                   key = "actual_top2_wagp")
    rmValues(c("adult_2013","actual_top2_wagp"))

###### Also for our data set we have 8 columns that use integer codes to represent categorical levels so we will coerce them to factor after the data read.

    for (j in c("COW", "SCHL", "MAR", "INDP", "RELP", "RAC1P", "SEX", "POBP")) {
      adult_2013_train[[j]] <- as.factor(adult_2013_train[[j]])
      adult_2013_test[[j]]  <- as.factor(adult_2013_test[[j]])
    }

### Fit a basic logistic regression model
    
    top2_wagp_glm_relp <- h2o.glm(x = "RELP", y = "TOP2_WAGP",
                                  training_frame = adult_2013_train,
                                  model_id  = "top2_wagp_glm_relp",
                                  family = "binomial",
                                  lambda = 0)
    top2_wagp_glm_relp

### Generate performance metrics from a single model

###### We can create detailed performance metrics using the `h2o.performance` function.

#### Binary classifier performance metrics
#####![](images/metrics.png)

    pred_top2_wagp_glm_relp <- h2o.predict(top2_wagp_glm_relp, adult_2013_test)
    pred_top2_wagp_glm_relp

    prob_top2_wagp_glm_relp <- h2o.assign(pred_top2_wagp_glm_relp[, 3L],
                                          key = "prob_top2_wagp_glm_relp")
    rmValues(c("adult_2013","actual_top2_wagp", "prob_top2_wagp_glm_relp"))

    f1_top2_wagp_glm_relp <- h2o.performance(top2_wagp_glm_relp,adult_2013_test)
    f1_top2_wagp_glm_relp
    class(f1_top2_wagp_glm_relp)
    getClassDef("H2OBinomialMetrics")

### Plot Receiver Operating Characteristic (ROC) curve and find its Area Under the Curve (AUC)

###### A receiver operating characteristic (ROC) curve is a graph of the true positive rate ( *recall* ) against the false positive rate (1 - *specificity*) for a binary classifier.

    plot(f1_top2_wagp_glm_relp)

###### We can extract the area under the ROC curve from our model object to see how close it is to the ideal of 1.

    h2o.auc(top2_wagp_glm_relp)
    h2o.auc(f1_top2_wagp_glm_relp)

###### The *Gini coefficient*, not to be confused with the *Gini impurity* splitting metric used in decision tree based algorithms such as `h2o.randomForest`, is a linear rescaling of AUC defined by $Gini = 2 * AUC - 1$.

    h2o.giniCoef(f1_top2_wagp_glm_relp)
    2 * h2o.auc(f1_top2_wagp_glm_relp) - 1

###### Another metric we could have used to determine a cutoff is the lift metric, which is a ratio of the probability of a positive prediction given the model and the baseline probability for the population, that is generated by the `h2o.gains` function, which was inspired by Craig Rolling's [gains](http://cran.r-project.org/web/packages/gains/index.html) package hosted on CRAN.

    #h2o.gains(actual_top2_wagp, prob_top2_wagp_glm_relp)

### Making class predictions based on a probability cutoff

###### Using the cutoff that results in the largest harmonic mean of precision and recall (F1 score), we can now classify our predictions on the test data set.

    BEST_CUTOFF <- 0.026703
    class_top2_wagp_glm_relp <- prob_top2_wagp_glm_relp > BEST_CUTOFF
    class_top2_wagp_glm_relp <- h2o.assign(class_top2_wagp_glm_relp,
                                           key = "class_top2_wagp_glm_relp")

### Fit an elastic net logistic regression model

###### Now that we are familiar with H2O model fitting in R, we can fit more sophisticated models involving a larger set of predictors.

    addpredset <- c("COW", "MAR", "INDP", "RAC1P", "SEX", "POBP", "AGEP",
                    "WKHP", "LOG_CAPGAIN", "LOG_CAPLOSS")

###### In the context of elastic net regularization, we need to search the parameter space defined by the mixing parameter `alpha` and the shrinkage parameter `lambda`. To aide us in this search H2O can produce a grid of models for all combinations of a discrete set of parameters.

###### We will use different methods for specifying the `alpha` and `lambda` values as they are dependent upon one another. For the `alpha` parameter, we will specify five values ranging from 0 (ridge) to 1 (lasso) by increments of 0.25. For `lambda`, we will turn on an automated `lambda` search by setting `lambda = TRUE` and specify the number of `lambda` values to 10 by setting `nlambda = 10`.

    top2_wagp_glm <- h2o.glm(x = c("RELP_SCHL", addpredset),
                                  y = "TOP2_WAGP",
                                  training_frame = adult_2013_train,
                                  model_id  = "top2_wagp_glm_grid",
                                  family = "binomial",
                                  lambda_search = TRUE,
                                  nlambda = 10,
                                  alpha = 0.5)

###### We now have an object of class `H2OBinomialModel`.

    class(top2_wagp_glm)
    getClassDef("H2OBinomialModel")

###### Plot the ROC curve from the model on the training data

    plot(top2_wagp_glm)

###### Plot the ROC curve from the model on the test data

    plot(h2o.performance(top2_wagp_glm, adult_2013_test))

###### We can find the number of non-zero coefficients in our best logistic regression model to gain some understanding about its overall complexity.

    table(h2o.coef(top2_wagp_glm) != 0)
    nzcoefs <- h2o.coef(top2_wagp_glm)
    nzcoefs <- names(nzcoefs)[nzcoefs != 0]
    nzcoefs <- unique(sub("\\..*$", "", nzcoefs))
    setdiff(c("RELP_SCHL", addpredset), nzcoefs) # all preds had non-zero coefs

### Fit a gradient boosting machine binomial regression model

###### Given that not all relationships can be reduced to a linear combination or terms, we can compare the GLM results with that of a gradient (tree) boosting machine. As with the final GLM exploration, we will fit a grid of GBM models by varying the number of trees and the shrinkage rate and select the best model with respect to the test data set.

    top2_wagp_gbm <- h2o.gbm(x = c("RELP", "SCHL", addpredset),
                                  y = "TOP2_WAGP",
                                  training_frame = adult_2013_train,
                                  key  = "top2_wagp_gbm_grid",
                                  distribution = "bernoulli",
                                  n.trees = 1000,
                                  shrinkage = 0.1,
                                  validation = adult_2013_test,
                                  importance = TRUE)
    top2_wagp_gbm

    class(top2_wagp_gbm)
    slotNames(top2_wagp_gbm)

    h2o.performance(top2_wagp_gbm, adult_2013_test)

### Fit a random forest classifier

###### We will fit a single random forest model with 200 trees of maximum depth 10.

    top2_wagp_forest <- h2o.randomForest(x = c("RELP", "SCHL", addpredset),
                                         y = "TOP2_WAGP",
                                         training_frame = adult_2013_train,
                                         model_id  = "top2_wagp_forest",
                                         depth = 20,
                                         ntree = 500,
                                         validation = adult_2013_test,
                                         seed = 8675309)
    top2_wagp_forest

### Fit a deep learning classifier

###### Lastly we will fit a single Deep Learning model with default settings (more details about Deep Learning follow later) and compare the errors across the four model types.

    top2_wagp_dl <- h2o::h2o.deeplearning(x = c("RELP", "SCHL", addpredset),
                                     y = "TOP2_WAGP",
                                     training_frame = adult_2013_train,
                                     model_id  = "top2_wagp_dl",
                                     validation = adult_2013_test)
    top2_wagp_dl


    h2o.mse(h2o.performance(top2_wagp_glm, adult_2013_test))
    h2o.mse(h2o.performance(top2_wagp_gbm, adult_2013_test))
    h2o.mse(h2o.performance(top2_wagp_forest, adult_2013_test))
    h2o.mse(h2o.performance(top2_wagp_dl, adult_2013_test))

