Ensure decimal place stability for MKI

[Damonamajor]

As found in issue ccao-data/data-architecture#737, there is some small (.00000000000x) non-deterministic aspect of the calculations. This is likely due to how decimals are handled in the calculations.

[Damonamajor]

Workflow

• Re-run the ratio_stats python script using DBT. This produces a new output in the dev environment.
• Run a query to identify values where key metrics do not match. This can be queried using the following.

SELECT 
    r.year,
    r.triad,
    r.geography_type,
    r.property_group,
    r.assessment_stage,
    r.geography_id,
    r.sale_year,
    ROUND(r.prb, 3) AS r_prb,
    ROUND(d.prb, 3) AS d_prb,
    CASE 
      WHEN ROUND(r.prb, 3) <> ROUND(d.prb, 3) THEN 'Mismatch' 
      ELSE 'Match' 
    END AS prb_status,
    ROUND(r.cod, 3) AS r_cod,
    ROUND(d.cod, 3) AS d_cod,
    CASE 
      WHEN ROUND(r.cod, 3) <> ROUND(d.cod, 3) THEN 'Mismatch' 
      ELSE 'Match' 
    END AS cod_status,
    ROUND(r.mki, 3) AS r_mki,
    ROUND(d.mki, 3) AS d_mki,
    CASE 
      WHEN ROUND(r.mki, 3) <> ROUND(d.mki, 3) THEN 'Mismatch' 
      ELSE 'Match' 
    END AS mki_status,
    ROUND(r.prd, 3) AS r_prd,
    ROUND(d.prd, 3) AS d_prd,
    CASE 
      WHEN ROUND(r.prd, 3) <> ROUND(d.prd, 3) THEN 'Mismatch' 
      ELSE 'Match' 
    END AS prd_status
FROM reporting.ratio_stats r 
JOIN z_dev_damajor_reporting.ratio_stats d
  ON r.year             = d.year
 AND r.triad            = d.triad
 AND r.geography_type   = d.geography_type
 AND r.property_group   = d.property_group
 AND r.assessment_stage = d.assessment_stage
 AND r.geography_id     = d.geography_id
 AND r.sale_year        = d.sale_year
WHERE (ROUND(r.prb, 3) <> ROUND(d.prb, 3))
   OR (ROUND(r.cod, 3) <> ROUND(d.cod, 3))
   OR (ROUND(r.mki, 3) <> ROUND(d.mki, 3))
   OR (ROUND(r.prd, 3) <> ROUND(d.prd, 3));

• Use the following query to identify mki values for a grouping which doesn't match (in reporting and z_dev_damajor_reporting).
• (0.8463302650269513 - 0.8462841130977653)

select * 
from z_dev_damajor_reporting.ratio_stats
where
year = '2021'
and triad = '1'
and geography_type = 'Town'
and property_group = 'MF'
and assessment_stage = 'assessor certified'
and geography_id = '71'
and sale_year = '2020'

• exported the attached csv to QC with via the following query

select * 
from reporting.ratio_stats_input
where
year = '2021'
and triad = '1'
and property_group = 'MF'
and assessment_stage = 'assessor certified'
and township_code = '71'
and sale_year = '2020'

Use the python test file with that csv where all outputs produce values of 0.8463093397231315.
This does not match either of the other dbt outputs.

@jeancochrane

test_data.csv

[jeancochrane]

I did a bunch more poking around at this and came to two different conclusions:

• MKI can be made fully reproducible by tweaking the Gini calculation to sort the input data by estimate in addition to sale_price, which suggests that it is indeed a tiebreaking problem. I'd like to read the Quintos paper to confirm that this is a reasonable approach before we make this change, however.

assesspy/assesspy/metrics.py

Line 215 in a0d0359

df.sort_values(by="sale_price", kind="mergesort", inplace=True)

• PRB/PRD/COD mismatches seem quite likely to be caused by performing repeated floating point arithmetic operations on small values. I was able to make PRB and COD reproducible by updating all of the arithmetic logic in metrics.py to use fixed-precision decimal math, but it's a big diff that changes some mission-critical code and it'll be tricky to unit test. We could either invest time in this change or instead ensure that our calling code always rounds to 6 decimal places, the threshold at which we've been able to consistently reproduce values.

[jeancochrane]

After reading both Quintos papers (2020 and 2021) and corresponding with the author, I'm feeling confident in a path forward for MKI:

Dr. Quintos advises we apply the same cumulative shares of property and value (the variables R and C_SP in the paper, respectively) to all sales with a matching rank. In the example she provided, which uses different but equivalent formulae to calculate the Gini coefficient of sale prices and the concentration index (CI) of assessments than this package does, this approach corresponds to the following observed behavior for a group of sales with the same price:

• The first sale in the group (chosen arbitrarily) receives share values SR and SP equivalent to the sum of those share values for all sales in the group
• The rest of the sales in the group (chosen arbitrarily) receive share values SR and SP of 0
• Because the cumulative shares R and C_SP correspond to the rolling sums of SR and SP, each sale in the group will have the same value of R and C_SP, which will corresponding to each sale getting the sum of SR and SP as applied to the group added to the rolling sum of shares of all sales with lower prices

Since we use different but equivalent formulae to calculate the Gini coefficient of sale prices and the CI of assessments, this method doesn't quite work for us out of the box. I think it might make sense for us to switch to using Dr. Quintos' formula so that we can hew as closely to the published methods as possible.

Here's some sample Python code that does this and replicates the values in the Excel doc shared with us:

--- a/assesspy/metrics.py
+++ b/assesspy/metrics.py
@@ -1,6 +1,7 @@
 import math
 from typing import Union

+import numpy as np
 import pandas as pd
 import statsmodels.api as sm

@@ -217,13 +218,35 @@ def _calculate_gini(
     a_sorted, sp_sorted = df["estimate"], df["sale_price"]
     n: int = a_sorted.size

-    assessed_sum: float = sum(a_sorted[i] * (i + 1) for i in range(n))
-    g_assessed: float = 2 * assessed_sum / a_sorted.sum() - (n + 1)
-    gini_assessed: float = g_assessed / float(n)
-
-    sale_price_sum: float = sum(sp_sorted[i] * (i + 1) for i in range(n))
-    g_sale_price: float = 2 * sale_price_sum / sp_sorted.sum() - (n + 1)
-    gini_sale_price: float = g_sale_price / float(n)
+    # Compute competition rank for sale prices in order to properly
+    # group sales with the same price
+    def _frequency(arr):
+        rank = pd.Series(arr).rank(method="dense").to_numpy()
+
+        frequency = []
+        _, first_occurrence = np.unique(rank, return_index=True)
+        for value_count in np.diff(np.append(first_occurrence, len(rank))):
+            frequency.append(int(value_count))
+            if value_count > 1:
+                frequency.extend([0] * (value_count - 1))
+        return np.array(frequency, dtype=int)
+
+    # Calculate frequencies
+    freq_sp = _frequency(sp_sorted)
+    freq_a = np.ones_like(a_sorted)
+
+    # Calculate shares
+    sr = freq_sp / np.sum(freq_sp)
+    sp = sp_sorted * freq_sp / np.sum(sp_sorted)
+    sa = a_sorted * freq_a / np.sum(a_sorted)
+
+    # Calculate cumulative share and covariance
+    r = np.cumsum(sr)
+
+    # numpy.cov() returns a matrix with all possible covariances, so we need
+    # to select the one that represents the covariance between `r` and `sp`
+    gini_assessed = (2 * np.cov(r, sa, bias=True)[0][1]) / np.mean(sa)
+    gini_sale_price = (2 * np.cov(r, sp, bias=True)[0][1]) / np.mean(sp)

     return gini_assessed, gini_sale_price

Let's switch to this method for both Python and R and use the Excel values Dr. Quintos shared with us as a new source of fixture data to test the results. We should also confirm that this method exactly matches the sample data from the original paper, which is already included in the test suite but does not have ordering problems.

[jeancochrane]

I'm going to put this on the backburner for now, since it's not urgent.

[claireboyd]

@jeancochrane hi from the NYC assessor's office! We also encountered this bug and talked at length (with Dr. Quintos) about the best way to approach this.

After a lot of deliberation, we decided the best way forward is to assume the "worst case scenario" in terms of MKI/KI metrics by sorting the data first by the ascending actual value (sale price) and then by the descending predicted value (modeled result). Not saying this has to be your solution, but wanted to share our thinking if helpful.

I couldn't find this equivalent issue in R - so sorry for commenting in R code in a python package - but here is how we decided to approach it.

#' Generate Inequality Metrics
#'
#' This function computes four commonly used inequality metrics—Gini, Concentration Index (CI),
#' Kakwani Index (KI), and Modified Kakwani Index (MKI)—to evaluate disparities between assessed
#' (predicted) values and actual prices. These metrics are widely used in health equity research,
#' including work by Kakwani et al., to assess health outcomes in relation to income or wealth.
#' We apply this same thinking to the property assessment application.
#'
#' The function supports two equivalent formulations: one based on covariance, and another based
#' on regression. While both produce the same inequality estimates, we use the regression-based
#' approach here to additionally obtain p-values for statistical significance.
#'
#' Special care is taken to ensure that the results are not affected by random variation due to
#' tied price values. In practice, we observed that the Concentration Index (CI) and Gini Index
#' can vary depending on the order of the assessed values when price values are repeated. The approach
#' proposed by Chen and Roy 2009—grouping tied values—does not meet our needs, as we aim to evaluate
#' variation in assessed values within each price level, not simply across prices.
#'
#' Therefore, we adopt a "worst-case scenario" strategy: we sort the data by ascending price and
#' descending predicted value. This means that when multiple properties have the same price,
#' the highest assessed values are ranked first. As a result, the generated metrics will be as
#' regressive as possible given the data. We believe this conservative choice is preferable, as
#' it removes arbitrary randomness from the results and ensures consistency across runs.
#'
#' #' ## Interpretation Guidelines for KI and MKI
#' To assess vertical equity using the KI and MKI indices, two sets of standards are commonly used:
#'
#' **1. Based on p-values:**
#' - For *homogenous* districts, a 5% significance level (p < 0.05) is used.
#' - For *non-homogenous* districts or small samples, a 10% significance level (p < 0.10) is appropriate.
#'
#' **2. Based on metric ranges (especially when sample assumptions are violated):**
#' - **KI**: should fall within (-0.01, 0.01) for homogenous districts, and (-0.02, 0.02) for non-homogenous or small samples.
#' - **MKI**: should fall within (0.95, 1.05) for homogenous districts, and (0.90, 1.10) for non-homogenous or small samples.
#'
#' These ranges provide a practical rule of thumb when p-value assumptions—such as large sample size,
#' clean data, and properly specified models—do not hold.
#'
#' Note that KI and MKI are derived from different regressions and capture distinct forms of vertical equity.
#' It is possible to reject vertical equity based on one but not the other. In such cases, there may be
#' underlying issues such as unresolved horizontal inequity or significant unmodeled heterogeneity.
#'
#' @param predicted A numeric vector of assessed or model-predicted values.
#' @param actual A numeric vector of actual prices (must be the same length as `predicted`).
#'
#' @return A named list containing:
#' \describe{
#'   \item{gini}{Gini index}
#'   \item{ci}{Concentration index}
#'   \item{ki}{Kakwani index}
#'   \item{mki}{Modified Kakwani index}
#'   \item{gini_p_value}{p-value for the Gini index}
#'   \item{ci_p_value}{p-value for the Concentration index}
#'   \item{ki_p_value}{p-value for the Kakwani index}
#'   \item{mki_p_value}{p-value for the Modified Kakwani index}
#' }
#'
#' @examples
#' data(quintos_sample)
#' generate_inequality_metrics(
#'     predicted = quintos_sample$assessment1,
#'     actual = quintos_sample$price
#' )
#'
#' @export
#' @importFrom stats lm
#' @import dplyr
#'
generate_inequality_metrics = function(predicted, actual){

    if(length(predicted) != length(actual)){
        stop("Length of predicted and actual inputs are not equal.")
    }

    if (any(is.na(predicted)) || any(is.na(actual))) {
        stop("Inputs contain NA values. Please remove or impute missing data.")
    }

    n = length(predicted)

    # Create a dataframe with computed ranks
    dataframe <- data.frame(
        predicted = predicted,
        actual = actual
    ) %>%
        dplyr::arrange(actual, desc(predicted)) %>%
        dplyr::mutate(rank = rank(actual, ties.method = "min")) %>%
        dplyr::mutate(
            SR = 1/n(),
            SP = actual*1/sum(actual),
            SA = predicted*1/sum(predicted))%>%
        dplyr::mutate(
            R = cumsum(SR),
            C_SP = cumsum(SP),
            C_SA = cumsum(SA)
        ) %>%
        dplyr::arrange(rank) %>%
        dplyr::mutate(
            R = cumsum(SR),
            C_SP = cumsum(SP),
            C_SA = cumsum(SA),
            y_gini = 2*var(R)*(SP/mean(SP)),
            y_ci = 2*var(R)*(SA/mean(SA)),
            y_ki = 2*var(R)*(SA/mean(SA) - SP/mean(SP)),
            y_mki_1 = SP/mean(SP),
            y_mki_2 = SA/mean(SA)
        )

    # covariance method - same results but no p-value, so depreciated
    # gini = 2 * cov(dataframe$R, dataframe$SP) / mean(dataframe$SP)
    # ci = 2 * cov(dataframe$R, dataframe$SA) / mean(dataframe$SA)
    #
    # covariance_results = list(
    #     gini_covariance = gini,
    #     ci_covariance = ci,
    #     ki_covariance = ci - gini,
    #     mki_covariance = ci / gini
    # )

    # gini regression
    gini_model = summary(stats::lm(y_gini ~ R, data=dataframe))
    gini = gini_model$coefficients['R','Estimate']
    gini_p_value = gini_model$coefficients['R','Pr(>|t|)']

    # ci regression
    ci_model = summary(stats::lm(y_ci ~ R, data=dataframe))
    ci = ci_model$coefficients['R','Estimate']
    ci_p_value = ci_model$coefficients['R','Pr(>|t|)']

    # ki regression
    ki_model = summary(stats::lm(y_ki ~ R, data=dataframe))
    ki = ki_model$coefficients['R','Estimate']
    ki_p_value = ki_model$coefficients['R','Pr(>|t|)']

    # mki stage 1 regression
    mki_model_1 = stats::lm(y_mki_1 ~ R, data=dataframe)
    dataframe$mki_predicted_values = predict(mki_model_1)

    # mki stage 2 regression
    mki_model_2 = summary(stats::lm(y_mki_2 ~ mki_predicted_values, data=dataframe))
    mki = mki_model_2$coefficients['mki_predicted_values','Estimate']
    mki_p_value = mki_model_2$coefficients['mki_predicted_values','Pr(>|t|)']

    regression_results = list(
        gini = gini,
        ci = ci,
        ki = ki,
        mki = mki,
        gini_p_value = round(as.numeric(format(gini_p_value, scientific = FALSE)),digits=7),
        ci_p_value = round(as.numeric(format(ci_p_value, scientific = FALSE)),digits=7),
        ki_p_value = round(as.numeric(format(ki_p_value, scientific = FALSE)),digits=7),
        mki_p_value = round(as.numeric(format(mki_p_value, scientific = FALSE)),digits=7)
    )

    return(regression_results)
}

[jeancochrane]

Thanks @claireboyd, super helpful to hear how your team resolved this issue! Switching to your worst-case approach would be a much simpler refactor for us, so we will likely follow your lead. Does your team happen to have any test cases for generate_inequality_metrics()?

[claireboyd]

sure do! here is the sample data set we use that has three tied values for price at $89,000, with three different assessment columns that sort the modeled values differently. here is our test suite for the function -- the second to last of which ensures that the results of the function are consistent across the three different sorts:

test_that("Output structure is as expected", {
    predicted <- quintos_sample$assessment1
    actual <- quintos_sample$price
    result <- generate_inequality_metrics(predicted, actual)

    expected_names <- c(
        "gini", "ci", "ki", "mki",
        "gini_p_value", "ci_p_value", "ki_p_value", "mki_p_value"
    )
    expect_type(result, "list")
    expect_named(result, expected_names)
    lapply(result, function(x) expect_type(x, "double"))
})

test_that("Handles unequal input lengths", {
    predicted <- c(1000, 123, 1234235)
    actual <- c(1000, 123, 1234235, 1245, 124352)
    expect_error(
        generate_inequality_metrics(predicted, actual),
        "Length of predicted and actual inputs are not equal"
    )
})

test_that("Returns reasonable values on known pattern", {
    predicted <- quintos_sample$assessment1
    actual <- quintos_sample$price
    result <- generate_inequality_metrics(predicted, actual)

    expect_true(is.numeric(result$gini))
    expect_true(is.numeric(result$ci))
    expect_true(is.numeric(result$ki))
    expect_true(is.numeric(result$mki))
})

test_that("Returns equivalent results for differently ordered data", {
    predicted1 <- quintos_sample$assessment1
    predicted2 <- quintos_sample$assessment2
    predicted3 <- quintos_sample$assessment3
    actual <- quintos_sample$price

    result1 <- generate_inequality_metrics(predicted1, actual)
    result2 <- generate_inequality_metrics(predicted2, actual)
    result3 <- generate_inequality_metrics(predicted3, actual)

    expect_equal(result1, result2)
    expect_equal(result2, result3)
    expect_equal(result1, result3)
})

test_that("Handles NA values appropriately", {
    predicted <- c(1000, 123, 1234235, NA, NA)
    actual <- c(1000, 123, 1234235, 1245, 124352)
    expect_error(
        generate_inequality_metrics(predicted, actual),
        "Inputs contain NA values"
    )
})