openavmkit.utilities.stats

Statistics engine for ratio studies, equity, and modeling diagnostics.

The numeric core behind OpenAVMKit's IAAO-aligned analysis. Implements:

• Ratio-study statistics — COD, PRD, PRB, with bootstrap confidence intervals (calc_cod, calc_prd, calc_prb, calc_ratio_stats_bootstrap)
• Horizontal equity — CHD (Coefficient of Horizontal Dispersion) per cluster
• Outlier trimming — trim_outlier_ratios, trim_outliers_mask
• Variable-selection diagnostics — correlation, VIF, p-values, t-values, ENR coefficients, adjusted R²
• Model fit metrics — MSE, R², adjusted R², MAPE
• ConfidenceStat — value plus bootstrap-derived confidence interval

Functions here are pure-numeric and have no settings or DataFrame-shape assumptions beyond what each function declares.

ConfidenceStat

ConfidenceStat(value, confidence_interval, low, high)

Any statistic along with it's confidence interval upper and lower bounds, and whether it is statistically significant

Attributes:

Name	Type	Description
value	float	The base value of the statistic
confidence_interval	float	The % value of the confidence interval (e.g. 0.95 for 95% confidence interval)
low	float	The lower bound of the confidence interval
high	float	The upper bound of the confidence interval

Source code in openavmkit/utilities/stats.py

def __init__(
    self,
    value : float,
    confidence_interval : float,
    low : float,
    high : float
):
    self.value = value
    self.confidence_interval = confidence_interval
    self.low = low
    self.high = high

calc_chds

calc_chds(df_in, field_cluster, field_value)

Calculate the Coefficient of Horizontal Dispersion (CHD) for each cluster in a DataFrame.

CHD is the same statistic as COD, the Coefficient of Dispersion, but calculated for horizontal equity clusters and used to measure horizontal dispersion, on the theory that similar properties in similar locations should have similar valuations. The use of the name "CHD" is chosen to avoid confusion because assessors strongly associate "COD" with sales ratio studies.

This function computes the CHD for each unique cluster in the input DataFrame based on the values in the specified field.

Parameters:

Name	Type	Description	Default
df_in	DataFrame	Input DataFrame	required
field_cluster	str	Name of the column representing cluster identifiers.	required
field_value	str	Name of the column containing the values for COD calculation	required

Returns:

Type	Description
A Series of COD values for each row, aligned with df_in

Source code in openavmkit/utilities/stats.py

def calc_chds(df_in: pd.DataFrame, field_cluster: str, field_value: str) -> pd.Series:
    """Calculate the Coefficient of Horizontal Dispersion (CHD) for each cluster in a
    DataFrame.

    CHD is the same statistic as COD, the Coefficient of Dispersion, but
    calculated for horizontal equity clusters and used to measure horizontal dispersion, on
    the theory that similar properties in similar locations should have similar valuations.
    The use of the name "CHD" is chosen to avoid confusion because assessors strongly
    associate "COD" with sales ratio studies.

    This function computes the CHD for each unique cluster in the input DataFrame based on
    the values in the specified field.

    Parameters
    ----------
    df_in : pd.DataFrame
        Input DataFrame
    field_cluster : str
        Name of the column representing cluster identifiers.
    field_value : str
        Name of the column containing the values for COD calculation

    Returns
    -------
    A Series of COD values for each row, aligned with df_in
    """
    # Create a dataframe matching the index of df_in with only the cluster id.
    df = df_in[[field_cluster]].copy()
    df["chd"] = 0.0

    clusters = df[field_cluster].unique()

    for cluster in clusters:
        df_cluster = df[df[field_cluster].eq(cluster)]
        # exclude negative and null values:
        df_cluster = df_cluster[
            ~pd.isna(df_cluster[field_value]) & df_cluster[field_value].gt(0)
        ]

        chd = calc_cod(df_cluster[field_value].values)
        df.loc[df[field_cluster].eq(cluster), "chd"] = chd

    return df["chd"]

calc_cod

calc_cod(values)

Calculate the Coefficient of Dispersion (COD) for an array of values.

COD is defined as the average absolute deviation from the median, divided by the median, multiplied by 100. Special cases are handled if the median is zero.

Parameters:

Name	Type	Description	Default
values	ndarray	Array of numeric values.	required

Returns:

Type	Description
float	The COD percentage.

Source code in openavmkit/utilities/stats.py

def calc_cod(values: np.ndarray) -> float:
    """
    Calculate the Coefficient of Dispersion (COD) for an array of values.

    COD is defined as the average absolute deviation from the median, divided by the
    median, multiplied by 100. Special cases are handled if the median is zero.

    Parameters
    ----------
    values : numpy.ndarray
        Array of numeric values.

    Returns
    -------
    float
        The COD percentage.
    """
    if len(values) == 0:
        return float("nan")

    median_value = np.median(values)
    abs_delta_values = np.abs(values - median_value)
    avg_abs_deviation = np.sum(abs_delta_values) / len(values)
    if median_value == 0:
        # if every value is zero, the COD is zero:
        if np.all(values == 0):
            return 0.0
        else:
            # if the median is zero but not all values are zero, return infinity
            return float("inf")
    cod = avg_abs_deviation / median_value
    cod *= 100
    return cod

calc_cod_bootstrap

calc_cod_bootstrap(values, confidence_interval=0.95, iterations=10000, seed=777)

Calculate COD using bootstrapping.

This function bootstraps the input values (resampling with replacement) to generate a distribution of CODs, then returns the median COD along with the lower and upper bounds of the confidence interval.

Parameters:

Name	Type	Description	Default
values	ndarray	Array of numeric values.	required
confidence_interval	float	The desired confidence level. Defaults to 0.95.	0.95
iterations	int	Number of bootstrap iterations. Defaults to 10000.	10000
seed	int	Random seed for reproducibility. Defaults to 777.	777

Returns:

Type	Description
tuple[float, float, float]	A tuple containing the median COD, lower bound, and upper bound of the confidence interval.

Source code in openavmkit/utilities/stats.py

def calc_cod_bootstrap(
    values: np.ndarray,
    confidence_interval: float = 0.95,
    iterations: int = 10000,
    seed: int = 777
) -> tuple[float, float, float]:
    """
    Calculate COD using bootstrapping.

    This function bootstraps the input values (resampling with replacement) to
    generate a distribution of CODs, then returns the median COD along with the
    lower and upper bounds of the confidence interval.

    Parameters
    ----------
    values : numpy.ndarray
        Array of numeric values.
    confidence_interval : float, optional
        The desired confidence level. Defaults to 0.95.
    iterations : int, optional
        Number of bootstrap iterations. Defaults to 10000.
    seed : int, optional
        Random seed for reproducibility. Defaults to 777.

    Returns
    -------
    tuple[float, float, float]
        A tuple containing the median COD, lower bound, and upper bound of the confidence interval.
    """

    n = len(values)
    if n == 0:
        return float("nan"), float("nan"), float("nan")
    np.random.seed(seed)

    # Replace negative values with zero:
    values = np.where(values < 0, 0, values)

    median = np.median(values)
    samples = np.random.choice(values, size=(iterations, n), replace=True)
    abs_delta_values = np.abs(samples - median)
    bootstrap_cods = np.mean(abs_delta_values, axis=1) / median * 100
    alpha = (1.0 - confidence_interval) / 2
    lower_bound, upper_bound = np.quantile(bootstrap_cods, [alpha, 1.0 - alpha])
    median_cod = np.median(bootstrap_cods)
    return median_cod, lower_bound, upper_bound

calc_correlations

calc_correlations(X, threshold=0.1, do_plots=False, min_strength=None, overlap=True)

Calculate correlations and iteratively drop variables with low combined scores.

This function computes the correlation matrix of X, then calculates a combined score for each variable based on its correlation strength with the target variable and its average correlation with other variables. Variables whose scores fall below the specified threshold are removed.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input DataFrame containing the variables to evaluate.	required
threshold	float	Minimum acceptable combined score for variables. Variables with a score below this value will be dropped. Defaults to 0.1.	0.1
do_plots	bool	If True, plot the initial and final correlation heatmaps. Defaults to False.	False
min_strength	float or None	A first pass: variables whose correlation with the target (the first column) is weaker than this, in absolute value, are dropped before anything else, so they neither get a combined score nor lower the others' clarity. None (the default) skips it.	None
overlap	bool	Whether overlap counts: a variable's overlap with the others lowers its score (its clarity, squared), and one more correlated with a stronger variable than with the target is dropped. False: the score is its correlation with the target alone and nothing is dropped for overlapping; overlap is left to another test (VIF). Defaults to True.	True

Returns:

Type	Description
dict	Dictionary with keys: "initial": pandas.Series of combined scores from the first iteration. "final": pandas.Series of combined scores after dropping low-scoring variables. "bad_vars": variables dropped (null correlations, and overlaps). "weak": variables dropped by the first pass (min_strength), if any. "matrix": the correlation matrix of every variable, the target first.

Source code in openavmkit/utilities/stats.py

def calc_correlations(
    X: pd.DataFrame,
    threshold: float = 0.1,
    do_plots: bool = False,
    min_strength: float | None = None,
    overlap: bool = True,
) -> dict:
    """
    Calculate correlations and iteratively drop variables with low combined scores.

    This function computes the correlation matrix of `X`, then calculates a combined score
    for each variable based on its correlation strength with the target variable and its
    average correlation with other variables. Variables whose scores fall below the
    specified `threshold` are removed.

    Parameters
    ----------
    X : pandas.DataFrame
        Input DataFrame containing the variables to evaluate.
    threshold : float, optional
        Minimum acceptable combined score for variables. Variables with a score below this
        value will be dropped. Defaults to 0.1.
    do_plots : bool, optional
        If True, plot the initial and final correlation heatmaps. Defaults to False.
    min_strength : float or None, optional
        A first pass: variables whose correlation with the target (the first column) is
        weaker than this, in absolute value, are dropped before anything else, so they
        neither get a combined score nor lower the others' clarity. None (the default)
        skips it.
    overlap : bool, optional
        Whether overlap counts: a variable's overlap with the others lowers its score (its
        clarity, squared), and one more correlated with a stronger variable than with the
        target is dropped. False: the score is its correlation with the target alone and
        nothing is dropped for overlapping; overlap is left to another test (VIF). Defaults
        to True.

    Returns
    -------
    dict
        Dictionary with keys:

        - "initial": pandas.Series of combined scores from the first iteration.
        - "final": pandas.Series of combined scores after dropping low-scoring variables.
        - "bad_vars": variables dropped (null correlations, and overlaps).
        - "weak": variables dropped by the first pass (``min_strength``), if any.
        - "matrix": the correlation matrix of every variable, the target first.
    """
    X = X.copy()
    first_run = None

    bad_vars = _corr_problem_columns(X)
    if bad_vars:
        X = X.drop(columns=bad_vars, errors="ignore")
        print("")
        print(f"Found {len(bad_vars)} columns with null correlations, recommend dropping:")
        for col in bad_vars:
            print(col)

    # Normalize all numerical values prior to computation:
    for col in X.columns:
        if X[col].dtype != "object":
            X[col] = (X[col] - X[col].mean()) / X[col].std()

    corr_full = X.corr()
    remaining = list(X.columns)

    # First pass (opt-in): no correlation with the target at all, out before the scores are worked out.
    weak = []
    if min_strength is not None and len(remaining) > 1:
        to_target = corr_full[remaining[0]].abs()
        weak = [c for c in remaining[1:] if not (to_target[c] >= min_strength)]
        remaining = [c for c in remaining if c not in weak]
        if len(remaining) < 2:
            empty = pd.DataFrame(columns=["variable", "corr_strength", "corr_clarity", "corr_score"])
            return {"initial": empty, "final": empty, "bad_vars": bad_vars, "weak": weak, "matrix": corr_full}

    # Keep a copy of the *initial* ordered corr for plotting later
    naive_corr_for_plot = None  # <--- NEW

    while True:
        naive_corr = corr_full.loc[remaining, remaining]

        # --- NEW: compute one explicit order and apply to both axes
        order = _order_by_target_corr(naive_corr, target_col=naive_corr.columns[0])
        naive_corr = naive_corr.loc[order, order]
        # -----------------------------------------------

        if naive_corr_for_plot is None:  # <--- NEW
            naive_corr_for_plot = naive_corr.copy()

        naive_corr_sans_target = naive_corr.iloc[1:, 1:]

        strength = naive_corr.iloc[:, 0].abs().iloc[1:]

        if overlap:
            clarity = 1 - (
                (naive_corr_sans_target.abs().sum(axis=1) - 1.0)
                / (len(naive_corr_sans_target.columns) - 1)
            )
        else:  # judged on the relation to the target alone; overlap is left to another test (VIF)
            clarity = pd.Series(1.0, index=naive_corr_sans_target.index)

        score = strength * clarity * clarity

        # Guard against all-NA scores (too few samples to compute correlations)
        if score.isna().all():
            break
        min_score_idx = score.idxmin()
        try:
            min_score = score[min_score_idx]
        except KeyError:
            min_score = float('nan')

        data = {"corr_strength": strength, "corr_clarity": clarity, "corr_score": score}
        df_score = pd.DataFrame(data).reset_index().rename(columns={"index": "variable"})

        if first_run is None:
            first_run = df_score.sort_values("corr_score", ascending=False)

        if min_score < threshold:
            remaining.remove(min_score_idx)
        else:
            break

    # ------------------------------------------------------------------
    # NEW DROP LOGIC (second pass):
    # For each variable, starting with the strongest (most correlated with target),
    # drop all other variables that are more correlated with current_variable than
    # they are with the target.
    # ------------------------------------------------------------------
    target_col = remaining[0]

    # process in descending strength-to-target order (excluding target itself)
    ordered_by_target = (
        corr_full.loc[remaining, target_col]
        .abs()
        .drop(labels=[target_col], errors="ignore")
        .sort_values(ascending=False)
        .index
        .tolist()
    )

    all_dropped = []

    # Without overlap, a variable more related to another than to the target stays: VIF judges overlap.
    for current_variable in ordered_by_target if overlap else []:
        if current_variable not in remaining:
            continue

        # Compare every other variable's corr to current_variable vs corr to target
        others = [v for v in remaining if v not in (target_col, current_variable)]
        to_drop = []
        for v in others:
            corr_to_current = abs(corr_full.loc[v, current_variable])
            corr_to_target = abs(corr_full.loc[v, target_col])
            if corr_to_current > corr_to_target:
                to_drop.append(v)

        if to_drop:
            for v in to_drop:
                if v in remaining:
                    remaining.remove(v)
                    all_dropped.append(v)

    all_dropped = list(set(all_dropped))
    if len(all_dropped) > 0:
        bad_vars = list(set(bad_vars+all_dropped))

    # ------------------------------------------------------------------

    # Recompute scores for the *final* remaining set (same math as above)
    naive_corr = corr_full.loc[remaining, remaining]
    order = _order_by_target_corr(naive_corr, target_col=naive_corr.columns[0])
    naive_corr = naive_corr.loc[order, order]

    naive_corr_sans_target = naive_corr.iloc[1:, 1:]
    strength = naive_corr.iloc[:, 0].abs().iloc[1:]
    if overlap:
        clarity = 1 - (
            (naive_corr_sans_target.abs().sum(axis=1) - 1.0)
            / (len(naive_corr_sans_target.columns) - 1)
        )
    else:
        clarity = pd.Series(1.0, index=strength.index)
    score = strength * clarity * clarity

    data = {"corr_strength": strength, "corr_clarity": clarity, "corr_score": score}
    df_score = (
        pd.DataFrame(data)
        .reset_index()
        .rename(columns={"index": "variable"})
        .sort_values("corr_score", ascending=False)
    )

    if do_plots and naive_corr_for_plot is not None:
        plot_correlation(naive_corr_for_plot, "Correlation of Variables (initial)")

    # recompute correlation matrix for remaining vars
    final_corr = corr_full.loc[remaining, remaining]

    # --- NEW: order final_corr the same way (by target corr within remaining set)
    final_order = _order_by_target_corr(final_corr, target_col=final_corr.columns[0])
    final_corr = final_corr.loc[final_order, final_order]
    # -----------------------------------------------

    if do_plots:
        plot_correlation(final_corr, "Correlation of Variables (final)")

    return {"initial": first_run, "final": df_score, "bad_vars": bad_vars, "weak": weak, "matrix": corr_full}

calc_cross_validation_score

calc_cross_validation_score(X, y)

Calculate cross-validation score using negative mean squared error.

This function fits a LinearRegression model using 5-fold cross validation and returns the mean cross-validated MSE (positive value).

Parameters:

Name	Type	Description	Default
X	array - like or DataFrame	Input features for modeling.	required
y	array - like or Series	Target variable.	required

Returns:

Type	Description
float	The mean cross-validated mean squared error.

Source code in openavmkit/utilities/stats.py

def calc_cross_validation_score(
    X: pd.DataFrame | np.ndarray,
    y: pd.Series | np.ndarray
) -> float:
    """
    Calculate cross-validation score using negative mean squared error.

    This function fits a LinearRegression model using 5-fold cross validation and
    returns the mean cross-validated MSE (positive value).

    Parameters
    ----------
    X : array-like or pandas.DataFrame
        Input features for modeling.
    y : array-like or pandas.Series
        Target variable.

    Returns
    -------
    float
        The mean cross-validated mean squared error.
    """
    # Impute NaN with column medians, matching calc_elastic_net_regularization,
    # calc_p_values_recursive_drop, calc_t_values_recursive_drop and
    # calc_vif_recursive_drop.
    #
    # Without this, sklearn raises "Input X contains NaN" inside each affected fold.
    # cross_val_score does not propagate that as an exception -- it records NaN for the
    # failed folds and warns -- so the `except ValueError` below never fires and this
    # returns NaN. The caller in model_runner.get_variable_recommendations then does
    # `if cv_score < best_score`, and `nan < x` is always False, so best_variables is
    # never updated: the cross-validation refinement silently does nothing on any
    # dataset with a NaN anywhere in X. Real assessment data routinely has some.
    if isinstance(X, pd.DataFrame):
        if X.isnull().values.any():
            warnings.warn(
                f"calc_cross_validation_score: NaN detected in "
                f"{list(X.columns[X.isnull().any()])}. "
                "Imputing with column medians for the cross-validation step only.",
                UserWarning,
            )
            # A column that is entirely NaN has no median, so the fill leaves it NaN
            # and we would be right back to the silent-NaN failure above. It carries
            # no signal either way, so flatten it to 0 and let the fit proceed.
            X = X.fillna(X.median(numeric_only=True)).fillna(0.0)
    else:
        X = np.asarray(X, dtype=float)
        if np.isnan(X).any():
            warnings.warn(
                "calc_cross_validation_score: NaN detected in X. "
                "Imputing with column medians for the cross-validation step only.",
                UserWarning,
            )
            medians = np.nanmedian(X, axis=0)
            # An all-NaN column has no median; fall back to 0 so the fit can proceed.
            medians = np.where(np.isnan(medians), 0.0, medians)
            X = np.where(np.isnan(X), medians, X)

    model = LinearRegression()
    # Use negative MSE and negate it to return positive MSE
    try:
        scores = cross_val_score(model, X, y, cv=5, scoring="neg_mean_squared_error")
    except ValueError:
        return float("nan")

    return -scores.mean()  # Convert negative MSE to positive

calc_elastic_net_regularization

calc_elastic_net_regularization(X, y, threshold_fraction=0.05)

Calculate Elastic Net regularization coefficients while iteratively dropping variables with low coefficients.

This function standardizes X, fits an Elastic Net model, and iteratively removes variables whose absolute coefficients fall below a fraction of the maximum coefficient.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input features DataFrame.	required
y	Series	Target variable series.	required
threshold_fraction	float	Fraction of the maximum coefficient below which variables are dropped. Defaults to 0.05.	0.05

Returns:

Type	Description
dict	Dictionary with two keys: "initial": pandas.Series of coefficients from the first Elastic Net fit. "final": pandas.Series of coefficients after dropping low-magnitude variables.

Source code in openavmkit/utilities/stats.py

def calc_elastic_net_regularization(
    X: pd.DataFrame,
    y: pd.Series,
    threshold_fraction: float = 0.05
) -> dict:
    """
    Calculate Elastic Net regularization coefficients while iteratively dropping variables with low coefficients.

    This function standardizes `X`, fits an Elastic Net model, and iteratively removes
    variables whose absolute coefficients fall below a fraction of the maximum coefficient.

    Parameters
    ----------
    X : pd.DataFrame
        Input features DataFrame.
    y : pd.Series
        Target variable series.
    threshold_fraction : float, optional
        Fraction of the maximum coefficient below which variables are dropped. Defaults to 0.05.

    Returns
    -------
    dict
        Dictionary with two keys:

        - "initial": pandas.Series of coefficients from the first Elastic Net fit.
        - "final": pandas.Series of coefficients after dropping low-magnitude variables.
    """

    X = X.copy()

    # Impute NaN with column medians before standardization.
    # ElasticNet (sklearn) does not accept NaN natively; median imputation is a neutral
    # choice for this variable-selection pre-pass — LightGBM training still sees real NaN.
    if X.isnull().values.any():
        import warnings
        warnings.warn(
            f"calc_elastic_net_regularization: NaN detected in "
            f"{list(X.columns[X.isnull().any()])}. "
            "Imputing with column medians for the ElasticNet variable-selection step only.",
            UserWarning,
        )
        X = X.fillna(X.median(numeric_only=True))

    # A column with no value at all has no median to fill it with: it can't be tested, so it's dropped (and said),
    # rather than failing the whole fit.
    empty = list(X.columns[X.isnull().all()])
    if empty:
        import warnings
        warnings.warn(
            f"calc_elastic_net_regularization: {empty} have no values at all; left out of the ElasticNet step.",
            UserWarning,
        )
        X = X.drop(columns=empty)

    # Standardize the features
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    first_run = None

    while True:

        # Apply Elastic Net regularization
        elastic_net = ElasticNet(alpha=1.0, l1_ratio=0.5)
        elastic_net.fit(X_scaled, y)

        # Calculate the absolute values of the coefficients
        abs_coefficients = np.abs(elastic_net.coef_)

        # Determine the threshold as a fraction of the largest coefficient
        max_coef = np.max(abs_coefficients)
        threshold = max_coef * threshold_fraction

        coefficients = elastic_net.coef_
        # align coefficients into a dataframe with variable names:
        coefficients = pd.DataFrame(
            {
                "variable": X.columns,
                "enr_coef": coefficients,
                "enr_coef_sign": np.sign(coefficients),
            }
        )
        coefficients = coefficients.sort_values(
            "enr_coef", ascending=False, key=lambda x: x.abs()
        )

        # identify worst variable:
        min_coef_idx = np.argmin(abs_coefficients)
        min_coef = abs_coefficients[min_coef_idx]

        if first_run is None:
            first_run = coefficients

        if min_coef < threshold:
            # remove the worst variable from X_scaled:
            X_scaled = np.delete(X_scaled, min_coef_idx, axis=1)
            # remove corresponding column from X:
            X = X.drop(X.columns[min_coef_idx], axis=1)
        else:
            break

    return {"initial": first_run, "final": coefficients}

calc_median_confidence_interval

calc_median_confidence_interval(values, confidence=0.9)

Confidence interval around a sample MEDIAN, per the IAAO exposure draft.

Implements Appendix D.2 of the Standard on Ratio Studies (Exposure Draft, May 2026): rank the values ascending, take the offset

j = z * sqrt(n) / 2        (n odd)
j = z * sqrt(n) / 2 + 0.5  (n even)

rounded up, then count j ranks up and down the array from the median.

This is a NONPARAMETRIC, order-statistic interval, not a t-interval around the mean. D.2 opens by saying why: "Unlike the mean, the median confidence interval does not depend on the assumption of a normal distribution, and it is not nearly as affected by outlier ratios." Ratio distributions are routinely skewed, so the distinction changes answers rather than decorating them.

Appendix D.4 gives an alternative binomial (Clapp, 1989) table for n <= 30, but the VEI does not use it: E.1 Step 4 points at D.2 specifically, and the worked example in E.4 confirms it -- with quartiles of 13 and 14 sales, the published bounds are the bare order statistics D.2 selects, not D.4's interpolated ones.

Parameters:

Name	Type	Description	Default
values	ndarray	Sample values (e.g. assessment ratios). NaNs are dropped.	required
confidence	float	Confidence level; 0.90 (the VEI's) or 0.95.	0.9

Returns:

Type	Description
tuple[float, float]	(lower, upper), both NaN for an empty sample. The bounds saturate at the smallest and largest observations when j reaches past the ends of the array, which is what a small sample can honestly support.

Raises:

Type	Description
ValueError	If confidence is neither 0.90 nor 0.95.

Source code in openavmkit/utilities/stats.py

def calc_median_confidence_interval(
    values: np.ndarray, confidence: float = 0.90
) -> tuple[float, float]:
    """Confidence interval around a sample MEDIAN, per the IAAO exposure draft.

    Implements Appendix D.2 of the Standard on Ratio Studies (Exposure Draft, May
    2026): rank the values ascending, take the offset

        j = z * sqrt(n) / 2        (n odd)
        j = z * sqrt(n) / 2 + 0.5  (n even)

    rounded up, then count j ranks up and down the array from the median.

    This is a NONPARAMETRIC, order-statistic interval, not a t-interval around the
    mean. D.2 opens by saying why: "Unlike the mean, the median confidence interval
    does not depend on the assumption of a normal distribution, and it is not nearly
    as affected by outlier ratios." Ratio distributions are routinely skewed, so the
    distinction changes answers rather than decorating them.

    Appendix D.4 gives an alternative binomial (Clapp, 1989) table for n <= 30, but
    the VEI does not use it: E.1 Step 4 points at D.2 specifically, and the worked
    example in E.4 confirms it -- with quartiles of 13 and 14 sales, the published
    bounds are the bare order statistics D.2 selects, not D.4's interpolated ones.

    Parameters
    ----------
    values : np.ndarray
        Sample values (e.g. assessment ratios). NaNs are dropped.
    confidence : float, optional
        Confidence level; 0.90 (the VEI's) or 0.95.

    Returns
    -------
    tuple[float, float]
        (lower, upper), both NaN for an empty sample. The bounds saturate at the
        smallest and largest observations when j reaches past the ends of the array,
        which is what a small sample can honestly support.

    Raises
    ------
    ValueError
        If `confidence` is neither 0.90 nor 0.95.
    """
    z = None
    for level, zz in _MEDIAN_CI_Z.items():
        if np.isclose(confidence, level):
            z = zz
            break
    if z is None:
        raise ValueError(
            f"calc_median_confidence_interval supports confidence levels "
            f"{sorted(_MEDIAN_CI_Z)}; got {confidence!r}"
        )

    v = np.asarray(values, dtype=float)
    v = v[~np.isnan(v)]
    n = v.size
    if n == 0:
        return float("nan"), float("nan")

    r = np.sort(v)  # r[0] is rank 1

    j = z * math.sqrt(n) / 2.0
    if n % 2 == 0:
        j += 0.5
    j = math.ceil(j)

    if n % 2 == 1:
        # Median is the single central rank; step j either side of it.
        median_rank = (n + 1) // 2
        lo_rank, hi_rank = median_rank - j, median_rank + j
    else:
        # Median falls between ranks n/2 and n/2 + 1. Counting j ranks outward from
        # each of those two central observations (inclusive) is what reproduces the
        # published example: n=14, j=4 gives ranks 4 and 11.
        lo_rank, hi_rank = (n // 2) - j + 1, (n // 2) + j

    lo_rank = max(1, min(n, lo_rank))
    hi_rank = max(1, min(n, hi_rank))
    return float(r[lo_rank - 1]), float(r[hi_rank - 1])

calc_mki

calc_mki(predictions, ground_truth)

Calculate the Modified Kakwani Index (MKI) of vertical equity.

The concentration index of assessments (their Gini coefficient with the parcels ordered by sale price) over the Gini coefficient of sale prices (IAAO Standard on Ratio Studies, exposure draft May 2026, Appendix F; Quintos 2020). 1 is vertical equity; under 1, regressive; over 1, progressive. Some consider 0.95 to 1.05 optimal for low-variability strata, 0.90 to 1.10 for higher-variability or small ones.

Parameters:

Name	Type	Description	Default
predictions	ndarray	Assessed (or predicted) values.	required
ground_truth	ndarray	Sale prices.	required

Returns:

Type	Description
float	The MKI; NaN with fewer than 2 observations or no spread in prices.

Source code in openavmkit/utilities/stats.py

def calc_mki(predictions: np.ndarray, ground_truth: np.ndarray) -> float:
    """
    Calculate the Modified Kakwani Index (MKI) of vertical equity.

    The concentration index of assessments (their Gini coefficient with the parcels ordered by sale price) over the
    Gini coefficient of sale prices (IAAO Standard on Ratio Studies, exposure draft May 2026, Appendix F; Quintos
    2020). 1 is vertical equity; under 1, regressive; over 1, progressive. Some consider 0.95 to 1.05 optimal for
    low-variability strata, 0.90 to 1.10 for higher-variability or small ones.

    Parameters
    ----------
    predictions : numpy.ndarray
        Assessed (or predicted) values.
    ground_truth : numpy.ndarray
        Sale prices.

    Returns
    -------
    float
        The MKI; NaN with fewer than 2 observations or no spread in prices.
    """
    a = np.asarray(predictions, dtype=float)
    p = np.asarray(ground_truth, dtype=float)
    ok = np.isfinite(a) & np.isfinite(p)
    a, p = a[ok], p[ok]
    n = len(p)
    if n < 2:
        return float("nan")
    order = np.argsort(p, kind="stable")
    a, p = a[order], p[order]
    rank = np.arange(1, n + 1)

    def gini(x: np.ndarray) -> float:
        total = np.sum(x)
        if total == 0:
            return float("nan")
        return (2 * np.sum(x * rank) / total - (n + 1)) / n

    g_price = gini(p)
    if not np.isfinite(g_price) or g_price == 0:
        return float("nan")
    return float(gini(a) / g_price)

calc_mse

calc_mse(prediction, ground_truth)

Calculate the Mean Squared Error (MSE) between predictions and ground truth.

Parameters:

Name	Type	Description	Default
prediction	ndarray	Array of predicted values.	required
ground_truth	ndarray	Array of true values.	required

Returns:

Type	Description
float	The MSE value.

Source code in openavmkit/utilities/stats.py

def calc_mse(
    prediction: np.ndarray,
    ground_truth: np.ndarray
) -> float:
    """
    Calculate the Mean Squared Error (MSE) between predictions and ground truth.

    Parameters
    ----------
    prediction : numpy.ndarray
        Array of predicted values.
    ground_truth : numpy.ndarray
        Array of true values.

    Returns
    -------
    float
        The MSE value.
    """
    return float(np.mean((prediction - ground_truth) ** 2))

calc_mse_r2_adj_r2

calc_mse_r2_adj_r2(predictions, ground_truth, num_vars)

Calculate the Mean Squared Error (MSE), r-squared, and adjusted r-squared

Parameters:

Name	Type	Description	Default
predictions	ndarray	Array of predicted values.	required
ground_truth	ndarray	Array of true values.	required
num_vars	int	Number of independent variables used to produce the predictions	required

Returns:

Type	Description
tuple[float, float, float]	A tuple containing three values: The MSE value The r-squared value The adjusted r-squared value

Source code in openavmkit/utilities/stats.py

def calc_mse_r2_adj_r2(
    predictions: np.ndarray, ground_truth: np.ndarray, num_vars: int
):
    """
    Calculate the Mean Squared Error (MSE), r-squared, and adjusted r-squared

    Parameters
    ----------
    predictions : numpy.ndarray
        Array of predicted values.
    ground_truth : numpy.ndarray
        Array of true values.
    num_vars : int
        Number of independent variables used to produce the predictions

    Returns
    -------
    tuple[float, float, float]

        A tuple containing three values:

        - The MSE value
        - The r-squared value
        - The adjusted r-squared value
    """

    mse = np.mean((ground_truth - predictions) ** 2)
    ss_res = np.sum((ground_truth - predictions) ** 2)
    ss_tot = np.sum((ground_truth - np.mean(ground_truth)) ** 2)

    r2 = 1 - (ss_res / ss_tot) if ss_tot != 0 else float("inf")

    n = len(predictions)
    k = num_vars
    divisor = n - k - 1
    if divisor == 0:
        adj_r2 = float("inf")
    else:
        adj_r2 = 1 - ((1 - r2) * (n - 1) / divisor)
    return mse, r2, adj_r2

calc_p_values_recursive_drop

calc_p_values_recursive_drop(X, y, sig_threshold=0.05)

Recursively drop variables with p-values above a specified significance threshold.

Fits an OLS model on X and iteratively removes the variable with the highest p-value until all remaining variables have p-values below the threshold.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input features DataFrame.	required
y	Series	Target variable series.	required
sig_threshold	float	Significance threshold for p-values. Variables with p-values above this threshold will be dropped. Defaults to 0.05.	0.05

Returns:

Type	Description
dict	Dictionary with two keys: "initial": pd.Series of p-values from the first OLS fit. "final": pd.Series of p-values after recursively dropping high-p-value variables.

Raises:

Type	Description
ValueError	If the OLS regression fails or no variables remain.

Source code in openavmkit/utilities/stats.py

def calc_p_values_recursive_drop(
    X: pd.DataFrame,
    y: pd.Series,
    sig_threshold: float = 0.05
) -> dict:
    """
    Recursively drop variables with p-values above a specified significance threshold.

    Fits an OLS model on `X` and iteratively removes the variable with the highest
    p-value until all remaining variables have p-values below the threshold.

    Parameters
    ----------
    X : pd.DataFrame
        Input features DataFrame.
    y : pd.Series
        Target variable series.
    sig_threshold : float, optional
        Significance threshold for p-values. Variables with p-values above this
        threshold will be dropped. Defaults to 0.05.

    Returns
    -------
    dict
        Dictionary with two keys:

        - "initial": pd.Series of p-values from the first OLS fit.
        - "final": pd.Series of p-values after recursively dropping high-p-value variables.

    Raises
    ------
    ValueError
        If the OLS regression fails or no variables remain.
    """

    X = X.copy()

    # Impute NaN with column medians — statsmodels OLS does not accept NaN natively.
    # This is only for the p-value variable-selection pre-pass; LightGBM training sees real NaN.
    if X.isnull().values.any():
        import warnings
        warnings.warn(
            f"calc_p_values_recursive_drop: NaN detected in "
            f"{list(X.columns[X.isnull().any()])}. "
            "Imputing with column medians for the OLS variable-selection step only.",
            UserWarning,
        )
        X = X.fillna(X.median(numeric_only=True))

    X = sm.add_constant(X, has_constant='add')
    X = X.astype(np.float64)

    model = None
    try:
        model = sm.OLS(y, X).fit()
    except MissingDataError as e:
        warnings.warn(f"OLS failed due to error: {e}")
        return None

    if model is None:
        return None

    first_run = None
    while True:
        max_p_value = model.pvalues.max()
        p_values = model.pvalues
        if first_run is None:
            first_run = p_values
        if max_p_value > sig_threshold:
            var_to_drop = p_values.idxmax()
            if var_to_drop == "const":
                # don't pick const, pick the next variable to drop:
                try:
                    var_to_drop = p_values.iloc[1:].idxmax()
                except ValueError as e:
                    break
            if pd.isna(var_to_drop):
                break
            X = X.drop(var_to_drop, axis=1)
            try:
                new_model = sm.OLS(y, X).fit()
            except ValueError as e:
                print(f"Error fitting model after dropping {var_to_drop}: {e}")
                break
            model = new_model
        else:
            break

    # align p_values into a dataframe with variable names:
    p_values = (
        pd.DataFrame({"p_value": model.pvalues})
        .sort_values("p_value", ascending=True)
        .reset_index()
        .rename(columns={"index": "variable"})
    )

    # do the same for "first_run":
    first_run = (
        pd.DataFrame({"p_value": first_run})
        .sort_values("p_value", ascending=True)
        .reset_index()
        .rename(columns={"index": "variable"})
    )

    return {
        "initial": first_run,
        "final": p_values,
    }

calc_prb

calc_prb(predictions, ground_truth, confidence_interval=0.95)

Calculate the Price Related Bias (PRB) metric using a regression-based approach.

This function fits an OLS model on the transformed ratios of predictions to ground truth, then returns the PRB value along with its lower and upper confidence bounds.

Parameters:

Name	Type	Description	Default
predictions	ndarray	Array of predicted values.	required
ground_truth	ndarray	Array of ground truth values.	required
confidence_interval	float	Desired confidence interval. Defaults to 0.95.	0.95

Returns:

Type	Description
tuple[float, float, float]	A tuple containing: PRB value Lower bound of the confidence interval Upper bound of the confidence interval

Raises:

Type	Description
ValueError	If predictions and ground_truth have different lengths.

Source code in openavmkit/utilities/stats.py

def calc_prb(
    predictions: np.ndarray,
    ground_truth: np.ndarray,
    confidence_interval: float = 0.95
) -> tuple[float, float, float]:
    """
    Calculate the Price Related Bias (PRB) metric using a regression-based approach.

    This function fits an OLS model on the transformed ratios of predictions to ground
    truth, then returns the PRB value along with its lower and upper confidence bounds.

    Parameters
    ----------
    predictions : np.ndarray
        Array of predicted values.
    ground_truth : np.ndarray
        Array of ground truth values.
    confidence_interval : float, optional
        Desired confidence interval. Defaults to 0.95.

    Returns
    -------
    tuple[float, float, float]
        A tuple containing:

        - PRB value
        - Lower bound of the confidence interval
        - Upper bound of the confidence interval

    Raises
    ------
    ValueError
        If `predictions` and `ground_truth` have different lengths.
    """
    # 1. Basic shape checks --------------------------------------------------
    predictions = np.asarray(predictions, dtype=float)
    ground_truth = np.asarray(ground_truth, dtype=float)
    if predictions.shape != ground_truth.shape:
        raise ValueError("predictions and ground_truth must have the same length/shape")

    # 2. Clean rows that cannot be used --------------------------------------
    mask = (
        ~np.isnan(predictions)
        & ~np.isnan(ground_truth)
        & (predictions > 0)          # cannot take log2 of non‑positive numbers
        & (ground_truth > 0)
    )
    n_ok = int(mask.sum())
    if n_ok < 3:                     # OLS needs at least 3 rows
        warnings.warn(
            f"Only {n_ok} valid observation(s) after cleaning – PRB not computed."
        )
        return np.nan, np.nan, np.nan

    preds = predictions[mask]
    truth = ground_truth[mask]

    # 3. Build transformed variables -----------------------------------------
    ratios = preds / truth
    median_ratio = np.median(ratios)

    left = (ratios - median_ratio) / median_ratio
    # IAAO value proxy: the average of the median-ratio-adjusted prediction and the
    # sale price. The 0.5 factor is a constant shift under log2, so it moves the
    # intercept only and leaves PRB (the slope) unchanged.
    right = np.log2((preds / median_ratio + truth) * 0.5)
    right = sm.add_constant(right, has_constant='add')   # adds intercept term

    # 4. Fit model + CI -------------------------------------------------------
    with np.errstate(all="ignore"):  # silence harmless internal numpy warnings
        model = sm.OLS(left, right).fit()

    # add_constant prepends the intercept, so params[0] is the constant and
    # params[1] is the slope. PRB is the SLOPE: the rate at which the percentage
    # deviation from the median ratio changes with value.
    # Guard against degenerate fit (rare but better to be explicit)
    if model.df_resid <= 0 or len(model.params) < 2 or not np.isfinite(model.params[1]):
        return np.nan, np.nan, np.nan

    prb = float(model.params[1])
    prb_lower, prb_upper = (
        model.conf_int(alpha=1.0 - confidence_interval)[1].tolist()
    )

    return prb, prb_lower, prb_upper

calc_prd

calc_prd(predictions, ground_truth)

Calculate the Price Related Differential (PRD).

PRD is computed as the ratio of the mean ratio to the weighted mean ratio of predictions to ground truth.

Parameters:

Name	Type	Description	Default
predictions	ndarray	Array of predicted values.	required
ground_truth	ndarray	Array of ground truth values.	required

Returns:

Type	Description
float	The PRD value.

Source code in openavmkit/utilities/stats.py

def calc_prd(predictions: np.ndarray, ground_truth: np.ndarray) -> float:
    """
    Calculate the Price Related Differential (PRD).

    PRD is computed as the ratio of the mean ratio to the weighted mean ratio of
    predictions to ground truth.

    Parameters
    ----------
    predictions : numpy.ndarray
        Array of predicted values.
    ground_truth : numpy.ndarray
        Array of ground truth values.

    Returns
    -------
    float
        The PRD value.
    """
    ratios = div_series_z_safe(predictions, ground_truth)
    if len(ratios) == 0:
        return float("nan")
    mean_ratio = np.mean(ratios)
    sum_ground_truth = np.sum(ground_truth)
    if sum_ground_truth == 0:
        return float("inf")
    weighted_mean_ratio = np.sum(predictions) / sum_ground_truth
    if weighted_mean_ratio == 0:
        return float("inf")
    prd = mean_ratio / weighted_mean_ratio
    return prd

calc_prd_bootstrap

calc_prd_bootstrap(predictions, ground_truth, confidence_interval=0.95, iterations=10000, seed=777)

Calculate PRD with bootstrapping.

This function bootstraps the prediction-to-ground_truth ratios to produce a distribution of PRD values and returns the lower bound, median, and upper bound of the confidence interval.

Parameters:

Name	Type	Description	Default
predictions	ndarray	Array of predicted values.	required
ground_truth	ndarray	Array of ground truth values.	required
confidence_interval	float	The desired confidence level. Defaults to 0.95.	0.95
iterations	int	Number of bootstrap iterations. Defaults to 10000.	10000
seed	int	Random seed for reproducibility. Defaults to 777.	777

Returns:

Type	Description
tuple[float, float, float]	A tuple containing median PRD, the lower bound, and upper bound of the confidence interval.

Source code in openavmkit/utilities/stats.py

def calc_prd_bootstrap(
    predictions: np.ndarray,
    ground_truth: np.ndarray,
    confidence_interval: float = 0.95,
    iterations: int = 10000,
    seed: int = 777,
) -> tuple[float, float, float]:
    """
    Calculate PRD with bootstrapping.

    This function bootstraps the prediction-to-ground_truth ratios to produce a
    distribution of PRD values and returns the lower bound, median, and upper bound of the
    confidence interval.

    Parameters
    ----------
    predictions : np.ndarray
        Array of predicted values.
    ground_truth : np.ndarray
        Array of ground truth values.
    confidence_interval : float, optional
        The desired confidence level. Defaults to 0.95.
    iterations : int, optional
        Number of bootstrap iterations. Defaults to 10000.
    seed : int, optional
        Random seed for reproducibility. Defaults to 777.

    Returns
    -------
    tuple[float, float, float]
        A tuple containing median PRD, the lower bound, and upper bound of the confidence interval.
    """

    np.random.seed(seed)
    n = len(predictions)
    ratios = predictions / ground_truth
    samples = np.random.choice(ratios, size=(iterations, n), replace=True)
    mean_ratios = np.mean(samples, axis=1)
    weighted_mean_ratios = np.sum(predictions) / np.sum(ground_truth)
    prds = mean_ratios / weighted_mean_ratios
    alpha = (1.0 - confidence_interval) / 2
    lower_bound, upper_bound = np.quantile(prds, [alpha, 1.0 - alpha])
    median_prd = np.median(prds)
    return median_prd, lower_bound, upper_bound

calc_r2

calc_r2(df, variables, y)

Calculate R² and adjusted R² values for each variable.

For each variable in the provided list, an OLS model is fit and the R², adjusted R², and the sign of the coefficient are recorded.

Parameters:

Name	Type	Description	Default
df	DataFrame	DataFrame containing the variables.	required
variables	list[str]	List of variable names to evaluate.	required
y	Series	Target variable series.	required

Returns:

Type	Description
DataFrame	DataFrame with columns for variable, R², adjusted R², and coefficient sign.

Source code in openavmkit/utilities/stats.py

def calc_r2(
    df: pd.DataFrame,
    variables: list[str],
    y: pd.Series
) -> pd.DataFrame:
    """
    Calculate R² and adjusted R² values for each variable.

    For each variable in the provided list, an OLS model is fit and the R², adjusted R²,
    and the sign of the coefficient are recorded.

    Parameters
    ----------
    df : pandas.DataFrame
        DataFrame containing the variables.
    variables : list[str]
        List of variable names to evaluate.
    y : pandas.Series
        Target variable series.

    Returns
    -------
    pandas.DataFrame
        DataFrame with columns for variable, R², adjusted R², and coefficient sign.
    """
    results = {"variable": [], "r2": [], "adj_r2": [], "coef_sign": []}

    for var in variables:
        # Build a joint frame so NA rows are dropped consistently for X and y
        data = pd.concat([df[var], y], axis=1).dropna()

        if len(data) < 3 or data[var].nunique() < 2:
            results["variable"].append(var)
            results["r2"].append(float("nan"))
            results["adj_r2"].append(float("nan"))
            results["coef_sign"].append(float("nan"))
            continue  # skip ill-posed models

        # Skip non-numeric columns (e.g. string categoricals); .astype(float) would raise ValueError
        if not pd.api.types.is_numeric_dtype(data[var]):
            results["variable"].append(var)
            results["r2"].append(float("nan"))
            results["adj_r2"].append(float("nan"))
            results["coef_sign"].append(float("nan"))
            continue

        X = sm.add_constant(data[var].astype(float), has_constant='add')

        # Align y with X using the same filtered rows
        if hasattr(y, "name") and y.name in data.columns:
            y_aligned = data[y.name]
        else:
            y_aligned = data.iloc[:, -1]

        try:
            model = sm.OLS(y_aligned, X).fit()
        except MissingDataError as e:
            print(f"Error fitting model for variable {var}: {e}")
            for v in variables:  # avoid shadowing 'var'
                if df[v].isna().any():
                    n = df[v].isna().sum()
                    print(f'Variable "{v}" has {n} missing values.')
            raise e

        results["variable"].append(var)
        results["r2"].append(model.rsquared)
        results["adj_r2"].append(model.rsquared_adj)
        results["coef_sign"].append(1 if model.params[var] >= 0 else -1)

    df_results = pd.DataFrame(data=results)
    return df_results

calc_ratio_stats_bootstrap

calc_ratio_stats_bootstrap(predictions, ground_truth, confidence_interval=0.95, iterations=10000, seed=777)

Calculate ratio study statistics (Median ratio, Mean ratio, COD, PRD) with bootstrap percentile confidence intervals, following IAAO definitions.

Parameters:

Name	Type	Description	Default
predictions	ndarray	An array of predicted values	required
ground_truth	ndarray	An array of corresponding ground truth (e.g. sale price) values	required
confidence_interval	float	The size of the confidence interval (e.g. 0.95 = 95% confidence)	0.95
iterations	int	The number of bootstrap iterations to perform. Defaults to 10,000.	10000
seed	int	Random seed, for reproducibility. Defaults to 777.	777

Returns:

Type	Description
dict	{ "median_ratio": ConfidenceStat, "mean_ratio": ConfidenceStat, "cod": ConfidenceStat, # COD = 100 * mean(|ri - median(r)|) / median(r) "prd": ConfidenceStat # PRD = mean(r) / weighted_mean(r) }

Source code in openavmkit/utilities/stats.py

def calc_ratio_stats_bootstrap(
    predictions: np.ndarray,
    ground_truth: np.ndarray,
    confidence_interval: float = 0.95,
    iterations: int = 10000,
    seed: int = 777
) -> dict:
    """
    Calculate ratio study statistics (Median ratio, Mean ratio, COD, PRD) with
    bootstrap percentile confidence intervals, following IAAO definitions.

    Parameters
    ----------
    predictions: np.ndarray
        An array of predicted values
    ground_truth: np.ndarray
        An array of corresponding ground truth (e.g. sale price) values
    confidence_interval: float
        The size of the confidence interval (e.g. 0.95 = 95% confidence)
    iterations: int, optional
        The number of bootstrap iterations to perform. Defaults to 10,000.
    seed: int, optional
        Random seed, for reproducibility. Defaults to 777.

    Returns
    -------
    dict
        {
          "median_ratio": ConfidenceStat,
          "mean_ratio": ConfidenceStat,
          "cod": ConfidenceStat,   # COD = 100 * mean(|ri - median(r)|) / median(r)
          "prd": ConfidenceStat    # PRD = mean(r) / weighted_mean(r)
        }
    """
    # --- input hygiene ---
    p = np.asarray(predictions, dtype=float).ravel()
    s = np.asarray(ground_truth,  dtype=float).ravel()
    if p.shape != s.shape:
        raise ValueError("predictions and ground_truth must have the same shape")

    mask = np.isfinite(p) & np.isfinite(s) & (s > 0)  # sales must be > 0
    p, s = p[mask], s[mask]
    n = p.size
    if n == 0:
        return None

    r = p / s

    # Point estimates from the original sample
    med_ratio_point = np.median(r)
    mean_ratio_point = np.mean(r)
    cod_point = (np.mean(np.abs(r - med_ratio_point)) / med_ratio_point) * 100.0
    weighted_mean_ratio_point = p.sum() / s.sum()
    prd_point = mean_ratio_point / weighted_mean_ratio_point

    # Bootstrap on indices of the (p, s) pairs
    rng = np.random.default_rng(seed)
    idx = rng.integers(0, n, size=(iterations, n))

    r_samp = (p[idx] / s[idx])                     # (B, n) ratios per replicate
    med_samp = np.median(r_samp, axis=1)           # (B,)
    mean_samp = np.mean(r_samp, axis=1)            # (B,)

    # per-replicate weighted mean ratio: sum(p)/sum(s) for the bootstrap sample
    wmr_samp = (p[idx].sum(axis=1) / s[idx].sum(axis=1))

    # per-replicate COD uses the replicate's own median
    cod_samp = (np.mean(np.abs(r_samp - med_samp[:, None]), axis=1) / med_samp) * 100.0

    # per-replicate PRD
    prd_samp = mean_samp / wmr_samp

    # percentile confidence intervals
    alpha = (1.0 - confidence_interval) / 2.0
    def ci(a):
        lo, hi = np.quantile(a, [alpha, 1.0 - alpha])
        return float(lo), float(hi)

    med_lo, med_hi   = ci(med_samp)
    mean_lo, mean_hi = ci(mean_samp)
    cod_lo, cod_hi   = ci(cod_samp)
    prd_lo, prd_hi   = ci(prd_samp)

    return {
        "median_ratio": ConfidenceStat(med_ratio_point,  confidence_interval, med_lo,  med_hi),
        "mean_ratio":   ConfidenceStat(mean_ratio_point, confidence_interval, mean_lo, mean_hi),
        "cod":          ConfidenceStat(cod_point,        confidence_interval, cod_lo,  cod_hi),
        "prd":          ConfidenceStat(prd_point,        confidence_interval, prd_lo,  prd_hi),
    }

calc_representation

calc_representation(X_sales, X_univ, do_plots=False)

Calculate how well-represented variables in the sales data frame are in the universe dataframe.

Representation is defined as the percentage of non-zero values, with NaNs treated as zero.

Also reports distribution stats (25th percentile, median, 75th percentile) for each variable in sales and universe.

Parameters:

Name	Type	Description	Default
X_sales	DataFrame	Input DataFrame with variables for the sales dataset	required
X_univ	DataFrame	Input DataFrame with variables for the universe dataset	required
do_plots	bool	If true, plots the results dataframe (Default = False)	False

Returns:

Type	Description
DataFrame	Columns: - name - sales_pct_gt0, univ_pct_gt0 - sales_mean, sales_median - univ_mean, univ_median

Source code in openavmkit/utilities/stats.py

def calc_representation(
    X_sales: pd.DataFrame,
    X_univ: pd.DataFrame,
    do_plots: bool = False
) -> pd.DataFrame:
    """
    Calculate how well-represented variables in the sales data frame are
    in the universe dataframe.

    Representation is defined as the percentage of non-zero values,
    with NaNs treated as zero.

    Also reports distribution stats (25th percentile, median, 75th percentile)
    for each variable in sales and universe.

    Parameters
    ----------
    X_sales : pandas.DataFrame
        Input DataFrame with variables for the sales dataset
    X_univ : pandas.DataFrame
        Input DataFrame with variables for the universe dataset
    do_plots : bool, optional
        If true, plots the results dataframe (Default = False)

    Returns
    -------
    pd.DataFrame
        Columns:
        - name
        - sales_pct_gt0, univ_pct_gt0
        - sales_mean, sales_median
        - univ_mean, univ_median
    """
    # Only compare shared columns
    common_cols = X_sales.columns.intersection(X_univ.columns)

    def _numeric_filled(df: pd.DataFrame) -> pd.DataFrame:
        # Coerce to numeric; treat non-numeric and NaNs as 0
        return df[common_cols].apply(pd.to_numeric, errors="coerce").fillna(0)

    sales_num = _numeric_filled(X_sales)
    univ_num = _numeric_filled(X_univ)

    # % non-zero
    sales_pct = (sales_num != 0).sum(axis=0) / len(sales_num)
    univ_pct = (univ_num != 0).sum(axis=0) / len(univ_num)

    result = (
        pd.DataFrame({
            "name": common_cols,
            "sales_pct_gt0": sales_pct.values,
            "univ_pct_gt0": univ_pct.values,
            "sales_mean": sales_num.mean().values,
            "sales_median": sales_num.median().values,
            "univ_mean": univ_num.mean().values,
            "univ_median": univ_num.median().values,
        })
        .sort_values("name")
        .reset_index(drop=True)
    )

    result = result.sort_values(by="univ_pct_gt0", ascending=False)
    # The shares as numbers, before they're formatted for display.
    raw = result[["name", "sales_pct_gt0", "univ_pct_gt0"]].copy()

    unrepresented_vars = result[
        result["sales_pct_gt0"].le(0) |
        result["univ_pct_gt0"].le(0)
    ]["name"].unique().tolist()

    for key in ["sales_pct_gt0", "univ_pct_gt0"]:
        result[key] = result[key].apply(lambda x: f"{x:.2%}")
    for key in ["sales_mean", "sales_median", "univ_mean", "univ_median"]:
        result[key] = result[key].apply(lambda x: _fancy_format(x))

    if len(unrepresented_vars) > 0:
        print("")
        print(f"Found {len(unrepresented_vars)} variables with 0 representation in either sales or universe dataset. You should NOT model with these:")
        for col in unrepresented_vars:
            print(col)
        print("")

    if do_plots:
        pd.set_option('display.max_rows', None)
        display(result)
        pd.set_option('display.max_rows', 15)

    return {
        "bad_vars": unrepresented_vars,
        "df": result,
        "raw": raw,
   }

calc_t_values

calc_t_values(X, y)

Calculate t-values for an OLS model.

Fits an ordinary least squares regression of y on X and returns the t-values of the estimated coefficients.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input features DataFrame (should include a constant term column).	required
y	Series	Target variable series.	required

Returns:

Type	Description
Series	Series of t-values corresponding to each coefficient in X.

Source code in openavmkit/utilities/stats.py

def calc_t_values(X: pd.DataFrame, y: pd.Series) -> pd.Series:
    """
    Calculate t-values for an OLS model.

    Fits an ordinary least squares regression of `y` on `X` and returns the t-values
    of the estimated coefficients.

    Parameters
    ----------
    X : pandas.DataFrame
        Input features DataFrame (should include a constant term column).
    y : pandas.Series
        Target variable series.

    Returns
    -------
    pandas.Series
        Series of t-values corresponding to each coefficient in `X`.
    """
    linear_model = sm.OLS(y, X)
    fitted_model = linear_model.fit()
    return fitted_model.tvalues

calc_t_values_recursive_drop

calc_t_values_recursive_drop(X, y, threshold=2)

Recursively drop variables with t-values below a given threshold.

Fits an OLS model on X and iteratively removes the variable with the smallest absolute t-value until all remaining variables have |t-value| above the threshold.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input features DataFrame.	required
y	Series	Target variable series.	required
threshold	float	Minimum acceptable absolute t-value. Variables with |t-value| below this threshold will be dropped. Defaults to 2.	2

Returns:

Type	Description
dict	Dictionary with two keys: "initial": pandas.Series of t-values and their signs from the first OLS fit. "final": pandas.Series of t-values and their signs after recursive dropping.

Source code in openavmkit/utilities/stats.py

def calc_t_values_recursive_drop(
    X: pd.DataFrame,
    y: pd.Series,
    threshold: float = 2
) -> dict:
    """
    Recursively drop variables with t-values below a given threshold.

    Fits an OLS model on `X` and iteratively removes the variable with the smallest
    absolute t-value until all remaining variables have |t-value| above the threshold.

    Parameters
    ----------
    X : pandas.DataFrame
        Input features DataFrame.
    y : pandas.Series
        Target variable series.
    threshold : float, optional
        Minimum acceptable absolute t-value. Variables with |t-value| below this threshold
        will be dropped. Defaults to 2.

    Returns
    -------
    dict
        Dictionary with two keys:

        - "initial": pandas.Series of t-values and their signs from the first OLS fit.
        - "final": pandas.Series of t-values and their signs after recursive dropping.
    """

    X = X.copy()

    # Impute NaN with column medians — statsmodels OLS does not accept NaN natively.
    # This is only for the t-value variable-selection pre-pass; LightGBM training sees real NaN.
    if X.isnull().values.any():
        import warnings
        warnings.warn(
            f"calc_t_values_recursive_drop: NaN detected in "
            f"{list(X.columns[X.isnull().any()])}. "
            "Imputing with column medians for the OLS variable-selection step only.",
            UserWarning,
        )
        X = X.fillna(X.median(numeric_only=True))

    X = sm.add_constant(X, has_constant='add')
    X = X.astype(np.float64)

    first_run = None
    i = 0
    while True:
        i += 1
        try:
            t_values = calc_t_values(X, y)
        except (ValueError, MissingDataError) as e:
            t_values = pd.Series([float("nan")] * X.shape[1], index=X.columns)
        if first_run is None:
            first_run = t_values

        if t_values.isna().all():
            min_t_var = float("nan")
        else:
            min_t_var = t_values.abs().idxmin()

        if pd.isna(min_t_var):
            min_t_var = 0

        if len(t_values) > 0:
            min_t_val = float("nan")
        else:
            min_t_val = t_values[min_t_var]

        if min_t_val < threshold:
            X = X.drop(min_t_var, axis=1)
        else:
            break

    # align t_values into a dataframe with variable names:
    t_values = (
        pd.DataFrame({"t_value": t_values, "t_value_sign": np.sign(t_values)})
        .sort_values("t_value", ascending=False, key=lambda x: x.abs())
        .reset_index()
        .rename(columns={"index": "variable"})
    )

    # do the same for "first_run":
    first_run = (
        pd.DataFrame({"t_value": first_run, "t_value_sign": np.sign(first_run)})
        .sort_values("t_value", ascending=False, key=lambda x: x.abs())
        .reset_index()
        .rename(columns={"index": "variable"})
    )

    return {"initial": first_run, "final": t_values}

calc_vif

calc_vif(X)

Calculate the Variance Inflation Factor (VIF) for each variable in a DataFrame.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input features DataFrame.	required

Returns:

Type	Description
DataFrame	DataFrame with columns: "variable": Name of each feature in X. "VIF": Variance Inflation Factor value for that feature.

Source code in openavmkit/utilities/stats.py

def calc_vif(X: pd.DataFrame) -> pd.DataFrame:
    """
    Calculate the Variance Inflation Factor (VIF) for each variable in a DataFrame.

    Parameters
    ----------
    X : pandas.DataFrame
        Input features DataFrame.

    Returns
    -------
    pandas.DataFrame
        DataFrame with columns:

        - "variable": Name of each feature in `X`.
        - "VIF": Variance Inflation Factor value for that feature.
    """
    vif_data = pd.DataFrame()
    vif_data["variable"] = X.columns

    if len(X.values) < 5:
        warnings.warn("Can't calculate VIF for less than 5 samples")
        vif_data["vif"] = [float("nan")] * len(X.columns)
        return vif_data

    if len(X.columns) == 1:
        warnings.warn("Can't calculate VIF for one column")
        vif_data["vif"] = [float("nan")] * len(X.columns)
        return vif_data

    # Calculate VIF for each column
    vif_data["vif"] = [
        variance_inflation_factor(X.values, i) for i in range(X.shape[1])
    ]

    return vif_data

calc_vif_recursive_drop

calc_vif_recursive_drop(X, threshold=10.0, settings=None)

Recursively drop variables with a Variance Inflation Factor (VIF) exceeding the threshold.

Parameters:

Name	Type	Description	Default
X	DataFrame	Input features DataFrame.	required
threshold	float	Maximum acceptable VIF. Variables with VIF above this threshold will be removed. Defaults to 10.0.	10.0
settings	dict	Settings dictionary containing field classifications, if needed for VIF computation. Defaults to None.	None

Returns:

Type	Description
dict	Dictionary with two keys: "initial": pandas.DataFrame of VIF values before dropping variables. "final": pandas.DataFrame of VIF values after recursively dropping high-VIF variables.

Raises:

Type	Description
ValueError	If no columns remain for VIF calculation.

Source code in openavmkit/utilities/stats.py

def calc_vif_recursive_drop(
    X: pd.DataFrame,
    threshold: float = 10.0,
    settings: dict = None
) -> dict:
    """
    Recursively drop variables with a Variance Inflation Factor (VIF) exceeding the threshold.

    Parameters
    ----------
    X : pandas.DataFrame
        Input features DataFrame.
    threshold : float, optional
        Maximum acceptable VIF. Variables with VIF above this threshold will be removed.
        Defaults to 10.0.
    settings : dict, optional
        Settings dictionary containing field classifications, if needed for VIF computation.
        Defaults to None.

    Returns
    -------
    dict
        Dictionary with two keys:

        - "initial": pandas.DataFrame of VIF values before dropping variables.
        - "final": pandas.DataFrame of VIF values after recursively dropping high-VIF variables.

    Raises
    ------
    ValueError
        If no columns remain for VIF calculation.
    """
    X = X.copy()

    # Impute NaN with column medians — VIF (OLS-based) does not accept NaN natively.
    # This is only for the VIF variable-selection pre-pass; LightGBM training sees real NaN.
    if X.isnull().values.any():
        import warnings
        warnings.warn(
            f"calc_vif_recursive_drop: NaN detected in "
            f"{list(X.columns[X.isnull().any()])}. "
            "Imputing with column medians for the VIF variable-selection step only.",
            UserWarning,
        )
        X = X.fillna(X.median(numeric_only=True))

    X = X.astype(np.float64)

    # Get boolean and categorical variables from settings if provided
    bool_fields = get_fields_boolean(settings, X)
    cat_fields = get_fields_categorical(settings, X, include_boolean=False)
    exclude_vars = bool_fields + cat_fields

    # Remove boolean and categorical variables
    X = X.drop(
        columns=[col for col in X.columns if col in exclude_vars], errors="ignore"
    )

    # Drop constant columns (VIF cannot be calculated for constant columns)
    X = X.loc[:, X.nunique() > 1]  # Keep only columns with more than one unique value

    # If no columns are left after removing constant columns or dropping NaN values, raise an error
    if X.shape[1] == 0:
        raise ValueError(
            "All columns are constant or have missing values; VIF cannot be computed."
        )
    first_run = None
    while True:
        vif_data = calc_vif(X)
        if first_run is None:
            first_run = vif_data
        if vif_data["vif"].max() > threshold:
            max_vif_idx = vif_data["vif"].idxmax()
            X = X.drop(X.columns[max_vif_idx], axis=1)
        else:
            break
    return {"initial": first_run, "final": vif_data}

plot_correlation

plot_correlation(corr, title='Correlation of Variables')

Plot a heatmap of a correlation matrix with adaptive formatting: - Shrinks axis label fonts as the matrix grows - Hides cell annotations when the plot becomes too dense

Source code in openavmkit/utilities/stats.py

def plot_correlation(corr: pd.DataFrame, title: str = "Correlation of Variables") -> None:
    """
    Plot a heatmap of a correlation matrix with adaptive formatting:
    - Shrinks axis label fonts as the matrix grows
    - Hides cell annotations when the plot becomes too dense
    """

    n = corr.shape[0]

    # Heuristics for readability
    if n <= 10:
        annot = True
        annot_size = 8
        label_size = 10
        fig_size = (10, 8)
    elif n <= 20:
        annot = True
        annot_size = 6
        label_size = 8
        fig_size = (12, 10)
    else:
        annot = False          # Hide numbers in cells
        annot_size = None
        label_size = 6
        fig_size = (14, 12)

    cmap = sns.diverging_palette(220, 10, as_cmap=True).reversed()

    plt.figure(figsize=fig_size)

    sns.heatmap(
        corr,
        annot=annot,
        fmt=".1f",
        cbar=True,
        cmap=cmap,
        vmax=1.0,
        vmin=-1.0,
        xticklabels=corr.columns.tolist(),
        yticklabels=corr.index.tolist(),
        annot_kws={"size": annot_size} if annot else None,
    )

    plt.title(title)

    # Axis label formatting
    plt.xticks(rotation=45, ha="right", fontsize=label_size)
    plt.yticks(rotation=0, fontsize=label_size)

    plt.tight_layout(pad=2)
    plt.show()

quick_median_chd

quick_median_chd(df, field_value, field_cluster)

Calculate the median CHD across the clusters of a DataFrame.

This function filters out missing and non-positive values for the given value field, groups the rest by the specified cluster field, computes each cluster's COD, and returns the median of those CODs. Rows with a missing cluster ID form a cluster of their own.

Vectorised: every cluster's median and mean absolute deviation come from one groupby, rather than a loop over the clusters.

Parameters:

Name	Type	Description	Default
df	DataFrame	Input DataFrame (only field_value and field_cluster are read).	required
field_value	str	Name of the field containing values for COD calculation.	required
field_cluster	str	Name of the field to group by for computing COD.	required

Returns:

Type	Description
float	The median COD value across all clusters (NaN if no value is left).

Source code in openavmkit/utilities/stats.py

def quick_median_chd(df: pd.DataFrame, field_value: str, field_cluster: str) -> float:
    """
    Calculate the median CHD across the clusters of a DataFrame.

    This function filters out missing and non-positive values for the given value field,
    groups the rest by the specified cluster field, computes each cluster's COD, and
    returns the median of those CODs. Rows with a missing cluster ID form a cluster of
    their own.

    Vectorised: every cluster's median and mean absolute deviation come from one groupby,
    rather than a loop over the clusters.

    Parameters
    ----------
    df : pd.DataFrame
        Input DataFrame (only ``field_value`` and ``field_cluster`` are read).
    field_value : str
        Name of the field containing values for COD calculation.
    field_cluster : str
        Name of the field to group by for computing COD.

    Returns
    -------
    float
        The median COD value across all clusters (NaN if no value is left).
    """
    sub = df[[field_cluster, field_value]]
    sub = sub[sub[field_value].notna() & sub[field_value].gt(0)]
    if sub.empty:
        return float("nan")

    def by_cluster(frame: pd.DataFrame):
        return frame.groupby(field_cluster, sort=False, dropna=False, observed=True)

    # Each value's distance from its cluster's median; COD = mean distance / median * 100.
    # Values are all positive here, so no cluster's median is zero (calc_cod's special cases).
    median = by_cluster(sub)[field_value].transform("median")
    sub = sub.assign(_abs_dev=(sub[field_value] - median).abs())
    per_cluster = by_cluster(sub).agg(
        _mean_abs_dev=("_abs_dev", "mean"), _median=(field_value, "median")
    )
    chds = per_cluster["_mean_abs_dev"] / per_cluster["_median"] * 100
    return float(np.median(chds.to_numpy()))

trim_outliers

trim_outliers(values, max_percent=0.1, iqr_factor=1.5)

Trim outliers using IQR fences per IAAO guidance, with a max trim cap. Fails immediately if NaNs are detected.

1) Compute Q1, Q3, IQR = Q3 - Q1. 2) Trim values outside [Q1 - iqr_factorIQR, Q3 + iqr_factorIQR]. 3) If more than max_percent would be removed, instead trim by symmetric quantile cut so total trimmed <= max_percent.

Parameters:

Name	Type	Description	Default
values	ndarray	1D numeric array with no NaNs allowed.	required
max_percent	float	Maximum fraction to remove (e.g., 0.10 = 10%).	0.1
iqr_factor	float	1.5 for standard outliers, 3.0 for extreme outliers.	1.5

Returns:

Type	Description
ndarray	Trimmed array according to the above rules.

Raises:

Type	Description
ValueError	If any NaN is detected in values.

Source code in openavmkit/utilities/stats.py

def trim_outliers(
    values: np.ndarray,
    max_percent: float = 0.10,
    iqr_factor: float = 1.5
) -> np.ndarray:
    """
    Trim outliers using IQR fences per IAAO guidance, with a max trim cap.
    Fails immediately if NaNs are detected.

    1) Compute Q1, Q3, IQR = Q3 - Q1.
    2) Trim values outside [Q1 - iqr_factor*IQR, Q3 + iqr_factor*IQR].
    3) If more than max_percent would be removed, instead trim by symmetric
       quantile cut so total trimmed <= max_percent.

    Parameters
    ----------
    values : np.ndarray
        1D numeric array with no NaNs allowed.
    max_percent : float, optional
        Maximum fraction to remove (e.g., 0.10 = 10%).
    iqr_factor : float, optional
        1.5 for standard outliers, 3.0 for extreme outliers.

    Returns
    -------
    np.ndarray
        Trimmed array according to the above rules.

    Raises
    ------
    ValueError
        If any NaN is detected in `values`.
    """
    if np.isnan(values).any():
        raise ValueError("NaN values detected — remove or impute them before trimming.")

    if values.size == 0:
        return values

    # IQR fences
    q1, q3 = np.quantile(values, [0.25, 0.75])
    iqr = q3 - q1
    lower_fence = q1 - iqr_factor * iqr
    upper_fence = q3 + iqr_factor * iqr

    within_fences = (values >= lower_fence) & (values <= upper_fence)
    trimmed_iqr = values[within_fences]

    n = values.size
    max_trim = int(np.floor(n * max_percent))

    # If IQR-based trimming exceeds cap, fall back to symmetric quantile trim
    trimmed_count = n - trimmed_iqr.size
    if trimmed_count > max_trim:
        tail_fraction = max_percent / 2.0
        low_q, high_q = tail_fraction, 1.0 - tail_fraction
        lo, hi = np.quantile(values, [low_q, high_q])
        keep = (values >= lo) & (values <= hi)
        trimmed = values[keep]
    else:
        trimmed = trimmed_iqr

    return trimmed

trim_outliers_mask

trim_outliers_mask(values, max_percent=0.1, iqr_factor=1.5)

Trim outliers using IQR fences per IAAO guidance, with a max trim cap. Fails immediately if NaNs are detected.

1) Compute Q1, Q3, IQR = Q3 - Q1. 2) Trim values outside [Q1 - iqr_factorIQR, Q3 + iqr_factorIQR]. 3) If more than max_percent would be removed, instead trim by symmetric quantile cut so total trimmed <= max_percent.

Parameters:

Name	Type	Description	Default
values	ndarray	1D numeric array with no NaNs allowed.	required
max_percent	float	Maximum fraction to remove (e.g., 0.10 = 10%).	0.1
iqr_factor	float	1.5 for standard outliers, 3.0 for extreme outliers.	1.5

Returns:

Type	Description
ndarray	Boolean array where True indicates values within the quantile bounds.

Raises:

Type	Description
ValueError	If any NaN is detected in values.

Source code in openavmkit/utilities/stats.py

def trim_outliers_mask(
    values: np.ndarray,
    max_percent: float = 0.10,
    iqr_factor: float = 1.5,
) -> np.ndarray:
    """
    Trim outliers using IQR fences per IAAO guidance, with a max trim cap.
    Fails immediately if NaNs are detected.

    1) Compute Q1, Q3, IQR = Q3 - Q1.
    2) Trim values outside [Q1 - iqr_factor*IQR, Q3 + iqr_factor*IQR].
    3) If more than max_percent would be removed, instead trim by symmetric
       quantile cut so total trimmed <= max_percent.

    Parameters
    ----------
    values : np.ndarray
        1D numeric array with no NaNs allowed.
    max_percent : float, optional
        Maximum fraction to remove (e.g., 0.10 = 10%).
    iqr_factor : float, optional
        1.5 for standard outliers, 3.0 for extreme outliers.

    Returns
    -------
    np.ndarray
        Boolean array where `True` indicates values within the quantile bounds.

    Raises
    ------
    ValueError
        If any NaN is detected in `values`.
    """
    if np.isnan(values).any():
        raise ValueError("NaN values detected — remove or impute them before trimming.")

    if values.size == 0:
        return values

    # IQR fences
    q1, q3 = np.quantile(values, [0.25, 0.75])
    iqr = q3 - q1
    lower_fence = q1 - iqr_factor * iqr
    upper_fence = q3 + iqr_factor * iqr

    within_fences = (values >= lower_fence) & (values <= upper_fence)
    trimmed_iqr = values[within_fences]

    n = values.size
    max_trim = int(np.floor(n * max_percent))

    # If IQR-based trimming exceeds cap, fall back to symmetric quantile trim
    trimmed_count = n - trimmed_iqr.size
    if trimmed_count > max_trim:
        tail_fraction = max_percent / 2.0
        low_q, high_q = tail_fraction, 1.0 - tail_fraction
        lo, hi = np.quantile(values, [low_q, high_q])
        return (values >= lo) & (values <= hi)
    else:
        return within_fences