diff --git a/src/model_drift/data/mgb_data.py b/src/model_drift/data/mgb_data.py new file mode 100644 index 0000000..0eb3793 --- /dev/null +++ b/src/model_drift/data/mgb_data.py @@ -0,0 +1,42 @@ +import datetime + + + +TRAIN_DATE_END = datetime.datetime(year=2019, month=10, day=1) +VAL_DATE_END = datetime.datetime(year=2020, month=1, day=1) +# shorter timeframe to speed up validation during training, measeurements start (year=2019, month=7, day=3) +#VAL_DATE_END = datetime.datetime(year=2019, month=10, day=21) + +RAW_LABELS = [ + 'No Finding', + 'Enlarged Cardiomediastinum', + 'Cardiomegaly', + 'Lung Lesion', + 'Lung Opacity', + 'Edema', + 'Consolidation', + 'Pneumonia', + 'Atelectasis', + 'Pneumothorax', + 'Pleural Effusion', + 'Pleural Other', + 'Fracture', + 'Support Devices', +] + +LABEL_GROUPINGS = { + 'Atelectasis': ['Atelectasis'], + 'Cardiomegaly': ['Cardiomegaly'], + #'Consolidation': ['Consolidation'], + 'Edema': ['Edema'], + #'Lung Lesion': ['Lung Lesion'], + # 'No Finding': ['No Finding'], + 'Lung Opacity': ['Lung Opacity', 'Pneumonia', 'Consolidation', 'Lung Lesion', 'Atelectasis', 'Edema'], + 'Pleural Other': ['Pleural Other'], + 'Pleural Effusion': ['Pleural Effusion'], + 'Pneumonia': ['Pneumonia'], + 'Pneumothorax': ['Pneumothorax'], + 'Support Devices': ['Support Devices'], + #'Enlarged Cardiomediastinum': ['Enlarged Cardiomediastinum', 'Cardiomegaly'], + #'Fracture': ['Fracture'], +} \ No newline at end of file diff --git a/src/model_drift/drift/__init__.py b/src/model_drift/drift/__init__.py index 9caded5..b1d2db9 100644 --- a/src/model_drift/drift/__init__.py +++ b/src/model_drift/drift/__init__.py @@ -3,10 +3,10 @@ # Licensed under the MIT License (MIT). See LICENSE in the repo root for license information. # ------------------------------------------------------------------------------------------ from .base import BaseDriftCalculator # noqa -from .categorical import ChiSqDriftCalculator # noqa +from .categorical import ChiSqDriftCalculator, ChiSqDriftCalculatorJackKnife, HellingerDriftCalculatorJackKnife # noqa from .histogram import HistIntersectionCalculator, KdeHistPlotCalculator # noqa from .collection import DriftCollectionCalculator # noqa -from .numeric import KSDriftCalculator, BasicDriftCalculator # noqa +from .numeric import KSDriftCalculator, BasicDriftCalculator, KSDriftCalculatorJackKnife, EMDDriftCalculatorJackKnife_1D, EMDDriftCalculatorJackKnife_woRef # noqa from .performance import AUROCCalculator, ClassificationReportCalculator # noqa from .sampler import Sampler # noqa -from .tabular import TabularDriftCalculator # noqa +from .tabular import TabularDriftCalculator # noqa \ No newline at end of file diff --git a/src/model_drift/drift/categorical.py b/src/model_drift/drift/categorical.py index 37beb98..f80010a 100644 --- a/src/model_drift/drift/categorical.py +++ b/src/model_drift/drift/categorical.py @@ -3,13 +3,13 @@ # Licensed under the MIT License (MIT). See LICENSE in the repo root for license information. # ------------------------------------------------------------------------------------------ from collections import Counter +import math import numpy as np from scipy.stats import chi2_contingency, chi2 from model_drift.drift.base import BaseDriftCalculator - def merge_freqs(ref_counts, sample): sample_counts = Counter(sample) keys = set().union(ref_counts, sample_counts) @@ -17,6 +17,11 @@ def merge_freqs(ref_counts, sample): obs = np.array([sample_counts.get(k, 0) for k in keys]) return exp, keys, obs +def hellinger_distance(p, q): + """Hellinger distance between two discrete distributions. + """ + return math.sqrt(sum([ (math.sqrt(p_i) - math.sqrt(q_i))**2 for p_i, q_i in zip(p, q) ]) / 2) + class ChiSqDriftCalculator(BaseDriftCalculator): name = "chi2" @@ -53,3 +58,122 @@ def _predict(self, sample): out['critical_diff'] = out['distance'] - out['critical_value'] return out + + + +class ChiSqDriftCalculatorJackKnife(BaseDriftCalculator): + name = "chi2_jackknife" + + def __init__(self, q_val=0.1, correction=True, lambda_=None, use_freq=False, include_critical_values=False, + **kwargs): + super().__init__(**kwargs) + self.q_val = q_val + self.correction = correction + self.lambda_ = lambda_ + self.use_freq = use_freq + self.include_critical_values = include_critical_values + + def convert(self, arg): + return arg.apply(str) + + def prepare(self, ref, **kwargs): + self._ref_counts = Counter(ref) + super().prepare(ref) + + def _predict(self, sample): + + nref = len(self._ref) + nobs = len(sample) + + ref1 = np.random.choice(self._ref, nobs) + ref2 = np.random.choice(self._ref, nobs) + + ref1_counts = Counter(ref1) + ref2_counts = Counter(ref2) + + exp_ref1_sam, keys, obs_ref1_sam = merge_freqs(ref1_counts, sample) + + exp_ref1_ref2, keys, obs_ref1_ref2 = merge_freqs(ref1_counts, ref2_counts) + + if self.use_freq: + exp_ref1_sam = exp_ref1_sam / exp_ref1_sam.sum() + obs_ref1_sam = obs_ref1_sam / obs_ref1_sam.sum() + + exp_ref1_ref2 = exp_ref1_ref2 / exp_ref1_ref2.sum() + obs_ref1_ref2 = obs_ref1_ref2 / obs_ref1_ref2.sum() + + + out = {} + + dist1, _, _, _ = chi2_contingency(np.vstack([exp_ref1_sam, obs_ref1_sam]), + correction=self.correction, + lambda_=self.lambda_) + + dist2, _, _, _ = chi2_contingency(np.vstack([exp_ref1_ref2, obs_ref1_ref2]), + correction=self.correction, + lambda_=self.lambda_) + + out["distance"] = max(dist1 - dist2, 0.0) + out['pval'] = float("NaN") + + + if self.include_critical_values: + raise NotImplementedError("Critical value not implemented for jackknife") + + return out + +class HellingerDriftCalculatorJackKnife(BaseDriftCalculator): + name = "hellinger_jackknife" + + def __init__(self, q_val=0.1, correction=True, lambda_=None, use_freq=True, include_critical_values=False, + **kwargs): + super().__init__(**kwargs) + self.q_val = q_val + self.correction = correction + self.lambda_ = lambda_ + self.use_freq = use_freq + self.include_critical_values = include_critical_values + + def convert(self, arg): + return arg.apply(str) + + def prepare(self, ref, **kwargs): + self._ref_counts = Counter(ref) + super().prepare(ref) + + def _predict(self, sample): + + nref = len(self._ref) + nobs = len(sample) + + ref1 = np.random.choice(self._ref, nobs) + ref2 = np.random.choice(self._ref, nobs) + + ref1_counts = Counter(ref1) + ref2_counts = Counter(ref2) + + exp_ref1_sam, keys, obs_ref1_sam = merge_freqs(ref1_counts, sample) + + exp_ref1_ref2, keys, obs_ref1_ref2 = merge_freqs(ref1_counts, ref2_counts) + + if self.use_freq: + exp_ref1_sam = exp_ref1_sam / exp_ref1_sam.sum() + obs_ref1_sam = obs_ref1_sam / obs_ref1_sam.sum() + + exp_ref1_ref2 = exp_ref1_ref2 / exp_ref1_ref2.sum() + obs_ref1_ref2 = obs_ref1_ref2 / obs_ref1_ref2.sum() + + out = {} + + dist1 = hellinger_distance(exp_ref1_sam, obs_ref1_sam) + dist2 = hellinger_distance(exp_ref1_ref2, obs_ref1_ref2) + + + out["distance"] = max(dist1 - dist2, 0.0) + out['pval'] = float("NaN") + + + if self.include_critical_values: + raise NotImplementedError("Critical value not implemented for jackknife") + + return out \ No newline at end of file diff --git a/src/model_drift/drift/numeric.py b/src/model_drift/drift/numeric.py index c10d871..53bf5c1 100644 --- a/src/model_drift/drift/numeric.py +++ b/src/model_drift/drift/numeric.py @@ -5,7 +5,8 @@ import numpy as np import pandas as pd from scipy.special import kolmogi -from scipy.stats import ks_2samp +from scipy.stats import ks_2samp, wasserstein_distance +import ot from model_drift.drift.base import BaseDriftCalculator @@ -48,6 +49,131 @@ def calc_critical_value(n1, n2, q=.01): return kolmogi(q) * np.sqrt((n1 + n2) / (n1 * n2)) +class KSDriftCalculatorJackKnife(NumericBaseDriftCalculator): + name = "ks_jackknife" + + def __init__(self, q_val=0.1, alternative='two-sided', mode='asymp', average='macro', include_critical_value=False, + **kwargs): + super().__init__(**kwargs) + self.q_val = q_val + self.alternative = alternative + self.mode = mode + self.average = average + self.include_critical_value = include_critical_value + + def _predict(self, sample): + nref = len(self._ref) + nobs = len(sample) + + ref1 = np.random.choice(self._ref, nobs) + ref2 = np.random.choice(self._ref, nobs) + out = {} + try: + dist1, _ = ks_2samp(ref1, sample, alternative=self.alternative, mode=self.mode) + dist2, _ = ks_2samp(ref1, ref2, alternative=self.alternative, mode=self.mode) + + out["distance"] = max(dist1 - dist2, 0.0) + out['pval'] = float("NaN") + + except TypeError: + out["distance"], out['pval'] = float("NaN"), float("NaN") + + if self.include_critical_value: + raise NotImplementedError("Critical value not implemented for jackknife") + + return out + + +class EMDDriftCalculatorJackKnife_woRef(NumericBaseDriftCalculator): + name = "emd_jackknife" + + def __init__(self, include_critical_value=False, **kwargs): + super().__init__(**kwargs) + self.include_critical_value = include_critical_value + + def convert(self, arg): + return arg + + def _predict(self, sample): + + def _emd_distance(tar, ref): + a = np.ones((len(ref))) / len(ref) # Uniform weights for the reference set + b = np.ones((len(tar))) / len(tar) # Uniform weights for the target set + M = ot.dist(ref, tar) + G0 = ot.emd(a, b, M, numItermax=1000000) + em_distance = np.sum(M * G0) + return em_distance + + nref = len(self._ref) + nobs = len(sample) + + sample_tuples = sample.apply(tuple) + ref_tuples = self._ref.apply(tuple) + ref_tuples_exclusive = ref_tuples[~ref_tuples.isin(sample_tuples)] + ref_lists_exclusive = ref_tuples_exclusive.apply(list) + + ref1 = np.random.choice(ref_lists_exclusive, nobs) + ref2 = np.random.choice(ref_lists_exclusive, nobs) + + sample_arr = np.array(sample.tolist()) + sample_arr = np.nan_to_num(sample_arr, nan=0.0, posinf=0.0, neginf=0.0) + + ref1_arr = np.array(ref1.tolist()) + ref1_arr = np.nan_to_num(ref1_arr, nan=0.0, posinf=0.0, neginf=0.0) + + ref2_arr = np.array(ref2.tolist()) + ref2_arr = np.nan_to_num(ref2_arr, nan=0.0, posinf=0.0, neginf=0.0) + + out = {} + try: + dist1 = _emd_distance(ref1_arr, sample_arr) + dist2 = _emd_distance(ref1_arr, ref2_arr) + + out["distance"] = max(dist1 - dist2, 0.0) + out['pval'] = float("NaN") + + except TypeError: + out["distance"], out['pval'] = float("NaN"), float("NaN") + + if self.include_critical_value: + raise NotImplementedError("Critical value not implemented for jackknife") + + return out + +class EMDDriftCalculatorJackKnife_1D(NumericBaseDriftCalculator): + name = "emd_jackknife_1d" + + def __init__(self, include_critical_value=False, **kwargs): + super().__init__(**kwargs) + self.include_critical_value = include_critical_value + + def _predict(self, sample): + nref = len(self._ref) + nobs = len(sample) + + ref1 = np.random.choice(self._ref, nobs) + ref2 = np.random.choice(self._ref, nobs) + out = {} + # drop NaNs from the arrays before computing the EMD + ref1 = ref1[~np.isnan(ref1)] + ref2 = ref2[~np.isnan(ref2)] + sample = sample[~np.isnan(sample)] + try: + dist1 = wasserstein_distance(ref1, sample) + dist2 = wasserstein_distance(ref1, ref2) + + out["distance"] = max(dist1 - dist2, 0.0) + out['pval'] = float("NaN") + + except TypeError: + out["distance"], out['pval'] = float("NaN"), float("NaN") + + if self.include_critical_value: + raise NotImplementedError("Critical value not implemented for jackknife") + + return out + + class BasicDriftCalculator(NumericBaseDriftCalculator): name = "stats" diff --git a/src/model_drift/drift/performance.py b/src/model_drift/drift/performance.py index 0c10668..46a42a2 100644 --- a/src/model_drift/drift/performance.py +++ b/src/model_drift/drift/performance.py @@ -47,14 +47,42 @@ def macro_auc(scores, labels, skip_missing=True): def micro_auc(scores, labels): return float(auroc(torch.tensor(scores), torch.tensor(labels).long(), average='micro').numpy()) +def youden_point(df, target_names): + """ + Calculate the Youden Index for each disease and return that as optimal operating point for the disease. + """ + operating_points = {} + for disease in target_names: + fpr, tpr, thresholds = metrics.roc_curve(df[f'label.{disease}'], df[f'activation.{disease}']) -def classification_report(scores, labels, target_names=None, th=0.5, ): + # Calculate Youden Index + idx = np.argmax(tpr - fpr) + best_cutoff = thresholds[idx] + + operating_points[disease] = format(best_cutoff, ".4f") + + return operating_points + + +def classification_report(scores, labels, target_names=None, th=None): keeps = (labels.sum(axis=0) > 0) if target_names is None: target_names = [str(i) for i in range(scores.shape[1])] target_names = np.array(target_names) - output = metrics.classification_report(labels, scores >= th, target_names=target_names, output_dict=True) + + if not isinstance(th, dict): + raise ValueError("Thresholds (th) must be provided as a dictionary with target names as keys.") + + #binarize scores according to their thresholds + binary_scores = np.zeros_like(scores, dtype=bool) + for i, target in enumerate(target_names): + if target in th: + binary_scores[:, i] = scores[:, i] >= float(th[target]) + else: + raise KeyError(f"No threshold provided for target '{target}'.") + + output = metrics.classification_report(labels, binary_scores, target_names=target_names, output_dict=True) for i, k in enumerate(target_names): if keeps[i] == 0: continue @@ -99,12 +127,31 @@ def _predict(self, sample): class ClassificationReportCalculator(BaseDriftCalculator): name = "class_report" - def __init__(self, label_col=None, score_col=None, target_names=None, th=0.5): + def __init__(self, label_col=None, score_col=None, target_names=None, th=None): super().__init__() self.label_col = label_col self.score_col = score_col self.target_names = target_names self.th = th + + def prepare(self, ref): + self._ref = self.convert(ref) + if 'activation' in self._ref and 'label' in self._ref: + + activations_df = pd.DataFrame(self._ref['activation'].tolist(), index=self._ref.index) + labels_df = pd.DataFrame(self._ref['label'].tolist(), index=self._ref.index) + combined_df = pd.concat([activations_df, labels_df], axis=1) + + combined_df.columns = ['activation.' + str(k) for i, k in enumerate(self.target_names)] + \ + ['label.' + str(k) for i, k in enumerate(self.target_names)] + + self._ref_df = combined_df + self._is_prepared = True + else: + print("Error: 'ref' does not contain required 'activation' and 'label' columns.") + + self.th = youden_point(self._ref_df, self.target_names) + self._is_prepared = True def convert(self, arg): if not isinstance(arg, pd.DataFrame): diff --git a/src/model_drift/drift/sampler.py b/src/model_drift/drift/sampler.py index e8bf7e6..1682002 100644 --- a/src/model_drift/drift/sampler.py +++ b/src/model_drift/drift/sampler.py @@ -24,3 +24,23 @@ def sample(self, sample, stratify=None): def sample_iterator(self, sample, n_samples=1, stratify=None): for _ in range(n_samples): yield self.sample(sample, stratify=stratify) + +class DummySampler(object): + """ + This is a dummy sampler that mimics the behavior of the Sampler class, but does + not perform any samples and instead returns the indicies unchanged. This is used + for the JackKnife resampling, in which the reference dataframe is resampled in the + metric itself and the sliding window sample is not resampled at all. In this case + we still want to create multiple copies of the sample window. + """ + def __init__(self, sample_size=0, replacement=True, random_state=None): + # These variables are all dummy variables, just kept for future + # compatibility. + self.sample_size = sample_size + self.replacement = replacement + self.random_state = random_state + + def sample_iterator(self, indices, n_samples=1, stratify=None): + + for _ in range(n_samples): + yield indices diff --git a/src/scripts/analysis/analysis_utils.py b/src/scripts/analysis/analysis_utils.py new file mode 100644 index 0000000..b900fbb --- /dev/null +++ b/src/scripts/analysis/analysis_utils.py @@ -0,0 +1,721 @@ +import json +import os +from pathlib import Path +from datetime import datetime +import logging + +import numpy as np +import pandas as pd +import matplotlib.dates as mdates +import matplotlib.pyplot as plt +import matplotlib as mpl +import plotly.graph_objs as go +import seaborn as sns +from sklearn.metrics import roc_auc_score, roc_curve + +from model_drift.data import mgb_data + +date_format = mdates.DateFormatter('%Y-%m') +month_locator = mdates.MonthLocator(interval=3) + +# Plotting parameters +plt.rcParams['svg.fonttype'] = 'none' +mpl.rcParams['svg.fonttype'] = 'none' +mpl.rcParams['font.size'] = 12 +mpl.rcParams['xtick.labelsize'] = 10 +mpl.rcParams['ytick.labelsize'] = 10 +plt.rcParams.update({ + 'svg.fonttype': 'none', + 'font.size': 12, + 'axes.titlesize': 14, + 'axes.labelsize': 12, + 'xtick.labelsize': 10, + 'ytick.labelsize': 10, + 'legend.fontsize': 10, + 'figure.titlesize': 16, + 'axes.grid': True, + 'grid.alpha': 0.3, + 'grid.color': '#cccccc' +}) + +logger = logging.getLogger(__name__) + +def create_performance_plots(df: pd.DataFrame, output_dir: Path, ref_start: str, ref_end: str): + + date_col = tuple(f'Unnamed: 0_level_{i}' for i in range(4)) + target_names = tuple(mgb_data.LABEL_GROUPINGS) + target_names = target_names + ('micro avg', 'macro avg') + num_cols = 2 + num_rows = (len(target_names) + num_cols - 1) // num_cols + + # Create a figure and flatten axs for indexing + fig, axs = plt.subplots(num_rows, num_cols, figsize=(10, num_rows * 5)) + axs = axs.flatten() + + # Define the date range to limit to reference and monitoring period + start_date = pd.to_datetime('2019-11-01') + end_date = pd.to_datetime('2021-07-01') + + # Convert date_col to datetime and filter the DataFrame + df['date'] = pd.to_datetime(df[date_col]) + df_filtered = df[(df['date'] >= start_date) & (df['date'] <= end_date)] + + # Loop over each label name and plot data + for i, name in enumerate(target_names): + # Extract the required series from the DataFrame + if not isinstance(df_filtered.index, pd.DatetimeIndex): + df_filtered.set_index(pd.to_datetime(df_filtered[date_col]), inplace=True) + auroc_series = df_filtered[('performance', name, 'auroc', 'mean')] + f1_score_series = df_filtered[('performance', name, 'f1-score', 'mean')] + date_series = df_filtered[date_col] + + # Calculate the average AUROC during the monitoring period + avg_auroc = df_filtered[('performance', name, 'auroc', 'mean')].mean() + + _normalized_auroc, avg_auroc_ref, std_auroc_ref = normalize_series(auroc_series, ref_start, ref_end, return_mean_std=True) + + std_3_upper = avg_auroc_ref + 3 * std_auroc_ref + std_3_lower = avg_auroc_ref - 3 * std_auroc_ref + + + # Plot each metric on its subplot + mask = (date_series >= ref_start) & (date_series <= ref_end) + axs[i].plot(date_series[~mask], auroc_series[~mask], label='AUROC', color='blue') + axs[i].plot(date_series[mask], auroc_series[mask], label='AUROC during reference period', color='cornflowerblue') + axs[i].axhline(y=avg_auroc_ref, color='gray', linestyle='--', label='Avg AUROC during Reference Window') + axs[i].fill_between(date_series, std_3_upper, std_3_lower, color='gray', alpha=0.35, label='3 Std Range of Reference Window') + #axs[i].plot(date_series, f1_score_series, label='F1-Score') + axs[i].set_title(name) + axs[i].legend() + #axs[i].grid(False) + axs[i].set_xlim(pd.to_datetime('2019-11-01').date(), pd.to_datetime('2021-07-01').date()) + axs[i].tick_params(axis='x', rotation=45) + axs[i].set_ylim(0.7, 1) + axs[i].spines['top'].set_visible(False) + axs[i].spines['right'].set_visible(False) + axs[i].yaxis.set_major_locator(plt.MultipleLocator(0.1)) + + # Save each plot individually + fig2, ax2 = plt.subplots(figsize=(10, 5)) + mask = (date_series >= ref_start) & (date_series <= ref_end) + ax2.plot(date_series[~mask], auroc_series[~mask], label='AUROC', color='blue') + ax2.plot(date_series[mask], auroc_series[mask], label='AUROC during reference period', color='cornflowerblue') + ax2.axhline(y=avg_auroc_ref, color='gray', linestyle='--', label='Avg AUROC during Reference Window') + ax2.fill_between(date_series, std_3_upper, std_3_lower, color='gray', alpha=0.35, label='3 Std Range of Reference Window') + # Add vertical line on Junary 1st, 2020 and March 10th + ax2.axvline(x=pd.to_datetime('2020-01-01'), color='darkblue', linestyle='--', linewidth=1) + ax2.axvline(x=pd.to_datetime('2020-03-10'), color='#5088A1', linestyle='--', linewidth=1) + ax2.set_title(name) + ax2.legend() + #ax2.grid(False) + ax2.set_xlim(pd.to_datetime('2019-11-01').date(), pd.to_datetime('2021-07-01').date()) + ax2.tick_params(axis='x', rotation=45) + ax2.set_ylim(0.7, 1) + ax2.spines['top'].set_visible(False) + ax2.spines['right'].set_visible(False) + ax2.yaxis.set_major_locator(plt.MultipleLocator(0.1)) + plt.tight_layout() + + + fig2.savefig(os.path.join(output_dir, f'performance_{name}.png'), dpi=600) + fig2.savefig(os.path.join(output_dir, f'performance_{name}.svg'), format='svg', bbox_inches='tight') + # Save performance data as CSV + performance_data = pd.DataFrame({ + 'date': date_series, + 'auroc': auroc_series, + 'avg_auroc_ref': avg_auroc_ref, + 'std_3_upper': std_3_upper, + 'std_3_lower': std_3_lower + }) + performance_data.to_csv(os.path.join(output_dir, f'performance_{name}.csv'), index=False) + + plt.close(fig2) + + # Hide unused axes if there are any + for j in range(i + 1, num_cols * num_rows): + axs[j].axis('off') + + plt.tight_layout() + plt.savefig(os.path.join(output_dir, 'performance_combined.png')) + plt.savefig(os.path.join(output_dir, 'performance_combined.svg'), format='svg', bbox_inches='tight') + plt.close() + +def normalize_series(data, ref_start, ref_end, return_mean_std=False): + """ Normalize the series data based on reference period mean and std """ + reference_period = data.loc[ref_start:ref_end] + mean = reference_period.mean() + std = reference_period.std() + norm_result = (data - mean) / std + + if return_mean_std: + return norm_result, mean, std + else: + return norm_result + +def create_normalized_performance_plots( + df: pd.DataFrame, output_dir: Path, ref_start: str, ref_end: str, + plot_start_date=pd.to_datetime('2019-11-01'), plot_end_date=pd.to_datetime('2021-07-01') +): + date_col = tuple(f'Unnamed: 0_level_{i}' for i in range(4)) + target_names = tuple(mgb_data.LABEL_GROUPINGS) + ('micro avg', 'macro avg') + num_cols = 2 + num_rows = (len(target_names) + num_cols - 1) // num_cols + + fig, axs = plt.subplots(num_rows, num_cols, figsize=(10, num_rows * 5)) + axs = axs.flatten() + + for i, name in enumerate(target_names): + # Extract and normalize the series data + if not isinstance(df.index, pd.DatetimeIndex): + df.set_index(pd.to_datetime(df[date_col]), inplace=True) + auroc_series = df[('performance', name, 'auroc', 'mean')] + #f1_score_series = df[('performance', name, 'f1-score', 'mean')] + date_series = df[date_col] + + normalized_auroc = normalize_series(auroc_series, ref_start, ref_end) + #normalized_f1_score = normalize_series(f1_score_series, ref_start, ref_end) + + axs[i].plot(date_series, normalized_auroc, label='Normalized AUROC') + #axs[i].plot(date_series, normalized_f1_score, label='Normalized F1-Score') + axs[i].axhline(y=0, color='grey', linestyle='--', linewidth=1) # + axs[i].set_title(name) + axs[i].legend() + #axs[i].grid(True) + axs[i].tick_params(axis='x', rotation=45) + axs[i].set_xlim(plot_start_date.date(), plot_end_date.date()) + axs[i].set_ylim(-3, 3) # Adjust the y-axis limits for normalized data + axs[i].yaxis.set_major_locator(plt.MultipleLocator(0.5)) + + for j in range(i + 1, num_cols * num_rows): + axs[j].axis('off') + + plt.tight_layout() + plt.savefig(os.path.join(output_dir, 'normalized_performance_combined.png')) + plt.savefig(os.path.join(output_dir, 'normalized_performance_combined.svg'), format='svg', bbox_inches='tight') + plt.close() + +def create_normalized_performance_plots_w_mmc( + df: pd.DataFrame, output_dir: Path, ref_start: str, ref_end: str, + mmc_df: pd.DataFrame, plot_start_date=pd.to_datetime('2019-11-01'), + plot_end_date=pd.to_datetime('2021-07-01') +): + date_col = tuple(f'Unnamed: 0_level_{i}' for i in range(4)) + target_names = tuple(mgb_data.LABEL_GROUPINGS) + ('micro avg', 'macro avg') + num_cols = 2 + num_rows = (len(target_names) + num_cols - 1) // num_cols + + fig, axs = plt.subplots(num_rows, num_cols, figsize=(10, num_rows * 5)) + axs = axs.flatten() + + for i, name in enumerate(target_names): + # Extract and normalize the series data + if not isinstance(df.index, pd.DatetimeIndex): + df.set_index(pd.to_datetime(df[date_col]), inplace=True) + auroc_series = df[('performance', name, 'auroc', 'mean')] + date_series = df.index + + normalized_auroc = normalize_series(auroc_series, ref_start, ref_end) + + ax1 = axs[i] + ax2 = ax1.twinx() # Create a second y-axis + + ax1.plot(date_series, normalized_auroc, label='Normalized AUROC', color='blue') + ax1.axhline(y=0, color='grey', linestyle='--', linewidth=1) + ax1.set_title(name) + #ax1.grid(True) + ax1.tick_params(axis='x', rotation=45) + ax1.set_xlim(plot_start_date.date(), plot_end_date.date()) + ax1.set_ylim(-3, 3) + ax1.yaxis.set_major_locator(plt.MultipleLocator(0.5)) + ax1.set_ylabel('Normalized AUROC') + + ax2.plot(date_series, mmc_df['mmc'], label='MMC', linestyle='--', color='orange') + ax2.set_ylim(min(mmc_df['mmc']), max(mmc_df['mmc'])) + ax2.set_ylabel('MMC') + + # Combine legends + lines1, labels1 = ax1.get_legend_handles_labels() + lines2, labels2 = ax2.get_legend_handles_labels() + ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper left') + + for j in range(i + 1, num_cols * num_rows): + axs[j].axis('off') + + plt.tight_layout() + plt.savefig(os.path.join(output_dir, 'normalized_performance_with_mmc_combined.png')) + plt.savefig(os.path.join(output_dir, 'normalized_performance_with_mmc_combined.svg'), format='svg', bbox_inches='tight') + plt.close() + +def create_joint_scatter_density_plots(df: pd.DataFrame, output_dir: Path, ref_start: str, ref_end: str, mmc_df: pd.DataFrame): + date_col = tuple(f'Unnamed: 0_level_{i}' for i in range(4)) + target_names = tuple(mgb_data.LABEL_GROUPINGS) + ('micro avg', 'macro avg') + + plt.rc('axes', titlesize=16) # Title font size + plt.rc('axes', labelsize=12) # Axis label font size + plt.rc('xtick', labelsize=5) # X-tick label font size + plt.rc('ytick', labelsize=5) # Y-tick label font size + + num_plots = len(target_names) + num_cols = 2 + num_rows = (num_plots + num_cols - 1) // num_cols + + fig, axs = plt.subplots(num_rows, num_cols, figsize=(10, num_rows * 5)) + axs = axs.flatten() + + # Initialize the final combined data with the MMC column + final_combined_data = mmc_df.loc[ref_end:, ['mmc']] + + for i, name in enumerate(target_names): + if not isinstance(df.index, pd.DatetimeIndex): + df.set_index(pd.to_datetime(df[date_col]), inplace=True) + + auroc_series = df[('performance', name, 'auroc', 'mean')] + normalized_auroc = normalize_series(auroc_series, ref_start, ref_end) + + # Exclude all the dates before ref_end + normalized_auroc = normalized_auroc.loc[ref_end:] + + # Add the normalized AUROC for this target to the final combined data + final_combined_data[f'normalized_auroc_{name}'] = normalized_auroc + + # Create combined_data for this iteration (for plotting) + combined_data = pd.concat([final_combined_data['mmc'], normalized_auroc], axis=1) + combined_data.columns = ['mmc', f'normalized_auroc_{name}'] + + combined_data.index = pd.to_datetime(combined_data.index) + combined_data['date_category'] = combined_data.index >= pd.Timestamp('2020-03-10') + + g = sns.jointplot( + data=combined_data, + x='mmc', y=f'normalized_auroc_{name}', hue='date_category', + kind="scatter", + legend=False, + joint_kws={'edgecolor':'w', 'linewidth':0.2}, + ) + #g.ax_joint.set_title(name, loc='left') + g.ax_joint.set_xlabel('MMC') + g.ax_joint.set_ylabel('Normalized AUROC') + + # Add horizontal lines at +3 and -3 + g.ax_joint.axhline(y=3, color='gray', linestyle='--') + g.ax_joint.axhline(y=-3, color='gray', linestyle='--') + + # Add vertical line at x=10 + g.ax_joint.axvline(x=10, color='gray', linestyle='-.') + + if False: #name == 'cardiomegaly': + y_min, y_max = df[auroc_col].min(), df[auroc_col].max() + max_limit = max(abs(y_min), abs(y_max)) + g.ax_joint.set_ylim([-max_limit, max_limit]) + else: + g.ax_joint.set_ylim([-16, 16]) + + handles, labels = g.ax_joint.get_legend_handles_labels() + labels = ['Before March 10, 2020', 'After March 10, 2020'] + + if i == 0: + g.ax_joint.legend(labels, title='Date', loc='lower right') + + # set grid lines false + g.ax_joint.grid(False) + + g.fig.suptitle(name, x=0.5, y=0.95, ha='center', fontsize=16) + g.fig.subplots_adjust(top=0.9) # Adjust to make room for the title + + g.savefig(os.path.join(output_dir, f'{name}_KDE.svg')) + g.savefig(os.path.join(output_dir, f'{name}_KDE.png'), dpi=600) + + + # Save the final combined data to CSV + final_combined_data.to_csv(os.path.join(output_dir, 'weighted_mmc_vs_performance_combined.csv')) + + # Create Table for the out of spec MMC values + proportions = {} + for name in target_names: + df_prop = final_combined_data[[f'normalized_auroc_{name}', 'mmc']] + + # select rows where mmc is smaller than 10 + df_prop_smaller = df_prop[df_prop['mmc'] < 10] + # calculate proportion where normalized_auroc is within [-3, 3] for mmc < 10 + if len(df_prop_smaller) > 0: + proportion_smaller = len(df_prop_smaller[(df_prop_smaller[f'normalized_auroc_{name}'] > -3) & (df_prop_smaller[f'normalized_auroc_{name}'] < 3)]) / len(df_prop_smaller) + else: + proportion_smaller = 0 + + # select rows where mmc is larger than or equal to 10 + df_prop_larger = df_prop[df_prop['mmc'] >= 10] + # calculate proportion where normalized_auroc is within [-3, 3] for mmc >= 10 + if len(df_prop_larger) > 0: + proportion_larger = len(df_prop_larger[(df_prop_larger[f'normalized_auroc_{name}'] > -3) & (df_prop_larger[f'normalized_auroc_{name}'] < 3)]) / len(df_prop_larger) + else: + proportion_larger = 0 + proportions[name] = (proportion_smaller, proportion_larger) + + df_proportions = pd.DataFrame(proportions).T + + # rename columns + df_proportions.columns = ['Proportion MMC < 10', 'Proportion MMC >= 10'] + df_proportions = df_proportions.round(3) + df_proportions.to_csv(os.path.join(output_dir, 'weighted_mmc_vs_performance_combined_proportions.csv')) + + +def create_mmc_performance_roc_plots( + df: pd.DataFrame, + output_dir: Path, + ref_start: str, + ref_end: str, + mmc_df: pd.DataFrame, + emergency_date: str = '2020-03-10', + delays_days=(0, 7, 15, 30), +): + """ROC curves for MMC+ separating windows before vs. after the state of emergency. + + The ground-truth label is purely date-based: windows dated on/after the + cutoff are labelled positive, earlier windows negative. One ROC curve is + evaluated per cutoff, where each cutoff is ``emergency_date`` plus a delay + from ``delays_days`` (in days). The MMC+ score is used as the classifier + score, so higher MMC+ should indicate a post-cutoff window. Reports the AUC + and the Youden-J operating point (the MMC+ threshold maximizing + sensitivity + specificity - 1) for each cutoff. + """ + plot_data = mmc_df[['mmc']].copy() + plot_data.index = pd.to_datetime(plot_data.index) + plot_data = plot_data.dropna().sort_index() + + y_score = plot_data['mmc'].to_numpy() + roc_summaries = [] + + fig, ax = plt.subplots(figsize=(6, 6)) + colors = plt.cm.viridis(np.linspace(0, 0.85, len(delays_days))) + + for delay, color in zip(delays_days, colors): + cutoff = pd.Timestamp(emergency_date) + pd.Timedelta(days=delay) + y_true = (plot_data.index >= cutoff).astype(int) + + if len(np.unique(y_true)) < 2: + logger.warning( + f"Skipping MMC date ROC for delay {delay}d: only one date class present around {cutoff.date()}." + ) + continue + + fpr, tpr, thresholds = roc_curve(y_true, y_score) + auc = roc_auc_score(y_true, y_score) + + # Youden's J: threshold maximizing (sensitivity + specificity - 1) = tpr - fpr + youden_j = tpr - fpr + best_idx = int(np.argmax(youden_j)) + best_threshold = thresholds[best_idx] + best_fpr = fpr[best_idx] + best_tpr = tpr[best_idx] + + roc_summaries.append({ + 'delay_days': delay, + 'cutoff_date': str(cutoff.date()), + 'auc': auc, + 'n_windows': int(len(plot_data)), + 'n_before': int((1 - y_true).sum()), + 'n_after': int(y_true.sum()), + 'youden_j': float(youden_j[best_idx]), + 'youden_mmc_threshold': float(best_threshold), + 'sensitivity_at_youden': float(best_tpr), + 'specificity_at_youden': float(1 - best_fpr), + }) + + ax.plot( + fpr, + tpr, + color=color, + linewidth=2, + label=f'+{delay}d (AUC = {auc:.3f}, Youden MMC+ = {best_threshold:.2f})', + ) + ax.scatter(best_fpr, best_tpr, color=color, s=60, zorder=5) + + plot_data[f'after_cutoff_{delay}d'] = y_true + + if not roc_summaries: + plt.close(fig) + logger.warning("Skipping MMC date ROC plot: no valid cutoff produced two classes.") + return + + ax.plot([0, 1], [0, 1], color='gray', linestyle='--', linewidth=1, label='Chance') + # Small margins so curves lying on the axis edges are not clipped + ax.set_xlim(-0.02, 1.02) + ax.set_ylim(-0.02, 1.02) + ax.set_xlabel('False Positive Rate (1 - Specificity)') + ax.set_ylabel('True Positive Rate (Sensitivity)') + ax.set_title('MMC+ discrimination of pre- versus post-state-of-emergency windows') + ax.legend(loc='lower right', fontsize=8) + ax.grid(True, alpha=0.3) + plt.tight_layout() + + fig.savefig(os.path.join(output_dir, 'mmc_date_roc.svg'), bbox_inches='tight') + fig.savefig(os.path.join(output_dir, 'mmc_date_roc.png'), dpi=600, bbox_inches='tight') + plt.close(fig) + + plot_data.to_csv(os.path.join(output_dir, 'mmc_date_roc_data.csv')) + + pd.DataFrame(roc_summaries).round(3).to_csv( + os.path.join(output_dir, 'mmc_date_roc_summary.csv'), index=False + ) + + +def create_mmc_plot(df, date_col, output_dir, title, col_plot='MMC', mmc_lower=None, + mmc_upper=None, plot_start_date=pd.to_datetime('2019-11-01'), + plot_end_date=pd.to_datetime('2021-07-01')): + + col_plot_display = 'MMC+' if col_plot.lower() == 'mmc' else col_plot + col_name = col_plot.lower() + + # Create complete date range and merge to introduce NaN values for missing dates + date_range = pd.date_range(start=plot_start_date, end=plot_end_date, freq='D') + date_df = pd.DataFrame({'date': date_range}) + + df = pd.DataFrame({ + 'date': df[date_col], + col_name: df[col_name], + }) + df = df.merge(date_df, on='date', how='right') + df.sort_values(by='date', inplace=True) + + # Check if there are NaN values in the 'mmc' columns and count them + nan_count = df[col_name].isna().sum() + if nan_count > 0: + logger.warning(f"Warning: There are {nan_count} NaN values in the '{col_name}' column.") + + # Interpolate NaN values only if there are less than 3 days of gap + df[col_name] = df[col_name].interpolate(method='linear', limit=2, limit_direction='both') + + # Check if there are still NaN values after interpolation + remaining_nan = df[col_name].isna().sum() + if remaining_nan > 0: + logger.warning(f"Warning: There are still {remaining_nan} NaN values in the '{col_name}' column after interpolation.") + logger.warning("These NaN values represent gaps of 3 or more days and were not interpolated.") + + # Create the figure with the original size + fig, ax = plt.subplots(figsize=(4.8, 2.4), facecolor='white') + + ax.plot(df['date'], df[col_name], label=col_plot_display, color='r') + + if mmc_lower is not None and mmc_upper is not None: + + mmc_lower = pd.DataFrame({ + 'date': mmc_lower[date_col], + 'mmc': mmc_lower['mmc'], + }) + mmc_upper = pd.DataFrame({ + 'date': mmc_upper[date_col], + 'mmc': mmc_upper['mmc'], + }) + mmc_lower = mmc_lower.merge(date_df, on='date', how='right') + mmc_upper = mmc_upper.merge(date_df, on='date', how='right') + + mmc_lower.sort_values(by='date', inplace=True) + mmc_upper.sort_values(by='date', inplace=True) + + # Interpolate NaN values for mmc_lower and mmc_upper if provided + nan_count_lower = mmc_lower['mmc'].isna().sum() + if nan_count_lower > 0: + logger.warning(f"Warning: There are {nan_count_lower} NaN values in the 'mmc' column of mmc_lower.") + mmc_lower['mmc'] = mmc_lower['mmc'].interpolate(method='linear', limit=2, limit_direction='both') + remaining_nan_lower = mmc_lower['mmc'].isna().sum() + if remaining_nan_lower > 0: + logger.warning(f"Warning: There are still {remaining_nan_lower} NaN values in the 'mmc' column of mmc_lower after interpolation.") + logger.warning("These NaN values represent gaps of 3 or more days and were not interpolated.") + + nan_count_upper = mmc_upper['mmc'].isna().sum() + if nan_count_upper > 0: + logger.warning(f"Warning: There are {nan_count_upper} NaN values in the 'mmc' column of mmc_upper.") + mmc_upper['mmc'] = mmc_upper['mmc'].interpolate(method='linear', limit=2, limit_direction='both') + remaining_nan_upper = mmc_upper['mmc'].isna().sum() + if remaining_nan_upper > 0: + logger.warning(f"Warning: There are still {remaining_nan_upper} NaN values in the 'mmc' column of mmc_upper after interpolation.") + logger.warning("These NaN values represent gaps of 3 or more days and were not interpolated.") + + + ax.fill_between(df['date'], mmc_lower['mmc'], mmc_upper['mmc'], + alpha=0.5, label='Mean ± 3 Std', color='gray') + + # Add vertical line on Junary 1st, 2020 and March 10th + ax.axvline(x=pd.to_datetime('2020-01-01'), color='darkblue', linestyle='--', linewidth=1) + ax.axvline(x=pd.to_datetime('2020-03-10'), color='#5088A1', linestyle='--', linewidth=1) + + ax.set_title(title, fontsize=8) + ax.set_xlabel('Date', fontsize=8) + ax.set_ylabel(col_plot_display, fontsize=8) + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.legend(fontsize=10) + ax.set_xlim(plot_start_date.date(), plot_end_date.date()) + + # TODO: Only for ER and WAC2 + #ax.set_ylim(-10, 85) + ## Set y-axis ticks + #ax.set_yticks([0, 10, 20, 30, 40, 50, 60, 70, 80]) + #ax.set_yticklabels(['0', '10', '20', '30', '40', '50', '60', '70', '80'], fontsize=6) + + # Standardize tick sizes and rotation + ax.tick_params(axis='both', which='major', labelsize=6) + ax.tick_params(axis='x', rotation=45) + + # Adjust layout to prevent cut-off labels + plt.tight_layout() + + fig = plt.gcf() + fig_width, fig_height = fig.get_size_inches() + + # Save the plot + plt.savefig(output_dir / f'{title.lower().replace(" ", "_")}.png', dpi=600, bbox_inches='tight') + fig.set_size_inches(fig_width, fig_height) + plt.savefig(output_dir / f'{title.lower().replace(" ", "_")}.svg', format='svg', bbox_inches='tight') + + # Resize the figure + fig.set_size_inches(10, 6) # New size + + # Save the plot in the new size + plt.savefig(output_dir / f'{title.lower().replace(" ", "_")}_large.png', dpi=600, bbox_inches='tight') + plt.savefig(output_dir / f'{title.lower().replace(" ", "_")}_large.svg', format='svg', bbox_inches='tight') + + # Save the plot data as a CSV + plot_data = pd.DataFrame({ + 'date': df['date'], + col_name: df[col_name], + }) + + if mmc_lower is not None and mmc_upper is not None: + plot_data['mmc_lower'] = mmc_lower['mmc'] + plot_data['mmc_upper'] = mmc_upper['mmc'] + + plot_data.to_csv(output_dir / f'{title.lower().replace(" ", "_")}.csv', index=False) + plt.close() + +def parse_date(filename): + try: + return datetime.strptime(filename.split('.')[0], '%Y-%m-%d') + except ValueError as e: + pass + +def plot_hist_feature(feature, basepath, output_dir): + dates = os.listdir(basepath) + + #select only the 1st and 15th of each month + dates_parsed = [(date, parse_date(date)) for date in dates if parse_date(date) is not None] + sorted_dates_parsed = sorted(dates_parsed, key=lambda x: x[1]) + sorted_dates_filtered = [date for date, date_obj in sorted_dates_parsed if date_obj.day in {1, 15}] + + all_categories = set() + for date in sorted_dates_filtered: + date_json = os.path.join(basepath, date) + with open(date_json, 'r') as f: + data = json.load(f) + if "histogram" in data["drilldowns"][feature]: + all_categories.update(data["drilldowns"][feature]["histogram"]["x"]) + + else: + break + # Convert set to sorted list to maintain order + all_categories = sorted(all_categories) + + + # Make a figure and add the traces for histogram and kde + fig = go.Figure() + max_y_value = 0 + trace_counter = 0 + for i, date in enumerate(sorted_dates_filtered): + date_json = os.path.join(basepath, date) + + with open(date_json, 'r') as f: + data = json.load(f) + + if "kdehistplot" in data["drilldowns"][feature] and "kde_x" in data["drilldowns"][feature]["kdehistplot"]: + # Extract the necessary data + hist = data["drilldowns"][feature]["kdehistplot"]["hist"] + edges = data["drilldowns"][feature]["kdehistplot"]["plot_edges"] + centers = data["drilldowns"][feature]["kdehistplot"]["plot_centers"] + kde_x = data["drilldowns"][feature]["kdehistplot"]["kde_x"] + kde = data["drilldowns"][feature]["kdehistplot"]["kde"] + + # Add traces for each date + fig.add_trace( + go.Bar(x=centers, y=hist, marker=dict(color='blue'), name=f'Histogram', opacity=0.75, visible=(i == 0)) + ) + fig.add_trace( + go.Scatter(x=kde_x, y=kde, mode='lines', line=dict(color='red'), name=f'KDE', visible=(i == 0)) + ) + trace_counter += 2 + max_hist_value = max(hist) + max_kde_value = max(kde) + max_y_value = max(max_y_value, max_hist_value, max_kde_value) + + elif "kdehistplot" in data["drilldowns"][feature]: + hist = data["drilldowns"][feature]["kdehistplot"]["hist"] + edges = data["drilldowns"][feature]["kdehistplot"]["plot_edges"] + centers = data["drilldowns"][feature]["kdehistplot"]["plot_centers"] + + # Add traces for each date + fig.add_trace( + go.Bar(x=centers, y=hist, marker=dict(color='blue'), name=f'Histogram', opacity=0.75, visible=(i == 0)) + ) + + trace_counter += 1 + max_hist_value = max(hist) + max_y_value = max(max_y_value, max_hist_value) + + else: + # Categorical data processing + probability = data["drilldowns"][feature]["histogram"]["probability"] + category_data = {cat: 0 for cat in all_categories} # Initialize all categories with 0 + for cat, prob in zip(data["drilldowns"][feature]["histogram"]["x"], probability): + category_data[cat] = prob + + fig.add_trace( + go.Bar(x=list(category_data.keys()), y=list(category_data.values()), marker=dict(color='blue'), name=f'Category Probability {date}', visible=(i == 0)) + ) + trace_counter += 1 + max_y_value = max(max_y_value, max(probability)) + + + # Create steps for the slider + steps = [] + visibility_array = [False] * trace_counter + + current_trace_index = 0 + for i, date in enumerate(sorted_dates_filtered): + visible = visibility_array[:] + data = json.load(open(os.path.join(basepath, sorted_dates_filtered[i]), 'r')) + if "kdehistplot" in data["drilldowns"][feature] and "kde_x" in data["drilldowns"][feature]["kdehistplot"]: + visible[current_trace_index] = True + visible[current_trace_index + 1] = True + current_trace_index += 2 + + elif "kdehistplot" in data["drilldowns"][feature]: + visible[current_trace_index] = True + current_trace_index += 1 + else: + visible[current_trace_index] = True + current_trace_index += 1 + + steps.append({ + 'method': 'update', + 'args': [{'visible': visible}, {'title': f"Histogram for {feature} on {date.split('.')[0]}"}], + 'label': date.split('.')[0] + }) + # Create and add slider + sliders = [dict( + active=0, + currentvalue={"prefix": "Date: "}, + pad={"t": 80}, + steps=steps + )] + + fig.update_layout( + sliders=sliders, + title_text=f"Histogram for {feature} on " + sorted_dates_filtered[0], + height=600, + width=1000, + title_x=0.5, + title_y=0.9, + ) + fig.update_yaxes(range=[0, max_y_value]) + + fig.write_html(os.path.join(output_dir,f'{feature}_histogram_interactive.html')) + #fig.show() \ No newline at end of file diff --git a/src/scripts/analysis/basic_performance_plots.py b/src/scripts/analysis/basic_performance_plots.py new file mode 100644 index 0000000..77b9901 --- /dev/null +++ b/src/scripts/analysis/basic_performance_plots.py @@ -0,0 +1,334 @@ +import json +import logging +import os +import sys +from datetime import datetime +from pathlib import Path + +import click +import matplotlib.dates as mdates +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from pycrumbs import tracked +from tqdm import tqdm + +from scripts.analysis import analysis_utils + +date_format = mdates.DateFormatter('%Y-%m') +month_locator = mdates.MonthLocator(interval=3) +plt.rcParams['svg.fonttype'] = 'none' + + +@click.command() +@click.argument('drift-csv-path', type=Path) +@click.argument('output-dir', type=Path) +@click.option('--window-length', type=str, default='30D') +@click.option( + '--equal-weights', + type=bool, + default=True, + help=( + 'If true, the VAE, score, and metadata metrics are weighted equally. ' + 'The correlation weights are only used to weight within the different ' + 'metadata values.' + ) +) +@click.option( + '--skip-histograms', + is_flag=True, + default=False, + help='Skip creating interactive drilldown histogram plots.', +) +@tracked(directory_parameter='output_dir') +def basic_performance_plots( + drift_csv_path: Path, + output_dir: Path, + window_length: str = '30D', + equal_weights: bool = True, + skip_histograms: bool = False, +): + """Makes some basic performance against time plots from a drift CSV.""" + + # Setup output directory + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Setup logging + log_file_path = os.path.join(output_dir, 'plotting_log.log') + logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s', handlers=[logging.FileHandler(log_file_path), logging.StreamHandler(sys.stdout)]) + logging.info(f"Saving plots to {output_dir}") + + + df = pd.read_csv(drift_csv_path, header=[0, 1, 2, 3]) + logging.info(f"Loaded drift CSV from: {drift_csv_path}") + + # Try to load the raw file with all exams in reference window to get the start and end dates + ref_csv_path = str(drift_csv_path).replace('output.csv', 'ref.csv') + if os.path.exists(ref_csv_path): + try: + ref_csv = pd.read_csv(ref_csv_path) + logging.info(f"Loaded reference CSV from: {ref_csv_path}") + ref_window_start_str = ref_csv["StudyDate"].min() + ref_window_end_str = ref_csv["StudyDate"].max() + + except Exception as e: + logging.error(f"Error reading reference CSV: {e}") + raise + else: + logging.warning(f"Reference CSV not found at {ref_csv_path}. Using standard dates for the reference window: 2019-10-01 to 2019-12-31.") + ref_window_start_str = '2019-10-01' + ref_window_end_str = '2019-12-31' + + ref_window_start = datetime.strptime(ref_window_start_str, "%Y-%m-%d") + ref_window_end = datetime.strptime(ref_window_end_str, "%Y-%m-%d") + + # Convert window length string in days to month float + try: + window_length_days = int(window_length.rstrip('D')) + except ValueError as e: + raise ValueError(f"{e} Note: only integer days are supported for window length.") + except Exception as e: + logging.error(f"An unexpected error occurred: {e}") + raise + + # add window_length to the ref_window_start to account for the overlap with the period before (1 window length into + # the reference window is the first day where only days that are actually within the reference window are included) + ref_window_start = ref_window_start + pd.DateOffset(days=window_length_days) + + + # The date column gets read in with a strange name + date_col = tuple(f'Unnamed: 0_level_{i}' for i in range(4)) + performance_col = ('performance', 'micro avg', 'auroc', 'mean') + + + df[date_col] = pd.to_datetime(df[date_col]) + + + analysis_utils.create_performance_plots(df, output_dir, ref_window_start, ref_window_end) + analysis_utils.create_normalized_performance_plots(df, output_dir, ref_window_start, ref_window_end) + + + # Unweighted MMC + mmc_cols = [ + col for col in df.columns + if not col[0].startswith('performance') + and col[2] == 'distance' + and col[3] == 'mean' + ] + mmc_cols_std = [ + col for col in df.columns + if not col[0].startswith('performance') + and col[2] == 'distance' + and col[3] == 'std' + ] + + vae_cols = [ + col for col in df.columns + if col[0].startswith('mu') | col[0].startswith('full_mu') + and col[2] == 'distance' + and col[3] == 'mean' + ] + + score_cols = [ + col for col in df.columns + if col[0].startswith('activation') + and col[2] == 'distance' + and col[3] == 'mean' + ] + + metadata_cols = [ + col for col in df.columns + if not col[0].startswith('performance') + if not col[0].startswith('mu') | col[0].startswith('full_mu') + if not col[0].startswith('activation') + and col[2] == 'distance' + and col[3] == 'mean' + ] + + + mmc_df = df[mmc_cols + [date_col]].copy() + mmc_df_lower = df[mmc_cols_std + [date_col]].copy() + mmc_df_upper = df[mmc_cols_std + [date_col]].copy() + + ref_df = mmc_df[(mmc_df[date_col] >= ref_window_start) & (mmc_df[date_col] <= ref_window_end)].copy() + + mmc_df_weights = df[mmc_cols + [date_col]+ [performance_col]].copy() + ref_df_weights = mmc_df_weights[(mmc_df_weights[date_col] >= ref_window_start) & (mmc_df_weights[date_col] <= ref_window_end)].copy() + + + # Normalize columns by mean and std of reference data + for c in mmc_cols: + mmc_df[c] = (mmc_df[c] - ref_df[c].mean()) / (ref_df[c].std() + 1e-6) + # replace mean word in c with std + c_list = list(c) + c_list[-1] = 'std' + c_std = tuple(c_list) + mmc_df_lower[c_std] = (df[c] - 3 * mmc_df_lower[c_std] - ref_df[c].mean()) / (ref_df[c].std() + 1e-6) + mmc_df_upper[c_std] = (df[c] + 3 * mmc_df_upper[c_std] - ref_df[c].mean()) / (ref_df[c].std() + 1e-6) + + mmc_df['mmc'] = mmc_df.mean(axis=1, numeric_only=True) + mmc_df_lower['mmc'] = mmc_df_lower.mean(axis=1, numeric_only=True) + mmc_df_upper['mmc'] = mmc_df_upper.mean(axis=1, numeric_only=True) + + + analysis_utils.create_mmc_plot(mmc_df, date_col, output_dir, title='Unweighted MMC+ with Range', mmc_lower=mmc_df_lower, mmc_upper=mmc_df_upper) + analysis_utils.create_mmc_plot(mmc_df, date_col, output_dir, title='Unweighted MMC+') + + #TODO: These plots currently use the unweighted MMC, but are not used in the paper + analysis_utils.create_normalized_performance_plots_w_mmc(df, output_dir, ref_window_start, ref_window_end, mmc_df) + + if equal_weights: + # For runs in the paper, we will weight each drift compoment (metadata, vae, activations) as 1/3. Within the metadata we are still weighting according to correlation with performance + correlation_matrix = ref_df_weights[metadata_cols + [performance_col]].corr() + performance_correlation = correlation_matrix[performance_col] + performance_correlation_df = pd.DataFrame(performance_correlation) + plot_df = performance_correlation_df.reset_index() + plot_df.columns = ['_'.join(col).strip() if isinstance(col, tuple) else col for col in plot_df.columns] + + plot_df.drop(columns=["level_1___", "level_2___", "level_3___"], inplace=True) + plot_df.columns = ["Metric", "Avg_AUROC_mean"] + + # drop row where Metric is performance + plot_df = plot_df[plot_df["Metric"] != "performance"] + + weights_raw = pd.Series(plot_df.Avg_AUROC_mean.values, index=plot_df.Metric).to_dict() + weights = {metric: abs(weight) for metric, weight in weights_raw.items()} + + #replace nan values with 0 for next step + weights = {metric: (0 if np.isnan(weight) else weight) for metric, weight in weights.items()} + + # normalize and take negative value. Note: Here the weights should sum to 1/3, as we will be adding the vae and scores each with 1/3 as well + weights = {metric: (1/3) * weight / (sum(weights.values())) for metric, weight in weights.items()} + + # add weight for vae and score + weights['full_mu'] = 1/3 + weights['activation'] = 1/3 + + logging.info(f'Equal weighting was used. Weights: {weights}') + + else: + # Weighted MMC + correlation_matrix = ref_df_weights.corr() + + # To get correlation with the performance column specifically + performance_correlation = correlation_matrix[performance_col] + performance_correlation_df = pd.DataFrame(performance_correlation) + plot_df = performance_correlation_df.reset_index() + plot_df.columns = ['_'.join(col).strip() if isinstance(col, tuple) else col for col in plot_df.columns] + + plot_df.drop(columns=["level_1___", "level_2___", "level_3___"], inplace=True) + plot_df.columns = ["Metric", "Avg_AUROC_mean"] + + # drop row where Metric is performance + plot_df = plot_df[plot_df["Metric"] != "performance"] + + weights_raw = pd.Series(plot_df.Avg_AUROC_mean.values, index=plot_df.Metric).to_dict() + weights = {metric: abs(weight) for metric, weight in weights_raw.items()} + + #replace nan values with 0 for next step + weights = {metric: (0 if np.isnan(weight) else weight) for metric, weight in weights.items()} + + # normalize and take negative value -> should that be applied before? + weights = {metric: (1) * weight / sum(weights.values()) for metric, weight in weights.items()} + + logging.info(f'Correlation weighting was used for all metrics. Weights: {weights}') + + + # save weights for future reference + with open(os.path.join(output_dir, 'correlation_weights.json'), 'w') as f: + json.dump(weights, f) + + # Create weighted MMC dataframes + mmc_df_weighted = mmc_df.copy() + mmc_df_lower_weighted = mmc_df_lower.copy() + mmc_df_upper_weighted = mmc_df_upper.copy() + + mmc_df_weighted.drop(columns=["mmc"], inplace=True) + mmc_df_lower_weighted.drop(columns=["mmc"], inplace=True) + mmc_df_upper_weighted.drop(columns=["mmc"], inplace=True) + + + # Apply weights to the MMC dataframes + for col in mmc_df_weighted.columns: + metric_name = [metric for metric in col if metric in weights] + if metric_name: + mmc_df_weighted[col] = mmc_df_weighted[col] * weights[metric_name[0]] + + else: + logging.warning(f"Column {col} does not match any metric name in the weights dictionary") + mmc_df_weighted['mmc'] = mmc_df_weighted[mmc_cols].sum(axis=1) + + for col in mmc_df_lower_weighted.columns: + metric_name = [metric for metric in col if metric in weights] + if metric_name: + mmc_df_lower_weighted[col] = mmc_df_lower_weighted[col] * weights[metric_name[0]] + + else: + logging.warning(f"Column {col} does not match any metric name in the weights dictionary") + mmc_df_lower_weighted['mmc'] = mmc_df_lower_weighted[mmc_cols_std].sum(axis=1) + + for col in mmc_df_upper_weighted.columns: + metric_name = [metric for metric in col if metric in weights] + if metric_name: + mmc_df_upper_weighted[col] = mmc_df_upper_weighted[col] * weights[metric_name[0]] + + else: + logging.warning(f"Column {col} does not match any metric name in the weights dictionary") + mmc_df_upper_weighted['mmc'] = mmc_df_upper_weighted[mmc_cols_std].sum(axis=1) + + # Create plots for weighted MMC + analysis_utils.create_mmc_plot(mmc_df_weighted, date_col, output_dir, title='Weighted MMC+ with Range', mmc_lower=mmc_df_lower_weighted, mmc_upper=mmc_df_upper_weighted) + analysis_utils.create_mmc_plot(mmc_df_weighted, date_col, output_dir, title='Weighted MMC+') + analysis_utils.create_joint_scatter_density_plots(df, output_dir, ref_window_start, ref_window_end, mmc_df_weighted) + analysis_utils.create_mmc_performance_roc_plots(df, output_dir, ref_window_start, ref_window_end, mmc_df_weighted) + + # Create plots for VAE features alone, using the weighted values + vae_df = mmc_df_weighted[vae_cols + [date_col]].copy() + vae_df['mean_vae_distance'] = vae_df[vae_cols].sum(axis=1) + + analysis_utils.create_mmc_plot(vae_df, date_col, output_dir, title='VAE', col_plot='mean_vae_distance') + + # Create plots for activation features alone, using the weighted values + score_df = mmc_df_weighted[score_cols + [date_col]].copy() + score_df['mean_activation_distance'] = score_df[score_cols].sum(axis=1) + + analysis_utils.create_mmc_plot(score_df, date_col, output_dir, title='Score', col_plot='mean_activation_distance') + + # Create plots for Metadata features alone, using the weighted values + metadata_df = mmc_df_weighted[metadata_cols + [date_col]].copy() + metadata_df['mean_metadata_distance'] = metadata_df[metadata_cols].sum(axis=1) + + analysis_utils.create_mmc_plot(metadata_df, date_col, output_dir, title='Metadata', col_plot='mean_metadata_distance') + + + # Create Histograms for the drilldown features if present + dirname = os.path.dirname(drift_csv_path) + base_path_drilldown = os.path.join(dirname, 'history') + + if skip_histograms: + logging.info("Skipping drilldown histogram plots (--skip-histograms).") + elif not os.path.isdir(base_path_drilldown): + logging.info(f"No drilldown history directory found at {base_path_drilldown}. Skipping histogram plots.") + else: + # Load one example date json to get the keys, needs to be adjusted if using a different dataset + date_json = os.path.join(base_path_drilldown, '2019-10-10.json') + if not os.path.exists(date_json): + logging.info(f"No example drilldown JSON found at {date_json}. Skipping histogram plots.") + else: + with open(date_json, 'r') as f: + data = json.load(f) + + keys = data['drilldowns'].keys() + + if keys: + output_dir_hist = Path(os.path.join(output_dir, 'histograms')) + output_dir_hist.mkdir(parents=True, exist_ok=True) + for feature in tqdm(keys, desc="Creating Histograms"): + analysis_utils.plot_hist_feature(feature, basepath=base_path_drilldown, output_dir=output_dir_hist) + else: + logging.info("There is no drilldown data present so no histograms could be created") + + +if __name__ == "__main__": + basic_performance_plots() \ No newline at end of file