Repository: rishi-kulkarni/hierarch (language: python)

Issue:
Allow to only receive resampling indices (make hierarch agnostic to specific statistical tests).
I am currently working with data from the Human Connectome Project. The dataset contains a lot of siblings, sometimes even dizygotic or even monozygotic twins. Basically for all my statistical tests, I would have to consider this nested structure when computing p-values (which I would like to derive from permutation test) or confidence intervals (which I would like to obtain from bootstrapping tests). This is how the demographic data would look like:

![image](https://user-images.githubusercontent.com/48254290/217231196-869663cd-5ab9-4ed6-a1a9-0480f3398b20.png)


The problem is, that I would like to conduct multivariate tests (i.e. CCA to be more specific). Which means that besides from the demographic information that could be used to conduct a hierarchical resampling procedure, I also have matrices  `X`  and `y` that contain the "actual data" and that I would use to perform my statistical test. When looking at the [Hypothesis Testing section](https://hierarch.readthedocs.io/en/stable/user/hypothesis.html#), I can see that hierarch currently seems to be restricted to simple univariate tests like t-tests and ANOVA? Do you think it would be feasible to allow users to only obtain resampled indices so they could apply their own tests?

This is my current approach which is quite tedious:
1.) Use [FSL-Palm](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM) to derive [exchangability blocks](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM/ExchangeabilityBlocks) using the [hcp2blocks.m](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM/ExchangeabilityBlocks#EBs_for_data_of_the_Human_Connectome_Project) function 
2.) Use a [python implementation](https://github.com/jameschapman19/scikit-palm) of FSL-Palm's quicker function (`skpalm.permutations.quickperms`) to derive a matrix of reampled indices.

Allow to only receive resampling indices (make hierarch agnostic to specific statistical tests).
I am currently working with data from the Human Connectome Project. The dataset contains a lot of siblings, sometimes even dizygotic or even monozygotic twins. Basically for all my statistical tests, I would have to consider this nested structure when computing p-values (which I would like to derive from permutation test) or confidence intervals (which I would like to obtain from bootstrapping tests). This is how the demographic data would look like:

![image](https://user-images.githubusercontent.com/48254290/217231196-869663cd-5ab9-4ed6-a1a9-0480f3398b20.png)


The problem is, that I would like to conduct multivariate tests (i.e. CCA to be more specific). Which means that besides from the demographic information that could be used to conduct a hierarchical resampling procedure, I also have matrices  `X`  and `y` that contain the "actual data" and that I would use to perform my statistical test. When looking at the [Hypothesis Testing section](https://hierarch.readthedocs.io/en/stable/user/hypothesis.html#), I can see that hierarch currently seems to be restricted to simple univariate tests like t-tests and ANOVA? Do you think it would be feasible to allow users to only obtain resampled indices so they could apply their own tests?

This is my current approach which is quite tedious:
1.) Use [FSL-Palm](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM) to derive [exchangability blocks](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM/ExchangeabilityBlocks) using the [hcp2blocks.m](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/PALM/ExchangeabilityBlocks#EBs_for_data_of_the_Human_Connectome_Project) function 
2.) Use a [python implementation](https://github.com/jameschapman19/scikit-palm) of FSL-Palm's quicker function (`skpalm.permutations.quickperms`) to derive a matrix of reampled indices.



Discussion / hints:
Many thanks for your note @JohannesWiesner!

You raise an interesting and quite useful question, which would ideally be very straightforward to implement - I could expose the helper function the `hypothesis_test` function uses to generate resampled datasets. I assume with that, you should be able to perform your desired analysis, right? 
Yup, that should do the trick! :) 
I'll take a crack at this over the next few weekends - hopefully it is going to be as simple as I hope. 
Any update on this? Still could use the bootstrapping implementation from hierarch :) AFAIK FSL Palm does not offer drawing with replacement.
@rishi-kulkarni thought I'd just comment to say I work on similar problems to @JohannesWiesner and agree there would be a lot of value to exposing the function to generate the indices to the extent it adds helpful motivation :)

Nice work on this package!
Many thanks for your note @JohannesWiesner!

You raise an interesting and quite useful question, which would ideally be very straightforward to implement - I could expose the helper function the `hypothesis_test` function uses to generate resampled datasets. I assume with that, you should be able to perform your desired analysis, right? 
Yup, that should do the trick! :) 
I'll take a crack at this over the next few weekends - hopefully it is going to be as simple as I hope. 
Any update on this? Still could use the bootstrapping implementation from hierarch :) AFAIK FSL Palm does not offer drawing with replacement.
@rishi-kulkarni thought I'd just comment to say I work on similar problems to @JohannesWiesner and agree there would be a lot of value to exposing the function to generate the indices to the extent it adds helpful motivation :)

Nice work on this package!

Fix patch:
diff --git a/README.md b/README.md
index 6b57514..0669f0b 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 ## A Hierarchical Resampling Package for Python
 
-Version 1.1.3
+Version 1.1.6
 
 hierarch is a package for hierarchical resampling (bootstrapping, permutation) of datasets in Python. Because for loops are ultimately intrinsic to cluster-aware resampling, hierarch uses Numba to accelerate many of its key functions.
 
@@ -15,12 +15,11 @@ hierarch has several functions to assist in performing resampling-based (and the
 3. [Documentation](#documentation)
 4. [Citation](#citation)
 
-
 <a name="introduction"></a>
-## Introduction 
 
-Design-based randomization tests represents the platinum standard for significance analyses [[1, 2, 3]](#1) - that is, they produce probability statements that depend only on the experimental design, not at all on less-than-verifiable assumptions about the probability distributions of the data-generating process. Researchers can use hierarch to quickly perform automated design-based randomization tests for experiments with arbitrary levels of hierarchy.
+## Introduction
 
+Design-based randomization tests represents the platinum standard for significance analyses [[1, 2, 3]](#1) - that is, they produce probability statements that depend only on the experimental design, not at all on less-than-verifiable assumptions about the probability distributions of the data-generating process. Researchers can use hierarch to quickly perform automated design-based randomization tests for experiments with arbitrary levels of hierarchy.
 
 <a id="1">[1]</a> Tukey, J.W. (1993). Tightening the Clinical Trial. Controlled Clinical Trials, 14(4), 266-285.
 
@@ -28,35 +27,39 @@ Design-based randomization tests represents the platinum standard for significan
 
 <a id="1">[3]</a> Berger, V.W. (2000). Pros and cons of permutation tests in clinical trials. Statistics in Medicine, 19(10), 1319-1328.
 
-
 <a name="setup"></a>
-## Setup 
+
+## Setup
 
 ### Dependencies
-* numpy
-* pandas (for importing data)
-* numba
-* scipy (for power analysis)
+
+- numpy
+- pandas (for importing data)
+- numba
+- scipy (for power analysis)
 
 ### Installation
 
-The easiest way to install hierarch is via PyPi. 
+The easiest way to install hierarch is via PyPi.
 
-```pip install hierarch```
+`pip install hierarch`
 
 Alternatively, you can install from Anaconda.
 
-```conda install -c rkulk111 hierarch```
-
+`conda install -c rkulk111 hierarch`
 
 <a name="documentation"></a>
+
 ## Documentation
+
 Check out our user guide at [readthedocs](https://hierarch.readthedocs.io/).
 
 <a name="citation"></a>
+
 ## Citation
+
 If hierarch helps you analyze your data, please consider citing it. The manuscript also
 contains a set of simulations validating hierarchical randomization tests in a variety of
 conditions.
 
-Kulkarni RU, Wang CL, Bertozzi CR (2022) Analyzing nested experimental designs—A user-friendly resampling method to determine experimental significance. PLoS Comput Biol 18(5): e1010061. https://doi.org/10.1371/journal.pcbi.1010061
\ No newline at end of file
+Kulkarni RU, Wang CL, Bertozzi CR (2022) Analyzing nested experimental designs—A user-friendly resampling method to determine experimental significance. PLoS Comput Biol 18(5): e1010061. https://doi.org/10.1371/journal.pcbi.1010061
diff --git a/docs/reference/_autosummary/hierarch.stats.hierarchical_randomization.rst b/docs/reference/_autosummary/hierarch.stats.hierarchical_randomization.rst
new file mode 100644
index 0000000..9459331
--- /dev/null
+++ b/docs/reference/_autosummary/hierarch.stats.hierarchical_randomization.rst
@@ -0,0 +1,6 @@
+hierarch.stats.hierarchical\_randomization
+==========================================
+
+.. currentmodule:: hierarch.stats
+
+.. autofunction:: hierarchical_randomization
\ No newline at end of file
diff --git a/docs/reference/_autosummary/hierarch.stats.rst b/docs/reference/_autosummary/hierarch.stats.rst
index c9c2917..2ee1908 100644
--- a/docs/reference/_autosummary/hierarch.stats.rst
+++ b/docs/reference/_autosummary/hierarch.stats.rst
@@ -16,6 +16,7 @@
       :nosignatures:
    
       confidence_interval
+      hierarchical_randomization
       hypothesis_test
       multi_sample_test
       studentized_covariance
diff --git a/hierarch/stats.py b/hierarch/stats.py
index 37dbde2..332c498 100644
--- a/hierarch/stats.py
+++ b/hierarch/stats.py
@@ -1,15 +1,18 @@
-import numpy as np
-from numba import jit
 import math
+from functools import lru_cache
 from itertools import combinations
+from typing import Collection, Generator, Optional, Union
+from warnings import simplefilter, warn
+
+import numpy as np
 import pandas as pd
+from numba import jit
+
 from hierarch.internal_functions import (
     GroupbyMean,
     bivar_central_moment,
 )
 from hierarch.resampling import Bootstrapper, Permuter
-from warnings import warn, simplefilter
-from functools import lru_cache
 
 
 def _preprocess_data(data):
@@ -1218,3 +1221,102 @@ def _cov_std_error(x, y):
     # correction of n - root(3) is applied
     denom_3 = ((n - 2) * (bivar_central_moment(x, y, pow=1, ddof=1.75) ** 2)) / (n - 1)
     return ((1 / (n - 1.5)) * (denom_1 + denom_2 - denom_3)) ** 0.5
+
+
+def hierarchical_randomization(
+    data_array: Union[np.ndarray, pd.DataFrame],
+    treatment_col: Union[int, str],
+    skip: Optional[Collection[int]] = None,
+    bootstraps: int = 100,
+    permutations: int = 1000,
+    random_state: Optional[Union[int, np.random.Generator]] = None,
+) -> Generator[np.ndarray, None, None]:
+    """Yields permuted datasets for a hierarchical randomization test.
+
+    Parameters
+    ----------
+    data_array : 2D numpy array or pandas DataFrame
+        Array-like containing both the independent and dependent variables to
+        be analyzed. It's assumed that the final (rightmost) column
+        contains the dependent variable values.
+    treatment_col : int or str
+        The index number of the column containing "N samples" to be compared.
+        Indexing starts at 0. If input data is a pandas DataFrame, this can be
+        the column name.
+    skip : list of ints, optional
+        Columns to skip in the bootstrap. Skip columns that were sampled
+        without replacement from the prior column, by default None
+    bootstraps : int, optional
+        Number of bootstraps to perform, by default 100. Can be set to 1 for a
+        permutation test without any bootstrapping.
+    permutations : int, optional
+        Number of permutations to perform PER bootstrap sample. "all"
+        for exact test (only works if there are only two treatments), by default 1000
+    random_state : int or numpy random Generator, optional
+        Seedable for reproducibility., by default None
+
+    Yields
+    ------
+    Generator[np.ndarray, None, None]
+        Permuted data for a hierarchical randomization test.
+
+    """
+    # turns the input array or dataframe into a float64 array
+    if isinstance(data_array, (np.ndarray, pd.DataFrame)):
+        if isinstance(data_array, pd.DataFrame) and isinstance(treatment_col, str):
+            treatment_col = int(data_array.columns.get_loc(treatment_col))
+        data = _preprocess_data(data_array)
+    else:
+        raise TypeError("Input data must be ndarray or DataFrame.")
+
+    assert isinstance(treatment_col, int)  # mypy
+
+    # set random state
+    rng = np.random.default_rng(random_state)
+
+    # enforce lower bound on skip
+    if skip is not None:
+        skip = list(skip)
+        for v in reversed(skip):
+            if v <= treatment_col + 1:
+                warn("No need to include columns before treated columns in skip.")
+                skip.remove(v)
+    else:
+        skip = []
+
+    # enforce bounds on bootstraps and permutations
+    if not isinstance(bootstraps, int) or bootstraps < 1:
+        raise TypeError("bootstraps must be an integer greater than 0")
+    if isinstance(permutations, str):
+        if permutations != "all":
+            raise TypeError("permutations must be 'all' or an integer greater than 0")
+    elif not isinstance(permutations, int) or permutations < 1:
+        raise TypeError("permutations must be 'all' or an integer greater than 0")
+
+    # initialize and fit the bootstrapper to the data
+    bootstrapper = Bootstrapper(random_state=rng, kind="indexes")
+    bootstrapper.fit(data, skip=skip)
+
+    for i in range(bootstraps):
+        # get a bootstrap sample
+        bootstrapped_sample = bootstrapper.transform(data, start=treatment_col + 2)
+
+        # initialize and fit the permuter to the aggregated data
+        # don't need to seed this, as numba's PRNG state is shared
+        permuter = Permuter()
+
+        if permutations == "all":
+            permuter.fit(bootstrapped_sample, treatment_col, exact=True)
+
+            # in the exact case, determine and set the total number of
+            # possible permutations
+            counts = np.unique(bootstrapped_sample[:, 0], return_counts=True)[1]
+            permutations = _binomial(counts.sum(), counts[0])
+
+        else:
+            # just fit the permuter if this is a randomized test
+            permuter.fit(bootstrapped_sample, treatment_col)
+
+        for j in range(permutations):
+            # yield a permuted sample
+            yield permuter.transform(bootstrapped_sample)
diff --git a/pyproject.toml b/pyproject.toml
index 9b97d23..704102e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "hierarch"
-version = "1.1.5"
+version = "1.1.6"
 description = "Hierarchical hypothesis testing library"
 authors = ["Rishi Kulkarni <rishi@kulkarni.science>"]
 homepage = "https://github.com/rishi-kulkarni/hierarch"


Test patch:
diff --git a/tests/test_stats.py b/tests/test_stats.py
index d9a3e7b..12ef80f 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -205,5 +205,17 @@ def test_conf(self):
         self.assertGreater(interval_95[1], interval_68[1])
 
 
+class TestHierarchicalRandomization(unittest.TestCase):
+    paramlist = [[0, 2], [stats.norm], [stats.norm]]
+    hierarchy = [2, 4, 3]
+    datagen = DataSimulator(paramlist, random_state=2)
+    datagen.fit(hierarchy)
+    data = datagen.generate()
+
+    def test_randomize(self):
+        for indices in hierarch.stats.hierarchical_randomization(self.data, 0):
+            self.assertEqual(len(indices), len(self.data))
+
+
 if __name__ == "__main__":
     unittest.main()