This R Notebook is the complement to my blog post Predicting And Mapping Arrest Types in San Francisco with LightGBM, R, ggplot2.
This notebook is licensed under the MIT License. If you use the code or data visualization designs contained within this notebook, it would be greatly appreciated if proper attribution is given back to this notebook and/or myself. Thanks! :)
Setup the R packages.
library(lightgbm)
library(Matrix)
library(caret)
library(viridis)
library(ggmap)
library(randomcoloR)
source("Rstart.R")
sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.3
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils
[6] datasets methods base
other attached packages:
[1] stringr_1.1.0 digest_0.6.11 RColorBrewer_1.1-2
[4] scales_0.4.1 extrafont_0.17 dplyr_0.5.0
[7] readr_1.0.0 ggmap_2.7 viridis_0.3.4
[10] caret_6.0-73 ggplot2_2.2.1 lattice_0.20-34
[13] Matrix_1.2-7.1 lightgbm_0.1 R6_2.2.0
loaded via a namespace (and not attached):
[1] reshape2_1.4.2 splines_3.3.2 colorspace_1.3-2
[4] stats4_3.3.2 mgcv_1.8-16 ModelMetrics_1.1.0
[7] nloptr_1.0.4 DBI_0.5-1 sp_1.2-4
[10] jpeg_0.1-8 foreach_1.4.3 plyr_1.8.4
[13] MatrixModels_0.4-1 munsell_0.4.3 gtable_0.2.0
[16] RgoogleMaps_1.4.1 mapproj_1.2-4 codetools_0.2-15
[19] knitr_1.15.1 SparseM_1.74 quantreg_5.29
[22] pbkrtest_0.4-6 parallel_3.3.2 Rttf2pt1_1.3.4
[25] proto_1.0.0 Rcpp_0.12.8 geosphere_1.5-5
[28] lme4_1.1-12 gridExtra_2.2.1 rjson_0.2.15
[31] png_0.1-7 stringi_1.1.2 tools_3.3.2
[34] bitops_1.0-6 magrittr_1.5 maps_3.1.1
[37] lazyeval_0.2.0 tibble_1.2 car_2.1-4
[40] extrafontdb_1.0 MASS_7.3-45 data.table_1.10.0
[43] assertthat_0.1 minqa_1.2.4 iterators_1.0.8
[46] nnet_7.3-12 nlme_3.1-128
Import data, and only keep relevant columns. Filter on Arrests only.
The data must be randomized for lightgbm
to give unbiased scores. Can do with dplyr’s sample_frac
.
# seed for sample_frac()
set.seed(123)
df <- df %>% sample_frac()
df %>% head()
There are 634,299 arrests in this dataset.
Engineer features for lightgbm
.
Year is # years since the lowest year (in this case, 2003, as noted in the dataset title)
df <- df %>%
mutate(month = factor(substring(Date, 1, 2)),
hour = factor(substring(Time, 1, 2)),
year = as.numeric(substring(Date, 7, 10)))
There were 50 or more warnings (use warnings() to see the first 50)
df %>% select(month, hour, year) %>% head()
Change DayOfWeek to Factor.
Since column become encoded as numeric instead of categorical, encode order of numerals such that Saturday and Sunday are adjacent for proper lte
/gte
behavior.
dow_order <- c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")
df <- df %>%
mutate(DayOfWeek = factor(DayOfWeek, levels=dow_order))
df %>% select(DayOfWeek) %>% head()
Map the category to an index. Labels must be zero-indexed.
df <- df %>%
mutate(category_index = as.numeric(factor(Category)) - 1)
df %>% select(category_index, Category) %>% head()
Use LightGBM’s categorical data feature for optimial performance.
Use caret
for train/test splitting since createDataPartition
ensures balanced distribution of categories between train and test.
# declare categorical feature names
categoricals <- NULL
# proportion of data to train on
split <- 0.7
set.seed(123)
trainIndex <- createDataPartition(df$category_index, p = split, list = FALSE, times = 1)
dtrain <- lgb.Dataset((df %>% select(X, Y, hour, month, year, DayOfWeek) %>% data.matrix())[trainIndex,],
colnames = c("X", "Y", "hour", "month", "year", "DayOfWeek"),
categorical_feature = categoricals,
label = df$category_index[trainIndex], free_raw_data=T)
dtest <- lgb.Dataset.create.valid(dtrain,
(df %>% select(X, Y, hour, month, year, DayOfWeek) %>% data.matrix())[-trainIndex,],
label = df$category_index[-trainIndex])
params <- list(objective = "multiclass", metric = "multi_logloss")
valids <- list(test=dtest)
num_classes <- length(unique(df$category_index))
# preformat sizes for use in data visualizations later
train_size_format <- length(trainIndex) %>% format(big.mark=",")
test_size_format <- (df %>% nrow() - length(trainIndex)) %>% format(big.mark=",")
The size of the training set is 444,011 and the size of the test set is 190,288.
# determine elapsed runtime
system.time(
# training output not printed to notebook since spammy. (verbose = 0 + record = T)
bst <- lgb.train(params,
dtrain,
nrounds = 500,
valids,
num_threads = 4,
num_class = num_classes,
verbose = 0,
record = T,
early_stopping_rounds = 5,
categorical_feature = categoricals
)
)[3]
elapsed
219.019
# multilogloss of final iteration on test set
paste("# Rounds:", bst$current_iter())
[1] "# Rounds: 147"
paste("Multilogloss of best model:", bst$record_evals$test$multi_logloss$eval %>% unlist() %>% tail(1))
[1] "Multilogloss of best model: 1.97979254519616"
Calculate variable importance. (note: takes awhile since single-threaded)
df_imp <- tbl_df(lgb.importance(bst, percentage = TRUE))
Group 1 summed to more than type 'integer' can hold so the result has been coerced to 'numeric' automatically, for convenience.
df_imp
preds
is a 1D vector of probabilities for each vector, of nrows x nclasses. Reshape accordingly and iterate through for the predicted label (label with the largest probability) and the corresponding probability.
test <- (df %>% select(X, Y, hour, month, year, DayOfWeek) %>% data.matrix())[-trainIndex,]
preds_matrix <- predict(bst, test, reshape=T)
preds_cor <- cor(preds_matrix)
preds_matrix[1:2,]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.0013608454 0.08942618 5.662892e-05 0.0004694346 0.056512186 0.03050866 0.00103404
[2,] 0.0008612873 0.13593373 2.485057e-05 0.0019040408 0.005572041 0.00338940 0.02117221
[,8] [,9] [,10] [,11] [,12] [,13]
[1,] 0.06979163 0.005465058 7.687412e-05 4.486161e-05 0.0002551846 0.004679881
[2,] 0.10628417 0.017717549 4.029097e-04 8.352999e-05 0.0009054910 0.004290717
[,14] [,15] [,16] [,17] [,18] [,19]
[1,] 0.004730996 1.015853e-04 0.001967209 0.030442149 0.0007290578 0.0125450954
[2,] 0.004323221 5.753834e-05 0.007895822 0.005553753 0.0027784636 0.0004657167
[,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 0.002873131 0.002140685 0.2488873 8.046594e-05 0.0025070133 6.230313e-05 0.01176494
[2,] 0.003402504 0.004404593 0.4764304 4.563271e-05 0.0004264386 2.261202e-04 0.06097371
[,27] [,28] [,29] [,30] [,31] [,32]
[1,] 2.049359e-05 0.035559700 0.01084553 2.862660e-05 0.032270854 4.121095e-05
[2,] 2.255887e-05 0.007379504 0.00127614 9.229268e-05 0.007677787 5.485632e-05
[,33] [,34] [,35] [,36] [,37] [,38] [,39]
[1,] 0.001107487 2.752509e-05 0.034151278 0.02489355 0.019603917 0.25434873 0.008587664
[2,] 0.001217908 1.290809e-05 0.004451146 0.01504172 0.009998015 0.07428322 0.012966080
# likely not most efficient method
results <- t(apply(preds_matrix, 1, function (x) {
max_index = which(x==max(x))
return (c(max_index-1, x[max_index]))
}))
df_results <- data.frame(results, label_act = df$category_index[-trainIndex]) %>%
tbl_df() %>%
transmute(label_pred = X1, prod_pred = X2, label_act)
df_results %>% arrange(desc(prod_pred)) %>% head(20)
rm(preds_matrix)
Confusion matrix:
cm <- confusionMatrix(df_results$label_pred, df_results$label_act)
longer object length is not a multiple of shorter object lengthLevels are not in the same order for reference and data. Refactoring data to match.
data.frame(cm$overall)
df_imp$Feature <- factor(df_imp$Feature, levels=rev(df_imp$Feature))
plot <- ggplot(df_imp, aes(x=Feature, y=Gain)) +
geom_bar(stat="identity", fill="#34495e", alpha=0.9) +
geom_text(aes(label=sprintf("%0.1f%%", Gain*100)), color="#34495e", hjust=-0.25, family="Open Sans Condensed Bold", size=2.5) +
fte_theme() +
coord_flip() +
scale_y_continuous(limits = c(0, 0.4), labels=percent) +
theme(plot.title=element_text(hjust=0.5), axis.title.y=element_blank()) +
labs(title="Feature Importance for SF Arrest Type Model", y="% of Total Gain in LightGBM Model")
max_save(plot, "imp", "SF Open Data", h=2)
Plot the confusion matrix. Fortunately, matrix is already in long format.
df_cm <- tbl_df(data.frame(cm$table))
df_cm %>% head(100)
Map the labels to the indices.
# create mapping df
df_labels <- df %>%
select(category_index, Category) %>%
group_by(category_index, Category) %>%
summarize() %>%
ungroup() %>%
mutate(category_index = factor(category_index))
df_cm <- df_cm %>%
left_join(df_labels, by = c("Prediction" = "category_index")) %>%
left_join(df_labels, by = c("Reference" = "category_index")) %>%
rename(label_pred = Category.x, label_act = Category.y)
df_cm %>% head(100)
Plot the confusion matrix. Since 39 labels, confusion matrix will be large to fit all labels. Will also need to log-scale.
# create a data frame of "correct values" to annotate
df_correct <- df_cm %>% filter(label_pred == label_act)
plot <- ggplot(df_cm, aes(x=label_act, y=label_pred, fill = Freq)) +
geom_tile() +
geom_point(data=df_correct, color="white", size=0.8) +
fte_theme() +
coord_equal() +
scale_x_discrete() +
scale_y_discrete() +
theme(legend.title = element_text(size=7, family="Open Sans Condensed Bold"), legend.position="top", legend.direction="horizontal", legend.key.width=unit(1.25, "cm"), legend.key.height=unit(0.25, "cm"), legend.margin=unit(0,"cm"), axis.text.x=element_text(angle=-90, size=6, vjust=0.5, hjust=0), axis.text.y=element_text(size=6), plot.title = element_text(hjust=1)) +
scale_fill_viridis(name="# of Preds", labels=comma, breaks=10^(0:4), trans="log10") +
labs(title = sprintf("Confusion Matrix between %s Predicted SFPD Arrest Labels and Actual", test_size_format),
x = "Actual Label of Arrest",
y = "Predicted Label of Arrest")
`legend.margin` must be specified using `margin()`. For the old behavior use legend.spacing
max_save(plot, "confusionMatrix", "SF Open Data", h=6, w=5, tall=T)
Covert the preds_cor
matrix into long (adapted from http://stackoverflow.com/a/26838774)
Requires reordering correlations for cleaner chart: http://www.sthda.com/english/wiki/ggplot2-quick-correlation-matrix-heatmap-r-software-and-data-visualization)
dd <- as.dist((1-preds_cor)/2)
hc <- hclust(dd, "centroid")
label_order <- hc$order
preds_cor_reorder <- preds_cor[label_order, label_order]
df_corr <- tbl_df(data.frame(Var1=c(row(preds_cor_reorder))-1, Var2=c(col(preds_cor_reorder))-1, value = c(preds_cor_reorder))) %>%
filter(Var1 <= Var2) %>%
mutate(Var1 = factor(Var1), Var2=factor(Var2))
df_corr %>% head(100)
Plot similar chart to confusion matrix.
df_corr <- df_corr %>%
left_join(df_labels, by = c("Var1" = "category_index")) %>%
left_join(df_labels, by = c("Var2" = "category_index")) %>%
mutate(label1 = factor(Category.x), label2 = factor(Category.y))
# fix the label order to the reordered order from the hclust
levels(df_corr$label1) <- levels(df_corr$label1)[label_order]
levels(df_corr$label2) <- levels(df_corr$label2)[label_order]
plot <- ggplot(df_corr, aes(x=label1, y=label2, fill=value)) +
geom_tile() +
fte_theme() +
scale_x_discrete() +
scale_y_discrete() +
coord_fixed() +
theme(legend.title = element_text(size=7, family="Open Sans Condensed Bold"), legend.position="top", legend.direction="horizontal", legend.key.width=unit(1.25, "cm"), legend.key.height=unit(0.25, "cm"), legend.margin=unit(0,"cm"), panel.margin=element_blank(), axis.text.x=element_text(angle=-90, vjust=0.5, hjust=0), axis.title.y=element_blank(), axis.title.x=element_blank(), plot.title=element_text(hjust=1, size=6)) +
scale_fill_gradient2(high = "#2ecc71", low = "#e74c3c", mid = "white",
midpoint = 0, limit = c(-0.5,0.5),
name="Pearson\nCorrelation", breaks=pretty_breaks(8)) +
labs(title = sprintf("Correlations between Predicted Multiclass Probabilities of %s SFPD Arrest Category Labels", test_size_format))
`panel.margin` is deprecated. Please use `panel.spacing` property instead`legend.margin` must be specified using `margin()`. For the old behavior use legend.spacing
max_save(plot, "correlationMatrix", "SF Open Data", h=6, w=5, tall=T)
Reusing my mapping previous SF code: https://github.com/minimaxir/sf-arrests-when-where/blob/master/crime_data_sf.ipynb
bbox = c(-122.516441,37.702072,-122.37276,37.811818)
map <- get_map(location = bbox, source = "stamen", maptype = "toner-lite")
Source : http://tile.stamen.com/toner-lite/13/1308/3165.png
Source : http://tile.stamen.com/toner-lite/13/1309/3165.png
Source : http://tile.stamen.com/toner-lite/13/1310/3165.png
Source : http://tile.stamen.com/toner-lite/13/1311/3165.png
Source : http://tile.stamen.com/toner-lite/13/1308/3166.png
Source : http://tile.stamen.com/toner-lite/13/1309/3166.png
Source : http://tile.stamen.com/toner-lite/13/1310/3166.png
Source : http://tile.stamen.com/toner-lite/13/1311/3166.png
Source : http://tile.stamen.com/toner-lite/13/1308/3167.png
Source : http://tile.stamen.com/toner-lite/13/1309/3167.png
Source : http://tile.stamen.com/toner-lite/13/1310/3167.png
Source : http://tile.stamen.com/toner-lite/13/1311/3167.png
Source : http://tile.stamen.com/toner-lite/13/1308/3168.png
Source : http://tile.stamen.com/toner-lite/13/1309/3168.png
Source : http://tile.stamen.com/toner-lite/13/1310/3168.png
Source : http://tile.stamen.com/toner-lite/13/1311/3168.png
Create 40000 million latitude/longitude points in SF to simulate locations (200 points on x axis, 2000 points on y axis)
grid_size = 200
df_points <- data.frame(expand.grid(X=seq(bbox[1], bbox[3], length.out=grid_size),
Y=seq(bbox[2], bbox[4], length.out=grid_size)
)
)
df_points %>% head()
df_points %>% nrow()
[1] 40000
Predict arrest types at each point on April 15th, 2017, at 8 PM.
Populate data with same format of data (i.e. add month, hour, year, DayOfWeek). Does not require much customization. (DayOfWeek is a Factor, however)
date_target <- as.POSIXct("2017-04-15 20:00:00")
df_points <- df_points %>%
mutate(hour = format(date_target, "%H"),
month = format(date_target, "%m"),
year = format(date_target, "%Y"),
DayOfWeek = which(levels(df$DayOfWeek) == format(date_target, "%A"))) %>%
data.matrix()
df_points %>% head()
X Y hour month year DayOfWeek
[1,] -122.5164 37.70207 20 4 2017 6
[2,] -122.5157 37.70207 20 4 2017 6
[3,] -122.5150 37.70207 20 4 2017 6
[4,] -122.5143 37.70207 20 4 2017 6
[5,] -122.5136 37.70207 20 4 2017 6
[6,] -122.5128 37.70207 20 4 2017 6
preds_matrix <- matrix(predict(bst, df_points), byrow=T, nrow(df_points), num_classes)
results <- t(apply(preds_matrix, 1, function (x) {
max_index = which(x==max(x))
return (c(max_index-1, x[max_index]))
}))
rm(preds_matrix)
df_results <- data.frame(X=df_points[,1], Y=df_points[,2], label=factor(results[,1]), prob=results[,2]) %>%
tbl_df() %>%
left_join(df_labels, by=c("label" = "category_index")) %>%
mutate(Category = factor(Category))
joining factors with different levels, coercing to character vector
df_results %>% head(20)
plot <- ggmap(map) +
geom_raster(data = df_results %>% filter(Category != "Other Offenses"), aes(x=X, y=Y, fill=Category), alpha=0.8, size=0) +
coord_cartesian() +
fte_theme() +
scale_fill_brewer(palette = "Dark2") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.title = element_text(size=7, family="Open Sans Condensed Bold"), legend.position="right", legend.key.width=unit(0.5, "cm"), legend.key.height=unit(2, "cm"), legend.margin=margin(1,0,1,0), plot.title=element_text(hjust=0, size=11)) +
labs(title = sprintf("Locations of Predicted Types of Arrests in San Francisco on %s",
format(date_target, '%B %d, %Y at%l %p')))
Ignoring unknown parameters: size
max_save(plot, sprintf("crime-%s", format(date_target, '%Y-%m-%d-%H')), "SF Open Data", w = 6, h = 6, tall=T)
Set each label such that is has a consistent color.
set.seed(123)
cols <- distinctColorPalette(num_classes)
names(cols) <- df_labels$Category
cols
Arson Assault
"#E44D3E" "#62E2B5"
Bad Checks Bribery
"#71EA49" "#947468"
Burglary Disorderly Conduct
"#DE99E3" "#7957DB"
Driving Under The Influence Drug/Narcotic
"#6DC0E4" "#E0BE43"
Drunkenness Embezzlement
"#D986A7" "#D9B97C"
Extortion Family Offenses
"#DFE5B0" "#CA77E5"
Forgery/Counterfeiting Fraud
"#5F9BDE" "#815493"
Gambling Kidnapping
"#66E4E2" "#7A9E53"
Larceny/Theft Liquor Laws
"#E08F47" "#61E386"
Loitering Missing Person
"#D7EA40" "#E14582"
Non-Criminal Other Offenses
"#ADE79E" "#7E35E8"
Pornography/Obscene Mat Prostitution
"#E348C4" "#AEDFE3"
Recovered Vehicle Robbery
"#E16DBC" "#E4E9DD"
Runaway Secondary Codes
"#6F7BDF" "#DA6B6D"
Sex Offenses, Forcible Sex Offenses, Non Forcible
"#D542E5" "#AEA3E2"
Stolen Property Suicide
"#E3B5D9" "#D9C1B7"
Suspicious Occ Trea
"#AEE8C8" "#CFCDE6"
Trespass Vandalism
"#DEE37D" "#E6A393"
Vehicle Theft Warrants
"#7C8DA6" "#A4E165"
Weapon Laws
"#70A596"
Plot the map for 24 hours (convert to a GIF using external tools). Reuse code above to generate a map for given date/time + hour delta.
system("mkdir -p map_ani")
create_arrest_map <- function(hour_delta, date) {
date_target <- date + hour_delta*60*60
grid_size <- 200
df_points <- data.frame(expand.grid(X=seq(bbox[1], bbox[3], length.out=grid_size),
Y=seq(bbox[2], bbox[4], length.out=grid_size)
)
)
df_points <- df_points %>%
mutate(hour = format(date_target, "%H"),
month = format(date_target, "%m"),
year = format(date_target, "%Y"),
DayOfWeek = which(levels(df$DayOfWeek) == format(date_target, "%A"))) %>%
data.matrix()
preds_matrix <- matrix(predict(bst, df_points), byrow=T, nrow(df_points), num_classes)
results <- t(apply(preds_matrix, 1, function (x) {
max_index = which(x==max(x))
return (c(max_index-1, x[max_index]))
}))
rm(preds_matrix)
df_results <- data.frame(X=df_points[,1], Y=df_points[,2], label=factor(results[,1]), prob=results[,2]) %>%
tbl_df() %>%
left_join(df_labels, by=c("label" = "category_index")) %>%
mutate(Category = factor(Category))
plot <- ggmap(map) +
geom_raster(data = df_results %>% filter(Category != "Other Offenses"), aes(x=X, y=Y, fill=Category), alpha=0.9, size=0) +
coord_cartesian() +
fte_theme() +
scale_fill_manual(values=cols) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.title = element_text(family="Open Sans Condensed Bold"), legend.position="right", legend.key.width=unit(0.5, "cm"), legend.margin=margin(0,0,0,0, "cm"), legend.key.height=unit(1, "cm"), plot.title=element_text(hjust=0, size=11), legend.text.align=0) +
labs(title = sprintf("Locations of Predicted Types of Arrests in San Francisco on %s",
format(date_target, '%B %d, %Y at %l %p')))
max_save(plot, sprintf("map_ani/crime-%s", format(date_target, '%Y-%m-%d-%H')), "SF Open Data", w = 6, h = 6, tall=T)
}
base_date <- as.POSIXct("2017-03-14 06:00:00")
hour_deltas <- 0:23
x <- lapply(hour_deltas, create_arrest_map, base_date)
joining factors with different levels, coercing to character vectorIgnoring unknown parameters: size
The categorical approach using LightGBM is better. Here is the former code using OHE.
Categorical Features must be factors for one-hot encoding.
Convert the factor variables into dummy variables: model.matrix()
can do this in R natively. (via Stack Overflow)
# model.matrix() adds an Intercept column: the "-1" removes it.
# Matrix converts the dense matrix to sparse (reduces memory footprint to 25%).
train <- Matrix(model.matrix(~ X + Y + hour + month + year + DayOfWeek - 1, df))
num_classes <- length(unique(df$category_index))
num_rows <- nrow(train)
train[1:10,]
The objective is multi_logloss
since there are many classes. The multiclass
objective returns a probability for each class.
Demo: https://github.com/Microsoft/LightGBM/blob/master/R-package/tests/testthat/test_basic.R#L29
preds <- predict(bst, train[1:2,])
preds
length(preds)
The MIT License (MIT)
Copyright (c) 2017 Max Woolf
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.