<jupyter_start><jupyter_text>Vantage-AI case projectFor this project we'll be using both R and python.R is used for more exploratory analyses and less automated analysis work. Once we understand the data and havea clear idea of which rules to implement for data cleaning, then we can continue working with Python to set up pipelines with the `scikit-learn` library. The goal of this case is to predict the need of repair of water pumps, effectively a fairly standard classification problem.What makes it hard is that there are many features that consists of an intractable amount of labels, meaning we can't just use them as categories since many categories will only have 1 or 2 observations; there won't be enough observations per 'category' to train a model Reticulate packageThe `reticulate` package is used to run python code. This way we can do some basic computation steps by writing python code within our R script, making it easier to port the resulting code to dedicated Python scripts to automate the analysis from raw data to the final result.After a LOT of trial and error the reticulate package now works pretty cleanly with `conda`. You can use `conda` as a package manager for both R and Python, as opposed to R package manager `packrat` or the widely used `virtualenv` for python. Project directory structureThe directory structure was taken from the BigData Republic GitHub! By doing this I hope to learn more about the preferred project structure at BDR.The cookiecutter-data-science repository can be found here https://github.com/BigDataRepublic/cookiecutter-data-science. Now, what makes this truly come full circle is that...    - The Tanzania Water Pump data was originally hosted by **DrivenData** for a data science competition       * Hosted here https://www.drivendata.org/competitions/7/pump-it-up-data-mining-the-water-table/       * Original data made available by [Taarifa](http://taarifa.org) and [Tanzania Ministry of Water](http://maji.go.tz/)   - The Github repository from BDR is actually a fork from https://github.com/drivendata/cookiecutter-data-science !       * Changes include getting rid of `virtualenv` , thus only using `conda`       * Source code directory named after the project instead of just `src`       While the directory structure is completely overkill for a project like this, it still serves its purpose as a standardized data science project directory. Even if project are small, finding the code isn't usually straightforward... Unless you already know the file structure, then it becomes a piece of cake!          Exploratory analysisBefore we move to modelling and actual data cleaning, let's explore the data to get a feeling for the data.<jupyter_code># Load libraries
library(tidyverse)
library(broom)
library(ggthemes)
library(DT)
# Library to use python code using our conda environment
library(reticulate)
py_config()<jupyter_output><empty_output><jupyter_text>Predicting Tanzania water pump maintenance The dataThe case presented here is about water pumps in Tanzania. The data set originates from the Tanzanian Ministry of Water. Currently, this ministry maintains its pumps based on a maintenance schedule or, of course, when they break down. We feel that the maintenance of the Tanzanian water pumps could improve in both the cost of maintenance and the prevention of break downs by introducing machine learning to predict if a water pump is in need of repair or even the moment of failure of each water pump.  ObjectiveOur objective is to develop a reproducible model that can predict which pumps will fail in the future, either the moment or just whether they fail. The resulting model is to be used by the Tanzanian government to effectively maintain their water pumps.Normally for this type of data you would expect that a single water pump was measuredover-time and using this data we could predict the need for repairs for a given pump. In this case however, at every practically every timepoint a different water point was measured! So we have to make some assumptions about thedata and the approach to how we create a model. Assumptions**We should know in advance which water point should be predicted**, there must be something like a maintenance schedule at least 3 days ahead. We make this assumption because we have to able to takeinto account the measured variables about the water pumps.!**Measured variable should all be assumed to be known in advance**, otherwise we should chooseto leave out some variables that may have only become known at the time of determining the functioning of a water pump.Variables like the `funder` and `district_code` are clearly known ahead of maintenance, but other variables like `water_quality` and `amount_tsh` (total amount water available to waterpoint) may not be known before determining the state of the pump. Getting a feeling for the dataBefore we start modelling the data we should understand the data. In this datasetwe find a total of 40 measured variables to be used as features in our models. In thissection we will aim to understand these features and make motivated choices to keep, remove or transform thesefeatures from our prediction model. A closer look at our featuresOur data consists of,  - `r nrow(df)` observations   - `r ncol(df)` features  - `r length(unique(labels[,2]))` types of labels to predict  Not all of the features are necessarily useful for predicting the moment of failure of water pumps, for instance the measurement `id` won't be of much use at it is unique.First we'll take a look at the raw data to fully understand the structure of our data,especially how exactly the data was measured over time. We have `r nrow(df)` observations. The raw dataThe raw data is not considered big-data in the sense that it will fit on a single machine; we do not have to worry about batch processing because of sheer file size with 'only' `r nrow(df)` observations. The number of features is also a mere `r ncol(df)`! Enough to make a interpretable predictive model and this will allow us to dive deeper into theindividual variables.<jupyter_code># Read local data
#   Expect to be synced with data hosten on Amazon S3 
data <- read.csv("../../data/raw/water_pump_set.csv")
labels <- read.csv("../../data/raw/water_pump_labels.csv")
features <- read.csv("../../data/raw/water_pump_features.csv")
# Merge labels and features into single dataset by id
df <- merge(data, labels, by="id")
head(df)
# 39 features,to use for modelling a unique id and the status_group labels
# Total of 59400 observations
dim(df)
# Show the feature descriptions
#   This excludes the 'status'group'
head(features)
# Let's add some more info to the feature info using tidyverse piping
feature_desc <- features %>% filter(Feature %in% colnames(data)) %>%
  mutate(example = sapply(data[1, -1], as.character), # First row of the data without the 'id' column
         unique = map_int(data[,-1],function(x) length(unique(x))),
         type = map_chr(data[,-1], function(x) class(unlist(x))), # 
         nonzero = map_int(data[,-1], function(x) sum(as.character(x) != "0")), # Nonzero values, assuming 0 usually means NA
         max_single = map_int(data[,-1], function(x) max(table(x))),      # Observations of most occuring value
         max_occurence = map_chr(data[,-1], function(x) names(which.max(table(x))))) # Highest occuring value
                                 
# Let's glance at the numeric values in the dataset                               
feature_desc %>% filter(type %in% c("numeric", "integer"))                
feature_desc %>% filter(type %in% c("numeric", "integer")) %>% dim()<jupyter_output><empty_output><jupyter_text>The number of unique values for each numerical feature is surprising in some cases. Continuous values like, - `gps_height` - `population` - `amount_tsh`  all have a very low number of unique values. So, even though we have not seen any `NA` values in the data, but there must be some missing values. These values are simply set to `0`, it is unclear to me whether we can distinguish true `0` values from missing values.   Non-numeric summaries<jupyter_code># Show non-numeric value summary
feature_desc %>% filter(!type %in% c("numeric", "integer"))
feature_desc %>% filter(!type %in% c("numeric", "integer")) %>% dim()<jupyter_output><empty_output><jupyter_text>Handling feature groupsThe variables in the data are far from independent, some variables practically mean the same thing.Some feature groups are so painfully obvious that we should clean these manually before moving on tomore generic analyses.In the `Description` column, we can see some overlap in values. These features are likely to represent the samevariable, so we can reduce the dimensionality of the dataset just by smartly combining some variables.Let's start with common sense and group some of the variables manually to prepare the datafor more automated pipelines. Here are some logical groups we can check,  - scheme_name, scheme_management, installer, funder  - extraction_type, extraction_type_group, extraction_type_class  - management, management_group  - scheme_management, scheme_name  - payment, payment_type  - water_quality, quality_group  - quantity_quantity_group  - source, source_class  - waterpoint_type, waterpoint_type_group    Scheme_name, scheme_management, installer, funder<jupyter_code>data %>% 
    select(scheme_name, scheme_management, installer, funder) %>%
    mutate_all(as.character) %>%
    mutate_all(tolower) %>% # To lowercase since there are many typing 'errors' or inconsistencies
    group_by(scheme_name, scheme_management, installer, funder) %>%
    tally() %>% 
    arrange(desc(n)) %>% 
    head(10)
# Scheme names
data %>% select(scheme_name) %>% table() %>% sort(decreasing = T) %>% head()
# Half the data is missing data, and the rest is extremely scattered
# Furthermore there are many unambiguous names and spelling errors like 'Adra', 'ADRA', 'WOULD BANK', etc..
df %>%
    filter(scheme_name != "") %>%
    group_by(scheme_name, status_group) %>%
    tally() %>% head(20)
# Based on this we remove the scheme_name entirely
# Scheme management
data %>% select(scheme_management) %>% table() %>% sort(decreasing = T) %>% head()
# Plot label per scheme_management
df %>% select(scheme_management, status_group) %>%
    ggplot(aes(x = scheme_management, fill = status_group)) +
    geom_bar() + coord_flip() +
    ggtitle("Water pomp status voor management_scheme") +
    theme_fivethirtyeight()

# Most of the water pumps are from a handful of installers. 
# We could consider grouping all the low occurences into a group 'other'
# Plot label per scheme_management
scheme_management_groups = df %>% select(scheme_management, status_group) %>%
    group_by(scheme_management) %>%
    tally() 

cat(sprintf("Scheme mangement in 10 groups accounts for %1.2f%% of the data. And just VWC is %1.2f%%", 
            100 * sum(scheme_management_groups$n) / nrow(df),
            100 * max(scheme_management_groups$n) / nrow(df)))

# installer
data %>% select(installer) %>% table() %>% sort(decreasing = T) %>% head()
length(unique(data$installer))

# Most of the water pumps are from a handful of installers. 
# We could consider grouping all the low occurences into a group 'other'
# Plot label per scheme_management
installer_groups = df %>% select(installer, status_group) %>%
    group_by(installer) %>%
    tally() %>%
    top_n(20, n)

# Plot installer versus status_group to verify 
df %>% select(installer, status_group) %>%
    mutate(installer = ifelse(installer %in% installer_groups$installer, as.character(installer), "other"),
           installer = ifelse(installer %in% c("", "0"), "unknown", installer)) %>%
    ggplot(aes(x = installer, fill = status_group)) +
    geom_bar() + coord_flip() +
    ggtitle("Water pomp status voor installer", subtitle = "Kleine groepen samengevoegd in 'other'") +
    theme_fivethirtyeight() 

# Show counts of installer labes, only top20
# Notice that many of these have about the same number of observations and grouping all the others into
# 'other' is quite arbitrary since the 'other' group will still be half the data. 
# The number of 'larger' groups and the number of spelling errors in the data require more detailed feature engineering
# like some kind of string matching- which we will not spend our time on.
# THIS FEATURE WILL BE REMOVED FROM THE MODELLING DATA
# funder
data %>% select(funder) %>% table() %>% sort(decreasing = T) %>% head()

cat(sprintf("Of the %i unique funder names, %i overlap with the installer names!",
            length(unique(data$funder)),
            length(intersect(tolower(as.character(data$funder)), tolower(as.character(data$installer))))))
# Due to the high overlap between funder and installer, we can assume that these two are heavily related and 
# will probably give us the same difficulty as the funder names by themselves.
# WE REMOVE THIS FEATURE FROM THE DATASET!
# Groups for each extraction type to determine which has the most information
df %>% select(extraction_type, extraction_type_group, extraction_type_class) %>%
  apply(2, function(x) length(unique(x)))

# Count label co-occurences to potentially remove some features as we have a 
# tractable number of features (not 100s this time!)
df %>%
    group_by(extraction_type, extraction_type_group, extraction_type_class) %>%
    tally()<jupyter_output><empty_output><jupyter_text>By far the most prevalent type of waterpump is the gravity water pump type. This type of pump is eco-friendly as it does not use energy other than gravity.Notice that there are some very small groups in the data,  mostly due to the inclusion of `extraction_type`. The `extraction_type` has some values that we could combine with other categories to reduce the dimensionality.  - **india mark ii** and **india mark iii** can be combined to just **india mark**  - **cemo** and **climax** are just two types of motorpumps, let's not differentiate these two  - **other - mkulima/shinyanga** is very specific, we can just classify this as handpump together with **play pump** and **walimi**  - Lastly the **swn 81** and **swn 80** are both handpumps, we can classifiy these as **swn**      * The SWN 80 and SWN 81 type of handle pumps only differ in their pump head so that the SWN 81 can reach deeper wells. [Read more about handpumps here:)](https://www.ircwash.org/sites/default/files/232.2-13249.pdf)      * Fun fact, the SWN stands for 'Sociale Werkplaats Nunspeet'! The SWN family of pumps was introduced partly by a workshop for the handicapped in the Netherlands sine 1976.Let's manually combine some of the feature groups into related groups.<jupyter_code># Temporary data about the extraction type
df_extraction <- df %>% 
    mutate(
        extraction_type = as.character(extraction_type),
        extraction_type = ifelse(extraction_type %in% c("india mark ii", "india mark iii"), "india mark", extraction_type),
        extraction_type = ifelse(extraction_type %in% c("cemo", "climax"), "motorpump", extraction_type),
        extraction_type = ifelse(extraction_type %in% c("other - swn 81", "swn 80"), "swn", extraction_type),
        extraction_type = ifelse(extraction_type %in% c('other - mkulima/shinyanga','other - play pump', 'walimi'), "handpump", extraction_type)
           ) %>%
    select(id, extraction_type, extraction_type_class)

df_extraction %>% group_by(extraction_type, extraction_type_class) %>% tally()

# Make a new table of group counts
df_extraction %>% 
    mutate(extraction_group = paste(extraction_type, extraction_type_class, sep="-")) %>%
    group_by(extraction_group, extraction_type, extraction_type_class) %>%
    tally() %>%
    arrange(n) %>%
    ungroup() %>%
    mutate(extraction_type = factor(extraction_type, levels=extraction_type)) %>% # Relevel to reorder plot
    ggplot(aes(x = extraction_type, y = n, label = n, fill = extraction_type_class)) +
        geom_bar(stat="identity", position="dodge") + 
        ggtitle("Water pomp types na combineren van groepen") +
        geom_text() + coord_flip() + theme_fivethirtyeight()<jupyter_output><empty_output><jupyter_text>Using two out of three extraction types, after some adjustments, in the data we can classify all the data in a way that none of the groups is truly underrepresented.  management, management_groupThese two variables seem to mean the exact same thing.<jupyter_code># Cross table of only two variables
table(df %>% select(management, management_group))
# The management variable has more information than management_group<jupyter_output><empty_output><jupyter_text>We will select the column with *more* information, so we keep the `management` variable and drop the `management_group`. Payment and payment_typeThe payment and payment type columns also mean the exact same thing.<jupyter_code># Cross table of only two variables
df %>% group_by(payment, payment_type) %>% tally()
# The payments are exactly the same, we can drop either one.<jupyter_output><empty_output><jupyter_text>The payment_type variables are more concise so we will drop the `payment` column<jupyter_code># Cross table of only two variables
df %>% group_by(water_quality, quality_group) %>% tally()<jupyter_output><empty_output><jupyter_text>The water quality intuitively seems like it would be predictive of the water pump status.<jupyter_code># Is water quality correlated to the water pump status?
df %>% group_by(water_quality, status_group) %>%
    tally() %>%
    group_by(water_quality) %>% mutate(n = n / sum(n)) %>%
    filter(water_quality != "soft") %>%
    ggplot(aes(x = water_quality, y = n, fill = status_group)) +
    geom_bar(stat="identity") + coord_flip() + 
    ggtitle("Water kwaliteit versus pomp status ratios") +
    theme_fivethirtyeight()

# Is water amount available correlated to the water quality?
df %>% group_by(water_quality) %>%
    filter(amount_tsh > 100) %>%
    ggplot(aes(x = water_quality, y = amount_tsh, col = water_quality)) +
    geom_boxplot() + coord_flip() + 
    ggtitle("Water kwaliteit versus pomp status ratios") +
    theme_fivethirtyeight()<jupyter_output><empty_output><jupyter_text>The `water_quality` contains more information so we will drop the `quality_group`.The `water_quality` is also related to the `amount_tsh` in the data, even though we are removing `amount_tsh`(!) The only pumps where a significant amount of water is available is the ones where 'soft' water was found. So even though we are leaving out `amount_tsh`, we will still capture some of its information by including the `water_quality`.We could improve the visualisation above but it will likely not affect the choices we already made, so let's move on... Quantity and quantity_groupThe `quantity` and `quantity_group` sounds like they have the same relation as the water quality. Let's check!<jupyter_code># Cross table of only two variables
df %>% group_by(quantity, quantity_group) %>% tally()<jupyter_output><empty_output><jupyter_text>welp.. we checked, and we'll just leave out the `quantity_group` column.  Source, source_class and source_type<jupyter_code># Cross table of only two variables
df %>% group_by(source, source_class, source_type) %>% tally()<jupyter_output><empty_output><jupyter_text>The water `source` is more detailed than `source_class` and contains the exact same values as `source_type`. So let's keep only the `source` and `source_type`.However, the `source_class` contains only 3 categories and this more general distinction may in fact help us make predictions.We can calculate the feature importance later, so let's keep both variables in our data and verify whether inclusion of `source_class` improves the model at all.Alternatively,<jupyter_code># Is water quality correlated to the water pump status?
df %>% group_by(source, source_class, status_group) %>%
    tally() %>%
    group_by(source) %>% mutate(n = n / sum(n)) %>%
    ggplot(aes(x = source, y = n, fill = status_group)) +
    geom_bar(stat="identity") + coord_flip() + 
    ggtitle("Water bron versus pomp status ratios") +
    theme_fivethirtyeight()<jupyter_output><empty_output><jupyter_text>The source is definitely not a perfect prediction of water pump status, but at least a distinction between classes could be beneficial. It remains to be seen whether the water source is truly indicative of the water pump status, why would a pump near a lake be more likely to be non-functional than that near a spring? Unless the water source has dried up, but this would rather affect multiple pumps in the same area.If there was more time available we could find if there is just a cluster of water pumps at a single lake that are non-functional because of regional differences, for example. The location of a pump probably determines the water source, and since we will already take the gps-coordinates into account we might not be including any new information with the `source` feature. Waterpoint_type and waterpoint_type_groupFinally, this group of two variables may be reduced to a single variable.<jupyter_code># Cross table of only two variables
df %>% group_by(waterpoint_type, waterpoint_type_group) %>% tally()
# They are the exact same variable. So we remove the waterpoint_type_group from our data, later.<jupyter_output><empty_output><jupyter_text>Conclusion about dealing with grouped data variables**Remove 9 features**  - `scheme_name`, `scheme_management` ,    too many unique  - `extraction_type_group`,    overlaps with `extraction_type` and `extraction_class`  - `funder` and `installer`,   many missing values, unique values and data-entry inconsistencies  - `source_class`, `payment`, `waterpoint_type_group`, `management_group` because they are duplicates of other variables    **Edit 1 feature**  - `extraction_type`,  rename some similar categories with very low occurencesLet's make a note of these findings. We will cleanly integrate these steps into our data preprocessing pipeline.<jupyter_code>df_combined_groups <- df %>% 
    select(-c(scheme_name, scheme_management, extraction_type_group, funder,
              installer, source_class, payment, waterpoint_type_group, management_group)) %>%
    mutate(
        extraction_type = as.character(extraction_type),
        extraction_type = ifelse(extraction_type %in% c("india mark ii", "india mark iii"), "india mark", extraction_type),
        extraction_type = ifelse(extraction_type %in% c("cemo", "climax"), "motorpump", extraction_type),
        extraction_type = ifelse(extraction_type %in% c("other - swn 81", "swn 80"), "swn", extraction_type),
        extraction_type = ifelse(extraction_type %in% c('other - mkulima/shinyanga','other - play pump', 'walimi'), "handpump", extraction_type)
           )
dim(df_combined_groups)
# Write df_combined_groups to csv
#    NOTE THIS IS ACTUALLY HAPPENING IN DEDICATED SCRIPTS, NOT FROM THIS NOTEBOOK
#df %>% write_csv("../../data/interim/data_grouped_vars.csv")

# This part was rewritten into vantage-project/vantage/data/clean_data.py
# The exploration is done in R but final code execution in Python<jupyter_output><empty_output><jupyter_text>Missing valuesNow that we have removed some categorical values that we have checked because they seemed grouped with other variables, we still have to take a look at the remaining variables.While there are no true `NA` values in the data, we have already seen cases where missing values are indicated by empty values, zeroes, or labels like 'none'. Let us look at the feature descriptions again, and take note of the number of unique values and maximum occuring value.<jupyter_code>feature_desc %>% filter(Feature %in% colnames(df_combined_groups))<jupyter_output><empty_output><jupyter_text>amount_tshWe'll **remove** this column as it consists mostly of zeroes, and there is no way for us to distinguish whether these are missing values or just dried up water pumps. I expect that not every included water pump was actually measured as detailed as every other pump. By far the most values are zeroes.<jupyter_code>table(df_combined_groups$amount_tsh)[1:10]<jupyter_output><empty_output><jupyter_text>wpt_nameThe `wpt_name`, or waterpump name, many unique values and the most occuring values is 'none'.Now, even the pumps with the same exact name don't actually seem to guarantee the same pump.This is reason enough to **remove** this feature as unique labels like this will only aid to overfit the model.<jupyter_code># Show some values with wpt_name==Zahanati
# Notice that these are not the same waterpoint at all! 
# So wpt_name is not a variable we'll be able to use
df %>% filter(wpt_name == "Zahanati") %>% select(wpt_name, everything())<jupyter_output><empty_output><jupyter_text>num_privateThe `num_private`, which came without a detailed description, has a loooooottt of zero values.Again, we cannot differentiate between true 0 values and missing values, since 0 does seem like possible value for `num_private`. We remove this column from further analysis.<jupyter_code># Ratio of values in the data
table(df$num_private)[1:5] / nrow(df) * 100<jupyter_output><empty_output><jupyter_text>construction_yearMissing values of `construction_year` are indicated by zero-values. How to replace these is not necessarily straightforward since we treat this as a continuous variable! Since the difference between 1955 and 2000 is much smaller than 1955 and 0 we should be careful here! If the construction year is predictive of the water pump status then zero-values will likely introduce a bias.For missing values we should definitely impute some value. Traditionally the most straightforward way is to just impute the mean or median value. We may correct for region or pump type to improve this imputation. A quick glance at factors that may indicate a missing construction_year we don't find any. So let us use a very simple imputation method to save time.<jupyter_code>df %>% filter(construction_year > 0) %>%
    ggplot(aes(x = construction_year, color = construction_year)) + geom_bar() + 
    ggtitle("Distribution of non-missing construction years",
           sprintf("Number of excluded missing construction years: %i", sum(df$construction_year == 0))) +
    theme_fivethirtyeight()

median(df$construction_year[df$construction_year > 0])

# Replace construction_year with median of non-zero construction years
df_combined_groups %>%
    mutate(construction_year = replace(construction_year, 
                                       construction_year==0,
                                       median(construction_year[construction_year>0])))<jupyter_output><empty_output><jupyter_text>GPS coordinatesThere are three types of gps coordinates in the data,  - gps_height  - longitude  - latitudeWe expect these to be continuous values.<jupyter_code># See which gps coordinates occur more than once to reveal missing values
# Maybe some duplicate gps coordinates?? Guess not..
df %>% group_by(gps_height, longitude, latitude) %>%
  tally() %>% filter(n > 2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08092399955576944, 'fraction_numerical': 0.00433124791766927, 'mean_word_length': 4.823237766760077, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 185, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '62090', 'license': 'Apache-2.0', 'n_tokens_mistral': 7550, 'n_tokens_neox': 7098, 'n_words': 3234, 'path': 'notebooks/dev/data-exploration-r.ipynb', 'repo_name': 'Lodewic/vantage-project', 'size': '27013'}