final-project.qmd

---
title: "Final Project"
author: "Ziwen Lu, Xuyan Xiu, Doris Yan"
format: 
  html:
    embed-resources: true
    code-fold: true
    theme: united
editor: visual
execute: 
  warning: false
  cache: true
---

```{r load packages, output = FALSE}

library(tidyverse)
library(tidymodels)
library(ggplot2)
library(parsnip)
library(doFuture)
library(doParallel)
library(xgboost)
library(vip)
library(haven)
library(here)
library(icpsrdata)
library(plotly)
dir.create("data")

```

# Data Wrangling

## ICPSR credential 

```{r}

options("icpsr_email" = "input email address", 
        "icpsr_password" = "input password")
```

## Reproducible data retrieval


```{r}
# crime predictive data takes 1-2 mins

icpsr_download(
  file_id = 38649,
  download_dir = here("data")
)

# crime implementation crime data takes 1-2 mins

icpsr_download(
  file_id = 37059,
  download_dir = here("data")
)

# Ses data takes 9-10 mins

icpsr_download(
  file_id = 38528,
  download_dir = here("data")
)

# school counts takes 1-2 mins

icpsr_download(
  file_id = 38569,
  download_dir = here("data")
)

# street connectivity takes 1-2 mins

icpsr_download(
  file_id = 38580,
  download_dir = here("data")
)

# pollution takes 1 min

icpsr_download(
  file_id = 38597,
  download_dir = here("data")
)

### land cover takes 1-2 mins

icpsr_download(
  file_id = 38598,
  download_dir = here("data")
)

```

## Read data

```{r}

# crime modeling data

crime_predictive <- 
  read_dta(
    here("data",
         "ICPSR_38649",
         "DS0001",
         "38649-0001-Data.dta"))

# crime implementation data

crime_implementation <- 
  read_dta(
    here("data",
         "ICPSR_37059",
         "DS0004",
         "37059-0004-Data.dta"))

# ses predictive data and implementation data

ses08_17 <-
  read_dta(
    here("data",
         "ICPSR_38528",
         "DS0002",
         "38528-0002-Data.dta"))

# school counts data and implementation data

schoolct <- 
  read_dta(
    here("data",
         "ICPSR_38569",
         "DS0001",
         "38569-0001-Data.dta"))

# street connectivity predictive data

streetcon_predictive <- 
  read_dta(
    here("data",
         "ICPSR_38580",
         "DS0001",
         "38580-0001-Data.dta"))

# street connectivity implementation data

streetcon_implementation <- 
  read_dta(
    here("data",
         "ICPSR_38580",
         "DS0003",
         "38580-0003-Data.dta"))

# pollution predictive data and implementation data

pollution <- 
  read_dta(
    here("data",
         "ICPSR_38597",
         "DS0001",
         "38597-0001-Data.dta"))

# land cover predictive data and implementation data

landcover <- 
  read_dta(
    here("data",
         "ICPSR_38598",
         "DS0001",
         "38598-0001-Data.dta"))

```

## Clean data

### crime predictive

We use annual measure of crime.

```{r}

crime_2010 <- 
  crime_predictive|>
  rename_all(tolower)|>
  rename(county = stcofips)|>
  select(county, year, viol, property, cpopcrim)|>
  filter(year == 2010)|>
  
  # generate crime rate per 100,000 people and the log transform
  
  mutate(viol_rate = (viol / cpopcrim) * 100000,
         property_rate = (property / cpopcrim) * 100000)|>
  mutate(log_viol = log(viol_rate),
       log_property = log(property_rate))

```

### crime implementation

```{r}

crime_2016 <- 
  crime_implementation|>
  rename_all(tolower)|>
  select(fips_st, fips_cty, cpopcrim, viol, property)|>
  
  # concatenate fips_st and fips_cty to make county fips
  
  mutate(
    fips_st = ifelse(fips_st >= 1 & fips_st <= 9, 
                          paste0("0", fips_st), 
                          as.character(fips_st)) 
         )|>
  mutate(
    fips_cty = str_pad(fips_cty, 
                       width = 3, 
                       pad = "0")
         )|>
  mutate(county = str_c(fips_st, fips_cty))|>
  select(-fips_cty, -fips_st)|>
  
  # generate crime rate per 100,000 people and the log transform
  
  mutate(viol_rate = (viol / cpopcrim) * 100000,
         property_rate = (property / cpopcrim) * 100000)|>
  mutate(log_viol = log(viol_rate),
       log_property = log(property_rate))|>
  select(county, everything())
  
```

```{r}

# create a function for effective data cleaning

clean_data<-
  function(data) {
    cleaned_data <- data |> 
      rename_all(tolower)|>
      mutate(
        county = substr((tract_fips10), 1, 5))|>
      select(-tract_fips10)

  return(cleaned_data)
  }

```

### ses

SES predictive data of 2010 is from ACS five-year estimate from 2008 to 2012. SES implementation data of 2016 is from ACS five-years estimate from 2013 to 2016.

```{r}

# ACS 2008 to 2012 for year 2010

ses_2010 <- 
  clean_data(ses08_17)|>
  group_by(county)|>
  summarise_all(sum, na.rm = TRUE)|>
  pivot_longer(
  cols = !county, 
  names_to = "indicator",
  values_to = "value"
  )|>
  filter(str_sub(indicator, -2) == "12")|>
  
  # pivot SES indicators back as variables
  
  pivot_wider(
    names_from = indicator,
    values_from = value
  )|>
  
  # cpopcrim and totpop08_12 are all total population estimate. 
  # cpopcrim is from crime dataset, totpop08_12 is from ses dataset
  
  # since we are predicting crime rate, it would be more appropriate to use 
  # the population variable in the crime dataset, thus we will use cpopcrim
  # as population estimate

  select(-totpop08_12)
   
# ACS 2013 to 2017 for year 2016

ses_2016 <- 
  clean_data(ses08_17)|>
  group_by(county)|>
  summarise_all(sum, na.rm = TRUE)|>
  pivot_longer(
  cols = !county, 
  names_to = "indicator",
  values_to = "value"
  )|>
  filter(str_sub(indicator, -2) == "17")|>
  
  # pivot SES indicators back as variables
  
  pivot_wider(
    names_from = indicator,
    values_from = value
  )|>
  select(-totpop13_17)

```

### school counts

We use annual school counts.

```{r}

# school counts at 2016

schoolct_2010 <- 
  clean_data(schoolct)|>
  group_by(county,year)|>
  summarise_all(sum, na.rm = TRUE)|>
  filter(year == 2010)

schoolct_2016 <- 
  clean_data(schoolct)|>
  group_by(county,year)|>
  summarise_all(sum, na.rm = TRUE)|>
  filter(year == 2016)

```

### street connectivity

We use street connectivity data at 2010 for predictive data, and street connectivity data at 2020 for implementation data.

```{r}

# from 2010 school count data

streetcon_2010 <- 
  clean_data(streetcon_predictive)|>
  group_by(county)|>
  summarise_all(sum, na.rm = TRUE)

# use 2020 data to represent 2016 school counts

streetcon_2016 <- 
  streetcon_implementation|>
  rename_all(tolower)|>
  mutate(county = substr((tract_fips20), 1, 5))|>
  select(-tract_fips20)|>
  group_by(county)|>
  summarise_all(sum, na.rm = TRUE)

```

### pollution

We use annual measure of pollution.

```{r}

# pollution in 2010

pollution_2010 <-
  clean_data(pollution)|>
  group_by(county,year)|>
  summarise_all(sum, na.rm = TRUE)|>
  filter(year == 2010)

# pollution in 2016

pollution_2016 <-
  clean_data(pollution)|>
  group_by(county,year)|>
  summarise_all(sum, na.rm = TRUE)|>
  filter(year == 2016)

```

### land cover

We use annual interpolation of land cover.

```{r}

# land cover at 2010

landcover_2010 <- 
  landcover|>
  rename_all(tolower)|>
  mutate(county = substr((tract_fips), 1, 5))|>
  select(-tract_fips, -intp_flag)|>
  group_by(county, year_intp)|>
  summarise_all(sum, na.rm = TRUE) |>
  filter(year_intp == 2010)

# land cover at 2016

landcover_2016 <- landcover|>
  rename_all(tolower)|>
  mutate(county = substr((tract_fips), 1, 5))|>
  select(-tract_fips, -intp_flag)|>
  group_by(county, year_intp)|>
  summarise_all(sum, na.rm = TRUE) |>
  filter(year_intp == 2016)

```

## Combine All Datasets

### predictive data

```{r}

# combine the data sets with the by "county" column

crime_all_2010 <- 
  full_join(
  crime_2010,
  full_join(
    ses_2010,
   full_join(
     schoolct_2010,
    full_join(
      streetcon_2010,
      full_join(
        pollution_2010, landcover_2010, by = "county"
        ),
      by = "county"
      ),
    by = "county"
    ),
   by = "county"
   ),
  by = "county"
  )

crime_all_2010<-crime_all_2010|>
  select(-year.y, -year.x, -year_intp)

```

### implementation data

```{r}

# combine the data sets with the by "county" and "year" columns

crime_all_2016 <- 
  full_join(
  crime_2016,
  full_join(
    ses_2016,
   full_join(
     schoolct_2016,
    full_join(
      streetcon_2016,
      full_join(
        pollution_2016, landcover_2010, by = "county"
        ),
      by = "county"
      ),
    by = "county"
    ),
   by = "county"
   ),
  by = "county"
  )

crime_all_2016 <- 
  crime_all_2016|>
  select(-year.y, -year.x, -year_intp)|>
  select(county, everything())|>
  mutate(year = 2016)

```

## 2010 predictive data EDA and data cleaning

We then removed STATA label and attributes in outcome variables.

```{r}

# remove STATA format and attributes

attr(crime_all_2010$viol, "format.stata") <- NULL
attr(crime_all_2010$property, "format.stata") <- NULL
attr(crime_all_2010$cpopcrim, "format.stata") <- NULL
attr(crime_all_2010$viol_rate, "format.stata") <- NULL
attr(crime_all_2010$property_rate, "format.stata") <- NULL
attr(crime_all_2010$log_viol, "format.stata") <- NULL
attr(crime_all_2010$log_property, "format.stata") <- NULL


attr(crime_all_2010$viol, "label") <- NULL
attr(crime_all_2010$property, "label") <- NULL
attr(crime_all_2010$cpopcrim, "label") <- NULL
attr(crime_all_2010$viol_rate, "label") <- NULL
attr(crime_all_2010$property_rate, "label") <- NULL
attr(crime_all_2010$log_viol, "label") <- NULL
attr(crime_all_2010$log_property, "label") <- NULL
 
str(crime_all_2010)

# the suffix for SES indicator 2008-2012 and 2013-2017 is different, we remove the suffix so that the variable names in predictor and implementation data are the same
remove_string_2010 <- 
  function(varname){
  sub("08_12","",varname)
}

crime_all_2010_cleaned <- 
  crime_all_2010|>
  
  # eliminate missing values and negative infinity for outcome variables
  
  filter(!is.na(log_viol) | !is.na(log_property)) |>
  filter(!(is.infinite(log_viol) | is.infinite(log_property)))|>
  rename_with(remove_string_2010, contains("08_12"))


# explore missing values

missing_values_2010 <-
  colSums(is.na(crime_all_2010_cleaned))

missing_values_2010 <-
  as.data.frame(missing_values_2010)

# look at the variables with the most missing values, 
# no variable has over 50% missing values

missing_values_percent_2010<-missing_values_2010|>
  mutate(missing_percent = missing_values_2010/nrow(crime_all_2010_cleaned))|>
  arrange(desc(missing_percent))|>
  top_n(10)

```

## 2016 implementation data EDA and data cleaning

Repeat the above data cleaning process for implementation data in 2016.

```{r}
attr(crime_all_2016$viol, "format.stata") <- NULL
attr(crime_all_2016$property, "format.stata") <- NULL
attr(crime_all_2016$cpopcrim, "format.stata") <- NULL
attr(crime_all_2016$viol_rate, "format.stata") <- NULL
attr(crime_all_2016$property_rate, "format.stata") <- NULL
attr(crime_all_2016$log_viol, "format.stata") <- NULL
attr(crime_all_2016$log_property, "format.stata") <- NULL


attr(crime_all_2016$viol, "label") <- NULL
attr(crime_all_2016$property, "label") <- NULL
attr(crime_all_2016$cpopcrim, "label") <- NULL
attr(crime_all_2016$viol_rate, "label") <- NULL
attr(crime_all_2016$property_rate, "label") <- NULL
attr(crime_all_2016$log_viol, "label") <- NULL
attr(crime_all_2016$log_property, "label") <- NULL
 
str(crime_all_2016)

remove_string_2016 <-
  function(varname){
  sub("13_17","",varname)
}

crime_all_2016_cleaned <-
  crime_all_2016|>
  
  # eliminate missing values for outcome variables
  
  filter(!is.na(log_viol) | !is.na(log_property)) |>
  filter(!(is.infinite(log_viol) | is.infinite(log_property)))|>
  rename_with(remove_string_2016, contains("13_17"))


# explore missing values

missing_values_2016 <- 
  colSums(is.na(crime_all_2016_cleaned))

missing_values_2016 <-
  as.data.frame(missing_values_2016)

# look at the variables with the most missing values, 
# no variable has over 50% missing values

missing_values_percent_2016 <-
  missing_values_2016|>
  mutate(missing_percent = missing_values_2016/nrow(crime_all_2010_cleaned))|>
  arrange(desc(missing_percent))|>
  top_n(10)

```

# Crime Rate Maps

Before starting the machine learning process, we used package `plotly` to map the average violent and property crime rates across the U.S. from 2002 to 2014, displayed on a state-by-state basis for better understanding.

## Build dataset

We generate crime rate as the number of crimes per 100,000 people.

$$
Rate = \left( \frac{Crime}{Population} \right)10 ^ 5
$$ {#eq-binomial}

```{r}

crime_02_14 <- 
  crime_predictive|>
  rename_all(tolower)|>
  rename(county=stcofips)|>
  mutate(viol = murder+rape+robbery+agasslt)|>
  mutate(property = burglry+larceny+mvtheft)|>
  
  select(county, year, viol, property, cpopcrim) |>
  mutate(viol_rate = (viol / cpopcrim) * 100000) |>
  mutate(property_rate = (property / cpopcrim) * 100000) |>
  
  filter(!is.na(viol_rate) | !is.na(property_rate)) |>
  filter(!(is.infinite(viol_rate) | is.infinite(property_rate)))

```

## Extract FIPS codes

We then converted county units to states based on the five-digit FIPS code to more visually compare differences between regions.

```{r}

crime_02_14$state_fips <- 
  substr(crime_02_14$county, 
         1, 2)

state_fips_to_abb <- 
  tibble(
  state_fips = c("01", "02", "04", "05", "06", 
                 "08", "09", "10", "11", "12", 
                 "13", "15", "16", "17", "18", 
                 "19", "20", "21", "22", "23", 
                 "24", "25", "26", "27", "28", 
                 "29", "30", "31", "32", "33", 
                 "34", "35", "36", "37", "38", 
                 "39", "40", "41", "42", "44", 
                 "45", "46", "47", "48", "49", 
                 "50", "51", "53", "54", "55", "56"),
  state_abb = c("AL", "AK", "AZ", "AR", "CA", "CO", 
                "CT", "DE", "DC", "FL", "GA", "HI", 
                "ID", "IL", "IN", "IA", "KS", "KY", 
                "LA", "ME", "MD", "MA", "MI", "MN", 
                "MS", "MO", "MT", "NE", "NV", "NH", 
                "NJ", "NM", "NY", "NC", "ND", "OH", 
                "OK", "OR", "PA", "RI", "SC", "SD", 
                "TN", "TX", "UT", "VT", "VA", "WA", 
                "WV", "WI", "WY")
)

```

## Choropleth map

Finally, we chose the choropleth map to illustrate the average violent crime and property crime rate in the U.S., 2002 - 2014. To keep the visualization informative, we used `plotly` package to add hover text and so on for interactivity. There are also other interactive settings for readers, like zoom in and out for details.

```{r}

# convert county-level data to state-level data, aggregate violent crime data

crime_data_map <- crime_02_14 |>
  group_by(state_fips) |>
  summarise(
    avg_viol_rate = mean(viol_rate, na.rm = TRUE),
    avg_property_rate = mean(property_rate, na.rm = TRUE)
  ) |>
  
  left_join(state_fips_to_abb, 
            by = 'state_fips') |>
  
  mutate(
    hover_violence = paste(state_abb, 
                           '<br>', 
                           "Violent Crimes Rate: ", 
                           round(avg_viol_rate, 2)),
    hover_property = paste(state_abb, 
                           '<br>', 
                           "Property Crimes Rate: ", 
                           round(avg_property_rate, 2))  
  ) |>
  
  pivot_longer(
    cols = c(avg_viol_rate, 
             avg_property_rate),
    names_to = "crime_type",
    values_to = "crime_count"
  )

```

## violent crime rate

```{r}

# create and output the map for violent crime

violence_map <- 
  plot_geo(crime_data_map |>
             filter(crime_type == "avg_viol_rate"),
           locationmode = 'USA-states') |>
  
  add_trace(
    z = ~crime_count, 
    text = ~hover_violence,  
    hoverinfo = 'text',     
    locations = ~state_abb,
    color = ~crime_count,
    zmin = 0, 
    zmax = 1000,
    colors = c("#1a9641", "#ffffbf", "#fdae61", "#d7191c"),
    colorbar = list(title = "Violent Crimes Rate",
                    tickvals = c(250, 500, 750))
  ) |>
  
  layout(
    title = 'Average US Crimes Rate by State, 2002 - 2014',
    geo = list(
      scope = 'usa',
      projection = list(type = 'albers usa'),
      showlakes = TRUE,
      lakecolor = toRGB('white')
    )
  )

```

## property crime rate

```{r}
# create and output the map for property crime rate

property_map <- 
  plot_geo(crime_data_map |>
             filter(crime_type == "avg_property_rate"),
           locationmode = 'USA-states') |>
  
  add_trace(
    z = ~crime_count, 
    text = ~hover_property,   
    hoverinfo = 'text',       
    locations = ~state_abb,
    color = ~crime_count,
    zmin = 0, 
    zmax = 5000,
    colors = c("#2b83ba", "#ffffbf", "#fdae61", "#d7191c"),
    colorbar = list(title = "Property Crimes Rate",
                    tickvals = c(1000, 2000, 3000, 4000))
  ) |>
  
  layout(
    title = 'Avergae US Crimes Rate by State, 2002 - 2014',
    geo = list(
      scope = 'usa',
      projection = list(type = 'albers usa'),
      showlakes = TRUE,
      lakecolor = toRGB('white')
    )
  )

```

## Combine the maps

```{r}

fig <- 
  plotly::subplot(violence_map,
                  property_map,
                  nrows = 2,
                  margin = 0.05)

# print the figure
fig

```

Areas with high crime rates are indicated by a darker red color. It can be seen that crime rate is higher in the south and near the border. There are also some interesting findings, such as the fact that Texas, which is usually considered to have a high rate of gun ownership, boasts a violent crime rate of less than 300, which could be influenced by population.

# Splitting Data

```{r}

set.seed(904)

crime_split <- 
  initial_split(crime_all_2010_cleaned, prop = 0.8)

crime_train <- 
  training(crime_split)
crime_test <- 
  testing(crime_split)

```

# Cross-validation

```{r}

# set up 5-fold cross validation

set.seed(904)
crime_folds <- 
  vfold_cv(data = crime_train, v = 5)

```

# Violent Crime

## Ridge Regression

### feature and target engineering

```{r}

# create a recipe

ridge_v_rec <- 
  recipe(log_viol ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_property,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())|>
  step_normalize(all_predictors())
 
# check if feature engineering works to remove variables

ridge_v_summary <-
  summary(ridge_v_rec)

# model
ridge_v_mod <- 
  linear_reg(penalty = tune(), mixture = 0) |>
  set_mode(mode = "regression") |>
  set_engine(engine = "glmnet")  

ridge_v_grid <-
  grid_regular(penalty(),
               levels = 20)
# workflow
ridge_v_wf <- 
  workflow() |>
  add_recipe(recipe = ridge_v_rec) |>
  add_model(spec = ridge_v_mod)

```

### choose the best model

```{r}

ridge_v_resamples <-
  ridge_v_wf |>
  tune_grid(
    resamples = crime_folds,
    grid = ridge_v_grid,
    metrics = metric_set(rmse)
    )

ridge_v_resamples |>
  collect_metrics()

top_ridge_v_models <-
  ridge_v_resamples |>
  show_best() |>
  arrange(penalty) 

ridge_v_resamples|>
  select_best()

ridge_v_resamples|>
  autoplot()

```

### prediction

```{r}

set.seed(904)
ridge_v_pre_wf <- 
  ridge_v_wf |> 
  finalize_workflow(select_best(ridge_v_resamples))

ridge_v_pre_wf <-
  fit(ridge_v_pre_wf, 
      data = crime_train)

ridge_v_predictions <- 
  predict(object = ridge_v_pre_wf, 
          new_data = crime_test)$.pred

actual <- 
  crime_test$log_viol

county <- 
  crime_test$county

ridge_v_predictions_df <-
  tibble(county,actual,ridge_v_predictions)

rmse_ridge_v<-rmse(ridge_v_predictions_df, actual, ridge_v_predictions)

```

## Lasso

### feature and target engineering

```{r}

lasso_v_rec <- 
  recipe(log_viol ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_property,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())|>
  step_normalize(all_predictors())

# check if feature engineering works to remove variables

lasso_v_summary <- 
  summary(lasso_v_rec)

# create a model

lasso_v_mod <- 
  linear_reg(penalty = tune(), mixture = 1) |>
  set_mode(mode = "regression") |>
  set_engine(engine = "glmnet")  
 
lasso_v_grid <-
  grid_regular(penalty(),
               levels = 20)

lasso_v_wf <-
  workflow() |>
  add_recipe(recipe = lasso_v_rec) |>
  add_model(spec = lasso_v_mod)

```

### choose the best model

```{r}

set.seed(904)
lasso_v_resamples <- 
  lasso_v_wf|>
  tune_grid(
    resamples = crime_folds,
    grid = lasso_v_grid,
    metrics = metric_set(rmse)
    )

lasso_v_resamples |>
  collect_metrics()
  
top_lasso_v_models <-
  lasso_v_resamples |>
  show_best() |>
  arrange(penalty) 

lasso_v_resamples|>
  select_best()

lasso_v_resamples|>
  autoplot()

```

### prediction

```{r}

set.seed(904)
lasso_v_pre_wf <- 
  lasso_v_wf |> 
  finalize_workflow(select_best(lasso_v_resamples))

lasso_v_pre_wf <-
  fit(lasso_v_pre_wf, 
      data = crime_train)

lasso_v_predictions <- 
  predict(object = lasso_v_pre_wf, 
          new_data = crime_test)$.pred

actual <- crime_test$log_viol

county <- crime_test$county

lasso_v_predictions_df <-
  tibble(county,actual,lasso_v_predictions)

rmse_lasso_v <- 
  rmse(lasso_v_predictions_df, actual, lasso_v_predictions)

```

## Random Forest

### feature and target engineering

```{r}

rf_v_rec <- 
  recipe(log_viol ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_property,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())

rf_v_mod <- 
  rand_forest(
  mtry = tune(),
  min_n = tune(),
  trees = 200
  ) |>
  set_mode(mode = "regression") |>
  set_engine(
    engine = "ranger", 
    importance = "impurity",
    num.threads = 4
  )

rf_v_wf <- 
  workflow() |>
  add_recipe(rf_v_rec) |>
  add_model(rf_v_mod)

rf_v_grid <- grid_regular(
  mtry(range = c(1, 15)),
  min_n(range = c(1, 15)),
  levels = 5
)

```

### choose the best model

```{r}

set.seed(904)
rf_v_resamples <- 
  tune_grid(
  rf_v_wf,
  resamples = crime_folds,
  grid = rf_v_grid)

collect_metrics(rf_v_resamples)

show_best(rf_v_resamples, 
          metric = "rmse") 

select_best(rf_v_resamples)
autoplot(rf_v_resamples)

```

### predictions

```{r}

set.seed(904)
rf_v_pre_wf <- 
  rf_v_wf |> 
  finalize_workflow(select_best(rf_v_resamples))

rf_v_pre_wf<-
  fit(rf_v_pre_wf, 
      data = crime_train)

rf_v_predictions <- 
  predict(object = rf_v_pre_wf, 
          new_data = crime_test)$.pred

actual <- 
  crime_test$log_viol

county <- 
  crime_test$county

rf_v_predictions_df<-
  tibble(county,actual,rf_v_predictions)

rmse_rf_v <- 
  rmse(rf_v_predictions_df, actual, rf_v_predictions)

```

### selecting the top 15 important variables

```{r}

set.seed(904)
rf_v_final_wf <- 
  rf_v_wf |> 
  finalize_workflow(select_best(rf_v_resamples))

rf_v_fit <- 
  rf_v_final_wf |>
  fit(data = crime_all_2010_cleaned)

crime_labels <- c(
  pnhblack = "Proportion of people non-Hispanic Black",
  pin1b_ = "Proportion of families with Income less than 15K",
  disadvantage  = "Disadvantage1",
  ppubas = "Prop. of households with public assistance income",
  ppov = "Prop. ppl w/ income past 12 months below poverty level",
  value_52 = "Road covered by Shrub/Scrub (sqm)",
  pfhfam = "Proportion female-headed families w/ kids",
  disadvantage2_ = "Disadvantage2",
  prop_value_41 = "Prop. road covered by deciduous Forest (sqm)",
  pnhwhite = "Proportion of people non-Hispanic White",
  ped1_ = "Prop. w/ less than high school diploma",
  prop_value_52 = "Prop. road covered by Shrub/Scrub (sqm)",
  phispanic = "Proportion of people of Hispanic origin",
  pin2b_ = "Proportion of families with Income 15-30K",
  pnwhite = "Proportion of people non-Hispanic White",
  value_41 = "Road covered by deciduous Forest (sqm)",
  punemp = "Prop. 16+ civ labor force unemployed",
  prop_value_21 = "Prop. Developed, Open Space (sqm)",
  value_24 = "Road Developed, High Intensity (sqm)",
  prop_value_24 = "Prop road Developed, High Intensity (sqm)",
  prop_value_23 = "Prop. Road Developed, Medium Intensity (sqm)",
  prop_value_22 = "Prop. Road Developed, Low Intensity (sqm)",
  p18_29  = "Proportion of population 18-29 years of age",
  popden = "Persons per square mile"
)

rf_v_fit|>
  extract_fit_parsnip()|>
  vip(num_features = 15) %>%
  .$data |>
  mutate(
    Importance = Importance / max(Importance),
    Variable = fct_reorder(.x = Importance, .f = Variable)
  ) |>
  ggplot(aes(Importance, Variable)) +
  scale_y_discrete(labels = crime_labels) +
  labs(title = "Top 15 Variables for Violent Crime Prediction") +
  geom_col()

```

Based on the result from random forest model, the top 15 important variables for predicting violent crime are above. Most of the variables are from the socioeconomic data and few are from the land cover data.

Disadvantage1 is the mean of proportion of people Non-Hispanic Black, proportion female-headed families with kids, proportion of households with public assistance income, proportion people with income past 12 months below poverty level, and proportion 16+ civilian labor force unemployed.

Disadvantage2 is the mean of proportion female-headed families with kids, proportion of households with public assistance income, proportion people with income past 12 months below poverty level, and proportion 16+ civilian labor force unemployed.

## MARS

### feature and target engineering

```{r}

# recipe

mars_v_rec <- 
  recipe(log_viol ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_property,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())

mars_v_rec_sum <-
  summary(mars_v_rec)

# model

mars_v_mod <- 
  mars(
  num_terms = tune(), 
  prod_degree = tune()
) |>
  set_mode(mode = "regression") |>
  set_engine(engine = "earth")

# tuning grid

mars_v_grid <- 
  grid_regular(
  num_terms(range = c(10, 100)),
  prod_degree(),
  levels = 10
)

# workflow

mars_v_wf <- 
  workflow() |>
  add_recipe(recipe = mars_v_rec) |>
  add_model(spec = mars_v_mod)

```

### choose the best model

```{r}

# resample   

set.seed(904)
mars_v_resamples <- 
  mars_v_wf|>
  tune_grid(
  resamples = crime_folds,
  grid = mars_v_grid,
  metrics = metric_set(rmse)
  )

collect_metrics(mars_v_resamples)

top_mars_v_models <-
  mars_v_resamples |>
  show_best()

mars_v_resamples|>
  select_best()

```

### prediction

```{r}

set.seed(904)
mars_v_pre_wf <- 
  mars_v_wf |> 
  finalize_workflow(select_best(mars_v_resamples))

mars_v_pre_wf<-
  fit(mars_v_pre_wf, 
      data = crime_train)

mars_v_predictions <- 
  predict(object = mars_v_pre_wf, 
          new_data = crime_test)$.pred

actual<-crime_test$log_viol

county<-crime_test$county

mars_v_predictions_df<-
  tibble(county,actual,mars_v_predictions)

rmse_mars_v<-rmse(mars_v_predictions_df, actual, mars_v_predictions)

```

## XG Boost

### feature and target engineering

```{r}

xgb_v_rec <- 
  recipe(log_viol ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_property,
              new_role = "id")|>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())

all_cores <- 
  parallel::detectCores(logical = FALSE)

registerDoFuture()

cl <- 
  makeCluster(all_cores - 3L)

plan(cluster, workers = cl)

xgb_v_spec <- 
  boost_tree(
  trees = 200,
  tree_depth = tune(), 
  min_n = tune(),
  loss_reduction = tune(),                     
  sample_size = tune(), 
  mtry = tune(),         
  learn_rate = tune()                   
) |>
  set_engine("xgboost") |>
  set_mode("regression")

set.seed(904)
xgb_v_grid <- 
  grid_latin_hypercube(
  tree_depth(),
  min_n(range = c(2, 10)),
  loss_reduction(range = c(-5, -3)),
  sample_size = sample_prop(),
  mtry(range = c(1, 10)),
  learn_rate(range = c(-5, -1)),
  size = 30
)

xgb_v_wf <- workflow() |>
  add_recipe(xgb_v_rec) |>
  add_model(xgb_v_spec)

```

### choose the best model

```{r}

set.seed(904)
xgb_v_resamples <- 
  tune_grid(
  xgb_v_wf,
  resamples = crime_folds,
  grid = xgb_v_grid)

collect_metrics(xgb_v_resamples)
show_best(xgb_v_resamples, metric = "rmse") 
select_best(xgb_v_resamples, metric= "rmse")

```

### prediction

```{r}

set.seed(904)
xgb_v_pre_wf <- 
  xgb_v_wf |> 
  finalize_workflow(select_best(xgb_v_resamples))

xgb_v_pre_wf <-
  fit(xgb_v_pre_wf, 
      data = crime_train)

xgb_v_predictions <- 
  predict(object = xgb_v_pre_wf, 
          new_data = crime_test)$.pred

actual <- 
  crime_test$log_viol

county <-
  crime_test$county

xgb_v_predictions_df <-
  tibble(county,actual,xgb_v_predictions)

rmse_xgb_v <- 
  rmse(xgb_v_predictions_df, actual, xgb_v_predictions)

```

## Final comparison of rmse

```{r}

violence_comparison <-
  bind_rows(rmse_ridge_v, rmse_lasso_v, rmse_rf_v, rmse_mars_v, rmse_xgb_v)

violence_comparison|>
  mutate(
    model = c("Ridge Regression","Lasso Regression",
              "Random Forest", "MARS", "XG Boost"))|>
  select(model, everything())

```

Because XG Boost has the lowest RMSE value, we will use XG boost to predict violent crime rate in 2016.

## Prediction on implementation violent crime data in 2016

```{r}

set.seed(904)
xgb_v_pre_wf_2016 <- 
  xgb_v_wf |> 
  finalize_workflow(select_best(xgb_v_resamples))

xgb_v_pre_wf_2016 <-
  fit(xgb_v_pre_wf, 
      data = crime_all_2010_cleaned)

xgb_v_predictions_2016 <- 
  predict(object = xgb_v_pre_wf_2016, 
          new_data = crime_all_2016_cleaned)$.pred

actual_implement_v <- crime_all_2016_cleaned$log_viol

county_implement <- crime_all_2016_cleaned$county

xgb_v_predictions_2016_df <-
  tibble(county_implement,actual_implement_v,xgb_v_predictions_2016)

rmse_xgb_v_2016 <-
  rmse(xgb_v_predictions_2016_df, actual_implement_v, xgb_v_predictions_2016)

rmse_xgb_v_2016

```

Final prediction for violent crime rate in 2016 has a RMSE of 0.593.

# Property Crime

## Ridge Regression

### feature and target engineering

```{r}

# create a recipe

ridge_p_rec <- 
  recipe(log_property ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_viol,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())|>
  step_normalize(all_predictors())
 
# check if feature engineering works to remove variables

ridge_p_summary <-
  summary(ridge_p_rec)

# model

ridge_p_mod <- 
  linear_reg(penalty = tune(), mixture = 0) |>
  set_mode(mode = "regression") |>
  set_engine(engine = "glmnet")  

ridge_p_grid<-
  grid_regular(penalty(),
               levels = 20)
# workflow

ridge_p_wf <- 
  workflow() |>
  add_recipe(recipe = ridge_p_rec) |>
  add_model(spec = ridge_p_mod)
```

### choose the best model

```{r}

ridge_p_resamples <-
  ridge_p_wf |>
  tune_grid(
    resamples = crime_folds,
    grid = ridge_p_grid,
    metrics = metric_set(rmse)
    )

ridge_p_resamples |>
  collect_metrics()

top_ridge_p_models <-
  ridge_p_resamples |>
  show_best() |>
  arrange(penalty) 

ridge_p_resamples|>
  select_best()

ridge_p_resamples|>
  autoplot()

```

### prediction

```{r}

set.seed(904)
ridge_p_pre_wf <- 
  ridge_p_wf |> 
  finalize_workflow(select_best(ridge_p_resamples))

ridge_p_pre_wf <-
  fit(ridge_p_pre_wf, 
      data = crime_train)

ridge_p_predictions <- 
  predict(object = ridge_p_pre_wf, 
          new_data = crime_test)$.pred

actual_p <- 
  crime_test$log_property

county <- 
  crime_test$county

ridge_p_predictions_df <-
  tibble(county,actual_p,ridge_p_predictions)

rmse_ridge_p <- 
  rmse(ridge_p_predictions_df, actual_p, ridge_p_predictions)

```

## Lasso

### feature and target engineering

```{r}

lasso_p_rec <- 
  recipe(log_property ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_viol,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())|>
  step_normalize(all_predictors())

# check if feature engineering works to remove variables

lasso_p_summary <-
  summary(lasso_p_rec)

# create a model

lasso_p_mod <- 
  linear_reg(penalty = tune(), mixture = 1) |>
  set_mode(mode = "regression") |>
  set_engine(engine = "glmnet")  
 
lasso_p_grid <-
  grid_regular(penalty(),
               levels = 20)
lasso_p_wf <-
  workflow() |>
  add_recipe(recipe = lasso_p_rec) |>
  add_model(spec = lasso_p_mod)

```

### choose the best model

```{r}

set.seed(904)
lasso_p_resamples <- 
  lasso_p_wf|>
  tune_grid(
    resamples = crime_folds,
    grid = lasso_p_grid,
    metrics = metric_set(rmse)
    )

lasso_p_resamples |>
  collect_metrics()
  
top_lasso_p_models <-
  lasso_p_resamples |>
  show_best() |>
  arrange(penalty) 

lasso_p_resamples|>
  select_best()

lasso_p_resamples|>
  autoplot()

```

### prediction

```{r}

set.seed(904)
lasso_p_pre_wf <- 
  lasso_p_wf |> 
  finalize_workflow(select_best(lasso_p_resamples))

lasso_p_pre_wf <-
  fit(lasso_p_pre_wf, 
      data = crime_train)

lasso_p_predictions <- 
  predict(object = lasso_p_pre_wf, 
          new_data = crime_test)$.pred

actual_p <- 
  crime_test$log_property

county <- 
  crime_test$county

lasso_p_predictions_df <-
  tibble(county,actual_p,lasso_p_predictions)

rmse_lasso_p <- 
  rmse(lasso_p_predictions_df, actual_p, lasso_p_predictions)

```

## Random Forest

### feature and target engineering

```{r}

rf_p_rec <- 
  recipe(log_property ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_viol,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())

rf_p_mod <- 
  rand_forest(
  mtry = tune(),
  min_n = tune(),
  trees = 200
  ) |>
  set_mode(mode = "regression") |>
  set_engine(
    engine = "ranger", 
    importance = "impurity",
    num.threads = 4
  )

rf_p_wf <- 
  workflow() |>
  add_recipe(rf_p_rec) |>
  add_model(rf_p_mod)

rf_p_grid <- 
  grid_regular(
  mtry(range = c(1, 15)),
  min_n(range = c(1, 15)),
  levels = 5
)

```

### choose the best model

```{r}

set.seed(904)
rf_p_resamples <- 
  tune_grid(
  rf_p_wf,
  resamples = crime_folds,
  grid = rf_p_grid)

collect_metrics(rf_p_resamples)

show_best(rf_p_resamples, 
          metric = "rmse") 

select_best(rf_p_resamples)
autoplot(rf_p_resamples)

```

### predictions

```{r}

set.seed(904)
rf_p_pre_wf <- 
  rf_p_wf |> 
  finalize_workflow(select_best(rf_p_resamples))

rf_p_pre_wf <-
  fit(rf_p_pre_wf, 
      data = crime_train)

rf_p_predictions <- 
  predict(object = rf_p_pre_wf, 
          new_data = crime_test)$.pred

actual_p <-
  crime_test$log_property

county <- 
  crime_test$county

rf_p_predictions_df<-
  tibble(county,actual_p,rf_p_predictions)

rmse_rf_p <- 
  rmse(rf_p_predictions_df, actual_p, rf_p_predictions)

```

### selecting the top 15 important variables

```{r}

rf_p_final_wf <- 
  rf_p_wf |> 
  finalize_workflow(select_best(rf_p_resamples))

rf_p_fit <- 
  rf_p_final_wf |>
  fit(data = crime_all_2010_cleaned)

set.seed(904)
rf_p_fit|>
  extract_fit_parsnip()|>
  vip(num_features = 15) %>%
  .$data |>
  mutate(
    Importance = Importance / max(Importance),
    Variable = fct_reorder(.x = Importance, .f = Variable)
  ) |>
  ggplot(aes(Importance, Variable)) +
  scale_y_discrete(labels = crime_labels) +
  labs(title = "Top 15 Variables for Property Crime Prediction") +
  geom_col()

```

Based on the result from random forest model, the top 15 important variables for predicting property crime are shown above. Most of the variables are from the socioeconomic data and few are from the land cover data.

## MARS

### feature and target engineering

```{r}

# recipe

mars_p_rec <- 
  recipe(log_property ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_viol,
              new_role = "id") |>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())

mars_p_rec_sum <- summary(mars_p_rec)

# model

mars_p_mod <- 
  mars(
  num_terms = tune(), 
  prod_degree = tune()
) |>
  set_mode(mode = "regression") |>
  set_engine(engine = "earth")

# tuning grid

mars_p_grid <- 
  grid_regular(
  num_terms(range = c(10, 100)),
  prod_degree(),
  levels = 10
)

# workflow

mars_p_wf <- 
  workflow() |>
  add_recipe(recipe = mars_p_rec) |>
  add_model(spec = mars_p_mod)

```

### choose the best model

```{r}

# resample   

set.seed(904)
mars_p_resamples <- 
  mars_p_wf|>
  tune_grid(
  resamples = crime_folds,
  grid = mars_p_grid,
  metrics = metric_set(rmse)
  )

collect_metrics(mars_p_resamples)

top_mars_p_models <-
  mars_p_resamples |>
  show_best()

mars_p_resamples|>
  select_best()

```

### prediction

```{r}

set.seed(904)
mars_p_pre_wf <- 
  mars_p_wf |> 
  finalize_workflow(select_best(mars_p_resamples))

mars_p_pre_wf <-
  fit(mars_p_pre_wf, 
      data = crime_train)

mars_p_predictions <- 
  predict(object = mars_p_pre_wf, 
          new_data = crime_test)$.pred

actual_p <- 
  crime_test$log_property

county <- 
  crime_test$county

mars_p_predictions_df <-
  tibble(county,actual_p,mars_p_predictions)

rmse_mars_p <- 
  rmse(mars_p_predictions_df, actual_p, mars_p_predictions)

```

## XG Boost

### feature and target engineering

```{r}

xgb_p_rec <-
  recipe(log_property ~ .,
         data = crime_train) |>
  update_role(county,
              year,
              viol,
              viol_rate,
              property_rate,
              property,
              log_viol,
              new_role = "id")|>
  step_impute_median(all_predictors())|>
  step_nzv(all_predictors())

all_cores <- 
  parallel::detectCores(logical = FALSE)

registerDoFuture()

cl <- 
  makeCluster(all_cores - 3L)

plan(cluster, workers = cl)

xgb_p_spec <- 
  boost_tree(
  trees = 200,
  tree_depth = tune(), 
  min_n = tune(),
  loss_reduction = tune(),                     
  sample_size = tune(), 
  mtry = tune(),         
  learn_rate = tune()                   
) |>
  set_engine("xgboost") |>
  set_mode("regression")


set.seed(904)
xgb_p_grid <- 
  grid_latin_hypercube(
  tree_depth(),
  min_n(range = c(2, 10)),
  loss_reduction(range = c(-5, -3)),
  sample_size = sample_prop(),
  mtry(range = c(1, 10)),
  learn_rate(range = c(-5, -1)),
  size = 30
)

xgb_p_wf <- workflow() |>
  add_recipe(xgb_p_rec) |>
  add_model(xgb_p_spec)

```

### choose the best model

```{r}

set.seed(904)
xgb_p_resamples <- 
  tune_grid(
  xgb_p_wf,
  resamples = crime_folds,
  grid = xgb_p_grid)

collect_metrics(xgb_p_resamples)
show_best(xgb_p_resamples, metric = "rmse") 
select_best(xgb_p_resamples, metric= "rmse")

```

### prediction

```{r}

set.seed(904)
xgb_p_pre_wf <- 
  xgb_p_wf |> 
  finalize_workflow(select_best(xgb_p_resamples))

xgb_p_pre_wf <-
  fit(xgb_p_pre_wf, 
      data = crime_train)

xgb_p_predictions <- 
  predict(object = xgb_p_pre_wf, 
          new_data = crime_test)$.pred

actual_p <- 
  crime_test$log_property

county <- 
  crime_test$county

xgb_p_predictions_df <-
  tibble(county,actual_p,xgb_p_predictions)

rmse_xgb_p <- 
  rmse(xgb_p_predictions_df, actual_p, xgb_p_predictions)

```

## Final comparison of rmse

```{r}

property_comparison <-
  bind_rows(rmse_ridge_p, rmse_lasso_p, rmse_rf_p, rmse_mars_p, rmse_xgb_p)

property_comparison|>
  mutate(
    model = c("Ridge Regression","Lasso Regression", 
              "Random Forest", "MARS", "XG Boost"))|>
  select(model, everything())

```

Because random forest has the lowest RMSE value, we will use random forest to predict property crime rate in 2016.

## Prediction on implementation property crime data in 2016

```{r}

set.seed(904)
rf_p_pre_wf_2016 <- 
  rf_p_wf |> 
  finalize_workflow(select_best(rf_p_resamples))

rf_p_pre_wf_2016 <-
  fit(rf_p_pre_wf, 
      data = crime_all_2010_cleaned)

rf_p_predictions_2016 <- 
  predict(object = rf_p_pre_wf_2016, 
          new_data = crime_all_2016_cleaned)$.pred

actual_implement_p <- crime_all_2016_cleaned$log_property

county_implement <- crime_all_2016_cleaned$county

rf_p_predictions_2016_df <-
  tibble(county_implement,actual_implement_p,rf_p_predictions_2016)

rmse_rf_p_2016 <- 
  rmse(rf_p_predictions_2016_df, actual_implement_p, rf_p_predictions_2016)

rmse_rf_p_2016

```

Final prediction for property crime rate in 2016 has a RMSE of 0.459.