<jupyter_start><jupyter_text>INFO-F-422 -  Statistical Foundations of Machine Learning > Student 1 - __[<EMAIL>](mailto:<EMAIL>) - Student ID 000528423__ / Driven Data Id : JulesNT> Student 2 - __[<EMAIL>](mailto:<EMAIL>) - Student ID 000528422__ / Driven Data Id : LucasBaussay> Student 3 - __[<EMAIL>](mailto:<EMAIL>) - Student ID 000528421__ / Driven Data Id : mhavardsedeno> Video presentation: https://youtu.be/rOhgU5WBHB8> Driven Data team : The French Dowsers Pump it Up : Data Mining the Water Table<jupyter_code>library(tidyverse)
library(greybox)
library(rpart)
library(randomForest)
library(ramify)
library(nnet)
library(devtools)
library(dummies)
library(ggplot2)
library(dplyr)
library(caret)
library(xgboost)<jupyter_output><empty_output><jupyter_text>Generic tools for the execution<jupyter_code>showMissingValues <- function(data,categorical_names) {
    missingVec = c()
    ratioVec = c()
    names = c()
    for (name in categorical_names) {
        nbMissingValues = sum(data[,name]=="")
        ratio = nbMissingValues/nrow(data)*100
        tryCatch(
            expr = {
                if (nbMissingValues > 0) {
                    missingVec = c(missingVec,nbMissingValues)
                    ratioVec = c(ratioVec,ratio)
                    names = c(names,name)
                }
            },
            error = function(e){},
            warning = function(w){},
            finally = {}
        )
    }
    table.results <- data.frame(feature=names,missing.values=missingVec,percentage.missing=ratioVec)
    return (table.results)
}

cramerAssociation <- function(data,categorical_names,threshold) {
    options(warn=-1)
    result = c()
    names = c()
    for (name in categorical_names) {
        cramer_value <- cramer(data[,name],data[,"status_group"])$value
        if (cramer_value < threshold){
            result = c(result,toString(cramer(data[,name],data[,"status_group"])$value))
            names = c(names,name)
        }
    }
    options(warn=0)
    table.results <- data.frame(features = names, association.statusgroup = result)
    return(table.results)
}

crossTableAssociation <- function(data,categorical_names,treshold){
    for(i in 1:length(categorical_names)) {
        var_name_1 = categorical_names[i]
        if (i < length(categorical_names)) {
            for(j in (i+1):length(categorical_names)) {
                var_name_2 = categorical_names[j]
                value <- cramer(data[,var_name_1],data[,var_name_2])$value
                if (value > treshold) {
                    cat(paste("[",var_name_1,"] is associated with [",var_name_2,"] -> ",value,"\n"))
                }
            }
        }
    }
}<jupyter_output><empty_output><jupyter_text>Predicting functional state of water-pumps in Tanzania > Tanzania has been in a water crisis for several years now. Even though the country partly includes three of the Great African Lakes, a third of Tanzania is arid desert. According to the most recent survey, almost half of the population doesn’t have access to safe water. This leads to several problem, not only about health and hygiene. Since women and children have to walk longer to collect water, it also impacts the access to education, the economic status and the happiness in the country.> One of the answer to this crisis is the installation of water pump all across the country. Plenty of data about it is available in a large database. In order to acknowledge the status of theses water-pumps, a survey was conducted by the ministry of water. This resulted in the classification of the pumps in three categories : functional, requires repair work and not functional.> Our goal today is to analyze the dataset provided and use prediction models in order to predict the state of the water-pumps. If we are able to correctly predict the state of the water-pump, this could save precious time and money to the government.   Loading and exploring the data>The features of the dataset are documented on the page of the competition : [Features description](https://www.drivendata.org/competitions/7/pump-it-up-data-mining-the-water-table/page/25/)<jupyter_code>train <- read.csv("../data/TrainingSetValues.csv")
trainLabels <- read.csv("../data/TrainingSetLabels.csv")
test <- read.csv("../data/TestSetValues.csv")<jupyter_output><empty_output><jupyter_text>Dimensions of the dataset<jupyter_code>dim(train)
dim(test)<jupyter_output><empty_output><jupyter_text>Structure of the dataset Quick look at the features values>We can see a lot more categorical features than numeric ones.<jupyter_code>#head(train)<jupyter_output><empty_output><jupyter_text>Distribution of the target feature> We observe here a strong inbalance between the different outputs.<jupyter_code>table(trainLabels$status_group)
plot(trainLabels$status_group)<jupyter_output><empty_output><jupyter_text>Types of the features and number of factors>The use of the function str() shows that :>- most features are categorical features and some of them have a huge number of factors.>- some groups of features apparently describe the same information (e.g. region and region_code)<jupyter_code>#str(train)
#summary(train)<jupyter_output><empty_output><jupyter_text>Data preprocessing> Main objective : maintain a good accuracy while reducing computation time and overfitting>- Remove features and factors that don't bring new information to the model>- Add more relevent features> Here are the main parts of the preprocessing pipeline :>- Combining the training set and the test set>- Handling missing values>- Reducing levels of categorical features>- Format dates>- Remove obviously useless features>- Remove features little associated with target value>- Identify similar features and choose only one per group  Combining the training set and the test set>This way the training and test set have the same factors and we don't encounter problems when predicting the test labels with the model. If we don't combine the two sets and the test set contains new factors, the model won't be able to predict anything.<jupyter_code>dataSet <- rbind(train,test)

dim(dataSet)<jupyter_output><empty_output><jupyter_text>>The following will help targeting a certain type of feature<jupyter_code>categorical_columns <- c(3,4,6,9,11,12,13,16,17,19,20,21,22,23,c(25:ncol(dataSet)))
categorical_names <- names(dataSet)[categorical_columns]
missing_values_alias <- c("none","None","No","no","0","Unknown",
                          "unknown","Not known","Not Known","not known",
                          "missing","Missing")<jupyter_output><empty_output><jupyter_text>Handling missing values Harmonize the dataset>First we replace alias of missing values (e.g. none, no, unknown) by the value "" to harmonize the dataSet. Indeed, in the feature "funder" a NA was initially described by the values "", "none", "Unknown", "0", "Not Known" and "No". It's thus important to only have one specific value across all the dataset.<jupyter_code>for (var in categorical_columns[categorical_columns!=9 & categorical_columns!=38 & (categorical_columns < 28 | categorical_columns > 36)]) {
    for (alias in missing_values_alias) {
        dataSet[var][dataSet[var]==alias] <- ""
    }
}
dataSet <- droplevels(dataSet)<jupyter_output><empty_output><jupyter_text>Proportion of missing values> In this section we compute the percentage of missing values for each feature that contains missing values.<jupyter_code>table.results <- showMissingValues(dataSet,categorical_names)
sorted.results <- table.results[order(table.results$missing.values,decreasing=T),]
sorted.results
ggplot(sorted.results,aes(x=percentage.missing,reorder(feature,percentage.missing)))+geom_bar(stat ="identity") + theme_light()<jupyter_output><empty_output><jupyter_text>Removing features that contains more than 20% of missing values>scheme_name contains 47% of missing values, we remove it from the dataSet and see how many rows still contain missing values.<jupyter_code>data <- dataSet
data[data==""] <- NA
dim(na.omit(data))

dataSet <- dataSet[, -which(names(dataSet) == "scheme_name")]

data <- dataSet
data[data==""] <- NA
dim(na.omit(data))<jupyter_output><empty_output><jupyter_text>Conclusion on missing values>Because we have to give a prediction on each of the individuals present in the test set, we can't remove lines from the dataset.>Thus we have to choose between imputing missing values and treating them as a separate category for each variable. The two options have their downsides, mostly adding bias to the data.>We choose to treat them as a separate category.  Reducing levels of categorical features>Some categorical features have a lot of factors. A large number of factors leads to high computation times and overfitting. Thus, those factors need to be reduced before handling the data to the model.>Finding a maximum number of factors is like making a trade of between underfitting and overfitting.>In the following we plot the frequencies of the 50 most frequent factors in each variable. This helps visualize how many factors we need to keep. The factor selection possess two objectives. The idea here is to minimize the number of factors we keep while maximizing the number of frequent factors we keep.>The method to reduce factors is the following : we observed the frequencies of the most frequent factors and selected the minimum number of factors that really stand out from the others. To sum up, we got rid of all the factors that were too little frequent. We used plots to visualise.>In the code section below, we show too different types of plots. The first shows a lot of useless factors that are going to be removed, the second one shows much less useless factors.<jupyter_code># Very flat distribution from 5 to 50
plot(summary(dataSet$installer)[1:50])
# Fairly steep decreasing slope
plot(summary(dataSet$subvillage)[1:50])<jupyter_output><empty_output><jupyter_text>>According to the plots we saw, we set a maximum number of factors for each feature. All the ones that are being removed are replaced by a factor "Other"<jupyter_code>factors_names = c("installer","funder","wpt_name",
                  "subvillage","lga","ward","extraction_type_class",
                  "extraction_type","extraction_type_group",
                  "quantity","quantity_group","waterpoint_type","waterpoint_type_group")
nb_levels = c(5,5,6,9,9,9,4,8,7,4,4,4,4) # + 1

names = list(
    installer_names = names(summary(dataSet$installer)[1:nb_levels[1]]),
    funder_names <- names(summary(dataSet$funder)[1:nb_levels[2]]),
    wpt_name_names <- names(summary(dataSet$wpt_name)[1:nb_levels[3]]),
    subvillage_names <- names(summary(dataSet$subvillage)[1:nb_levels[4]]),
    lga_names <- names(summary(dataSet$lga)[1:nb_levels[5]]),
    ward_names <- names(summary(dataSet$ward)[1:nb_levels[6]]),
    extraction_type_class_names <- names(summary(dataSet$extraction_type_class)[1:nb_levels[7]]),
    extraction_type_names <- names(summary(dataSet$extraction_type)[1:nb_levels[8]]),
    extraction_type_group_names <- names(summary(dataSet$extraction_type_group)[1:nb_levels[9]]),
    t1 <- names(summary(dataSet$quantity)[1:nb_levels[10]]),
    t2 <- names(summary(dataSet$quantity_group)[1:nb_levels[11]]),
    t3 <- names(summary(dataSet$waterpoint_type)[1:nb_levels[12]]),
    t4 <- names(summary(dataSet$waterpoint_type_group)[1:nb_levels[13]])
)

index = 1
for (var in factors_names){
    variable <- factor(dataSet[,var], levels=c(unlist(unname(names[index])), "Other"))
    variable[is.na(variable)] <- "Other"
    dataSet[,var] <- variable
    index = index + 1
}<jupyter_output><empty_output><jupyter_text>Format dates>In the dataset, date_recorded is a categorical feature in the format : "year-month-day". This format can't truly be interpreted by a model. Indeed, the model can't extract, from this format, the useful informations like year, month, day. We thus need to remove this variable format and replace it by a new one.>There are at least two ways to modify the format of the date :>- We can add tree new features, year, month and day. This way the model can exactly know when the line was recorded.>- We can add a new feature describing the number of days that have passed between the date and the oldest date in the dataset (we set the lower bound to 2001-01-01)>Each method brings information that the other doesn't. For instance, the first one explicitly differenciate time entities like months and years. If a month is highly correlated with a status "non_functional", the model will quickly know about it. In the other hand, it is easier with the second method to know precisely how old is the date recording.>We choose to create all the four new features and wait for the futures analysis to tell us if some features need to be removed.>We leave the feature contruction_year unchanged because having an integer to describe the year of the pump construction is already a good thing.<jupyter_code># Computing the new features
days_since_2001 <- as.numeric(as.Date(dataSet$date_recorded) - as.Date("2001-01-01"))
year_recorded <- as.numeric(format(as.Date(dataSet$date_recorded), "%Y"))
month_recorded <- as.numeric(format(as.Date(dataSet$date_recorded), "%m"))
day_recorded <- as.numeric(format(as.Date(dataSet$date_recorded), "%d"))
# Removing the old feature
dataSet <- dataSet[, -which(names(dataSet) == "date_recorded")]
# Adding the new features to the dataset
dataSet <- cbind(dataSet, year_recorded)
dataSet <- cbind(dataSet, month_recorded)
dataSet <- cbind(dataSet, day_recorded)
dataSet <- cbind(dataSet, days_since_2001)<jupyter_output><empty_output><jupyter_text>Remove obviously useless features>Here we remove features that are clearly useless for the learning process :>- num_private has not documentation and is full of zeros.>- recorded_by is constant> IMPORTANT : we won't use the id feature when running the model, but we need it on the dataSet to identify lines. We remove id just before running a model, inside the learning functions.<jupyter_code>dataSet <- dataSet[, -which(names(dataSet) == "num_private")]
dataSet <- dataSet[, -which(names(dataSet) == "recorded_by")]<jupyter_output><empty_output><jupyter_text>Remove features little associated with target value>In this section we use the cramer V measure of association between two categorical features. The function cramerAssociation(data,features_studied,threshold), defined at the beginning of the notebook, computes the cramer V association between the features and the target value, and prints all the features that have an association value smaller than 0.2.>The cramer V is not designed to measure association between continuous and categorical features, but we didn't find another measure that worked for us. Point biserial correlation is only taking binary features instead of categorical features, and Pearson’s correlation coefficient is only taking continous features. Thus, the results on the association between continous features like amout_tsh and the target value are to be studied carefully.>The Cramer's V measure of association is computed as follows :\begin{align}V & = 𝑠𝑞𝑟𝑡(\frac{𝜒2}{(𝑛𝑞−1)}) \\q & = min(nrow,ncol) \\n & = grand.total.of.observations\end{align}<jupyter_code>train <- dataSet[1:nrow(train),]
test <- dataSet[(nrow(train)+1):nrow(dataSet),]

train <- merge(train,trainLabels)

table.results <- cramerAssociation(train,names(train),threshold=0.2)

head(table.results)

# See all the values
# table.results
# Remove features less associated with target value than 20%
for (name in table.results$features) {
    dataSet <- dataSet[, -which(names(dataSet) == name)]
}
dim(dataSet)<jupyter_output><empty_output><jupyter_text>Identify similar features and choose only one per group>In this section, we look at the association between paired features to identify groups of similar features. We set the association threshold to 0.75, meaning the function crossTableAssociation will signal similar features if their cramer V association is greater than 0.75.>We want to take only one variable per group.<jupyter_code>categorical_names = c("region","extraction_type_group",
                      "extraction_type_class","quantity",
                      "quantity_group","waterpoint_type_group")
crossTableAssociation(dataSet,categorical_names,0.75)<jupyter_output>[ extraction_type_group ] is associated with [ extraction_type_class ] ->  0.925697664597424 
[ quantity ] is associated with [ quantity_group ] ->  1<jupyter_text>>We thus keep only the best variable of each group of similar features. We keep the one having the highest association with the target value.<jupyter_code>table.results <- cramerAssociation(train,names(train),threshold=1)
table.results[table.results$features=="quantity" | table.results$features=="quantity_group" |
              table.results$features=="extraction_type_group" | table.results$features=="extraction_type_class",]<jupyter_output><empty_output><jupyter_text>>We choose **extraction_type_group** and **quantity**<jupyter_code>dataSet <- dataSet[, -which(names(dataSet) == "quantity_group")]
dataSet <- dataSet[, -which(names(dataSet) == "extraction_type_class")]<jupyter_output><empty_output><jupyter_text>>We also remove the variable region_code that is clearly a alias of region.<jupyter_code>dataSet <- dataSet[, -which(names(dataSet) == "region_code")]<jupyter_output><empty_output><jupyter_text>Last check on the data>Here we give a last look to the data, making sure that:>- features have few levels>- features have an association with status_group greater than 20%.<jupyter_code>train <- dataSet[1:nrow(train),]
test <- dataSet[(nrow(train)+1):nrow(dataSet),]

train <- merge(train,trainLabels)

head(train)

dim(train)
dim(test)

str(train)

cat(paste("\n"))

cramerAssociation(train,names(train[,-1]),1)<jupyter_output><empty_output><jupyter_text>The data is ready !  2) Model selection>In this section we estimate the accuracy of three different models, using the preprocessed data constructed in part 1.>The models we will study are :>- **single-hidden-layer neural network**>- **decision tree**>- **random forest**>We estimate the classification rate of each model by running a **10 fold cross validation** procedure, except for the random forest because of the high computation cost. We thus operate a **5 fold cross validation** to validate the random forest model, to avoid a crash from jupyter kernel.  Choose your execution parameter>You can set the parameter **EXECUTION** to a value between 1,2 and 3. One being the minimal run, and 3 the expensive one. By default the parameter is set to 2, where all the important algorithms are called.<jupyter_code>#EXECUTION <- 1 # MINIMAL : all the lines that take more than 1 minute of compilation won't run (VERY FAST)
EXECUTION <- 2 # NORMAL : all the main algorithms will run, but the hyperparameters tunings won't
#EXECUTION <- 3 # LONG : all the lines will run (VERY EXPENSIVE)<jupyter_output><empty_output><jupyter_text>Single-Hidden-Layer neural network>The single-hidden-layer neural network structure is explicit. The model is a neural network containing only one hidden layer. We choose to allow skip connections (parameter skip = TRUE) to minimize vanishing gradient.>Because the performance of this model highly depends on it's number of nodes, we choose to study the impact of the number of nodes on the classification rate. To limit computation time, we don't go further than 10 hidden nodes, testing only 3, 5 and 10 hidden nodes.>We obtain the best results (classification rate around 60%) with 5 hidden nodes. You can verify this by running the **neural_network_optimize_hidden_nodes(data,list_nodes)**, wich takes some time to execute.>We choose to use integer encoding for the categorical features. The nnet function automatically encode the categorical features using integer encoding. Thus we don't need to modify our existing data before giving it to the nnet model.<jupyter_code>neural_network_fold_cross_validation <- function(trainingSet,k,nb_hidden_nodes) {
    
    # INITIALISATION
    
    start_time <- Sys.time()
    # removing column id in dataset
    trainingSet <- trainingSet[, -which(names(trainingSet) == "id")]
    # indexing status_group
    target_variable <- ncol(trainingSet)
    # vector storing the classification rates
    classification_rate <- array(0,k)
    # shuffle trainingSet
    train_idx <- sample(1:nrow(trainingSet))
    # split the trainingset into k groups
    max <- ceiling(nrow(trainingSet)/k)
    splits <- split(train_idx, ceiling(seq_along(train_idx)/max))
    
    # 10 FOLD CROSS VALIDATION

    for (i in 1:k){
        # separate the training set from the fold test set
        test_data <- trainingSet[splits[[i]],]
        train_data <- trainingSet[-splits[[i]],]
        
        # Fit the model to the training set
        model <- model <- nnet(status_group ~ .,
                                 data=train_data,
                                 size=nb_hidden_nodes,
                                 skip=TRUE,
                                 MaxNWts=10000,
                                 trace=F)

        # Predict the label of test_data
        Y_pred <- predict(model,test_data[,-target_variable])
        Y_hat <- ifelse(Y_pred[,1] >= Y_pred[,2] & Y_pred[,1] >= Y_pred[,3],
                     "functional",
                     ifelse(Y_pred[,2] >= Y_pred[,3], "functional needs repair", "non functional"))
        
        # True label
        Y <- test_data[,target_variable]

        # Compute confusion matrix 
        confusion_matrix <- table(Y_hat,Y)
        
        # Compute the classification rate
        if ((dim(confusion_matrix)[1]) == 1){
            classification_rate[i] = (confusion_matrix[1,1])/length(Y)
        } else if((dim(confusion_matrix)[1]) == 2) {
            classification_rate[i] = (confusion_matrix[1,1]+confusion_matrix[2,2])/length(Y)
        } else {
            classification_rate[i] = (confusion_matrix[1,1]+confusion_matrix[2,2]+confusion_matrix[3,3])/length(Y)
        }
    }
    
    end_time <- Sys.time()
    print(end_time - start_time)
    
    return(mean(classification_rate))
}
neural_network_optimize_hidden_nodes <- function(trainingSet,k,values_hidden_nodes) {
    classification_rate <- c()
    # Computing a 10 fold cross validation for each number of hidden nodes
    for (nb_hidden_nodes in values_hidden_nodes) {
        classification_rate <- c(classification_rate,neural_network_fold_cross_validation(trainingSet,k,nb_hidden_nodes))
    }
    cat(paste("\nClassification rates :\n"))
    print(classification_rate)
    plot(classification_rate)
    bestIndex <- argmax(mat(classification_rate,nrow=1))
    bestClassificationRate <- max(classification_rate)
    cat(paste("The neural network works best with",
              values_hidden_nodes[bestIndex],
              "hidden nodes, obtaining a classification rate of",bestClassificationRate))
    return(bestClassificationRate)
}
if (EXECUTION == 2) {
    accuracy_nn = neural_network_fold_cross_validation(train,10,5)
    accuracy_nn
}

if (EXECUTION == 3) {
    best_accuracy_nn = neural_network_optimize_hidden_nodes(train,10,c(3,5,10))
    best_accuracy_nn
}<jupyter_output>Time difference of 2.009154 mins<jupyter_text>Decision tree>This model is a simple single decision tree. We chose this model to better illustrate the improvement of the random forest over a single tree.>Tunning of hyperparameters : we studied the impact of the depth of the tree over the classification rate. We obtained best results with ??? (classification rate around ???%). You can replicate this experiment by running the function **decision_tree_optimize_depth(data,list_depths,verbose)**, which again is quite expensive in terms of computation times.<jupyter_code>decision_tree_fold_crossValidation <- function(trainingSet,k,depth,verbose) {
    
    # INITIALISATION
    
    start_time <- Sys.time()
    # removing column id in dataset
    trainingSet <- trainingSet[, -which(names(trainingSet) == "id")]
    # indexing status_group
    target_variable <- ncol(trainingSet)
    # vector storing the classification rates
    classification_rate <- array(0,k)
    # Shuffle trainingSet
    train_idx <- sample(1:nrow(trainingSet))
    # Split the dataset into k groups
    max <- ceiling(nrow(trainingSet)/k)
    splits <- split(train_idx, ceiling(seq_along(train_idx)/max))
    
    # 10 FOLD CROSS VALIDATION

    for (i in 1:k){
        # separate the training set from the fold test set
        test_data <- trainingSet[splits[[i]],]
        train_data <- trainingSet[-splits[[i]],]

        # Fit the model to the training set
        model <- rpart(status_group ~ ., method="class",data=train_data,control = list(maxdepth = depth))
        
        if (i == 1 & verbose) {
            plot(model, uniform=TRUE, main="Classification Tree - Rpart")
            text(model, use.n=TRUE, all=TRUE, cex=.8)
        }

        # Predict the label of test_data
        Y_pred <- predict(model,test_data[,-target_variable])
        Y_hat <- ifelse(Y_pred[,1] >= Y_pred[,2] & Y_pred[,1] >= Y_pred[,3],
                     "functional",
                     ifelse(Y_pred[,2] >= Y_pred[,3], "functional needs repair", "non functional")
                    )
        # True label
        Y <- test_data[,target_variable]
        
        # Compute confusion matrix 
        confusion_matrix <- table(Y_hat,Y)

        # Compute the classification rate
        if ((dim(confusion_matrix)[1]) == 1){
            classification_rate[i] = (confusion_matrix[1,1])/length(Y)
        } else if((dim(confusion_matrix)[1]) == 2) {
            classification_rate[i] = (confusion_matrix[1,1]+confusion_matrix[2,2])/length(Y)
        } else {
            classification_rate[i] = (confusion_matrix[1,1]+confusion_matrix[2,2]+confusion_matrix[3,3])/length(Y)
        }
    }
    
    end_time <- Sys.time()
    print(end_time - start_time)
    return(mean(classification_rate))
}
decision_tree_optimize_depth <- function(trainingSet,k,values_depth,verbose) {
    classification_rate <- c()
    # Computing a 10 fold cross validation for each depth
    for (depth in values_depth) {
        classification_rate <- c(classification_rate,decision_tree_fold_crossValidation(trainingSet,k,depth,verbose))
    }
    cat(paste("\nClassification rates :\n"))
    print(classification_rate)
    plot(classification_rate)
    bestIndex <- argmax(mat(classification_rate,nrow=1))
    bestClassificationRate <- max(classification_rate)
    cat(paste("The decision tree works best with a depth equal to",
              values_depth[bestIndex],
              ", obtaining a classification rate of",bestClassificationRate))
    return(bestClassificationRate)
}
if (EXECUTION == 2) {
    accuracy_dt <- decision_tree_fold_crossValidation(train,10,9,verbose=TRUE)
    cat(paste("Classification rate :",accuracy_dt))
}

if (EXECUTION == 3) {
    best_accuracy_dt <- decision_tree_optimize_depth(train,10,c(3,5,7,9,11,13,15,17,19,21),verbose=FALSE)
    best_accuracy_dt
}<jupyter_output>Time difference of 45.58932 secs
Classification rate : 0.524360269360269<jupyter_text>Random forest>The random forest builds several decision trees, using for each one a random number of features. The output of the model is obtained by passing the input inside each tree. The output that is the more frequent is considered the final output.>For this model, because of the big computation times, we don't dive into hyperparameters tuning, leaving the default values.<jupyter_code>random_forest_fold_cross_validation <- function(trainingSet,k,verbose) {
    
    # INITIALISATION
    
    start_time <- Sys.time()
    
    # removing column id in dataset
    trainingSet <- trainingSet[, -which(names(trainingSet) == "id")]
    # indexing status_group
    target_variable <- ncol(trainingSet)
    # vector storing the classification rates
    classification_rate <- array(0,k)
    # Shuffle trainingSet
    train_idx <- sample(1:nrow(trainingSet))
    # Split the dataset into k groups
    max <- ceiling(nrow(trainingSet)/k)
    splits <- split(train_idx, ceiling(seq_along(train_idx)/max))

    # 10 FOLD CROSS VALIDATION
    
    for (i in 1:k){
        # separate the training set from the fold test set
        test_data <- trainingSet[splits[[i]],]
        train_data <- trainingSet[-splits[[i]],]

        # Fit the model to the training set
        model <- randomForest(status_group ~ ., method="class",data=train_data)
        
        # get the importance of the features in predicting the target value
        importance = varImp(model)
        
        # print the importances (only for the first iteration)
        if (verbose & i == 1) {
            perm <- order(importance$Overall,decreasing=T)
            column_name = importance[perm,0]
            column_number = importance[perm,1]
            ordered_importance <- data.frame(column_name = column_name, column_number=column_number)
        }
        
        # Predict the label of test_data
        Y_pred <- predict(model,test_data[,-target_variable])
        
        # True label
        Y <- test_data[,target_variable]
        
        # Compute confusion matrix
        confusion_matrix <- table(Y_pred,Y)

        # Compute the classification rate
        if ((dim(confusion_matrix)[1]) == 1){
            classification_rate[i] = (confusion_matrix[1,1])/length(Y)
        } else if((dim(confusion_matrix)[1]) == 2) {
            classification_rate[i] = (confusion_matrix[1,1]+confusion_matrix[2,2])/length(Y)
        } else {
            classification_rate[i] = (confusion_matrix[1,1]+confusion_matrix[2,2]+confusion_matrix[3,3])/length(Y)
        }
    }
    
    end_time <- Sys.time()
    print(end_time - start_time)
    return(list(accuracy=mean(classification_rate),importance=ordered_importance))
}
if (EXECUTION >= 2) {
    results_rf <- random_forest_fold_cross_validation(train,5,verbose=TRUE)
    cat(paste("\nEstimated classification rate of the model\n"))
    results_rf$accuracy
}
if (EXECUTION >= 2) {
    cat(paste("Importance of the features"))
    results_rf$importance
}<jupyter_output>Importance of the features<jupyter_text>Discussion>The tree models accessed in part 2 obtain very different results.>First, the decision tree shows pretty bad results, predicting only 53% of correct target values. This behavior can be explain by the huge dimensionality reduction that has been done in part 1. The tree has not so many variables to work with and struggles in creating a powerful predictor.>Second, the single-hidden-layer neural network works better than the decision tree, having a estimated classification rate up to 60%. The simple structure of the neural network is probably the reason that more precise rules haven't been discovered.>At last, the random forest beats all records with a estimated classification rate equal to 79%. 20% above the neural network.>If we needed to select a model from those three, it would be the random forest.  3) Learning procedure with Gradient Boosting (XGBOOST)> In this chapter we will go back to another model validation procedure, this time acessing the classification rate of Gradient Boosting algorithm. We use the model implementation from the library **xgboost**.  One-hot encoding>Because the model xgboost only work on numeric features, we choose to encode the categorical features with one-hot encoding. Only the target value "status_group" will be encoded with integer encoding because xgboost doesn't allow multiple target features. We thus have 1 for functional, 2 for functional needs repair, and 3 for non functional.  Model validation : 10 fold cross validation>One more time, we choose to estimate our model classification rate by using a 10 fold cross validation method.<jupyter_code>xg_boost_fold_cross_validation <- function(trainingSet,k) {
    
    # INITIALISATION
    
    start_time <- Sys.time()
    
    # removing column id in dataset
    trainingSet <- trainingSet[, -which(names(trainingSet) == "id")]
    
    # encode factors into dummy variables (one-hot encoding)
    options(warn=-1)
    trainingSet <- cbind(trainingSet,dummy(trainingSet$region))
    trainingSet <- cbind(trainingSet,dummy(trainingSet$extraction_type_group))
    trainingSet <- cbind(trainingSet,dummy(trainingSet$quantity))
    trainingSet <- cbind(trainingSet,dummy(trainingSet$waterpoint_type_group))
    trainingSet <- trainingSet[,-c(4:8)]
    options(warn=0)
    
    # indexing status_group
    target_variable <- ncol(trainingSet)
    # vector storing the classification rates
    classification_rate <- array(0,k)
    # Shuffle trainingSet
    train_idx <- sample(1:nrow(trainingSet))
    # Split the dataset into k groups
    max <- ceiling(nrow(trainingSet)/k)
    splits <- split(train_idx, ceiling(seq_along(train_idx)/max))
    
    # 10 FOLD CROSS VALIDATION
    
    for (i in 1:k) {
        # separate the training set from the fold test set
        options(warn=-1)
        test_data <- trainingSet[splits[[i]],]
        label_test <- as.numeric(test_data$status_group)
        test_data <- test_data[, -which(names(test_data) == "status_group")]

        train_data <- trainingSet[-splits[[i]],]
        label_train <- as.numeric(train_data$status_group)
        train_data <- train_data[, -which(names(train_data) == "status_group")]
        options(warn=0)
        
        # make sure the data set is numeric
        test_data <- as.matrix(as.data.frame(lapply(test_data, as.numeric)))
        train_data <- as.matrix(as.data.frame(lapply(train_data, as.numeric)))
        
        # transform the dataset into a xgb.DMatrix for reasons of efficiency
        train.DMatrix <- xgb.DMatrix(data = train_data,label = label_train, missing = NA)
        
        # Fit the model to the training set
        model <- xgboost(data = train.DMatrix, objective = "multi:softmax", booster = "gbtree",
                         eval_metric="mlogloss", nrounds=10,num_class = 4,eta = .7,
                         max_depth = 14,colsample_bytree =1.,verbose=F)
        
        # print the importances of the features in predicting the target value (only for the first iteration)
        if (i == 1) {
            importance <- xgb.importance(feature_names = colnames(train_data), model = model)
            xgb.plot.importance(importance_matrix = importance)
        }
        
        # Predict the label of test_data
        Y_pred <- predict(model,test_data)
        Y_pred[Y_pred==1] <- "functional"
        Y_pred[Y_pred==2] <- "functional needs repair"
        Y_pred[Y_pred==3] <- "non functional"
        
        # True label
        label_test[label_test==1] <- "functional"
        label_test[label_test==2] <- "functional needs repair"
        label_test[label_test==3] <- "non functional"
        
        # Compute confusion matrix
        confusion_matrix <- table(Y_pred,label_test)
        
        # Compute the classification rate
        classification_rate[i] = (confusion_matrix[1,1]+
                                  confusion_matrix[2,2]+confusion_matrix[3,3])/length(label_test)
    }

    end_time <- Sys.time()
    print(end_time - start_time)
    
    return(mean(classification_rate))
}
if (EXECUTION >= 2) {
    accuracy_xgb <- xg_boost_fold_cross_validation(train,10)
    cat(paste("\nEstimated classification rate of the model\n"))
    accuracy_xgb
}<jupyter_output>Time difference of 41.9445 secs

Estimated classification rate of the model<jupyter_text>> We are satisfied with this accuracy. We can thus launch the learning process.  Learning phase>Now we can run the model on the entire training set and apply this model to the test set to predict the missing labels.<jupyter_code>learning_xg_boost <- function(dataSet,lenTrain,train_labels) {
    start_time <- Sys.time()
    
    # store the id column of the dataset
    id_column <- dataSet$id
    # remove the id column
    dataSet <- dataSet[, -which(names(dataSet) == "id")]
    
    # encode the dataset using one-hot encoding
    options(warn=-1)
    dataSet <- cbind(dataSet,dummy(dataSet$region))
    dataSet <- cbind(dataSet,dummy(dataSet$extraction_type_group))
    dataSet <- cbind(dataSet,dummy(dataSet$quantity))
    dataSet <- cbind(dataSet,dummy(dataSet$waterpoint_type_group))
    dataSet <- dataSet[,-c(4:8)]
    options(warn=0)
    
    # split the dataset into training and test set (actual ones this time)
    train_set <- dataSet[1:lenTrain,]
    test_set  <- dataSet[(lenTrain+1):nrow(dataSet),]
    
    # encode target value using integer encoding
    label <- as.numeric(train_labels$status_group)
    
    # making sure the datasets are numeric
    train_set <- as.matrix(as.data.frame(lapply(train_set, as.numeric)))
    test_set <- as.matrix(as.data.frame(lapply(test_set, as.numeric)))
    
    # transforming training data into xgb.DMatrix
    train.DMatrix <- xgb.DMatrix(data = train_set,label = label, missing = NA)
    
    # learning time for the xgboost model
    model <- xgboost(data = train.DMatrix, objective = "multi:softmax", booster = "gbtree",eval_metric="mlogloss",
                     nrounds=10,num_class = 4,eta = .7,max_depth = 14,colsample_bytree =1.,verbose=F)

    # predict the labels of the testset
    Y_pred <- predict(model,test_set)
    Y_pred[Y_pred==1] <- "functional"
    Y_pred[Y_pred==2] <- "functional needs repair"
    Y_pred[Y_pred==3] <- "non functional"
    
    # associate each prediction with its id
    result <- data.frame(id=id_column[(lenTrain+1):nrow(dataSet)],status_group=Y_pred)
    
    cat(paste("Time spent : ",Sys.time() - start_time))

    return(result)
}
result <- learning_xg_boost(dataSet,nrow(train),trainLabels)
write.csv(result,"../submissions/submisssion_9_vars_XGBOOST_whole_train.csv",quote=FALSE,row.names=FALSE)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09751201732607892, 'fraction_numerical': 0.01259838360361312, 'mean_word_length': 4.00899589892843, 'pattern_counts': {'":': 0, '<': 300, '<?xml version=': 0, '>': 215, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '313096', 'license': 'MIT', 'n_tokens_mistral': 11472, 'n_tokens_neox': 10683, 'n_words': 3404, 'path': 'src/notebook.ipynb', 'repo_name': 'jnicolasthouvenin/Predicting_functional_state_of_water_pumps', 'size': '37978'}