Part 2
Eric Langowski
1 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.
iaao <- cmfproperty::iaao_graphs(stats=stats, ratios=ratios, min_reporting_yr = min_yr, max_reporting_yr = max_yr, jurisdiction_name = 'Cook County')For 2023, the COD in Cook County was 19.24 which did not meet the IAAO standard for uniformity.
iaao[[2]]
In 2023, the PRD in Cook County, was 1.053 which does not meet the IAAO standard for vertical equity.
iaao[[4]]
In 2023, the PRB in Cook County was 0.004 which indicates that sales ratios increase by 0.4% when home values double. This meets the IAAO standard.
iaao[[6]]
bs <- cmfproperty::binned_scatter(ratios = ratios, min_reporting_yr = min_yr, max_reporting_yr = max_yr, jurisdiction_name = 'Cook County')In 2023, the most expensive homes (the top decile) were assessed at 77.1% of their value and the least expensive homes (the bottom decile) were assessed at 91.8%. In other words, the least expensive homes were assessed at 1.19 times the rate applied to the most expensive homes. Across our sample from 2021 to 2023, the most expensive homes were assessed at 77.2% of their value and the least expensive homes were assessed at 82.6%, which is 1.07 times the rate applied to the most expensive homes.
bs[[2]]
2 Classification Model of Overassessment (Part B)
ratios_chars <-
ratios %>%
left_join(characteristics %>% select(-township_code, -class) %>% distinct(), by=c('pin')) %>%
mutate(year = as.character(TAX_YEAR)) 2.1 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.
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.
2.2 Initial Model Specification
class ~ ASSESSED_VALUE + SALE_YEAR + char_yrblt + char_bldg_sf + char_land_sf + township_name
cor_output %>%
kableExtra::kable() %>%
kableExtra::kable_material(c("striped", "hover"))| term | ASSESSED_VALUE | SALE_YEAR | char_yrblt | char_bldg_sf | char_land_sf |
|---|---|---|---|---|---|
| ASSESSED_VALUE | .03 | .19 | .70 | .17 | |
| SALE_YEAR | .03 | -.01 | -.02 | .00 | |
| char_yrblt | .19 | -.01 | .33 | .14 | |
| char_bldg_sf | .70 | -.02 | .33 | .31 | |
| char_land_sf | .17 | .00 | .14 | .31 |
2.2.1 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.
2.3 Workflow Prep
Our tidymodels workflow requires a model, formula, and
recipe.
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.
2.4 Specification 1, Strict Class
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))## i Fold1: preprocessor 1/1## ✓ Fold1: preprocessor 1/1## i Fold1: preprocessor 1/1, model 1/1## ✓ Fold1: preprocessor 1/1, model 1/1## i Fold1: preprocessor 1/1, model 1/1 (extracts)## i Fold1: preprocessor 1/1, model 1/1 (predictions)## i Fold2: preprocessor 1/1## ✓ Fold2: preprocessor 1/1## i Fold2: preprocessor 1/1, model 1/1## ✓ Fold2: preprocessor 1/1, model 1/1## i Fold2: preprocessor 1/1, model 1/1 (extracts)## i Fold2: preprocessor 1/1, model 1/1 (predictions)## i Fold3: preprocessor 1/1## ✓ Fold3: preprocessor 1/1## i Fold3: preprocessor 1/1, model 1/1## ✓ Fold3: preprocessor 1/1, model 1/1## i Fold3: preprocessor 1/1, model 1/1 (extracts)## i Fold3: preprocessor 1/1, model 1/1 (predictions)## i Fold4: preprocessor 1/1## ✓ Fold4: preprocessor 1/1## i Fold4: preprocessor 1/1, model 1/1## ✓ Fold4: preprocessor 1/1, model 1/1## i Fold4: preprocessor 1/1, model 1/1 (extracts)## i Fold4: preprocessor 1/1, model 1/1 (predictions)## i Fold5: preprocessor 1/1## ✓ Fold5: preprocessor 1/1## i Fold5: preprocessor 1/1, model 1/1## ✓ Fold5: preprocessor 1/1, model 1/1## i Fold5: preprocessor 1/1, model 1/1 (extracts)## i Fold5: preprocessor 1/1, model 1/1 (predictions)#collect_metrics(model_fit, summarize=FALSE)
our_results <- collect_metrics(model_fit, summarize=TRUE)
our_results %>%
kableExtra::kable() %>%
kableExtra::kable_material(c("striped", "hover"))| .metric | .estimator | mean | n | std_err | .config |
|---|---|---|---|---|---|
| accuracy | binary | 0.8796092 | 5 | 0.0000774 | Preprocessor1_Model1 |
| roc_auc | binary | 0.6697350 | 5 | 0.0028675 | Preprocessor1_Model1 |
#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"))| .pred_class | n |
|---|---|
| no | 62121 |
| yes | 10 |
conf_mat(our_preds, estimate=.pred_class, truth=strict_def)## Truth
## Prediction no yes
## no 54646 7475
## yes 5 5Our model has pretty mediocre accuracy of 0.88.
2.5 Classifier Accuracy Metrics
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"))| .metric | .estimator | .estimate |
|---|---|---|
| recall | binary | 0.9999085 |
| specificity | binary | 0.0006684 |
| precision | binary | 0.8796703 |
| accuracy | binary | 0.8796092 |
| f_meas | binary | 0.9359435 |
Some initial views on our model:
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"))| bin | avg_sp | share_correct | share_over |
|---|---|---|---|
| 1 | $100,157 | 73% | 27% |
| 2 | $176,665 | 83% | 17% |
| 3 | $226,096 | 89% | 11% |
| 4 | $262,663 | 89% | 11% |
| 5 | $298,145 | 89% | 11% |
| 6 | $337,854 | 90% | 10% |
| 7 | $386,758 | 90% | 10% |
| 8 | $461,476 | 92% | 8% |
| 9 | $620,364 | 92% | 8% |
| 10 | $1,312,556 | 91% | 9% |
roc_curve(our_preds, strict_def, .pred_no) %>%
autoplot()
2.6 Discussion of Thresholds
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.
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')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
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.
j_index(our_preds, strict_def, .pred_class)Let’s vary across many thresholds replicating this analysis
threshold_data <- our_preds %>%
probably::threshold_perf(strict_def, .pred_no,
thresholds=seq(0.5, 1, by=0.0025))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.
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)## Truth
## Prediction no yes
## no 34242 2814
## yes 20409 4666We 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.
3 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.
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()## Rows: 46,598
## Columns: 43
## $ ASSESSED_VALUE <dbl> 12.11220, 12.50651, 12.10260, 14.24078, 12…
## $ SALE_YEAR <dbl> 2021, 2021, 2022, 2022, 2022, 2022, 2022, …
## $ char_yrblt <dbl> 1951, 1957, 1967, 1909, 1973, 1899, 1951, …
## $ char_bldg_sf <dbl> 864, 1308, 1154, 4605, 1339, 2084, 1410, 1…
## $ char_land_sf <dbl> 7100, 11340, 13200, 15900, 8750, 3000, 657…
## $ sp_log <dbl> 12.40901, 12.76569, 12.18587, 14.82274, 12…
## $ township_name_Berwyn <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Bloom <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ township_name_Bremen <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Calumet <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Cicero <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Elk.Grove <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Evanston <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Hanover <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, …
## $ township_name_Hyde.Park <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Jefferson <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, …
## $ township_name_Lake <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
## $ township_name_Lake.View <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Lemont <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, …
## $ township_name_Leyden <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Lyons <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Maine <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_New.Trier <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Niles <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_North.Chicago <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Northfield <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Norwood.Park <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Oak.Park <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Orland <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Palatine <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Palos <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Proviso <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Rich <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_River.Forest <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Riverside <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Rogers.Park <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Schaumburg <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_South.Chicago <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Stickney <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_Thornton <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ township_name_West.Chicago <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, …
## $ township_name_Wheeling <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, …
## $ township_name_Worth <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …reg_workflow <-
workflow() %>%
add_model(reg_model) %>%
add_recipe(reg_recipe)3.1 Model Evaluation
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"))| .metric | .estimator | .estimate |
|---|---|---|
| mape | standard | 1.5191941 |
| rmse | standard | 0.2573236 |
| rsq | standard | 0.8559879 |
Actual assessments:
multi_metric(our_preds %>% mutate(av_log = log(ASSESSED_VALUE)), truth=sp_log, estimate=av_log) %>%
kableExtra::kable() %>%
kableExtra::kable_material(c("striped", "hover"))| .metric | .estimator | .estimate |
|---|---|---|
| mape | standard | 2.7184020 |
| rmse | standard | 0.4380267 |
| rsq | standard | 0.8400065 |
ratios <- cmfproperty::reformat_data(
our_preds %>% mutate(av_pred = exp(.pred)),
'SALE_PRICE',
'av_pred',
'SALE_YEAR'
)## [1] "Renaming already present column 'ASSESSED_VALUE' to 'ASSESSED_VALUE_2'."
## [1] "Filtered out non-arm's length transactions"
## [1] "Inflation adjusted to 2021"bs <- cmfproperty::binned_scatter(ratios = ratios, min_reporting_yr = 2021, max_reporting_yr = 2022,
jurisdiction_name = 'Cook County')In 2022, the most expensive homes (the top decile) were assessed at 95.3% of their value and the least expensive homes (the bottom decile) were assessed at 111.4%. In other words, the least expensive homes were assessed at 1.17 times the rate applied to the most expensive homes. Across our sample from 2021 to 2022, the most expensive homes were assessed at 95.5% of their value and the least expensive homes were assessed at 111.3%, which is 1.17 times the rate applied to the most expensive homes.
bs[[2]]