Merge pull request #97 from liuzhenqi77/dominance-stats

[ENH] Adds dominance stats function

Author: rmarkello

docs/api.rst

@@ -111,6 +111,7 @@ Python Reference API
    permtest_1samp
    permtest_rel
    permtest_pearsonr
+   get_dominance_stats
 
 .. _ref_metrics:
 

netneurotools/stats.py

@@ -6,9 +6,12 @@
 import warnings
 
 import numpy as np
+from tqdm import tqdm
+from itertools import combinations
 from scipy import optimize, spatial, special, stats as sstats
 from scipy.stats.stats import _chk2_asarray
 from sklearn.utils.validation import check_random_state
+from sklearn.linear_model import LinearRegression
 
 from . import utils
 
@@ -784,3 +787,141 @@ def gen_spinsamples(coords, hemiid, n_rotate=1000, check_duplicates=True,
         return spinsamples, cost
 
     return spinsamples
+
+
+def get_dominance_stats(X, y, use_adjusted_r_sq=True, verbose=False):
+    """
+    Returns the dominance analysis statistics for multilinear regression.
+
+    This is a rewritten & simplified version of [DA1]_. It is briefly
+    tested against the original package, but still in early stages.
+    Please feel free to report any bugs.
+
+    Warning: Still work-in-progress. Parameters might change!
+
+    Parameters
+    ----------
+    X : (N, M) array_like
+        Input data
+    y : (N,) array_like
+        Target values
+    use_adjusted_r_sq : bool, optional
+        Whether to use adjusted r squares. Default: True
+    verbose : bool, optional
+        Whether to print debug messages. Default: False
+
+    Returns
+    -------
+    model_metrics : dict
+        The dominance metrics, currently containing `individual_dominance`,
+        `partial_dominance`, `total_dominance`, and `full_r_sq`.
+    model_r_sq : dict
+        Contains all model r squares
+
+    Notes
+    -----
+    Example usage
+
+    .. code:: python
+
+        from netneurotools.stats import get_dominance_stats
+        from sklearn.datasets import load_boston
+        X, y = load_boston(return_X_y=True)
+        model_metrics, model_r_sq = get_dominance_stats(X, y)
+
+    To compare with [DA1]_, use `use_adjusted_r_sq=False`
+
+    .. code:: python
+
+        from dominance_analysis import Dominance_Datasets
+        from dominance_analysis import Dominance
+        boston_dataset=Dominance_Datasets.get_boston()
+        dominance_regression=Dominance(data=boston_dataset,
+                                       target='House_Price',objective=1)
+        incr_variable_rsquare=dominance_regression.incremental_rsquare()
+        dominance_regression.dominance_stats()
+
+    References
+    ----------
+    .. [DA1] https://github.com/dominance-analysis/dominance-analysis
+
+    """
+
+    # this helps to remove one element from a tuple
+    def remove_ret(tpl, elem):
+        lst = list(tpl)
+        lst.remove(elem)
+        return tuple(lst)
+
+    # sklearn linear regression wrapper
+    def get_reg_r_sq(X, y):
+        lin_reg = LinearRegression()
+        lin_reg.fit(X, y)
+        yhat = lin_reg.predict(X)
+        SS_Residual = sum((y - yhat) ** 2)
+        SS_Total = sum((y - np.mean(y)) ** 2)
+        r_squared = 1 - (float(SS_Residual)) / SS_Total
+        adjusted_r_squared = 1 - (1 - r_squared) * \
+            (len(y) - 1) / (len(y) - X.shape[1] - 1)
+        if use_adjusted_r_sq:
+            return adjusted_r_squared
+        else:
+            return r_squared
+
+    # generate all predictor combinations in list (num of predictors) of lists
+    n_predictor = X.shape[-1]
+    # n_comb_len_group = n_predictor - 1
+    predictor_combs = [list(combinations(range(n_predictor), i))
+                       for i in range(1, n_predictor + 1)]
+    if verbose:
+        print(f"[Dominance analysis] Generated \
+              {len([v for i in predictor_combs for v in i])} combinations")
+
+    # get all r_sq's
+    model_r_sq = dict()
+    for len_group in tqdm(predictor_combs, desc='num-of-predictor loop',
+                          disable=not verbose):
+        for idx_tuple in tqdm(len_group, desc='insider loop',
+                              disable=not verbose):
+            r_sq = get_reg_r_sq(X[:, idx_tuple], y)
+            model_r_sq[idx_tuple] = r_sq
+    if verbose:
+        print(f"[Dominance analysis] Acquired {len(model_r_sq)} r^2's")
+
+    # getting all model metrics
+    model_metrics = dict([])
+
+    # individual dominance
+    individual_dominance = []
+    for i_pred in range(n_predictor):
+        individual_dominance.append(model_r_sq[(i_pred,)])
+    individual_dominance = np.array(individual_dominance).reshape(1, -1)
+    model_metrics["individual_dominance"] = individual_dominance
+
+    # partial dominance
+    partial_dominance = [[]] * (n_predictor - 1)
+    for i_len in range(n_predictor - 1):
+        i_len_combs = list(combinations(range(n_predictor), i_len + 2))
+        for j_node in range(n_predictor):
+            j_node_sel = [v for v in i_len_combs if j_node in v]
+            reduced_list = [remove_ret(comb, j_node) for comb in j_node_sel]
+            diff_values = [
+                model_r_sq[j_node_sel[i]] - model_r_sq[reduced_list[i]]
+                for i in range(len(reduced_list))]
+            partial_dominance[i_len].append(np.mean(diff_values))
+
+    # save partial dominance
+    partial_dominance = np.array(partial_dominance)
+    model_metrics["partial_dominance"] = partial_dominance
+    # get total dominance
+    total_dominance = np.mean(
+        np.r_[individual_dominance, partial_dominance], axis=0)
+    # test and save total dominance
+    assert np.allclose(total_dominance.sum(),
+                       model_r_sq[tuple(range(n_predictor))]), \
+           "Sum of total dominance is not equal to full r square!"
+    model_metrics["total_dominance"] = total_dominance
+    # save full r^2
+    model_metrics["full_r_sq"] = model_r_sq[tuple(range(n_predictor))]
+
+    return model_metrics, model_r_sq

requirements.txt

@@ -5,3 +5,4 @@ nilearn
 numpy>=1.16
 scikit-learn
 scipy>=1.4.0
+tqdm