Update interpret pipeline stage to optionally generate leaf node comps for all properties

.gitignore

@@ -26,3 +26,6 @@ cache/
 
 # Ignore scratch documents
 scratch*.*
+
+# Python files
+__pycache__

Dockerfile

@@ -14,9 +14,9 @@ ENV RENV_PATHS_CACHE /setup/cache
 RUN apt-get update && \
     apt-get install --no-install-recommends -y \
         libcurl4-openssl-dev libssl-dev libxml2-dev libgit2-dev git \
-        libudunits2-dev python3-dev python3-pip libgdal-dev libgeos-dev \
-        libproj-dev libfontconfig1-dev libharfbuzz-dev libfribidi-dev pandoc \
-        curl gdebi-core && \
+        libudunits2-dev python3-dev python3-pip python3-venv libgdal-dev \
+        libgeos-dev libproj-dev libfontconfig1-dev libharfbuzz-dev \
+        libfribidi-dev pandoc curl gdebi-core && \
     rm -rf /var/lib/apt/lists/*
 
 # Install Quarto
@@ -28,7 +28,7 @@ RUN gdebi -n quarto-linux-amd64.deb
 RUN pip install --no-cache-dir dvc[s3]
 
 # Copy R bootstrap files into the image
-COPY renv.lock .Rprofile DESCRIPTION ./
+COPY renv.lock .Rprofile DESCRIPTION requirements.txt ./
 COPY renv/profiles/reporting/renv.lock reporting-renv.lock
 COPY renv/ renv/
 

R/helpers.R

@@ -182,6 +182,24 @@ extract_num_iterations <- function(x) {
   length(evals)
 }
 
+# Extract weights for model features based on tree importance. Assumes that
+# the model was trained with the `valids` parameter set such that error metrics
+# are saved for each tree on the model$record_evals attribute. The output
+# weights are useful for computing comps using leaf node assignments
+extract_weights <- function(model, init_score, metric = "rmse") {
+  # Index into the errors list, and un-list so it is a flat/1dim list
+  record_evals <- model$record_evals
+  errors <- unlist(record_evals$tree_errors[[metric]]$eval)
+  errors <- c(init_score, errors)
+  diff_in_errors <- diff(errors, 1, 1)
+
+  # Take proportion of diff in errors over total diff in
+  # errors from all trees
+  weights <- diff_in_errors / sum(diff_in_errors)
+
+  return(weights)
+}
+
 # Given the result of a CV search, get the number of iterations from the
 # result set with the best performing hyperparameters
 select_iterations <- function(tune_results, metric, type = "mean") {

R/setup.R

@@ -72,6 +72,10 @@ shap_enable <- as.logical(Sys.getenv(
   "SHAP_ENABLE_OVERRIDE",
   unset = get(params_obj_name)$toggle$shap_enable
 ))
+comp_enable <- as.logical(Sys.getenv(
+  "COMP_ENABLE_OVERRIDE",
+  unset = get(params_obj_name)$toggle$comp_enable
+))
 upload_enable <- as.logical(Sys.getenv(
   "UPLOAD_ENABLE_OVERRIDE",
   unset = get(params_obj_name)$toggle$upload_enable

dvc.yaml

@@ -89,7 +89,7 @@ stages:
     - output/test_card/model_test_card.parquet
     - output/assessment_pin/model_assessment_pin.parquet
     params:
-    - assessment.data_year
+    - assessment
     - ratio_study
     outs:
     - output/performance/model_performance_test.parquet:
@@ -110,9 +110,13 @@ stages:
       importance metrics for each feature
     deps:
     - input/assessment_data.parquet
+    - input/training_data.parquet
+    - output/assessment_card/model_assessment_card.parquet
     - output/workflow/fit/model_workflow_fit.zip
     - output/workflow/recipe/model_workflow_recipe.rds
     params:
+    - toggle.shap_enable
+    - toggle.comp_enable
     - model.predictor.all
     outs:
     - output/shap/model_shap.parquet:
@@ -121,6 +125,8 @@ stages:
         cache: false
     - output/intermediate/timing/model_timing_interpret.parquet:
         cache: false
+    - output/comp/model_comp.parquet:
+        cache: false
 
   finalize:
     cmd: Rscript pipeline/05-finalize.R
@@ -171,6 +177,7 @@ stages:
     - output/performance/model_performance_assessment.parquet
     - output/performance_quantile/model_performance_quantile_assessment.parquet
     - output/shap/model_shap.parquet
+    - output/comp/model_comp.parquet
     - output/feature_importance/model_feature_importance.parquet
     - output/metadata/model_metadata.parquet
     - output/timing/model_timing.parquet

misc/file_dict.csv

@@ -23,8 +23,9 @@ output,performance_quantile_test_linear,3,evaluate,ccao-model-results-us-east-1,
 output,performance_quantile_assessment,3,evaluate,ccao-model-results-us-east-1,output/performance_quantile/model_performance_quantile_assessment.parquet,performance_quantile/year={year}/stage=assessment/{run_id}.parquet,performance_quantile,geography [by class] by quantile,"year, run_id, stage, geography_type, geography_id, by_class, class, quantile",Performance metrics by quantile within class and geography,Assessment set uses the prior year sales to compare to the assessed value
 output,shap,4,interpret,ccao-model-results-us-east-1,output/shap/model_shap.parquet,shap/,shap,card,"year, run_id, township_code, meta_pin, meta_card_num",SHAP values for each feature for each card in the assessment data,NOTE: Each run adds new partitions to S3 which must be added via a Glue crawler
 output,feature_importance,4,interpret,ccao-model-results-us-east-1,output/feature_importance/model_feature_importance.parquet,feature_importance/year={year}/{run_id}.parquet,feature_importance,predictor,"year, run_id, model_predictor_all_name","Feature importance values (gain, cover, and frequency) for the run",
+output,comp,4,interpret,ccao-model-results-us-east-1,output/comp/model_comp.parquet,comp/,comp,card,"year, run_id, meta_pin, meta_card_num",Comparables for each card (computed using leaf node assignments),
 output,report_performance,5,finalize,ccao-model-results-us-east-1,reports/performance/performance.html,report/year={year}/report_type=performance/{run_id}.html,,model run,,Rendered Quarto doc with model performance statistics,
 output,report_pin,5,finalize,ccao-model-results-us-east-1,reports/pin/,report/year={year}/report_type=pin/run_id={run_id}/,,model run,,Rendered Quarto doc for individual PINs,
 output,metadata,5,finalize,ccao-model-results-us-east-1,output/metadata/model_metadata.parquet,metadata/year={year}/{run_id}.parquet,metadata,model run,"year, run_id","Information about each run, including parameters, run ID, git info, etc.",
 intermediate,timing,,all,,output/intermediate/timing/,,,model stage,"year, msg",Parquet files for each stage containing the stage time elapsed,Converted into a one-row data frame in the finalize stage
-output,timing,,all,ccao-model-results-us-east-1,output/timing/model_timing.parquet,timing/year={year}/{run_id}.parquet,timing,model run,"year, run_id",Finalized time elapsed for each stage of the run,"Each row represents one run, while columns represent the stages"
+output,timing,,all,ccao-model-results-us-east-1,output/timing/model_timing.parquet,timing/year={year}/{run_id}.parquet,timing,model run,"year, run_id",Finalized time elapsed for each stage of the run,"Each row represents one run, while columns represent the stages"

params.yaml

@@ -27,6 +27,9 @@ toggle:
   # desirable to save time when testing many models
   shap_enable: FALSE
 
+  # Should comps be calculated for this run in the interpret stage?
+  comp_enable: FALSE
+
   # Upload all modeling artifacts and results to S3 in the upload stage. Set
   # to FALSE if you are not a CCAO employee
   upload_enable: TRUE

pipeline/01-train.R

@@ -145,7 +145,9 @@ lgbm_model <- parsnip::boost_tree(
     # using floor(log2(num_leaves)) + add_to_linked_depth. Useful since
     # otherwise Bayesian opt spends time exploring irrelevant parameter space
     link_max_depth = params$model$parameter$link_max_depth,
-
+    # Always initialize save_tree_error to false until we're ready to train
+    # the final model, since it's incompatible with CV
+    save_tree_error = FALSE,
 
     ### 4.1.2. Tuned Parameters ------------------------------------------------
 
@@ -346,12 +348,14 @@ if (cv_enable) {
 
 # Finalize the model specification by disabling early stopping, instead using
 # the maximum number of iterations used during the best cross-validation round
-# OR the default `num_iterations` if CV was not performed
+# OR the default `num_iterations` if CV was not performed. Also enable comps
+# if they're configured to run, since they're incompatible with CV
 lgbm_model_final <- lgbm_model %>%
   set_args(
     stop_iter = NULL,
     validation = 0,
-    trees = lgbm_final_params$num_iterations
+    trees = lgbm_final_params$num_iterations,
+    save_tree_error = comp_enable
   )
 
 # Fit the final model using the training data and our final hyperparameters

pipeline/04-interpret.R

@@ -24,8 +24,8 @@ message("Loading model fit and recipe")
 lgbm_final_full_fit <- lightsnip::lgbm_load(paths$output$workflow_fit$local)
 lgbm_final_full_recipe <- readRDS(paths$output$workflow_recipe$local)
 
-if (shap_enable) {
-  message("Loading assessment data for SHAP calculation")
+if (shap_enable || comp_enable) {
+  message("Loading assessment data for SHAP and comp calculation")
 
   # Load the input data used for assessment. This is the universe of CARDs (not
   # PINs) that need values. Will use the the trained model to calc SHAP values
@@ -39,6 +39,13 @@ if (shap_enable) {
   )
 }
 
+if (comp_enable) {
+  message("Loading predicted values for comp calculation")
+
+  assessment_card <- read_parquet(paths$output$assessment_card$local) %>%
+    as_tibble()
+}
+
 
 
 
@@ -105,6 +112,131 @@ lightgbm::lgb.importance(lgbm_final_full_fit$fit) %>%
   rename_with(~ paste0(.x, "_value"), gain:frequency) %>%
   write_parquet(paths$output$feature_importance$local)
 
+
+
+
+#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# 4. Find Comparables  ---------------------------------------------------------
+#- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+if (comp_enable) {
+  message("Finding comparables")
+
+  # Calculate the leaf node assignments for every predicted value.
+  # Due to integer overflow problems with leaf node assignment, we need to
+  # chunk our data such that they are strictly less than the limit of 1073742
+  # rows. More detail here: https://github.com/microsoft/LightGBM/issues/1884
+  chunk_size <- 500000
+  chunks <- split(
+    assessment_data_prepped,
+    ceiling(seq_along(assessment_data_prepped[[1]]) / chunk_size)
+  )
+  chunked_leaf_nodes <- chunks %>%
+    map(\(chunk) {
+      predict(
+        object = lgbm_final_full_fit$fit,
+        newdata = as.matrix(chunk),
+        type = "leaf",
+      )
+    })
+  # Prefer do.call(rbind, ...) over bind_rows() because the chunks are
+  # not guaranteed to have the same number of rows, and bind_rows() will raise
+  # an error in that case
+  leaf_nodes <- do.call(rbind, chunked_leaf_nodes) %>% as_tibble()
+
+  # Calculate weights representing feature importance, so that we can weight
+  # leaf node assignments based on the most important features.
+  # To do this, we need the training data so that we can compute the mean sale
+  # price and use it as the base model error
+  message("Extracting weights from training data")
+  training_data <- read_parquet(paths$input$training$local) %>%
+    filter(!ind_pin_is_multicard, !sv_is_outlier) %>%
+    as_tibble()
+
+  tree_weights <- extract_weights(
+    model = lgbm_final_full_fit$fit,
+    init_score = mean(training_data$meta_sale_price, na.rm = TRUE),
+    metric = params$model$objective
+  )
+  if (length(tree_weights) == 0 || sum(tree_weights) != 1) {
+    stop("Unable to extract tree weights: length 0 or weight sum != 1")
+  }
+
+  # Get predicted values and leaf node assignments for the training data
+  training_data_prepped <- recipes::bake(
+    object = lgbm_final_full_recipe,
+    new_data = training_data,
+    all_predictors()
+  )
+  training_leaf_nodes <- predict(
+    object = lgbm_final_full_fit$fit,
+    newdata = as.matrix(training_data_prepped),
+    type = "leaf"
+  ) %>%
+    as_tibble()
+  training_leaf_nodes$predicted_value <- predict(
+    object = lgbm_final_full_fit$fit,
+    newdata = as.matrix(training_data_prepped)
+  ) %>%
+    # Round predicted values down for binning
+    floor()
+
+  # Get predicted values for the assessment set, which we already have in
+  # the assessment card set
+  leaf_nodes$predicted_value <- assessment_data %>%
+    left_join(assessment_card, by = c("meta_pin", "meta_card_num")) %>%
+    # Round predicted values down for binning
+    mutate(pred_card_initial_fmv = floor(pred_card_initial_fmv)) %>%
+    dplyr::pull(pred_card_initial_fmv)
+
+  # Make sure that the leaf node tibbles are all integers, which is what
+  # the comps algorithm expects
+  leaf_nodes <- leaf_nodes %>%
+    mutate(across(everything(), ~ as.integer(.x)))
+  training_leaf_nodes <- training_leaf_nodes %>%
+    mutate(across(everything(), ~ as.integer(.x)))
+
+  # Do the comps calculation in Python because the code is simpler and faster
+  message("Calling out to python/comps.py to perform comps calculation")
+  comps_module <- import("python.comps")
+  tryCatch(
+    {
+      comps <- comps_module$get_comps(
+        leaf_nodes, training_leaf_nodes, tree_weights,
+        n = as.integer(20)
+      )
+    },
+    error = function(e) {
+      # Log the full Python traceback in case of an error
+      print(py_last_error())
+      stop("Encountered error in python/comps.py")
+    }
+  )
+  # Correct for the fact that Python is 0-indexed by incrementing the
+  # comp indexes by 1
+  comps[[1]] <- comps[[1]] + 1
+
+  # Translate comp indexes to PINs
+  comps[[1]] <- comps[[1]] %>%
+    mutate(across(everything(), \(idx_row) {
+      training_data[idx_row, ]$meta_pin
+    })) %>%
+    cbind(
+      pin = assessment_data$meta_pin,
+      card = assessment_data$meta_card_num
+    ) %>%
+    relocate(pin, card) %>%
+    rename_with(\(colname) gsub("comp_idx_", "comp_pin_", colname))
+
+  # Combine the comp indexes and scores into one dataframe and write to a file
+  cbind(comps[[1]], comps[[2]]) %>%
+    write_parquet(paths$output$comp$local)
+} else {
+  # If comp creation is disabled, we still need to write an empty stub file
+  # so DVC doesn't complain
+  arrow::write_parquet(data.frame(), paths$output$comp$local)
+}
+
 # End the stage timer and write the time elapsed to a temporary file
 tictoc::toc(log = TRUE)
 bind_rows(tictoc::tic.log(format = FALSE)) %>%

pipeline/05-finalize.R

@@ -89,6 +89,7 @@ metadata <- tibble::tibble(
   ratio_study_near_column = params$ratio_study$near_column,
   ratio_study_num_quantile = list(params$ratio_study$num_quantile),
   shap_enable = shap_enable,
+  comp_enable = comp_enable,
   cv_enable = cv_enable,
   cv_num_folds = params$cv$num_folds,
   cv_fold_overlap = params$cv$fold_overlap,

pipeline/06-upload.R

@@ -204,6 +204,20 @@ if (upload_enable) {
     relocate(run_id) %>%
     write_parquet(paths$output$feature_importance$s3)
 
+  # Upload comps
+  if (comp_enable) {
+    message("Uploading comps")
+    read_parquet(paths$output$comp$local) %>%
+      mutate(run_id = run_id, year = params$assessment$working_year) %>%
+      group_by(year, run_id) %>%
+      arrow::write_dataset(
+        path = paths$output$comp$s3,
+        format = "parquet",
+        hive_style = TRUE,
+        compression = "snappy"
+      )
+  }
+
 
   # 2.5. Finalize --------------------------------------------------------------
   message("Uploading run metadata, timings, and reports")