Fighting Crime with Tidyverse

This notebook illustrates the use of Tidyverse libraries (dplyr, tidyr, etc. ) for feature engineering. The data is used for this notebook the San Francisco Crime Classification competition on Kaggle. Let us first load all the libraries that we will use

# Load libraries
library(plyr)
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(purrr)
library(lubridate)
library(Matrix)
library(MatrixModels)
library(caret)

The libraries dplyr, tidyr, readr, stringr, and purrr are parts of the tidyverse collection. We load the data using the read_csv function of readr which is more efficient than the corressponding read.csv of base R.

train <- read_csv("./data/train.csv", col_types = "Tccccccnn")
test <- read_csv("./data/test.csv", col_types = "iTcccnn")

The variables test and train are tibble objects, a more efficient version of the base R data.frame. Let’s print train and test to explore what they look like

train

and

test

The print method of the tibble object never attempts to print the whole data.frame as opposed to base R, which is helpful especially in data exploration. The Category column is the target category which we need to predict for the test data. As can be seen, the test data does not contain the columns Resolution and Descript, so we remove them from train as well. It may be possible to use the information in these columns to engineer some features, but for now we ignore them for simplicity.

train$Descript <- NULL
train$Resolution <- NULL

Now, we are ready for some feature engineering!

Feature Engineeing

The data contains several useful features that are already suitable for using in a model. PdDistrict is a categorical variable that encodes which police district the recorded crime is associated with, and X and Y are the longitude and latitude of the location. One obvious attempt to generate new features would be to use the Dates column and create year, month, hour (day is already encoded in DayOfWeek). Not as straightforward is the engineering of features from the Address column. Let’s first look at how many distinct addresses there are in the dataset

length(unique(train$Address))

[1] 23228

length(unique(test$Address))

[1] 23184

In principle, one can construct a categorical Address variable with around 20,000 levels (or better perform feature hashing, see here for an example). This would lead to very long training times for models that we want to use, so we rather resort to engineering simpler features from Address. One idea is to construct numeric features from Address using the ratio of each crime in a given category to the total crimes (of all categories) recorded in a given Address. Namely, we use the following as a numeric feature

\[r_{A,c} = 1 + \frac{n_{A,c}}{\sum_{i=1}^{39} n_{A,i}}\] where \(n_{A,c}\) is the number of crimes of category \(c\) in a given address \(A\). There are 39 different crime categories, which explains the limits of the sum in the denominator, while the 1 is added to the ratio since we will compute the log of this feature eventually.

levels(as.factor(train$Category))

 [1] "ARSON"                       "ASSAULT"                     "BAD CHECKS"                 
 [4] "BRIBERY"                     "BURGLARY"                    "DISORDERLY CONDUCT"         
 [7] "DRIVING UNDER THE INFLUENCE" "DRUG/NARCOTIC"               "DRUNKENNESS"                
[10] "EMBEZZLEMENT"                "EXTORTION"                   "FAMILY OFFENSES"            
[13] "FORGERY/COUNTERFEITING"      "FRAUD"                       "GAMBLING"                   
[16] "KIDNAPPING"                  "LARCENY/THEFT"               "LIQUOR LAWS"                
[19] "LOITERING"                   "MISSING PERSON"              "NON-CRIMINAL"               
[22] "OTHER OFFENSES"              "PORNOGRAPHY/OBSCENE MAT"     "PROSTITUTION"               
[25] "RECOVERED VEHICLE"           "ROBBERY"                     "RUNAWAY"                    
[28] "SECONDARY CODES"             "SEX OFFENSES FORCIBLE"       "SEX OFFENSES NON FORCIBLE"  
[31] "STOLEN PROPERTY"             "SUICIDE"                     "SUSPICIOUS OCC"             
[34] "TREA"                        "TRESPASS"                    "VANDALISM"                  
[37] "VEHICLE THEFT"               "WARRANTS"                    "WEAPON LAWS"                

This idea is easy to implement (especially with tidyverse functionalities, as will be shown below), however requires some care. In a naive implementation, one would compute the crime by address ratios (defined above) for each Address in the train set. Then, these ratios can be merged with test by the Address column (i.e. by using left_join function of dplyr with by="Address"). Doing so, one will immediately realize (after obtaining predictions for Category) that this leads to overfitting. The reason is that we have used the target variable Category to construct the new crime by address ratio features. As a result, the trained model found much higher weights for these features, since they are highly correlated to Category by construction. Thus, the model memorizes the train data and cannot generalize properly when it encounters new data.

This looks pretty bad and it seems like constructing new features using the target Category is doomed. However, if one splits the train data into 2 pieces, and construct crime by address ratios from piece_1 and merge them with piece_2 (and repeat vice versa from piece_2 to piece_1) then the overfitting could be mitigated. The reason that this works is because the new features are constructed by using out-of-sample target values (Category) and so the crime by address ratios of each piece is not memorized.

As an illustration consider the following piece of example data:

Let’s split this data such that piece 1 contains rows (1,2,3,4,5,8) and piece 2 contains (6,7). Now, construct the crime by Address ratios (using the above formula) using piece 1, which would result in

and merge these with piece 2, which results in

The NAs are the result of the fact that Address E was not a part of piece 1. In a large dataset, such NA values would reflect the fact that in that address there is not much crime so one can impute them with the default value 1.0 (i.e. no crimes) as an approximation.

Now, it would be ideal to have more divisions of train so that each piece contains Address values that are present in the remaining pieces, so that a precise estimate of its crime by address feature can be constructed (thus fewer NAs would be encountered). This calls for the use of a k-folds divison of train. This is in fact a reminescent of the stacking predictions idea, which is used for combining predictions of different models, and to remove the bias on the best perfroming model.

Now it is time to construct some functions to implement this idea. The first one is the simple function for imputing NAs with a default value (1.0)

# Fill NA's with default ratio
fill_na <- function(x, default_ratio){
  x[is.na(x)] <- default_ratio
  x
}

The second one constructs k-folds division of train and computes the crime by address ratios. Notice the use of spread function from tidyr which is used to contsruct new columns from the rows. The cool aspect of Tidyverse is that all the data wrangling operations are extremely intuitive.

create_ratios_folds <- function(train, nfolds, seed = 1234){
  # Create folds
  set.seed(seed)
  folds <- createFolds(train$Category, k = nfolds)
  
  # Initiate a list where the holdout data frames from each fold will be stored
  list_out <- list()
  
  for (ifold in 1:length(folds)){
    train_in <- train[-folds[[ifold]], ]
    train_ho <- train[folds[[ifold]], ]
    
    # Group by Address & Category
    df <- train_in %>%
      mutate(Address = as.factor(Address)) %>%
      mutate(Category = as.factor(Category)) %>%
      group_by(Address, Category) %>% count() %>%
      rename(tot_address_category = n)
    
    # Group by Address
    df2 <- train_in %>%
      mutate(Address = as.factor(Address)) %>%
      mutate(Category = as.factor(Category)) %>%
      group_by(Address) %>% count() %>% 
      rename(tot_address = n)
    
    # Join and determine the ratio of each Category per Address
    df3 <- left_join(df, df2, by = "Address")
    df3 <- df3 %>%
      mutate(freq = ifelse(tot_address_category >= 5, 1.0 + tot_address_category / tot_address, 1.0))
    
    # Now spread the ratios to columns
    df3$tot_address_category <- NULL # Not needed
    df3$tot_address <- NULL          # Not needed
    df4 <- spread(df3, key = Category, value = freq, fill = 1.0) # Fill NA's with 1
    
    # Collapse the feature names (useful for future use in model Matrices)
    names(df4) <- str_replace_all(names(df4), "( )", "")
    names(df4) <- str_replace(names(df4), "/", "_")
    names(df4) <- str_replace(names(df4), "-", "_")
    
    # Join df4 with train_ho
    train_ho <- left_join(train_ho,
                           df4 %>% ungroup(Address) %>% mutate(Address = as.character(Address)),
                           by = "Address")
    
    # Impute NAs with 1.0 using fill_na function
    train_ho <- train_ho %>% mutate_at(names(df4)[-1], fill_na, default_ratio=1.0)
  
    # Add to list
    list_out[[ifold]] <- train_ho
  }
  
  # Return train_out
  bind_rows(list_out)
}

Now, we can create the crime by address ratios from the train set by a 5-fold division

# Create new fetures using 5 folds
train_new <- create_ratios_folds(train, 5)

OK, it is time to merge with the test set so that it has these newly engineered features. Below is the function that we use to perform this step

create_ratios_test <- function(train, test){
  # Features to be created
  log_feats <- names(train)[8:46]
  
  # At each address, the median of the determined Address features
  df <- train %>% select_(.dots = c(log_feats, "Address")) %>% distinct() %>%
    group_by(Address) %>% summarise_at(log_feats, median)
  
  # Join df4 with train_out
  test_out <- left_join(test,
                        df %>% ungroup(Address),
                        by = "Address")
  
  # Impute NAs with 1.0
  test_out <- test_out %>% mutate_at(log_feats, fill_na, default_ratio=1.0)
  
  # Return test_out
  test_out
}

Notice that we do not perform a simple merge by left_join. First, we create a table df which contains the crime by address ratios and the corressponding addresses. This table has duplicates (more than one crime happened in one address in the whole data) so we get rid of them by the distinct function. However, since we have used out-of-sample values for crime ratios, it happens that for a given address, there are multiple possible values of each crime ratio. A simple solution to obtain unique values is simply to use the median value in such cases, which is implemented by the summarize_at function above. The rest is simply merging test with df and returning the new data test_out that contains the newly engineered features!

# Cerate new features in test
test_new <- create_ratios_test(train_new, test)

           used  (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells  1617791  86.4    2637877  140.9   2637877  140.9
Vcells 82123966 626.6  163649973 1248.6 135650070 1035.0

Finally, we compute the log of the new crime by address ratios (which leads to less skewed features), and the mutate_at function is suitable for this task

# Compute log
log_feats <- names(train_new)[8:46]
train_new <- train_new %>% mutate_at(log_feats, log)
test_new <- test_new %>% mutate_at(log_feats, log)

Further Feature Engineering

Now that we have performed the most intricate part of feature engineering, let’s contsruct a few more based on Dates and the text features of the Address column. It is always a good idea to merge train and test first to perform these wrangling operations

# Combine test and train
test_ids <- test_new$Id # Save id's
test_new$Id <- NULL
test_new$Category <- NA # Fill test Category with NAs
full_data <- bind_rows(train_new, test_new) # Filter test by is.na(Category) == TRUE
rm(train_new, test_new) # train_new and test_new not needed anymore

Let’s construct hour, month and year from Dates using lubridate functions

# Date and time
full_data <- full_data %>%
  mutate(year = year(ymd_hms(Dates))) %>%
  mutate(month = month(ymd_hms(Dates))) %>%
  mutate(hour = hour(ymd_hms(Dates)))
full_data$Dates <- NULL

Now, let’s squeeze some more features from Address, based on the text it contains. We will look for the following features:

1. Has the crime occured in an intersection (look for “/”)?
2. Has the crime ocured in an address that contains the word “Block”?
3. What type of abbreviation is used in the address (e.g. ST, AV, BLD,..)?
4. Is the Block numbered 0 in the address (instead of a valid numeric value like 200)?
5. Is the street numbered (e.g. 23RD, 15TH)?

# Address feature: is it an intersection?
find_intersection <- function(x){ as.integer(str_detect(x, "/")) }
full_data <- full_data %>%
  mutate(is.intersection = find_intersection(Address))
# Address feature: is it block?
find_block <- function(x) { as.integer(str_detect(x, "Block")) }
full_data <- full_data %>% 
  mutate(is.block = find_block(Address))
# Address features: is it ST, AV, DR, WY, PZ, LN, BL, RD, zeroBlock, numberedST
find_type <- function(x, t){
  as.integer(str_detect(x, t))
}
full_data <- full_data %>% 
  mutate(is.AV = find_type(Address, "( AV)$|( AVE)$|( AV )| ( AVE )")) %>%
  mutate(is.ST = find_type(Address, "( ST)$|( ST )")) %>%
  mutate(is.DR = find_type(Address, "( DR)$| ( DR )")) %>%
  mutate(is.WY = find_type(Address, "( WY)$|( WAY)$|( WY )|( WAY )")) %>%
  mutate(is.PZ = find_type(Address, "( PZ)$|( PZ )")) %>%
  mutate(is.PL = find_type(Address, "( PL)$|( PL )")) %>%
  mutate(is.LN = find_type(Address, "( LN)$|( LN )")) %>%
  mutate(is.BL = find_type(Address, "( BL)$|( BLVD)$|( BL )|( BLVD )")) %>%
  mutate(is.RD = find_type(Address, "( RD)$|( RD )")) %>%
  mutate(is.zeroBlock = find_type(Address, "^(0 Block)")) %>%
  mutate(is.numberedST = find_type(Address, "[0-9](RD)|[0-9](TH)"))
# Remove Adress feature
full_data$Address <- NULL

The stringr function str_detect has been very useful to capture regular expressions for this step. Let’s convert year, month, DayOfWeek, Category and PdDistrict into categorical variables using mutate

# Convert Date features, PdDistrict into factor
full_data <- full_data %>% 
  mutate(Category = as.factor(Category)) %>%
  mutate(PdDistrict = as.factor(PdDistrict)) %>%
  mutate(year = as.factor(year)) %>%
  mutate(month = as.factor(month)) %>%
  mutate(DayOfWeek = as.factor(DayOfWeek))

We have chosen to keep hour numeric since there seems to be a continuous relationship between hour and the type of crime. One thing that we noticed is that there are some values of latitude and longitude which do not make sense, since thay are out of San Francisco city limits

summary(full_data$X)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -122.5  -122.4  -122.4  -122.4  -122.4  -120.5 

summary(full_data$Y)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  37.71   37.75   37.78   37.77   37.78   90.00 

Definitely, one cannot have a crime located at North Pole to be recorded in San Fransisco! Let’s remove these strange outliers by calculating the median values of X and Y by PdDistrict and then replace outliers with the median corressponding to their PdDistrict

# Medians by PdDistrict
med_table <- full_data %>% group_by(PdDistrict) %>% summarise(med_X = median(X), med_Y = median(Y))
full_data <- full_data %>%
  left_join(med_table, by = "PdDistrict")
# Replace outliers
outs <- full_data$Y > 40 # Outliers, this removes outliers in Y too 
full_data[outs, ]$X <- full_data[outs, ]$med_X
full_data[outs, ]$Y <- full_data[outs, ]$med_Y
full_data$med_X <- NULL
full_data$med_Y <- NULL

Finally, let’s standardize the numeric columns X, Y and hour

# Normalize X,Y,hour
normalize_feat <- function(x){
  (x - mean(x))/sd(x)
}
full_data <- full_data %>%
  mutate(X = normalize_feat(X)) %>%
  mutate(Y = normalize_feat(Y)) %>%
  mutate(hour = normalize_feat(hour))

Now our data is ready for modelling!

Modelling with XGBoost

By popular demand, the XGBoost implementation of gradient boosted tree algorithm is our choice here. First we split back into train and test and remove the variable full_data to clean up some space in the memory

# Split back into train and test
train <- full_data %>% filter(is.na(Category) == FALSE)
test <- full_data %>% filter(is.na(Category) == TRUE)
rm(full_data); gc()

           used  (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells  1602234  85.6    5103933  272.6   7717408  412.2
Vcells 96741344 738.1  236034262 1800.8 235989328 1800.5

Then, we construct a sparse model matrix to encode the categorical variables we have created above

Y <- train$Category
X <- model.Matrix(~., 
                  data = train %>% select(-Category), 
                  sparse = TRUE)[,-1] 

Let’s prepare the data for XGBoost

# XGBoost
library(xgboost)

# XGboost wants levels to start with 0
class.names <- levels(Y)
levels(Y) <- 0:38
Y <- as.numeric(levels(Y))[Y]

and split the original train into training and validation sets:

# Split into train and validation
library(caret)
set.seed(123456789)
idx_train <- createDataPartition(train$Category, p = 0.7)[[1]]
X_tr <- X[idx_train, ]; X_vld <- X[-idx_train, ]
Y_tr <- Y[idx_train]; Y_vld <- Y[-idx_train]
# XGB style matrices
dtrain <- xgb.DMatrix(data = X_tr, label = Y_tr)
dvalid <- xgb.DMatrix(data = X_vld, label = Y_vld)
watchlist <- list(train=dtrain, test=dvalid)

Now train (this takes a while)

### Train using standard params
# Parameters
params <- list(booster = "gbtree",
               eval_metric = "mlogloss",
               objective = "multi:softprob",
               subsample = 1,
               colsample_bytree = 1,
               eta = 0.3,  
               min_child_weight = 1,
               max_depth = 6,
               gamma = 0.0,
               seed = 123)
# Fit 
fit_tr <- xgb.train(params = params,
                    data = dtrain,
                    num_class = 39,
                    nrounds = 500,
                    watchlist = watchlist,
                    early_stopping_rounds = 10,
                    print_every_n = 10,
                    maximize = FALSE)

[1] train-mlogloss:2.988178 test-mlogloss:2.994273 
Multiple eval metrics are present. Will use test_mlogloss for early stopping.
Will train until test_mlogloss hasn't improved in 10 rounds.

[11]    train-mlogloss:2.371029 test-mlogloss:2.404007 
[21]    train-mlogloss:2.286314 test-mlogloss:2.342143 
[31]    train-mlogloss:2.254193 test-mlogloss:2.329931 
[41]    train-mlogloss:2.234521 test-mlogloss:2.326505 
[51]    train-mlogloss:2.218702 test-mlogloss:2.324803 
[61]    train-mlogloss:2.204673 test-mlogloss:2.323799 
[71]    train-mlogloss:2.192091 test-mlogloss:2.323144 
[81]    train-mlogloss:2.179863 test-mlogloss:2.322827 
[91]    train-mlogloss:2.167418 test-mlogloss:2.322321 
[101]   train-mlogloss:2.156822 test-mlogloss:2.322483 
Stopping. Best iteration:
[93]    train-mlogloss:2.165367 test-mlogloss:2.322304

A few points deserve more explanation here. The loss function that is being minimized here is the multi-class log-loss which is defined as \[ L = -\frac{1}{N} \sum_{i=1}^N\, \sum_{c=1}^{K}\, Y_{ic}\, \log(P_{ic}) \] where \(N\) is the number of training examples and \(K\) is the number of classes (crime categories here). \(Y_{ic}\) is the one-hot encoded matrix of true classes which assings the value 1 to \(Y_{ic}\) if the example \(i\) is of class \(c\) and 0 otherwise. \(P_{ic}\) is the matrix of predicted class probabilities from the model. This loss function awards high predicted probabilities for correct class and low predicted probabilities for others, so that it contains more information than a simple accuracy measure.

This post is not about XGBoost and how to tune its hyperparameters. However, a few observations from the above training log deserves some explanation. In this dataset, it is hard to perform a cross-validation sweep over the hyperparameters on a grid, since training is very slow. Thus, one would need to watch train/test loss by tweaking the hyperparameters by hand. As can be seen, early stopping kicked in just after about 100 rounds. This means that we started to overfit rather quickly. In this situation, it would be useful to reduce eta and increase min_child_weight. In addition, reducing subsample will help, but we have to make sure that there are enough nrounds to let the training loss to become low enough. The next step would be to try to increase max_depth to capture more complex patterns in each round, while lowering eta and increasing min_child_weight to prevent overfitting. It could also be useful to reduce colsample_bytree which would help reduce overfitting, since at each step a different subset of all features will be used to grow a tree. Again in such a case, one has to make sure nrounds is large enough. I would not expect gamma having much of a use in this case, since we do not see the test loss skyrocketing at any point.

Now let’s look at the importance of each feature (measured by the amount they reduce the loss function, see here for more details)

importance <- xgb.importance(feature_names = X_tr@Dimnames[[2]], model = fit_tr)
importance

and plot the first 25 most important features

xgb.plot.importance(importance_matrix = head(importance, 25))

As can be seen, the crime by address features are on the top of the importance list! The other features that look important are hour, is.intersection, X and Y, which are expected by intuition. The hour and location at which the crime happens correlate with which crime it is.

Final words

In this post we have shown how to wrangle columns in a dataset to engineer new features. We have also described a method to avoid pitfalls of overfitting when the target variable in the training data is incorporated in feature engineering. All of the data wrangling is made quite simple and intuitive by the Tidyverse collection of libraries. I strongly recommend the book R for data science to learn more about the these functionalities and further examples of usage.