Repository: ME-ICA/tedana

Issue:
New quality control measure: Fraction of accepted/rejected common variance
### Summary

I had a conversation with @afni-rickr about orthogonalizing the rejected time series from the accepted time series (i.e. the --tedort option). Even if this option is not selected, this might be the basis of a good QC measure.

### Additional Detail

Proposed QC measure:
1. Fit the time series of the rejected components to the accepted components
2. For each accepted component, calculate the % variance explained by the rejected components
3. For the accepted components, calculate the weighted sum of the (% variance explained by the rejected components) * (% of total variance explained by the accepted components).

| Accepted Component # | % total variance explained | % variance in component explained by rejected |
| --- | --- | --- |
| 1 | 3% | 5% |
| 2 | 2% | 10% |
| 12 | 0.1% | 100% |

QC measure would be 0.03*0.05 + 0.02+0.1 + 0.001*1 = 12.25%
That is, an accepted component that is 100% explained by rejected components, but is only 0.1% of the total variance is less weighted than a component that is only 5% explained by rejected components, but is 3% of the total variance.

### Next Steps

1. Do others agree this is useful
2. Someone volunteers to add this.

Hints from issue comments:

Solution patch:
diff --git a/docs/api.rst b/docs/api.rst
index 1722002..529e81a 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -197,9 +197,33 @@ API
    tedana.io.writeresults
    tedana.io.writeresults_echoes
 
-
 .. _api_stats_ref:
 
+******************************************
+:mod:`tedana.reporting`: Reporting functions
+******************************************
+
+.. automodule:: tedana.reporting
+   :no-members:
+   :no-inherited-members:
+
+.. currentmodule:: tedana.reporting
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   tedana.reporting.html_report.generate_report
+   tedana.reporting.html_report.get_elbow_val
+   tedana.reporting.quality_metrics.calculate_rejected_components_impact
+   tedana.reporting.static_figures.comp_figures
+   tedana.reporting.static_figures.pca_results
+   tedana.reporting.static_figures.plot_t2star_and_s0
+   tedana.reporting.static_figures.plot_rmse
+   tedana.reporting.static_figures.plot_adaptive_mask
+
+.. _api_reporting_ref:
+
 ******************************************
 :mod:`tedana.stats`: Statistical functions
 ******************************************
diff --git a/tedana/metrics/collect.py b/tedana/metrics/collect.py
index 0365752..02b4a7d 100644
--- a/tedana/metrics/collect.py
+++ b/tedana/metrics/collect.py
@@ -672,6 +672,19 @@ def get_metadata(component_table: pd.DataFrame) -> Dict:
                 -1: "Component is flipped prior to metric calculation.",
             },
         }
+    if "Var Exp of rejected to accepted" in component_table:
+        metric_metadata["Var Exp of rejected to accepted"] = {
+            "LongName": "100*R2 of fit of rejected to accepted mixing matrix time series",
+            "Description": (
+                "The time series for each component are not independent. "
+                "Time series from rejected components are regressed from the data. "
+                "For each accepted component, this is a calculation of how much "
+                "the variance of the accept component times series is explained "
+                "by the rejected component time series. "
+                "This is not used in decision trees but might be a useful quality metric."
+            ),
+            "Units": "percent",
+        }
 
     # There are always components in the component_table, definitionally
     metric_metadata["Component"] = {
diff --git a/tedana/reporting/__init__.py b/tedana/reporting/__init__.py
index bb52bca..f6ca226 100644
--- a/tedana/reporting/__init__.py
+++ b/tedana/reporting/__init__.py
@@ -1,6 +1,12 @@
 """Reporting code for tedana."""
 
 from tedana.reporting.html_report import generate_report
+from tedana.reporting.quality_metrics import calculate_rejected_components_impact
 from tedana.reporting.static_figures import comp_figures, pca_results
 
-__all__ = ["generate_report", "comp_figures", "pca_results"]
+__all__ = [
+    "calculate_rejected_components_impact",
+    "comp_figures",
+    "generate_report",
+    "pca_results",
+]
diff --git a/tedana/reporting/dynamic_figures.py b/tedana/reporting/dynamic_figures.py
index 99c1a6f..d95bdd2 100644
--- a/tedana/reporting/dynamic_figures.py
+++ b/tedana/reporting/dynamic_figures.py
@@ -87,6 +87,7 @@ def _create_data_struct(comptable_path, color_mapping=color_mapping):
 
     # remove space from column name
     df.rename(columns={"variance explained": "var_exp"}, inplace=True)
+    df.rename(columns={"Var Exp of rejected to accepted": "var_exp_rej"}, inplace=True)
 
     # For providing sizes based on Var Explained that are visible
     mm_scaler = MinMaxScaler(feature_range=(4, 20))
@@ -119,6 +120,7 @@ def _create_data_struct(comptable_path, color_mapping=color_mapping):
             kappa=df["kappa"],
             rho=df["rho"],
             varexp=df["var_exp"],
+            varexprej=df["var_exp_rej"],
             kappa_rank=df["kappa_rank"],
             rho_rank=df["rho_rank"],
             varexp_rank=df["var_exp_rank"],
@@ -159,6 +161,7 @@ def _create_kr_plt(comptable_cds, kappa_elbow=None, rho_elbow=None):
             ("Kappa", "@kappa{0.00}"),
             ("Rho", "@rho{0.00}"),
             ("Var. Expl.", "@varexp{0.00}%"),
+            ("Var. Expl. by Rej.", "@varexprej{0.00}%"),
             ("Tags", "@classtag"),
         ]
     )
@@ -278,6 +281,7 @@ def _create_sorted_plt(
             ("Kappa", "@kappa{0.00}"),
             ("Rho", "@rho{0.00}"),
             ("Var. Expl.", "@varexp{0.00}%"),
+            ("Var. Expl. by Rej.", "@varexprej{0.00}%"),
             ("Tags", "@classtag"),
         ]
     )
@@ -331,6 +335,7 @@ def _create_varexp_pie_plt(comptable_cds):
             ("Kappa", "@kappa{0.00}"),
             ("Rho", "@rho{0.00}"),
             ("Var. Expl.", "@varexp{0.00}%"),
+            ("Var. Expl. by Rej.", "@varexprej{0.00}%"),
             ("Tags", "@classtag"),
         ]
     )
diff --git a/tedana/reporting/quality_metrics.py b/tedana/reporting/quality_metrics.py
new file mode 100644
index 0000000..941a85a
--- /dev/null
+++ b/tedana/reporting/quality_metrics.py
@@ -0,0 +1,74 @@
+"""File to assess quality metrics for components."""
+
+import numpy as np
+
+from tedana.stats import fit_model
+
+
+def calculate_rejected_components_impact(selector, mixing):
+    """Calculate the % variance explained by the rejected components for accepted components.
+
+    The final metric is the weighted sum across accepted components of
+    the total variance explained by each accepted component times
+    the varianced explained (100*R^2) of each accepted component by the rejected comp time series.
+    This quantifies the impact of rejected components on the overall variance explained
+    by accepted components.
+
+    Parameters
+    ----------
+    selector : :obj: tedana.selection.component_selector.ComponentSelector
+        Uses `variance explained` and `classification` columns in ``selector.component_table_``.
+    mixing : (T [x C]) array_like
+        Mixing matrix for converting input data to component space, where `C`
+        is components and `T` is time
+
+    Returns
+    -------
+    None
+        Updates elements the selector object:
+        ``selector.component_table_`` has an added column for "Var Exp of rejected to accepted"
+        ``selector.cross_component_metrics_`` has an added value
+        `total_var_exp_rejected_components_on_accepted`
+        `variance explained` is a percentage ranging from 0-100% and
+        `total_var_exp_rejected_components_on_accepted` is also a percent with the same range.
+    """
+    component_table = selector.component_table_
+
+    rej = component_table.loc[component_table["classification"] == "rejected"].index
+    acc = component_table.loc[component_table["classification"] == "accepted"].index
+
+    if len(rej) == 0 or len(acc) == 0:
+        component_table["Var Exp of rejected to accepted"] = np.nan
+        selector.cross_component_metrics_["total_var_exp_rejected_components_on_accepted"] = np.nan
+        return
+
+    rej_arrs, acc_arrs = mixing[:, rej], mixing[:, acc]
+
+    # Calculate the R-squared values for each accepted component
+    _, sse, _ = fit_model(rej_arrs, acc_arrs)
+    # Note: Since the mixing matrix time series are typically mean centered with stdev 1,
+    #    ss_total will be the number of time points (+ rounding errors)
+    ss_total = np.sum(
+        (acc_arrs - np.tile(np.mean(acc_arrs, axis=0), [acc_arrs.shape[0], 1])) ** 2, axis=0
+    )
+    r2 = 1 - (sse / ss_total)
+
+    if "Var Exp of rejected to accepted" in component_table.columns:
+        # Reset to NaN before calculating for accepted components
+        component_table["Var Exp of rejected to accepted"] = np.nan
+    else:
+        # initialize new column in component table that's left of the classification & tag columns
+        num_columns = len(component_table.columns)
+        component_table.insert(num_columns - 2, "Var Exp of rejected to accepted", np.nan)
+
+    # Calculated r2 is assigned to accepted components and the rest remain NaN
+    component_table.loc[acc, "Var Exp of rejected to accepted"] = 100 * r2
+
+    # Update selector with the overall impact measure
+    selector.cross_component_metrics_["total_var_exp_rejected_components_on_accepted"] = (
+        np.sum(
+            component_table["Var Exp of rejected to accepted"][acc]
+            * component_table["variance explained"][acc]
+        )
+        / 100
+    )
diff --git a/tedana/selection/tedica.py b/tedana/selection/tedica.py
index 5ea6b79..ed88324 100644
--- a/tedana/selection/tedica.py
+++ b/tedana/selection/tedica.py
@@ -2,8 +2,6 @@
 
 import logging
 
-from tedana.metrics import collect
-
 LGR = logging.getLogger("GENERAL")
 RepLGR = logging.getLogger("REPORT")
 
@@ -63,6 +61,5 @@ def automatic_selection(component_table, selector, **kwargs):
     # TODO external_regressor_config=external_regressor_config, was here
     component_table["classification_tags"] = ""
     selector.select(component_table, cross_component_metrics=kwargs)
-    selector.metadata_ = collect.get_metadata(selector.component_table_)
 
     return selector
diff --git a/tedana/workflows/ica_reclassify.py b/tedana/workflows/ica_reclassify.py
index 736df94..6025983 100644
--- a/tedana/workflows/ica_reclassify.py
+++ b/tedana/workflows/ica_reclassify.py
@@ -436,6 +436,10 @@ def ica_reclassify_workflow(
         print(tc)
         io_generator.save_file(ioh.get_file_contents(tc), tc)
 
+    # calculate the fit of rejected to accepted components to use as a quality measure
+    # Note: This adds a column to component_table & needs to run before the table is saved
+    reporting.quality_metrics.calculate_rejected_components_impact(selector, mixing)
+
     # Save component selector and tree
     selector.to_files(io_generator)
 
diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py
index 367931c..a45ad44 100644
--- a/tedana/workflows/tedana.py
+++ b/tedana/workflows/tedana.py
@@ -936,6 +936,10 @@ def tedana_workflow(
     betas_oc = utils.unmask(computefeats2(data_optcom, mixing, mask_denoise), mask_denoise)
     io_generator.save_file(betas_oc, "z-scored ICA components img")
 
+    # calculate the fit of rejected to accepted components to use as a quality measure
+    # Note: This adds a column to component_table & needs to run before the table is saved
+    reporting.quality_metrics.calculate_rejected_components_impact(selector, mixing)
+
     # Save component selector and tree
     selector.to_files(io_generator)
     # Save metrics and metadata


Test patch:
diff --git a/tedana/tests/test_reporting.py b/tedana/tests/test_reporting.py
index 6d05325..53b4af1 100644
--- a/tedana/tests/test_reporting.py
+++ b/tedana/tests/test_reporting.py
@@ -3,6 +3,8 @@
 import numpy as np
 
 from tedana import reporting
+from tedana.tests.test_external_metrics import sample_mixing_matrix
+from tedana.tests.test_selection_utils import sample_selector
 
 
 def test_smoke_trim_edge_zeros():
@@ -11,4 +13,60 @@ def test_smoke_trim_edge_zeros():
     assert reporting.static_figures._trim_edge_zeros(arr) is not None
 
 
-# TODO: Test other functions in reporting?
+def test_calculate_rejected_components_impact():
+    selector = sample_selector()
+    mixing = sample_mixing_matrix()
+
+    component_table = selector.component_table_
+
+    reporting.calculate_rejected_components_impact(selector, mixing)
+
+    # Total matches expecting value in testing data
+    assert np.round(
+        selector.cross_component_metrics_["total_var_exp_rejected_components_on_accepted"], 4
+    ) == np.round(4.0927834, 4)
+
+    assert "Var Exp of rejected to accepted" in selector.component_table_.columns
+    rej = component_table[component_table["classification"] == "rejected"]
+    acc = component_table[component_table["classification"] == "accepted"]
+
+    assert rej["Var Exp of rejected to accepted"].isna().all()
+    assert not acc["Var Exp of rejected to accepted"].isna().all()
+
+
+def test_calculate_rejected_components_impact_no_rej():
+    selector = sample_selector()
+    mixing = sample_mixing_matrix()
+
+    component_table = selector.component_table_
+
+    component_table.drop(
+        component_table[component_table["classification"] == "rejected"].index, inplace=True
+    )
+
+    reporting.calculate_rejected_components_impact(selector, mixing)
+
+    assert "Var Exp of rejected to accepted" in selector.component_table_.columns
+    assert selector.component_table_["Var Exp of rejected to accepted"].isna().all()
+    assert np.isnan(
+        selector.cross_component_metrics_["total_var_exp_rejected_components_on_accepted"]
+    )
+
+
+def test_calculate_rejected_components_impact_no_acc():
+    selector = sample_selector()
+    mixing = sample_mixing_matrix()
+
+    component_table = selector.component_table_
+
+    component_table.drop(
+        component_table[component_table["classification"] == "accepted"].index, inplace=True
+    )
+
+    reporting.calculate_rejected_components_impact(selector, mixing)
+
+    assert "Var Exp of rejected to accepted" in selector.component_table_.columns
+    assert selector.component_table_["Var Exp of rejected to accepted"].isna().all()
+    assert np.isnan(
+        selector.cross_component_metrics_["total_var_exp_rejected_components_on_accepted"]
+    )