---
title: "Part 2"
author: "Eric Langowski"
output: 
  html_document:
    code_folding: hide
    df_print: paged
    theme: sandstone
    toc: true
    toc_float: true
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE)
library(tidyverse)
library(lubridate)
library(tidymodels)
theme_set(theme_bw())

if(file.exists("database/cook.sqlite")){
  f <- "database/cook.sqlite"
} else {
  f <- "static/solutions/database/cook.sqlite"
}


make_table <- function(inputdata){
    DT::datatable(
      inputdata %>% 
        mutate(across(where(is.numeric), ~round(.x, 3))), 
      extensions = 'Buttons', 
      rownames = FALSE,
      options = list(
        dom = 'Bfrtip', 
        buttons = c('excel', 'pdf'),
        pageLength = 20)
  )
}


con <- DBI::dbConnect(RSQLite::SQLite(), f)

assessments <- dplyr::tbl(con, 'assessments') %>% collect()
sales <- dplyr::tbl(con, 'sales') %>% collect() %>% mutate(sale_date = as_datetime(sale_date))
characteristics <- dplyr::tbl(con, 'characteristics') %>% collect()
geo_uni <- dplyr::tbl(con, 'geospatial_universe') %>% collect()

joined <- sales %>%
  mutate(year = year(ymd(sale_date))) %>%
  filter(sale_filter_less_than_10k == 0,
         sale_filter_deed_type == 0,
         is_multisale == 0) %>%
  left_join(assessments %>% select(-township_code, -class), 
            join_by('pin', 'year'=='tax_year'))

ratios <- cmfproperty::reformat_data(
  joined %>% mutate(assessed_value = if_else(year == 2023, 10*certified_tot, 10*board_tot)),
  'sale_price',
  'assessed_value',
  'year'
)

min_yr <- 2021
max_yr <- 2023


stats <- cmfproperty::calc_iaao_stats(ratios)

output <- cmfproperty::diagnostic_plots(ratios = ratios,
                                        stats = stats, 
                                        min_reporting_yr = min_yr,
                                        max_reporting_yr = max_yr)
```

# Introduction

Cook County, Illinois includes almost 2 million parcels ranging from suburban mansions to thousand feet skyscrapers. The system is administered by a series of countywide agencies and it is responsible for fairly dividing property tax levies (the amounts that taxing bodies such as school districts levy against properties in their districts) among all these properties.

In the past, assessments in Cook County were highly regressive, where lower valued properties paid a significantly higher tax rate than higher valued properties. A study by the Center for Municipal Finance at the University of Chicago showed that from 2011 to 2015 that these regressive assessments lead to over $2 billion of property tax liability being paid by lower valued homes instead of higher valued ones.

This study, alongside a public reckoning, resulted in the election of a new assessor in 2018.

This analysis focuses on single family homes Additionally, using the `cmfproperty` package, the IAAO arm's length standard was applied to the data. This will present a conservative picture of the housing market in Cook County. Note that homes in Cook County as supposed to be assessed at 10% of their fair market value. For clarity, assessed values here are converted to market values.

The sales ratio is a key unit for analyzing the accuracy of assessments. It is calculated by taking the ratio of a property's sale price to its assessment at the time of sale. The sales ratios can be evaluated using metrics from the International Association of Assessing Officers.


```{r}
iaao <- cmfproperty::iaao_graphs(stats=stats, ratios=ratios, min_reporting_yr = min_yr, max_reporting_yr = max_yr, jurisdiction_name = 'Cook County')
```

`r iaao[[1]]`

```{r}
iaao[[2]]
```

`r iaao[[3]]`


```{r}
iaao[[4]]
```

`r iaao[[5]]`

```{r}
iaao[[6]]
```

```{r}
bs <- cmfproperty::binned_scatter(ratios = ratios, min_reporting_yr = min_yr, max_reporting_yr = max_yr, jurisdiction_name = 'Cook County')
```

`r bs[[1]]`

```{r}
bs[[2]]
```

# Classification Model of Overassessment (Part B)

```{r}
ratios_chars <- 
  ratios %>%
  left_join(characteristics %>% select(-township_code, -class) %>% distinct(), by=c('pin')) %>%
  mutate(year = as.character(TAX_YEAR)) 
```

## Constructing a Classifier

Overassessment can have many reasonable definitions. Strictly, residential properties are supposed to be assessed at 100% of their value so properties which had sales ratios greater than 100% are technically overassessed. An alternate definition might be a property assessed at some rate greater than the average rate of assessment (e.g. higher than average). For this initial version, let's set the alternate definition as 20% higher than the median rate of assessment.

```{r}
part_b_data <- ratios_chars %>% filter(SALE_YEAR < 2023, !is.na(char_yrblt)) %>%
  mutate(strict_def = if_else(RATIO > 1, 'yes', 'no'),
         alternate_def = if_else(RATIO > 1.2 * median(RATIO), 'yes', 'no'))

part_b_data %>%
  count(strict_def, alternate_def) %>%
  make_table()

ggplot(part_b_data, aes(x=RATIO)) +
  geom_density() +
  geom_vline(aes(xintercept = 1), linetype='dashed') +
  geom_vline(aes(xintercept = 1.2 * median(RATIO)), linetype='dotted') +
  labs(title='Distribution of Ratios and Classifier Boundaries', 
       caption='Dashed line is strict defintion, dotted line is alternate defintion')
```

Note that the distribution of ratios is somewhat continuous---meaning that this will be a difficult classification problem! For example, a property which has a ratio of 0.99 and one which has a ratio of 1.01 will be classified differently even though there is not a significant difference in those rates of assessments.

## Initial Model Specification

`class ~ ASSESSED_VALUE + SALE_YEAR + char_yrblt + char_bldg_sf + char_land_sf + township_name`


```{r, include=FALSE}
cor_output <- corrr::correlate(part_b_data %>% 
                   select(ASSESSED_VALUE,
                          SALE_YEAR,
                          char_yrblt,
                          char_bldg_sf,
                          char_land_sf)) %>%
  corrr::fashion()
```


```{r}
cor_output %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
```

### recipe

```{}
recipe(strict_def ~ ASSESSED_VALUE + SALE_YEAR + 
                         char_yrblt + char_bldg_sf + 
                         char_land_sf + township_name, 
                       data=part_b_data) %>%
  step_log(ASSESSED_VALUE) %>%
  step_interact(~c(ASSESSED_VALUE, SALE_YEAR)) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  step_corr(all_predictors()) %>%
  prep()
```

Key steps highlighted:

- Logging assessed value, ensures we capture variability in assessed value.

## Workflow Prep

Our `tidymodels` workflow requires a model, formula, and recipe.

```{r}
folds <- rsample::vfold_cv(part_b_data, v = 5, strata = strict_def)

class_model <-
  rand_forest(trees = 100) %>%
  set_mode('classification') %>%
  set_engine('ranger')

class_recipe <- recipe(strict_def ~ ASSESSED_VALUE + SALE_YEAR + 
                         char_yrblt + char_bldg_sf + 
                         char_land_sf + township_name, 
                       data=part_b_data) %>%
  step_log(ASSESSED_VALUE) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
  step_interact(~char_bldg_sf:char_land_sf) %>%
  step_corr(all_predictors())

class_recipe %>% prep() %>% summary() %>% slice_sample(n=10) %>% make_table()

bake(class_recipe %>% prep(),
     part_b_data) %>% slice_sample(n=10) %>% make_table()
```

For this initial specification, we utilize 5-fold cross validation to create our training/testing data and a random forest set to have 100 trees.

## Specification 1, Strict Class

```{r}
first_workflow <-
  workflow() %>%
  add_model(class_model) %>%
  add_recipe(class_recipe)

model_fit <- first_workflow %>%
  fit_resamples(folds, control=control_resamples(save_pred=TRUE, verbose=TRUE))
```


```{r}
#collect_metrics(model_fit, summarize=FALSE)
our_results <- collect_metrics(model_fit, summarize=TRUE)
our_results %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
#collect_predictions(model_fit)
our_preds <- collect_predictions(model_fit, summarize=TRUE) 
our_preds %>%
  count(.pred_class) %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
conf_mat(our_preds, estimate=.pred_class, truth=strict_def)
```

Our model has pretty mediocre accuracy of `r str_glue("{our_results %>% filter(.metric=='accuracy') %>% pull(mean) %>% round(3)}")`.

## Classifier Accuracy Metrics

```{r}
our_preds <- collect_predictions(model_fit, summarize=TRUE) 
multi_metric <- metric_set(recall, specificity, precision, accuracy, f_meas)
multi_metric(our_preds, truth=strict_def, estimate=.pred_class) %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
```

Some initial views on our model:

```{r}
part_b_data %>%
  mutate(pred = our_preds$.pred_class,
         bin = ntile(SALE_PRICE, 10)) %>%
  group_by(bin) %>%
  summarize(avg_sp = dollar(mean(SALE_PRICE)),
            share_correct = percent(sum(strict_def == pred) / n()),
            share_over = percent(sum(strict_def == 'yes') / n()))  %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
```

```{r}
roc_curve(our_preds, strict_def, .pred_no) %>%
  autoplot()
```

## Discussion of Thresholds

```{r}
our_preds %>% 
  ggplot(aes(x=.pred_yes, fill=strict_def)) + 
  geom_density(alpha=0.5) + 
  geom_vline(xintercept = 0.5, linetype='dashed') +
  labs(x='Probability of Yes', y='Density', title='Probabilty Density by Actual Classification')
```

We can see from the above plot a little bit about why this model predicts so few yes's since probabilities of 0.5 or greater are needed based on the default threshold. We can see that, on average, our model gives true yes's a higher probability of yes than the no's.


```{r}
our_preds %>% 
  ggplot(aes(x=.pred_yes, fill=strict_def)) + 
  geom_histogram(alpha=0.5) + 
  geom_vline(xintercept = 0.5, linetype='dashed') +
  labs(x='Probability of Yes', y='Actual Counts', title='Counts by Actual Classification')
```

This is somewhat now the classic problem. We can increase the number of overassessed cases predicted correctly, but that will also increase incorrect predictions of "no" overassessed predicted as "yes."

So functionally, to "tune" the threshold, we need to recalculate our model's performance at a variety of thresholds. Fortunately, the `probably` package allows us to do this easily.

A way to measure this tradeoff is called the `j_index` which is calculated as adding the sensitivity (true positive rate) and specificity (true negative rate) and subtracting 1. This metric is 1 when there are no false positive or no false negatives.

```{r}
j_index(our_preds, strict_def, .pred_class)
```

Let's vary across many thresholds replicating [this analysis](https://probably.tidymodels.org/articles/where-to-use.html)

```{r}
threshold_data <- our_preds %>%
  probably::threshold_perf(strict_def, .pred_no,
                 thresholds=seq(0.5, 1, by=0.0025))
```

```{r}
ggplot(threshold_data,
       aes(x=.threshold, y=.estimate, color=.metric)) +
  geom_line() +
  geom_vline(xintercept = threshold_data %>%
  filter(.metric == "j_index") %>%
  filter(.estimate == max(.estimate)) %>%
  pull(.threshold)) +
  labs(x="'No' Threshold\n(probs above threshold are no)", y='Metric Estimate', title="Threshold to maximize J-Index")
```

The optimal threshold for maximizing j-index is a threshold of 0.875. This means that cases with a predicted probability of greater than 0.875 are classified as "no" and those less than are classified as "yes", or in other words a probability of 0.125 of "yes" for overassessment is sufficient to predict overassessment.

```{r}
hard_pred_0.875 <- our_preds %>%
  mutate(
    .pred = probably::make_two_class_pred(
      estimate = .pred_no,
      levels = levels(strict_def),
      threshold = 0.875
    )
  ) %>%
  select(strict_def, contains(".pred"))

conf_mat(hard_pred_0.875, estimate=.pred, truth=strict_def)
```

We now predict many more cases of overassessment and functionally we have optimized so that we predicted currently true positives and true negatives at a rate of about 63%.

This doesn't necessarily improve our model but is an instructive way to think about how to convert classification model probabilities to hard class predictions.

# Regression Model of Market Valuations (Part C)

Generating our own assessments (for 2023) is very similar to modeling overassessment. Initially, I will use the same recipe and formula except replacing class with sale price.  I will also demonstrate the xgboost (boosted decision tree) package. Boosted decision trees are similar to random forests except each tree is created iteratively in a process of continuous improvement.

Predictions will then be made for 2023 compared to the baseline of actual assessed values. Workflow is the same except that xgboost requires us to bake our data first.

```{r}
part_c_data <- ratios_chars %>% 
  filter(SALE_YEAR < 2023, !is.na(char_yrblt)) %>%
  mutate(sp_log = log(SALE_PRICE_ADJ))

split <- rsample::initial_split(part_c_data)
train <- training(split) 
test <- testing(split)

reg_model <- boost_tree(trees=200) %>%
  set_mode('regression') %>%
  set_engine('xgboost')

reg_recipe <- recipe(sp_log ~ ASSESSED_VALUE + SALE_YEAR + 
                         char_yrblt + char_bldg_sf + 
                         char_land_sf + township_name, 
                       data=train) %>%
  step_log(ASSESSED_VALUE) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_corr(all_predictors()) %>% prep()

reg_recipe %>% bake(train) %>% glimpse()

reg_workflow <-
  workflow() %>%
  add_model(reg_model) %>%
  add_recipe(reg_recipe)
```

## Model Evaluation

```{r}
model_fit_reg <- reg_workflow %>%
  fit(train)

our_preds <- model_fit_reg %>% augment(test)

multi_metric <- metric_set(mape, rmse, rsq)
multi_metric(our_preds, truth=sp_log, estimate=.pred) %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
```

Actual assessments:

```{r}
multi_metric(our_preds %>% mutate(av_log = log(ASSESSED_VALUE)), truth=sp_log, estimate=av_log) %>%
  kableExtra::kable() %>%
  kableExtra::kable_material(c("striped", "hover"))
```

```{r}
ratios <- cmfproperty::reformat_data(
  our_preds %>% mutate(av_pred = exp(.pred)),
  'SALE_PRICE',
  'av_pred',
  'SALE_YEAR'
)

bs <- cmfproperty::binned_scatter(ratios = ratios, min_reporting_yr = 2021, max_reporting_yr = 2022,
                                  jurisdiction_name = 'Cook County')
```

`r bs[[1]]`

```{r}
bs[[2]]
```