From 8f8d233af7e78a613e2d3443934dab4603a7dfd6 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 11 May 2023 16:09:38 -0700
Subject: [PATCH 01/50] intial peds-bootstrap

---
 q2_fmt/_peds.py        | 42 ++++++++++++++++++++++++++++++++++++++++++
 q2_fmt/plugin_setup.py | 41 ++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 82 insertions(+), 1 deletion(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index a89475e..d2b6cac 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -10,6 +10,7 @@
 import pandas as pd
 import numpy as np
 import warnings
+from scipy.stats import mannwhitneyu
 
 
 def sample_peds(table: pd.DataFrame, metadata: qiime2.Metadata,
@@ -329,3 +330,44 @@ def _create_masking(time_metadata, donor_df, recip_df, reference_column):
 def _mask_recipient(donor_mask, recip_df):
     maskedrecip = donor_mask & recip_df
     return maskedrecip
+
+
+def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
+                   time_column: str, reference_column: str,
+                   subject_column: str,
+                   filter_missing_references: bool = False,
+                   drop_incomplete_subjects: bool = False,
+                   bootstrap_replicates: int = 999):
+    fake_donor = []
+    for i in range(0, bootstrap_replicates):
+        if i == 0:
+            peds, = \
+             sample_peds(
+              table=table,
+              metadata=metadata,
+              time_column=time_column,
+              subject_column=subject_column,
+              reference_column=reference_column,
+              drop_incomplete_subjects=drop_incomplete_subjects,
+              filter_missing_references=filter_missing_references)
+            peds = peds.view(pd.DataFrame).set_index("id")
+            real_donor = peds["measure"].to_list()
+        else:
+            recipient = metadata[reference_column].dropna()
+            donor = metadata[metadata[reference_column].isna()]
+            shifted_list = recipient[reference_column].sample(frac=1).to_list()
+            recipient[reference_column] = shifted_list
+            metadata_df = pd.concat([recipient, donor])
+            metadata = qiime2.Metadata(metadata_df)
+            peds, = \
+                sample_peds(
+                 table=table,
+                 metadata=metadata,
+                 time_column=time_column,
+                 subject_column=subject_column,
+                 reference_column=reference_column,
+                 drop_incomplete_subjects=drop_incomplete_subjects,
+                 filter_missing_references=filter_missing_references)
+            peds = peds.view(pd.DataFrame).set_index("id")
+            fake_donor = fake_donor + peds["measure"].to_list()
+    s, p = mannwhitneyu(real_donor, fake_donor, alternative='greater')
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 7abf839..7cd64e9 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -9,7 +9,8 @@
 import importlib
 
 from qiime2.plugin import (Str, Plugin, Metadata, TypeMap,
-                           Bool, Choices, Visualization, Properties, Citations)
+                           Bool, Choices, Visualization, Properties, Citations,
+                           Int, Range)
 from q2_types.sample_data import SampleData, AlphaDiversity
 from q2_types.distance_matrix import DistanceMatrix
 
@@ -264,6 +265,44 @@
     }
 )
 
+plugin.methods.register_function(
+    function=q2_fmt.peds_bootstrap,
+    inputs={'table': FeatureTable[Frequency | RelativeFrequency |
+                                  PresenceAbsence]},
+    parameters={'metadata': Metadata, 'time_column': Str,
+                'reference_column': Str, 'subject_column': T_subject,
+                'filter_missing_references': Bool,
+                'drop_incomplete_subjects': Bool,
+                'bootstrap_replicates': Int % Range(999, None)},
+    outputs=[('peds_dists', GroupDist[Ordered, Matched] % Properties("peds"))],
+    parameter_descriptions={
+        'metadata': 'The sample metadata.',
+        'time_column': 'The column within the `metadata` that the'
+                       ' `table` should be grouped by. This column'
+                       ' should contain simple integer values.',
+        'reference_column': 'The column within the `metadata` that contains'
+                            ' the sample to use as a reference'
+                            ' for a given `table`.'
+                            ' For example, this may be the relevant donor'
+                            ' sample to compare against.',
+        'subject_column': 'The column within the `metadata` that contains the'
+                          ' subject ID to be tracked against timepoints.',
+        'filter_missing_references': 'Filter out references contained within'
+                                     ' the metadata that are not present'
+                                     ' in the table.'
+                                     ' Default behavior is to raise an error.',
+        'drop_incomplete_subjects': 'Filter out subjects that do not have'
+                                    ' a sample at every timepoint.'
+                                    ' Default behavior is to raise an error.',
+    },
+    output_descriptions={
+        # TODO: I dont know what it is going to output yet?
+    },
+    name='',
+    description='',
+    citations=[citations['aggarwala_precise_2021']]
+)
+
 plugin.visualizers.register_function(
     function=q2_fmt.plot_heatmap,
     inputs={'data': GroupDist[Ordered, Matched] % Properties('peds')},

From 2a3aa9dadc39e00dd06a6e67fec08622d1c60e40 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 8 May 2024 10:00:21 -0700
Subject: [PATCH 02/50] updates from EMC paper to peds-bootstrap

---
 q2_fmt/_peds.py                  | 42 +++++++++------------
 q2_fmt/tests/test_engraftment.py | 64 +++++++++++++++++++++++++++++++-
 2 files changed, 80 insertions(+), 26 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index d2b6cac..6ba0843 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -338,36 +338,28 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
                    filter_missing_references: bool = False,
                    drop_incomplete_subjects: bool = False,
                    bootstrap_replicates: int = 999):
+    metadata_df = metadata.to_dataframe
+    ## TODO: Grab Donor in a more logic way
+    donor = metadata_df.loc[metadata_df['Location'] == body_site]
+    recipient = metadata_df.loc[metadata_df['Location'] != body_site]
     fake_donor = []
-    for i in range(0, bootstrap_replicates):
+    for i in range(0, bootstrap_replicates+1):
         if i == 0:
-            peds, = \
-             sample_peds(
-              table=table,
-              metadata=metadata,
-              time_column=time_column,
-              subject_column=subject_column,
-              reference_column=reference_column,
-              drop_incomplete_subjects=drop_incomplete_subjects,
-              filter_missing_references=filter_missing_references)
+            peds, = sample_peds(table=table, metadata=metadata, 
+                                time_column=time_column,
+                                reference_column=reference_column,
+                                subject_column=subject_column)
             peds = peds.view(pd.DataFrame).set_index("id")
-            real_donor = peds["measure"].to_list()
+            real_temp = peds["measure"].to_list()
         else:
-            recipient = metadata[reference_column].dropna()
-            donor = metadata[metadata[reference_column].isna()]
             shifted_list = recipient[reference_column].sample(frac=1).to_list()
-            recipient[reference_column] = shifted_list
-            metadata_df = pd.concat([recipient, donor])
+            recipient.loc[:, reference_column] = shifted_list
+            metadata_df = pd.concat([donor, recipient])
             metadata = qiime2.Metadata(metadata_df)
-            peds, = \
-                sample_peds(
-                 table=table,
-                 metadata=metadata,
-                 time_column=time_column,
-                 subject_column=subject_column,
-                 reference_column=reference_column,
-                 drop_incomplete_subjects=drop_incomplete_subjects,
-                 filter_missing_references=filter_missing_references)
+            peds, = sample_peds(table=table, metadata=metadata,
+                                time_column=time_column,
+                                reference_column=reference_column,
+                                subject_column=subject_column)
             peds = peds.view(pd.DataFrame).set_index("id")
             fake_donor = fake_donor + peds["measure"].to_list()
-    s, p = mannwhitneyu(real_donor, fake_donor, alternative='greater')
+    s, p = mannwhitneyu(real_temp, fake_donor, alternative='greater')
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index b63a5dd..d29df21 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -18,7 +18,7 @@
                           _filter_associated_reference,
                           _check_reference_column, _check_for_time_column,
                           _check_subject_column, _check_column_type,
-                          feature_peds)
+                          feature_peds, peds_bootstrap)
 
 
 class TestBase(TestPluginBase):
@@ -1022,3 +1022,65 @@ def test_reference_column_name_is_ID(self):
                                     " not be the same as the index of"
                                     " metadata: `id`"):
             _check_reference_column(metadata_df, "id")
+
+
+class TestBoot(TestBase):
+    def test_high_donor_overlap(self):
+        metadata_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Ref': ['donor1', 'donor2', 'donor3', float("Nan"), float("Nan"),
+                    float("Nan")],
+            'subject': ['sub1', 'sub2', 'sub3', float("Nan"), float("Nan"),
+                        float("Nan")],
+            'group': [1, 1, 1, float("Nan"), float("Nan"),
+                      float("Nan")],
+            "Location": [float("Nan"), float("Nan"),
+                         float("Nan"), 'test', 'test',
+                         'test']}).set_index('id')
+
+        table_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Feature1': [1, 0, 0, 1, 0, 0],
+            'Feature2': [0, 1, 0, 0, 1, 0],
+            'Feature3': [0, 0, 1, 0, 0, 1]}).set_index('id')
+
+        p, real_mean, fake_mean = peds_bootstrap(metadata_df=metadata_df,
+                                                 table=table_df,
+                                                 time_column="group",
+                                                 reference_column="Ref",
+                                                 subject_column="subject",
+                                                 iters=999)
+        self.assertGreater(real_mean, fake_mean)
+
+    def test_low_donor_overlap(self):
+        metadata_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Ref': ['donor1', 'donor2', 'donor3', float("Nan"), float("Nan"),
+                    float("Nan")],
+            'subject': ['sub1', 'sub2', 'sub3', float("Nan"), float("Nan"),
+                        float("Nan")],
+            'group': [1, 1, 1, float("Nan"), float("Nan"),
+                      float("Nan")],
+            "Location": [float("Nan"), float("Nan"),
+                         float("Nan"), 'test', 'test',
+                         'test']}).set_index('id')
+
+        table_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Feature1': [1, 0, 0, 0, 0, 1],
+            'Feature2': [0, 1, 0, 1, 0, 0],
+            'Feature3': [0, 0, 1, 0, 1, 0]}).set_index('id')
+
+        p, real_median, fake_median = peds_bootstrap(metadata_df=metadata_df,
+                                                     table=table_df,
+                                                     time_column="group",
+                                                     reference_column="Ref",
+                                                     subject_column="subject",
+                                                     iters=999)
+        self.assertGreater(fake_median, real_median)
+
+

From 82a2e53154730f272c9d2a47466c289717241731 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 8 May 2024 10:47:20 -0700
Subject: [PATCH 03/50] lint from merging

---
 q2_fmt/_peds.py                  | 8 ++++----
 q2_fmt/tests/test_engraftment.py | 5 ++---
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 70af2f9..641e6c0 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -478,13 +478,13 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
                    drop_incomplete_subjects: bool = False,
                    bootstrap_replicates: int = 999):
     metadata_df = metadata.to_dataframe
-    ## TODO: Grab Donor in a more logic way
-    donor = metadata_df.loc[metadata_df['Location'] == body_site]
-    recipient = metadata_df.loc[metadata_df['Location'] != body_site]
+    # TODO: Grab Donor in a more logic way
+    # donor = metadata_df.loc[metadata_df['Location'] == body_site]
+    # recipient = metadata_df.loc[metadata_df['Location'] != body_site]
     fake_donor = []
     for i in range(0, bootstrap_replicates+1):
         if i == 0:
-            peds, = sample_peds(table=table, metadata=metadata, 
+            peds, = sample_peds(table=table, metadata=metadata,
                                 time_column=time_column,
                                 reference_column=reference_column,
                                 subject_column=subject_column)
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 98725a1..dac2324 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -19,8 +19,8 @@
                           _check_reference_column, _check_for_time_column,
                           _check_subject_column, _check_column_type,
                           _drop_incomplete_timepoints, feature_peds,
-                          _check_column_missing, _rename_features, peds_bootstrap)
-
+                          _check_column_missing, _rename_features,
+                          peds_bootstrap)
 
 
 class TestBase(TestPluginBase):
@@ -1189,4 +1189,3 @@ def test_low_donor_overlap(self):
                                                      subject_column="subject",
                                                      iters=999)
         self.assertGreater(fake_median, real_median)
-

From 8370ad89fc5e239fb008eceec6f62064117b4f5c Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 8 May 2024 10:51:04 -0700
Subject: [PATCH 04/50] add bootstrap to init

---
 q2_fmt/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/q2_fmt/__init__.py b/q2_fmt/__init__.py
index 87dc12f..4b8ce33 100644
--- a/q2_fmt/__init__.py
+++ b/q2_fmt/__init__.py
@@ -14,4 +14,4 @@
 del get_versions
 
 __all__ = ['engraftment', 'sample_peds', 'feature_peds',
-           'peds', 'peds_heatmap', 'group_timepoints']
+           'peds', 'peds_heatmap', 'group_timepoints', 'peds_bootstrap']

From 0be6145bb484079ade78218080049c506efec6f4 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Fri, 10 May 2024 08:53:09 -0700
Subject: [PATCH 05/50] updates

---
 q2_fmt/_peds.py        | 6 ++++--
 q2_fmt/plugin_setup.py | 5 +----
 2 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 641e6c0..633f420 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -479,8 +479,10 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
                    bootstrap_replicates: int = 999):
     metadata_df = metadata.to_dataframe
     # TODO: Grab Donor in a more logic way
-    # donor = metadata_df.loc[metadata_df['Location'] == body_site]
-    # recipient = metadata_df.loc[metadata_df['Location'] != body_site]
+    donor = metadata_df[metadata.index.isin(
+        _check_reference_column(metadata=metadata_df,
+                                reference_column=reference_column))]
+    recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
     fake_donor = []
     for i in range(0, bootstrap_replicates+1):
         if i == 0:
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 1a52bbb..063f838 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -299,7 +299,7 @@
     }
 )
 
-plugin.methods.register_function(
+plugin.visualizers.register_function(
     function=q2_fmt.peds_bootstrap,
     inputs={'table': FeatureTable[Frequency | RelativeFrequency |
                                   PresenceAbsence]},
@@ -329,9 +329,6 @@
                                     ' a sample at every timepoint.'
                                     ' Default behavior is to raise an error.',
     },
-    output_descriptions={
-        # TODO: I dont know what it is going to output yet?
-    },
     name='',
     description='',
     citations=[citations['aggarwala_precise_2021']]

From 89306772380f762f8e21da608646671bd4468c1c Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Mon, 13 May 2024 17:07:18 -0700
Subject: [PATCH 06/50] code dump

---
 q2_fmt/__init__.py       |   3 +-
 q2_fmt/_peds.py          | 112 +++++++++++++++++++++++++++++++--------
 q2_fmt/assets/index.html |  61 +++++++++++++++++----
 q2_fmt/plugin_setup.py   |  23 ++++----
 4 files changed, 151 insertions(+), 48 deletions(-)

diff --git a/q2_fmt/__init__.py b/q2_fmt/__init__.py
index 4b8ce33..7e18587 100644
--- a/q2_fmt/__init__.py
+++ b/q2_fmt/__init__.py
@@ -8,7 +8,8 @@
 
 from ._version import get_versions
 from ._engraftment import engraftment, group_timepoints
-from ._peds import sample_peds, feature_peds, peds, peds_heatmap
+from ._peds import (sample_peds, feature_peds, peds, peds_heatmap,
+                    peds_bootstrap)
 
 __version__ = get_versions()['version']
 del get_versions
diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 633f420..11c3f80 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -18,6 +18,7 @@
 import json
 
 from q2_fmt._util import json_replace
+from q2_stats._visualizer import _make_stats
 
 
 def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
@@ -31,6 +32,7 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
 
     if peds_metric == 'sample':
         sample_peds = ctx.get_action('fmt', 'sample_peds')
+        peds_bootstrap = ctx.get_action('fmt', 'peds_bootstrap')
         peds_dist = sample_peds(
             table=table, metadata=metadata, time_column=time_column,
             subject_column=subject_column, reference_column=reference_column,
@@ -38,6 +40,15 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
             drop_incomplete_timepoint=drop_incomplete_timepoint,
             filter_missing_references=filter_missing_references)
 
+        peds_stats = peds_bootstrap(
+            table=table, metadata=metadata, time_column=time_column,
+            subject_column=subject_column,
+            reference_column=reference_column,
+            drop_incomplete_subjects=drop_incomplete_subjects,
+            drop_incomplete_timepoint=drop_incomplete_timepoint,
+            filter_missing_references=filter_missing_references)
+        stat = peds_stats[0]
+
     else:
         if drop_incomplete_subjects or drop_incomplete_timepoint:
             warnings.warn('Feature PEDS was selected as the PEDS metric, which'
@@ -51,15 +62,18 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
             table=table, metadata=metadata, time_column=time_column,
             subject_column=subject_column, reference_column=reference_column,
             filter_missing_references=filter_missing_references)
-    results += peds_heatmap(data=peds_dist[0], level_delimiter=level_delimiter)
+        stat = None
+    results += peds_heatmap(data=peds_dist[0], stats=stat,
+                            level_delimiter=level_delimiter)
 
     return tuple(results)
 
 
 def peds_heatmap(output_dir: str, data: pd.DataFrame,
-                 level_delimiter: str = None):
+                 level_delimiter: str = None, stats: pd.DataFrame = None):
     _rename_features(data=data, level_delimiter=level_delimiter)
-
+    if stats is not None:
+        table1, stats = _make_stats(stats)
     J_ENV = jinja2.Environment(
         loader=jinja2.PackageLoader('q2_fmt', 'assets')
     )
@@ -89,9 +103,10 @@ def peds_heatmap(output_dir: str, data: pd.DataFrame,
                              title=title, measure=gradient,
                              measure_name=measure_name, order=order)
 
-    with open(os.path.join(output_dir, "index.html"), "w") as fh:
+    with open(os.path.join(output_dir, 'index.html'), 'w') as fh:
         spec_string = json.dumps(full_spec)
-        fh.write(index.render(spec=spec_string))
+        fh.write(index.render(spec=spec_string, stats=stats,
+                              table1=table1))
 
 
 def sample_peds(table: pd.DataFrame, metadata: qiime2.Metadata,
@@ -476,31 +491,82 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
                    subject_column: str,
                    filter_missing_references: bool = False,
                    drop_incomplete_subjects: bool = False,
-                   bootstrap_replicates: int = 999):
-    metadata_df = metadata.to_dataframe
-    # TODO: Grab Donor in a more logic way
-    donor = metadata_df[metadata.index.isin(
+                   drop_incomplete_timepoint: list = None,
+                   bootstrap_replicates: int = 999) -> (pd.DataFrame):
+
+    metadata_df = metadata.to_dataframe()
+    donor = metadata_df[metadata_df.index.isin(
         _check_reference_column(metadata=metadata_df,
                                 reference_column=reference_column))]
+
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
-    fake_donor = []
+    fake_donor = pd.DataFrame([])
     for i in range(0, bootstrap_replicates+1):
         if i == 0:
-            peds, = sample_peds(table=table, metadata=metadata,
-                                time_column=time_column,
-                                reference_column=reference_column,
-                                subject_column=subject_column)
-            peds = peds.view(pd.DataFrame).set_index("id")
-            real_temp = peds["measure"].to_list()
+            peds =\
+             sample_peds(table=table, metadata=metadata,
+                         time_column=time_column,
+                         reference_column=reference_column,
+                         subject_column=subject_column,
+                         filter_missing_references=filter_missing_references,
+                         drop_incomplete_subjects=drop_incomplete_subjects,
+                         drop_incomplete_timepoint=drop_incomplete_timepoint)
+            real_temp = peds["measure"]
         else:
             shifted_list = recipient[reference_column].sample(frac=1).to_list()
             recipient.loc[:, reference_column] = shifted_list
             metadata_df = pd.concat([donor, recipient])
             metadata = qiime2.Metadata(metadata_df)
-            peds, = sample_peds(table=table, metadata=metadata,
-                                time_column=time_column,
-                                reference_column=reference_column,
-                                subject_column=subject_column)
-            peds = peds.view(pd.DataFrame).set_index("id")
-            fake_donor = fake_donor + peds["measure"].to_list()
-    s, p = mannwhitneyu(real_temp, fake_donor, alternative='greater')
+            peds =\
+             sample_peds(table=table, metadata=metadata,
+                         time_column=time_column,
+                         reference_column=reference_column,
+                         subject_column=subject_column,
+                         filter_missing_references=filter_missing_references,
+                         drop_incomplete_subjects=drop_incomplete_subjects,
+                         drop_incomplete_timepoint=drop_incomplete_timepoint)
+            fake_donor[i] = peds["measure"]
+    fake_donor_series = fake_donor.median(axis=1)
+    # Common Langauge effect size calcs in prep for stats refactor.
+    # These are not being used currently but will be soon
+    fake_donors_CLES = fake_donor.median(axis=0)
+    agree = ((fake_donors_CLES <
+              real_temp.median()).sum())/bootstrap_replicates
+    disagree = 1 - agree
+    cles = agree - disagree
+    cles
+
+    s, p = mannwhitneyu(real_temp, fake_donor_series.to_list(),
+                        alternative='greater')
+
+    stats_df = pd.DataFrame([["real_values", real_temp.size,
+                              real_temp.median(), "fake_values",
+                              fake_donor_series.size,
+                              fake_donor_series.median(),
+                              fake_donor_series.size,
+                              s, p, np.nan]])
+    stats_df.columns = ['A:group', 'A:n', 'A:measure', 'B:group', 'B:n',
+                        'B:measure', 'n', 'test-statistic', 'p-value',
+                        'q-value']
+
+    stats_df['A:group'].attrs.update({'title': 'real_values',
+                                     'description': '...'})
+    stats_df['B:group'].attrs.update({'title': 'randomized_values',
+                                     'description': '...'})
+    n = {'title': 'count', 'description': '...'}
+    stats_df['A:n'].attrs.update(n)
+    stats_df['B:n'].attrs.update(n)
+    measure = {
+        'title': 'Median PEDS Value',
+        'description': '...'
+    }
+    stats_df['A:measure'].attrs.update(measure)
+    stats_df['B:measure'].attrs.update(measure)
+    stats_df['n'].attrs.update(dict(title='count', description='...'))
+    stats_df['test-statistic'].attrs.update(dict(title='Mann-Whitney U',
+                                                 description='...'))
+    stats_df['p-value'].attrs.update(dict(title='one-tail', description='...'))
+    stats_df['q-value'].attrs.update(
+        dict(title='Benjamini–Hochberg', description='...'))
+
+    return stats_df
diff --git a/q2_fmt/assets/index.html b/q2_fmt/assets/index.html
index 209dce4..8e2b81a 100644
--- a/q2_fmt/assets/index.html
+++ b/q2_fmt/assets/index.html
@@ -13,36 +13,77 @@
       width: 100%;
       box-sizing: border-box;
       margin: 0px;
-      padding: 10px;
-      /* padding: 50px 212px 0px 50px; */
+      padding: 50px 212px 0px 50px;
     }
     #viz {
-      margin: 100px auto;
+      margin: 0 auto;
+      display: block;
+      width: 100%;
       position: relative;
     }
     .vega-bindings {
-      top: 25px;
-      left: 25px;
-      display: flex;
+      top: 90px;
+      right: 50px;
+      width: 200px;
       position: fixed;
       font-family: monospace;
     }
     .vega-bind {
-    margin-right: 50px;
-    vertical-align: middle  
+      margin-right: 15px;
+    }
+    .vega-bind-name {
+      margin-right: 5px;
+    }
+    .vega-embed summary {
+      position: fixed !important;
+      top: 50px !important;
+      right: 222px !important;
+    }
+    .vega-embed .vega-actions {
+      position: fixed !important;
+      right: 212px !important;
+      top: 85px !important;
+    }
+    #stats table {
+      margin-top: 50px;
+      width: 100%;
+      text-align: right;
+      border-collapse: collapse;
+      border: 0px;
+      font-family: monospace;
     }
 
-    .vega-bind span {
-    display: flex;
+    #stats thead tr {
+      border: none;
     }
 
+    #stats tr th {
+      border: none;
+      padding-right: 10px;
+    }
+
+    #stats td {
+      border: none;
+      padding-right: 10px;
+    }
 
+    #stats tr {
+      border-top: 1px solid;
+    }
   </style>
 </head>
 <body>
 
 <div id="viz"></div>
 
+
+{% if stats %}
+<div id='stats'>
+{{ stats }}
+</div>
+<p><strong>Table 1.</strong> {{ table1 }}</p>
+{% endif %}
+
 <script type="text/javascript">
   var spec = JSON.parse(document.getElementById('spec').innerHTML);
   vegaEmbed('#viz', spec, {'renderer': 'svg'}).then(function(result) {
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 063f838..988e6a6 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -204,7 +204,8 @@
 
 plugin.visualizers.register_function(
     function=q2_fmt.peds_heatmap,
-    inputs={'data': GroupDist[Ordered, Matched] % Properties("peds")},
+    inputs={'data': GroupDist[Ordered, Matched] % Properties("peds"),
+            'stats': StatsTable[Pairwise]},
     parameters={'level_delimiter': Str},
     parameter_descriptions={},
     name='PEDS Heatmap',
@@ -299,7 +300,7 @@
     }
 )
 
-plugin.visualizers.register_function(
+plugin.methods.register_function(
     function=q2_fmt.peds_bootstrap,
     inputs={'table': FeatureTable[Frequency | RelativeFrequency |
                                   PresenceAbsence]},
@@ -307,8 +308,9 @@
                 'reference_column': Str, 'subject_column': T_subject,
                 'filter_missing_references': Bool,
                 'drop_incomplete_subjects': Bool,
+                'drop_incomplete_timepoint': List[Str],
                 'bootstrap_replicates': Int % Range(999, None)},
-    outputs=[('peds_dists', GroupDist[Ordered, Matched] % Properties("peds"))],
+    outputs=[('stats', StatsTable[Pairwise])],
     parameter_descriptions={
         'metadata': 'The sample metadata.',
         'time_column': 'The column within the `metadata` that the'
@@ -329,21 +331,14 @@
                                     ' a sample at every timepoint.'
                                     ' Default behavior is to raise an error.',
     },
+    output_descriptions={
+        'stats': 'stats table for comparing real PEDS values'
+        ' from fake Bootstrapped values'
+    },
     name='',
     description='',
     citations=[citations['aggarwala_precise_2021']]
 )
 
-plugin.visualizers.register_function(
-    function=q2_fmt.plot_heatmap,
-    inputs={'data': GroupDist[Ordered, Matched] % Properties('peds')},
-    parameters={},
-    name='Plot Heatmap',
-    description='',
-    examples={
-        'peds_heatmap': ex.peds_heatmap
-     }
-)
-
 
 importlib.import_module('q2_fmt._transformer')

From 12f6e02b66861135536e41571041a76fadcb2d70 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Mon, 13 May 2024 17:12:45 -0700
Subject: [PATCH 07/50] lint

---
 q2_fmt/_peds.py | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 11c3f80..c4a6280 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -503,28 +503,28 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
     fake_donor = pd.DataFrame([])
     for i in range(0, bootstrap_replicates+1):
         if i == 0:
-            peds =\
-             sample_peds(table=table, metadata=metadata,
-                         time_column=time_column,
-                         reference_column=reference_column,
-                         subject_column=subject_column,
-                         filter_missing_references=filter_missing_references,
-                         drop_incomplete_subjects=drop_incomplete_subjects,
-                         drop_incomplete_timepoint=drop_incomplete_timepoint)
+            peds = sample_peds(
+                table=table, metadata=metadata,
+                time_column=time_column,
+                reference_column=reference_column,
+                subject_column=subject_column,
+                filter_missing_references=filter_missing_references,
+                drop_incomplete_subjects=drop_incomplete_subjects,
+                drop_incomplete_timepoint=drop_incomplete_timepoint)
             real_temp = peds["measure"]
         else:
             shifted_list = recipient[reference_column].sample(frac=1).to_list()
             recipient.loc[:, reference_column] = shifted_list
             metadata_df = pd.concat([donor, recipient])
             metadata = qiime2.Metadata(metadata_df)
-            peds =\
-             sample_peds(table=table, metadata=metadata,
-                         time_column=time_column,
-                         reference_column=reference_column,
-                         subject_column=subject_column,
-                         filter_missing_references=filter_missing_references,
-                         drop_incomplete_subjects=drop_incomplete_subjects,
-                         drop_incomplete_timepoint=drop_incomplete_timepoint)
+            peds = sample_peds(
+                table=table, metadata=metadata,
+                time_column=time_column,
+                reference_column=reference_column,
+                subject_column=subject_column,
+                filter_missing_references=filter_missing_references,
+                drop_incomplete_subjects=drop_incomplete_subjects,
+                drop_incomplete_timepoint=drop_incomplete_timepoint)
             fake_donor[i] = peds["measure"]
     fake_donor_series = fake_donor.median(axis=1)
     # Common Langauge effect size calcs in prep for stats refactor.

From 30c0ae7e2ac7099dbe5e1c9caf671bc23ab8d465 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 14 May 2024 11:06:19 -0700
Subject: [PATCH 08/50] tests and examples

---
 q2_fmt/_examples.py              | 22 +++++++++++++++++
 q2_fmt/plugin_setup.py           |  5 +++-
 q2_fmt/tests/test_engraftment.py | 42 ++++++++++++++++++--------------
 3 files changed, 50 insertions(+), 19 deletions(-)

diff --git a/q2_fmt/_examples.py b/q2_fmt/_examples.py
index 331fe64..8ce2ff2 100644
--- a/q2_fmt/_examples.py
+++ b/q2_fmt/_examples.py
@@ -229,3 +229,25 @@ def peds_heatmap(use):
     )
 
     peds_heatmap.assert_output_type('Visualization')
+
+
+def bootstrap_peds_method(use):
+    md = use.init_metadata('md', peds_md_factory)
+    peds_table = use.init_artifact('peds_table', peds_feature_table_factory)
+
+    peds_stats, = use.action(
+        use.UsageAction('fmt', 'peds_bootstrap'),
+        use.UsageInputs(
+            table=peds_table,
+            metadata=md,
+            time_column='time_point',
+            reference_column='Donor',
+            subject_column='SubjectID'
+        ),
+        use.UsageOutputNames(
+            stats='stats'
+        )
+
+    )
+
+    peds_stats.assert_output_type("StatsTable[Pairwise]")
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 988e6a6..5231d1c 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -337,7 +337,10 @@
     },
     name='',
     description='',
-    citations=[citations['aggarwala_precise_2021']]
+    citations=[citations['aggarwala_precise_2021']],
+    examples={
+        'peds_methods': ex.bootstrap_peds_method
+    }
 )
 
 
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index dac2324..1169207 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1122,7 +1122,6 @@ def test_rename_features_with_blank_label(self):
                                        reference_column="Ref",
                                        subject_column="subject")
         _rename_features(data=feature_peds_df, level_delimiter=";")
-        print(feature_peds_df)
         Fs1 = feature_peds_df.set_index("id").at['Feature 1 __',
                                                  'subject']
         Fs2 = feature_peds_df.set_index("id").at['Feature 2',
@@ -1152,14 +1151,17 @@ def test_high_donor_overlap(self):
             'Feature1': [1, 0, 0, 1, 0, 0],
             'Feature2': [0, 1, 0, 0, 1, 0],
             'Feature3': [0, 0, 1, 0, 0, 1]}).set_index('id')
+        metadata = Metadata(metadata_df)
 
-        p, real_mean, fake_mean = peds_bootstrap(metadata_df=metadata_df,
-                                                 table=table_df,
-                                                 time_column="group",
-                                                 reference_column="Ref",
-                                                 subject_column="subject",
-                                                 iters=999)
-        self.assertGreater(real_mean, fake_mean)
+        stats = peds_bootstrap(metadata=metadata,
+                               table=table_df,
+                               time_column="group",
+                               reference_column="Ref",
+                               subject_column="subject",
+                               bootstrap_replicates=999)
+        real_median = stats["A:measure"].values
+        fake_median = stats["B:measure"].values
+        self.assertGreater(real_median, fake_median)
 
     def test_low_donor_overlap(self):
         metadata_df = pd.DataFrame({
@@ -1178,14 +1180,18 @@ def test_low_donor_overlap(self):
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
                    'donor1', 'donor2', 'donor3'],
-            'Feature1': [1, 0, 0, 0, 0, 1],
-            'Feature2': [0, 1, 0, 1, 0, 0],
-            'Feature3': [0, 0, 1, 0, 1, 0]}).set_index('id')
-
-        p, real_median, fake_median = peds_bootstrap(metadata_df=metadata_df,
-                                                     table=table_df,
-                                                     time_column="group",
-                                                     reference_column="Ref",
-                                                     subject_column="subject",
-                                                     iters=999)
+            'Feature1': [1, 0, 0, 0, 1, 1],
+            'Feature2': [0, 1, 0, 1, 0, 1],
+            'Feature3': [0, 0, 1, 1, 1, 0]}).set_index('id')
+        metadata = Metadata(metadata_df)
+
+        stats = peds_bootstrap(metadata=metadata,
+                               table=table_df,
+                               time_column="group",
+                               reference_column="Ref",
+                               subject_column="subject",
+                               bootstrap_replicates=999)
+
+        real_median = stats["A:measure"].values
+        fake_median = stats["B:measure"].values
         self.assertGreater(fake_median, real_median)

From e5cef3b80eec8a434d6825d67eca538ebc4f6ead Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 14 May 2024 11:37:56 -0700
Subject: [PATCH 09/50] merge py3.9 fixes

---
 merge py3.9 fixes | 6 ++++++
 1 file changed, 6 insertions(+)
 create mode 100644 merge py3.9 fixes

diff --git a/merge py3.9 fixes b/merge py3.9 fixes
new file mode 100644
index 0000000..66504fa
--- /dev/null
+++ b/merge py3.9 fixes	
@@ -0,0 +1,6 @@
+Merge branch 'dev' of https://github.com/qiime2/q2-fmt into peds-bootstrapping
+# Please enter a commit message to explain why this merge is necessary,
+# especially if it merges an updated upstream into a topic branch.
+#
+# Lines starting with '#' will be ignored, and an empty message aborts
+# the commit.

From cc45fd5c74f85839bd29398425d1bf06c1c797f3 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 14 May 2024 11:39:11 -0700
Subject: [PATCH 10/50] im alittle confused

---
 merge py3.9 fixes | 6 ------
 1 file changed, 6 deletions(-)
 delete mode 100644 merge py3.9 fixes

diff --git a/merge py3.9 fixes b/merge py3.9 fixes
deleted file mode 100644
index 66504fa..0000000
--- a/merge py3.9 fixes	
+++ /dev/null
@@ -1,6 +0,0 @@
-Merge branch 'dev' of https://github.com/qiime2/q2-fmt into peds-bootstrapping
-# Please enter a commit message to explain why this merge is necessary,
-# especially if it merges an updated upstream into a topic branch.
-#
-# Lines starting with '#' will be ignored, and an empty message aborts
-# the commit.

From 418fd4a2cec1f15ae730a5d61cea6a711890aba8 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 14:57:57 -0700
Subject: [PATCH 11/50] lint from conflict merge

---
 q2_fmt/plugin_setup.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index a80e4c6..18b6566 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -214,8 +214,7 @@
 
 plugin.visualizers.register_function(
     function=q2_fmt.peds_heatmap,
-
-   inputs={'data': Dist1D[Ordered, Matched] % Properties("peds"),
+    inputs={'data': Dist1D[Ordered, Matched] % Properties("peds"),
             'stats': StatsTable[Pairwise]},
     parameters={'level_delimiter': Str},
     parameter_descriptions={},

From fdec4e5b6216887c9360c9385614f2fea8376518 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 15:13:06 -0700
Subject: [PATCH 12/50] simulation not bootstrap

---
 q2_fmt/_peds.py        | 20 ++++++++++----------
 q2_fmt/plugin_setup.py | 18 +++++++++++++-----
 2 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 72d8d57..77a2c32 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -492,13 +492,13 @@ def _mask_recipient(donor_mask, recip_df):
     return maskedrecip
 
 
-def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
-                   time_column: str, reference_column: str,
-                   subject_column: str,
-                   filter_missing_references: bool = False,
-                   drop_incomplete_subjects: bool = False,
-                   drop_incomplete_timepoint: list = None,
-                   bootstrap_replicates: int = 999) -> (pd.DataFrame):
+def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
+                    time_column: str, reference_column: str,
+                    subject_column: str,
+                    filter_missing_references: bool = False,
+                    drop_incomplete_subjects: bool = False,
+                    drop_incomplete_timepoint: list = None,
+                    replicates: int = 999) -> (pd.DataFrame):
 
     metadata_df = metadata.to_dataframe()
     donor = metadata_df[metadata_df.index.isin(
@@ -507,7 +507,7 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
 
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
     fake_donor = pd.DataFrame([])
-    for i in range(0, bootstrap_replicates+1):
+    for i in range(0, replicates+1):
         if i == 0:
             peds = sample_peds(
                 table=table, metadata=metadata,
@@ -537,11 +537,11 @@ def peds_bootstrap(table: pd.DataFrame, metadata: qiime2.Metadata,
     # These are not being used currently but will be soon
     fake_donors_CLES = fake_donor.median(axis=0)
     agree = ((fake_donors_CLES <
-              real_temp.median()).sum())/bootstrap_replicates
+              real_temp.median()).sum())/replicates
     disagree = 1 - agree
     cles = agree - disagree
     cles
-
+    # TODO: REFACTOR STATISTICS
     s, p = mannwhitneyu(real_temp, fake_donor_series.to_list(),
                         alternative='greater')
 
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 18b6566..273ba3b 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -311,7 +311,7 @@
 )
 
 plugin.methods.register_function(
-    function=q2_fmt.peds_bootstrap,
+    function=q2_fmt.peds_simulation,
     inputs={'table': FeatureTable[Frequency | RelativeFrequency |
                                   PresenceAbsence]},
     parameters={'metadata': Metadata, 'time_column': Str,
@@ -319,7 +319,7 @@
                 'filter_missing_references': Bool,
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str],
-                'bootstrap_replicates': Int % Range(999, None)},
+                'replicates': Int % Range(999, None)},
     outputs=[('stats', StatsTable[Pairwise])],
     parameter_descriptions={
         'metadata': 'The sample metadata.',
@@ -340,13 +340,21 @@
         'drop_incomplete_subjects': 'Filter out subjects that do not have'
                                     ' a sample at every timepoint.'
                                     ' Default behavior is to raise an error.',
+        'drop_incomplete_timepoint': 'Filter out timepoints that do not have'
+                                    ' enough samples.'
+                                    ' Default behavior is to raise an error.',
+        'replicates': 'The number of replicates to run the Monte Carlo'
+                      ' simulation on'
     },
     output_descriptions={
-        'stats': 'stats table for comparing real PEDS values'
-        ' from fake Bootstrapped values'
+        'stats': 'stats table for comparing if the actual PEDS values'
+        ' to Shuffled simulated values'
     },
     name='',
-    description='',
+    description='A Monte Carlo Simulation that shuffles the recipient and'
+                ' donated microbiome in order to investigate if the true'
+                ' donated microbiome is the most similar to their true FMT'
+                ' recipient',
     citations=[citations['aggarwala_precise_2021']],
     examples={
         'peds_methods': ex.bootstrap_peds_method

From a0506de2d190b795f26378682a3baba6cc82ca10 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 15:56:23 -0700
Subject: [PATCH 13/50] simulation not bootstrap pt2

---
 q2_fmt/__init__.py               |  4 ++--
 q2_fmt/tests/test_engraftment.py | 28 ++++++++++++++--------------
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/q2_fmt/__init__.py b/q2_fmt/__init__.py
index 7e18587..fa5a3f1 100644
--- a/q2_fmt/__init__.py
+++ b/q2_fmt/__init__.py
@@ -9,10 +9,10 @@
 from ._version import get_versions
 from ._engraftment import engraftment, group_timepoints
 from ._peds import (sample_peds, feature_peds, peds, peds_heatmap,
-                    peds_bootstrap)
+                    peds_simulation)
 
 __version__ = get_versions()['version']
 del get_versions
 
 __all__ = ['engraftment', 'sample_peds', 'feature_peds',
-           'peds', 'peds_heatmap', 'group_timepoints', 'peds_bootstrap']
+           'peds', 'peds_heatmap', 'group_timepoints', 'peds_simulation']
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 70b729d..6d30bdb 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -20,7 +20,7 @@
                           _check_subject_column, _check_column_type,
                           _drop_incomplete_timepoints, feature_peds,
                           _check_column_missing, _rename_features,
-                          peds_bootstrap)
+                          peds_simulation)
 
 
 class TestBase(TestPluginBase):
@@ -1259,7 +1259,7 @@ def test_rename_features_with_blank_label(self):
         self.assertEqual("2", Fs2)
 
 
-class TestBoot(TestBase):
+class TestSim(TestBase):
     def test_high_donor_overlap(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
@@ -1282,12 +1282,12 @@ def test_high_donor_overlap(self):
             'Feature3': [0, 0, 1, 0, 0, 1]}).set_index('id')
         metadata = Metadata(metadata_df)
 
-        stats = peds_bootstrap(metadata=metadata,
-                               table=table_df,
-                               time_column="group",
-                               reference_column="Ref",
-                               subject_column="subject",
-                               bootstrap_replicates=999)
+        stats = peds_simulation(metadata=metadata,
+                                table=table_df,
+                                time_column="group",
+                                reference_column="Ref",
+                                subject_column="subject",
+                                bootstrap_replicates=999)
         real_median = stats["A:measure"].values
         fake_median = stats["B:measure"].values
         self.assertGreater(real_median, fake_median)
@@ -1314,12 +1314,12 @@ def test_low_donor_overlap(self):
             'Feature3': [0, 0, 1, 1, 1, 0]}).set_index('id')
         metadata = Metadata(metadata_df)
 
-        stats = peds_bootstrap(metadata=metadata,
-                               table=table_df,
-                               time_column="group",
-                               reference_column="Ref",
-                               subject_column="subject",
-                               bootstrap_replicates=999)
+        stats = peds_simulation(metadata=metadata,
+                                table=table_df,
+                                time_column="group",
+                                reference_column="Ref",
+                                subject_column="subject",
+                                replicates=999)
 
         real_median = stats["A:measure"].values
         fake_median = stats["B:measure"].values

From bf2fd5c10a0feed946dce745ae0adef18f0826d4 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 16:13:52 -0700
Subject: [PATCH 14/50] secret boots

---
 q2_fmt/tests/test_engraftment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 6d30bdb..0f1d52c 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1287,7 +1287,7 @@ def test_high_donor_overlap(self):
                                 time_column="group",
                                 reference_column="Ref",
                                 subject_column="subject",
-                                bootstrap_replicates=999)
+                                replicates=999)
         real_median = stats["A:measure"].values
         fake_median = stats["B:measure"].values
         self.assertGreater(real_median, fake_median)

From 75ddf76e59dcb6d4b82c1679738dbce34a40dfaa Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 16:25:58 -0700
Subject: [PATCH 15/50] ex

---
 q2_fmt/_examples.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/q2_fmt/_examples.py b/q2_fmt/_examples.py
index 0795150..7ac6960 100644
--- a/q2_fmt/_examples.py
+++ b/q2_fmt/_examples.py
@@ -236,12 +236,12 @@ def peds_heatmap(use):
     peds_heatmap.assert_output_type('Visualization')
 
 
-def bootstrap_peds_method(use):
+def simulation_peds_method(use):
     md = use.init_metadata('md', peds_md_factory)
     peds_table = use.init_artifact('peds_table', peds_feature_table_factory)
 
     peds_stats, = use.action(
-        use.UsageAction('fmt', 'peds_bootstrap'),
+        use.UsageAction('fmt', 'peds_simulation'),
         use.UsageInputs(
             table=peds_table,
             metadata=md,

From b1039533895ca9f6fa0b0ba43f0726e6069f7f28 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 16:35:48 -0700
Subject: [PATCH 16/50] imports

---
 q2_fmt/plugin_setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 273ba3b..c5a8772 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -357,7 +357,7 @@
                 ' recipient',
     citations=[citations['aggarwala_precise_2021']],
     examples={
-        'peds_methods': ex.bootstrap_peds_method
+        'peds_methods': ex.simulation_peds_method
     }
 )
 

From 27ed8223cdd626051cb5ce0c9a3054a32228a027 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 17:32:10 -0700
Subject: [PATCH 17/50] np.nan

---
 q2_fmt/_peds.py                  |  51 ++---
 q2_fmt/tests/test_engraftment.py | 314 +++++++++++++++----------------
 2 files changed, 183 insertions(+), 182 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 77a2c32..cebf876 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -507,31 +507,32 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
 
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
     fake_donor = pd.DataFrame([])
-    for i in range(0, replicates+1):
-        if i == 0:
-            peds = sample_peds(
-                table=table, metadata=metadata,
-                time_column=time_column,
-                reference_column=reference_column,
-                subject_column=subject_column,
-                filter_missing_references=filter_missing_references,
-                drop_incomplete_subjects=drop_incomplete_subjects,
-                drop_incomplete_timepoint=drop_incomplete_timepoint)
-            real_temp = peds["measure"]
-        else:
-            shifted_list = recipient[reference_column].sample(frac=1).to_list()
-            recipient.loc[:, reference_column] = shifted_list
-            metadata_df = pd.concat([donor, recipient])
-            metadata = qiime2.Metadata(metadata_df)
-            peds = sample_peds(
-                table=table, metadata=metadata,
-                time_column=time_column,
-                reference_column=reference_column,
-                subject_column=subject_column,
-                filter_missing_references=filter_missing_references,
-                drop_incomplete_subjects=drop_incomplete_subjects,
-                drop_incomplete_timepoint=drop_incomplete_timepoint)
-            fake_donor[i] = peds["measure"]
+    peds = sample_peds(
+           table=table, metadata=metadata,
+           time_column=time_column,
+           reference_column=reference_column,
+           subject_column=subject_column,
+           filter_missing_references=filter_missing_references,
+           drop_incomplete_subjects=drop_incomplete_subjects,
+           drop_incomplete_timepoint=drop_incomplete_timepoint)
+    real_temp = peds["measure"]
+    for i in range(replicates):
+        recipient[reference_column].sample(frac=1).reset_index(drop=True)
+        # TODO: Deprecated Below:
+        shuffled_list = recipient[reference_column].sample(frac=1).to_list()
+        recipient.loc[:, reference_column] = shuffled_list
+        metadata_df = pd.concat([donor, recipient])
+        # TODO: ^^ Deprecated
+        metadata = qiime2.Metadata(metadata_df)
+        peds = sample_peds(
+               table=table, metadata=metadata,
+               time_column=time_column,
+               reference_column=reference_column,
+               subject_column=subject_column,
+               filter_missing_references=filter_missing_references,
+               drop_incomplete_subjects=drop_incomplete_subjects,
+               drop_incomplete_timepoint=drop_incomplete_timepoint)
+        fake_donor[i] = peds["measure"]
     fake_donor_series = fake_donor.median(axis=1)
     # Common Langauge effect size calcs in prep for stats refactor.
     # These are not being used currently but will be soon
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 0f1d52c..7b7f8a4 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -631,12 +631,12 @@ def test_get_donor(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor2', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub2', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 1, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor2', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub2', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 1, 2, np.nan,
+                      np.nan]}).set_index('id')
         reference_series = metadata_df['Ref'].dropna()
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -662,12 +662,12 @@ def test_get_subject(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor2', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub2', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 1, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor2', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub2', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 1, 2, np.nan,
+                      np.nan]}).set_index('id')
         reference_series = metadata_df['Ref'].dropna()
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -693,12 +693,12 @@ def test_timepoint(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor2', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub2', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 1, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor2', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub2', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 1, 2, np.nan,
+                      np.nan]}).set_index('id')
         reference_series = metadata_df['Ref'].dropna()
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -724,12 +724,12 @@ def test_no_donors(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': [float("Nan"), float("Nan"), float("Nan"), float("Nan"),
-                    float("Nan"), float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub2', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 1, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': [np.nan, np.nan, np.nan, np.nan,
+                    np.nan, np.nan],
+            'subject': ['sub1', 'sub1', 'sub2', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 1, 2, np.nan,
+                      np.nan]}).set_index('id')
         reference_series = metadata_df['Ref']
         with self.assertRaisesRegex(KeyError, 'Missing references for'
                                     ' the associated sample data. Please make'
@@ -747,12 +747,12 @@ def test_incomplete_timepoints(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -777,12 +777,12 @@ def test_incomplete_timepoints_with_flag(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -811,12 +811,12 @@ def test_incorrect_reference_column_name(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         with self.assertRaisesRegex(KeyError, ".*the provided"
                                     " `--p-reference-column`: `R` in the"
                                     " metadata"):
@@ -826,12 +826,12 @@ def test_incorrect_group_column_name(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         with self.assertRaisesRegex(KeyError,
                                     ".*the provided `--p-time-column`: `time`"
                                     " in the metadata"):
@@ -841,12 +841,12 @@ def test_incorrect_subject_column_name(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         with self.assertRaisesRegex(KeyError, ".*the provided"
                                     " `--p-subject-column`: `sub` in the"
                                     " metadata"):
@@ -856,12 +856,12 @@ def test_no_feature_overlap(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -889,12 +889,12 @@ def test_feature_overlap(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -921,12 +921,12 @@ def test_peds_calc(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -954,12 +954,12 @@ def test_peds_calc(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -978,12 +978,12 @@ def test_unique_subjects_in_timepoints(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 2, 1, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 2, 1, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -1003,9 +1003,9 @@ def test_feature_peds_calc(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
                    'donor1'],
-            'Ref': ['donor1', 'donor1', 'donor1', float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', float("Nan")],
-            'group': [1, 1, 1, float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', np.nan],
+            'group': [1, 1, 1, np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
@@ -1028,12 +1028,12 @@ def test_sample_id_match(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 2, 1, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 2, 1, np.nan,
+                      np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
             'id': ['s1', 's2', 's3', 's4',
@@ -1053,12 +1053,12 @@ def test_column_type_nonnumeric(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': ["t1", "t2", "t3", "t2", float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': ["t1", "t2", "t3", "t2", np.nan,
+                      np.nan]}).set_index('id')
         metadata_obj = Metadata(metadata_df)
         column_properties = metadata_obj.columns
         with self.assertRaisesRegex(AssertionError, ".*Column with non-numeric"
@@ -1069,12 +1069,12 @@ def test_column_type_noncategorical(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': [1, 1, 1, 2, float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': [1, 1, 1, 2, np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata_obj = Metadata(metadata_df)
         column_properties = metadata_obj.columns
         with self.assertRaisesRegex(AssertionError, ".*Column with"
@@ -1087,12 +1087,12 @@ def test_reference_series_not_in_table(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['1', '1', '2', '2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub2', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 1, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['1', '1', '2', '2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub2', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 1, 2, np.nan,
+                      np.nan]}).set_index('id')
         reference_series = metadata_df['Ref'].dropna()
         table_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
@@ -1116,12 +1116,12 @@ def test_column_name_is_ID(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         with self.assertRaisesRegex(KeyError, ".*`--p-subject-column` can not"
                                     " be the same as the index of"
                                     " metadata: `id`"):
@@ -1131,12 +1131,12 @@ def test_drop_incomplete_timepoints(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata_df = _drop_incomplete_timepoints(metadata_df, "group", [3])
         self.assertEqual(metadata_df["group"].unique()[0], float(1))
         self.assertEqual(metadata_df["group"].unique()[1], float(2))
@@ -1145,12 +1145,12 @@ def test_drop_incomplete_timepoints_list(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3', 'sample4',
                    'donor1', 'donor2'],
-            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub1', 'sub1', 'sub2', float("Nan"),
-                        float("Nan")],
-            'group': [1, 2, 3, 2, float("Nan"),
-                      float("Nan")]}).set_index('id')
+            'Ref': ['donor1', 'donor1', 'donor1', 'donor2', np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub1', 'sub1', 'sub2', np.nan,
+                        np.nan],
+            'group': [1, 2, 3, 2, np.nan,
+                      np.nan]}).set_index('id')
         metadata_df = _drop_incomplete_timepoints(metadata_df, "group", [3, 2])
         self.assertEqual(metadata_df["group"].dropna().unique(), [float(1)])
 
@@ -1158,9 +1158,9 @@ def test_rename_features_with_delim(self):
         metadata_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
                        'donor1'],
-                'Ref': ['donor1', 'donor1', 'donor1', float("Nan")],
-                'subject': ['sub1', 'sub1', 'sub1', float("Nan")],
-                'group': [1, 1, 1, float("Nan")]}).set_index('id')
+                'Ref': ['donor1', 'donor1', 'donor1', np.nan],
+                'subject': ['sub1', 'sub1', 'sub1', np.nan],
+                'group': [1, 1, 1, np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
@@ -1184,9 +1184,9 @@ def test_rename_features_with_no_delim(self):
         metadata_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
                        'donor1'],
-                'Ref': ['donor1', 'donor1', 'donor1', float("Nan")],
-                'subject': ['sub1', 'sub1', 'sub1', float("Nan")],
-                'group': [1, 1, 1, float("Nan")]}).set_index('id')
+                'Ref': ['donor1', 'donor1', 'donor1', np.nan],
+                'subject': ['sub1', 'sub1', 'sub1', np.nan],
+                'group': [1, 1, 1, np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
@@ -1210,9 +1210,9 @@ def test_rename_features_with_wrong_delim(self):
         metadata_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
                        'donor1'],
-                'Ref': ['donor1', 'donor1', 'donor1', float("Nan")],
-                'subject': ['sub1', 'sub1', 'sub1', float("Nan")],
-                'group': [1, 1, 1, float("Nan")]}).set_index('id')
+                'Ref': ['donor1', 'donor1', 'donor1', np.nan],
+                'subject': ['sub1', 'sub1', 'sub1', np.nan],
+                'group': [1, 1, 1, np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
@@ -1236,9 +1236,9 @@ def test_rename_features_with_blank_label(self):
         metadata_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
                        'donor1'],
-                'Ref': ['donor1', 'donor1', 'donor1', float("Nan")],
-                'subject': ['sub1', 'sub1', 'sub1', float("Nan")],
-                'group': [1, 1, 1, float("Nan")]}).set_index('id')
+                'Ref': ['donor1', 'donor1', 'donor1', np.nan],
+                'subject': ['sub1', 'sub1', 'sub1', np.nan],
+                'group': [1, 1, 1, np.nan]}).set_index('id')
         metadata = Metadata(metadata_df)
         table_df = pd.DataFrame({
                 'id': ['sample1', 'sample2', 'sample3',
@@ -1264,14 +1264,14 @@ def test_high_donor_overlap(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
                    'donor1', 'donor2', 'donor3'],
-            'Ref': ['donor1', 'donor2', 'donor3', float("Nan"), float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub2', 'sub3', float("Nan"), float("Nan"),
-                        float("Nan")],
-            'group': [1, 1, 1, float("Nan"), float("Nan"),
-                      float("Nan")],
-            "Location": [float("Nan"), float("Nan"),
-                         float("Nan"), 'test', 'test',
+            'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan,
+                        np.nan],
+            'group': [1, 1, 1, np.nan, np.nan,
+                      np.nan],
+            "Location": [np.nan, np.nan,
+                         np.nan, 'test', 'test',
                          'test']}).set_index('id')
 
         table_df = pd.DataFrame({
@@ -1296,14 +1296,14 @@ def test_low_donor_overlap(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
                    'donor1', 'donor2', 'donor3'],
-            'Ref': ['donor1', 'donor2', 'donor3', float("Nan"), float("Nan"),
-                    float("Nan")],
-            'subject': ['sub1', 'sub2', 'sub3', float("Nan"), float("Nan"),
-                        float("Nan")],
-            'group': [1, 1, 1, float("Nan"), float("Nan"),
-                      float("Nan")],
-            "Location": [float("Nan"), float("Nan"),
-                         float("Nan"), 'test', 'test',
+            'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan,
+                        np.nan],
+            'group': [1, 1, 1, np.nan, np.nan,
+                      np.nan],
+            "Location": [np.nan, np.nan,
+                         np.nan, 'test', 'test',
                          'test']}).set_index('id')
 
         table_df = pd.DataFrame({

From 4280d8986ba85f8ecc0cc1a1f34e32a540ec9858 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 19:25:17 -0700
Subject: [PATCH 18/50] address small comments and start stats

---
 q2_fmt/_peds.py | 94 ++++++++++++++++++++++++++++++-------------------
 1 file changed, 57 insertions(+), 37 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index cebf876..997c167 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -506,7 +506,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                                 reference_column=reference_column))]
 
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
-    fake_donor = pd.DataFrame([])
+    shuffled_donor = pd.DataFrame([])
     peds = sample_peds(
            table=table, metadata=metadata,
            time_column=time_column,
@@ -515,15 +515,10 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
            filter_missing_references=filter_missing_references,
            drop_incomplete_subjects=drop_incomplete_subjects,
            drop_incomplete_timepoint=drop_incomplete_timepoint)
-    real_temp = peds["measure"]
+    actual_temp = peds["measure"]
     for i in range(replicates):
-        recipient[reference_column].sample(frac=1).reset_index(drop=True)
-        # TODO: Deprecated Below:
-        shuffled_list = recipient[reference_column].sample(frac=1).to_list()
-        recipient.loc[:, reference_column] = shuffled_list
-        metadata_df = pd.concat([donor, recipient])
-        # TODO: ^^ Deprecated
-        metadata = qiime2.Metadata(metadata_df)
+        metadata = _shuffle_donor_associations(recipient, reference_column,
+                                               donor)
         peds = sample_peds(
                table=table, metadata=metadata,
                time_column=time_column,
@@ -532,48 +527,73 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                filter_missing_references=filter_missing_references,
                drop_incomplete_subjects=drop_incomplete_subjects,
                drop_incomplete_timepoint=drop_incomplete_timepoint)
-        fake_donor[i] = peds["measure"]
-    fake_donor_series = fake_donor.median(axis=1)
+        shuffled_donor[i] = peds["measure"]
+    shuffled_donor_series = shuffled_donor.median(axis=1)
     # Common Langauge effect size calcs in prep for stats refactor.
     # These are not being used currently but will be soon
-    fake_donors_CLES = fake_donor.median(axis=0)
-    agree = ((fake_donors_CLES <
-              real_temp.median()).sum())/replicates
+
+    # Calculating per-subject p-values
+    actual_element_wise = actual_temp.values[:, None]
+    agree_df = shuffled_donor < actual_element_wise
+    agree = agree_df.sum(axis=1)
+    # adding 1 here because you can mathmatically can get p-value of 0 from a 
+    # Monte Carlo Simulation
     disagree = 1 - agree
+    per_subject_p = (agree + 1)/(replicates+1)
+    raise ValueError(per_subject_p)
+    
     cles = agree - disagree
-    cles
+    cles()
     # TODO: REFACTOR STATISTICS
-    s, p = mannwhitneyu(real_temp, fake_donor_series.to_list(),
+    s, p = mannwhitneyu(actual_temp, shuffled_donor_series.to_list(),
                         alternative='greater')
-
-    stats_df = pd.DataFrame([["real_values", real_temp.size,
-                              real_temp.median(), "fake_values",
-                              fake_donor_series.size,
-                              fake_donor_series.median(),
-                              fake_donor_series.size,
-                              s, p, np.nan]])
-    stats_df.columns = ['A:group', 'A:n', 'A:measure', 'B:group', 'B:n',
-                        'B:measure', 'n', 'test-statistic', 'p-value',
-                        'q-value']
-
-    stats_df['A:group'].attrs.update({'title': 'real_values',
-                                     'description': '...'})
-    stats_df['B:group'].attrs.update({'title': 'randomized_values',
-                                     'description': '...'})
-    n = {'title': 'count', 'description': '...'}
+    p = 50
+
+    stats_df = pd.DataFrame([["actual_values", actual_temp.size,
+                              actual_temp.median(), "shuffled_values",
+                              shuffled_donor_series.size,
+                              shuffled_donor_series.median(),
+                              shuffled_donor_series.size,
+                              s, p, np.nan]],
+                            columns=['A:group', 'A:n', 'A:measure',
+                                     'B:group', 'B:n', 'B:measure', 'n',
+                                     'test-statistic', 'p-value', 'q-value'])
+
+    stats_df['A:group'].attrs.update({'title': 'actual_values',
+                                     'description': 'PEDS values calculated'
+                                      ' with actual recipient and donated'
+                                      ' microbiome pairing'})
+    stats_df['B:group'].attrs.update({'title': 'shuffled_values',
+                                     'description': 'PEDS values calculated'
+                                      ' with shuffled recupient and donated'
+                                      ' microbiome pairings'})
+    n = {'title': 'count', 'description': 'Number of recipients and donated'
+         ' microbiome pairings'}
     stats_df['A:n'].attrs.update(n)
     stats_df['B:n'].attrs.update(n)
     measure = {
         'title': 'Median PEDS Value',
-        'description': '...'
+        'description': 'Median PEDS Value'
     }
     stats_df['A:measure'].attrs.update(measure)
     stats_df['B:measure'].attrs.update(measure)
-    stats_df['n'].attrs.update(dict(title='count', description='...'))
+    stats_df['n'].attrs.update(dict(title='count', description='Number of'
+                                    'comparisons'))
     stats_df['test-statistic'].attrs.update(dict(title='Mann-Whitney U',
-                                                 description='...'))
-    stats_df['p-value'].attrs.update(dict(title='one-tail', description='...'))
+                                                 description='Mann-Whitney U'
+                                                             'test statistic'))
+    stats_df['p-value'].attrs.update(dict(title='one-tail p-value',
+                                          description='one-tail p-value'))
     stats_df['q-value'].attrs.update(
-        dict(title='Benjamini–Hochberg', description='...'))
+        dict(title='Benjamini–Hochberg', description='FDR corrections using'
+             'Benjamini–Hochberg procedure'))
 
     return stats_df
+
+
+def _shuffle_donor_associations(recipient, reference_column, donor):
+    shuffled_list = recipient[reference_column].sample(frac=1).to_list()
+    recipient.loc[:, reference_column] = shuffled_list
+    metadata_df = pd.concat([donor, recipient])
+    metadata = qiime2.Metadata(metadata_df)
+    return metadata

From 69f842a5414a94ceee0e3139c06418c711ccfe2a Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 20:16:53 -0700
Subject: [PATCH 19/50] added per-subject-stats

---
 q2_fmt/_peds.py | 93 +++++++++++++++++++++++++------------------------
 1 file changed, 47 insertions(+), 46 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 997c167..77cf977 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -10,7 +10,7 @@
 import pandas as pd
 import numpy as np
 import warnings
-from scipy.stats import mannwhitneyu
+from scipy.stats import false_discovery_control, combine_pvalues
 from collections import Counter
 import os
 import pkg_resources
@@ -528,67 +528,68 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                drop_incomplete_subjects=drop_incomplete_subjects,
                drop_incomplete_timepoint=drop_incomplete_timepoint)
         shuffled_donor[i] = peds["measure"]
-    shuffled_donor_series = shuffled_donor.median(axis=1)
-    # Common Langauge effect size calcs in prep for stats refactor.
-    # These are not being used currently but will be soon
 
     # Calculating per-subject p-values
     actual_element_wise = actual_temp.values[:, None]
-    agree_df = shuffled_donor < actual_element_wise
-    agree = agree_df.sum(axis=1)
-    # adding 1 here because you can mathmatically can get p-value of 0 from a 
+    agree_df = shuffled_donor > actual_element_wise
+    agree_series = agree_df.sum(axis=1)
+    # adding 1 here because you can mathmatically can get p-value of 0 from a
     # Monte Carlo Simulation
-    disagree = 1 - agree
-    per_subject_p = (agree + 1)/(replicates+1)
-    raise ValueError(per_subject_p)
-    
-    cles = agree - disagree
-    cles()
-    # TODO: REFACTOR STATISTICS
-    s, p = mannwhitneyu(actual_temp, shuffled_donor_series.to_list(),
-                        alternative='greater')
-    p = 50
-
-    stats_df = pd.DataFrame([["actual_values", actual_temp.size,
-                              actual_temp.median(), "shuffled_values",
-                              shuffled_donor_series.size,
-                              shuffled_donor_series.median(),
-                              shuffled_donor_series.size,
-                              s, p, np.nan]],
-                            columns=['A:group', 'A:n', 'A:measure',
-                                     'B:group', 'B:n', 'B:measure', 'n',
-                                     'test-statistic', 'p-value', 'q-value'])
-
-    stats_df['A:group'].attrs.update({'title': 'actual_values',
-                                     'description': 'PEDS values calculated'
-                                      ' with actual recipient and donated'
-                                      ' microbiome pairing'})
-    stats_df['B:group'].attrs.update({'title': 'shuffled_values',
-                                     'description': 'PEDS values calculated'
-                                      ' with shuffled recupient and donated'
-                                      ' microbiome pairings'})
+    per_subject_p = (agree_series + 1)/(replicates+1)
+    per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
+    # Common Langauge effect size calcs in prep for stats refactor.
+    disagree = 1 - per_subject_p
+    rank_biserial = per_subject_p - disagree
+
+    per_sub_stats = pd.DataFrame({'A:group': "actual_values",
+                                 'A:n': actual_temp.size,
+                                  'A:measure': actual_temp,
+                                  'B:group': "shuffled_values",
+                                  'B:n': replicates,
+                                  'B:measure': shuffled_donor.median(axis=1),
+                                  'n': replicates*actual_temp.size,
+                                  'test-statistic': agree_series,
+                                  'p-value': per_subject_p,
+                                  'q-value': per_subject_q,
+                                  'disagree': disagree,
+                                  'rank-biserial-correlation': rank_biserial})
+    per_sub_stats['A:group'].attrs.update({'title': 'actual_values',
+                                          'description': 'PEDS values'
+                                           ' calculated with actual recipient'
+                                           ' and donated microbiome pairing'})
+    per_sub_stats['B:group'].attrs.update({'title': 'shuffled_values',
+                                          'description': 'PEDS values'
+                                           ' calculated with shuffled'
+                                           ' recipient and donated microbiome'
+                                           ' pairings'})
     n = {'title': 'count', 'description': 'Number of recipients and donated'
          ' microbiome pairings'}
-    stats_df['A:n'].attrs.update(n)
-    stats_df['B:n'].attrs.update(n)
+    per_sub_stats['A:n'].attrs.update(n)
+    per_sub_stats['B:n'].attrs.update(n)
     measure = {
         'title': 'Median PEDS Value',
         'description': 'Median PEDS Value'
     }
-    stats_df['A:measure'].attrs.update(measure)
-    stats_df['B:measure'].attrs.update(measure)
-    stats_df['n'].attrs.update(dict(title='count', description='Number of'
+    per_sub_stats['A:measure'].attrs.update(measure)
+    per_sub_stats['B:measure'].attrs.update(measure)
+    per_sub_stats['n'].attrs.update(dict(title='count', description='Number of'
                                     'comparisons'))
-    stats_df['test-statistic'].attrs.update(dict(title='Mann-Whitney U',
+    per_sub_stats['test-statistic'].attrs.update(dict(title='Mann-Whitney U',
                                                  description='Mann-Whitney U'
                                                              'test statistic'))
-    stats_df['p-value'].attrs.update(dict(title='one-tail p-value',
+    per_sub_stats['p-value'].attrs.update(dict(title='one-tail p-value',
                                           description='one-tail p-value'))
-    stats_df['q-value'].attrs.update(
+    per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
-
-    return stats_df
+    per_sub_stats['disagree'].attrs.update(
+        dict(title='values that agree with the null hypothesis',
+             description='values that agree with the null hypothesis'))
+    per_sub_stats['rank-biserial-correlation'].attrs.update(
+        dict(title='rank biserial correlation',
+             description='rank biserial correlation'))
+
+    return per_sub_stats
 
 
 def _shuffle_donor_associations(recipient, reference_column, donor):

From 18b88117ec6f5683265e7268f745b1db67347582 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 10 Jul 2024 21:22:09 -0700
Subject: [PATCH 20/50] stats refactor

---
 q2_fmt/_peds.py        | 71 +++++++++++++++++++++++++++++++++++-------
 q2_fmt/plugin_setup.py |  7 +++--
 2 files changed, 65 insertions(+), 13 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 77cf977..7f55b8b 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -498,14 +498,30 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                     filter_missing_references: bool = False,
                     drop_incomplete_subjects: bool = False,
                     drop_incomplete_timepoint: list = None,
-                    replicates: int = 999) -> (pd.DataFrame):
+                    replicates: int = 999) -> (pd.DataFrame, pd.DataFrame):
 
     metadata_df = metadata.to_dataframe()
     donor = metadata_df[metadata_df.index.isin(
         _check_reference_column(metadata=metadata_df,
                                 reference_column=reference_column))]
-
+    try:
+        donor.index.unique() != 1
+    except AssertionError as e:
+        raise AssertionError("There is only one donated microbiome in your"
+                             " data. A Monte Carlo simulation shuffles"
+                             " your donated microbiome and recipient pairing"
+                             "  and needs more than one donated microbiome"
+                             " to successfully shuffle") from e
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
+    try:
+        recipient.index.unique() != 1
+    except AssertionError as e:
+        raise AssertionError("There is only one recipient in your"
+                             " data. A Monte Carlo simulation shuffles"
+                             " your donated microbiome and recipient pairing"
+                             "  and needs more than one recipient"
+                             " to successfully shuffle") from e
+
     shuffled_donor = pd.DataFrame([])
     peds = sample_peds(
            table=table, metadata=metadata,
@@ -528,6 +544,21 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                drop_incomplete_subjects=drop_incomplete_subjects,
                drop_incomplete_timepoint=drop_incomplete_timepoint)
         shuffled_donor[i] = peds["measure"]
+        per_sub_stats = _per_subject_stats(actual_temp, shuffled_donor,
+                                           replicates)
+        global_stats = _global_stats(per_sub_stats['p-value'])
+    return per_sub_stats, global_stats
+
+
+def _shuffle_donor_associations(recipient, reference_column, donor):
+    shuffled_list = recipient[reference_column].sample(frac=1).to_list()
+    recipient.loc[:, reference_column] = shuffled_list
+    metadata_df = pd.concat([donor, recipient])
+    metadata = qiime2.Metadata(metadata_df)
+    return metadata
+
+
+def _per_subject_stats(actual_temp, shuffled_donor, replicates): 
 
     # Calculating per-subject p-values
     actual_element_wise = actual_temp.values[:, None]
@@ -574,8 +605,9 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     per_sub_stats['B:measure'].attrs.update(measure)
     per_sub_stats['n'].attrs.update(dict(title='count', description='Number of'
                                     'comparisons'))
-    per_sub_stats['test-statistic'].attrs.update(dict(title='Mann-Whitney U',
-                                                 description='Mann-Whitney U'
+    # TODO: ask about this
+    per_sub_stats['test-statistic'].attrs.update(dict(title='...',
+                                                 description='....'
                                                              'test statistic'))
     per_sub_stats['p-value'].attrs.update(dict(title='one-tail p-value',
                                           description='one-tail p-value'))
@@ -588,13 +620,30 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     per_sub_stats['rank-biserial-correlation'].attrs.update(
         dict(title='rank biserial correlation',
              description='rank biserial correlation'))
-
     return per_sub_stats
 
 
-def _shuffle_donor_associations(recipient, reference_column, donor):
-    shuffled_list = recipient[reference_column].sample(frac=1).to_list()
-    recipient.loc[:, reference_column] = shuffled_list
-    metadata_df = pd.concat([donor, recipient])
-    metadata = qiime2.Metadata(metadata_df)
-    return metadata
+def _global_stats(p_series):
+    stats, p = combine_pvalues(p_series, method="stouffer")
+    global_stats = pd.DataFrame([["p-values", p_series.size, stats, p,
+                                  np.nan]],
+                                columns=['Measure', 'n', 'test-statistic',
+                                         'p-value', 'q-value'])
+    global_stats['Measure'].attrs.update({'title': ('p-value'),
+                                          'description': ("p-value")})
+    global_stats['test-statistic'].attrs.update(dict(title="Stouffer's",
+                                                description="Stouffer's Z"
+                                                " score method"))
+    global_stats['p-value'].attrs.update(dict(title='one-tail p-value',
+                                         description='one-tail p-value'))
+    global_stats['q-value'].attrs.update(
+        dict(title='Benjamini–Hochberg', description='FDR corrections using'
+             'Benjamini–Hochberg procedure'))
+    return global_stats
+    # TODO: Do I need?
+    # per_sub_stats['disagree'].attrs.update(
+    #   dict(title='values that agree with the null hypothesis',
+    #         description='values that agree with the null hypothesis'))
+    # per_sub_stats['rank-biserial-correlation'].attrs.update(
+    #    dict(title='rank biserial correlation',
+    #         description='rank biserial correlation'))
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index c5a8772..4c3926d 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -320,7 +320,8 @@
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str],
                 'replicates': Int % Range(999, None)},
-    outputs=[('stats', StatsTable[Pairwise])],
+    outputs=[('per-subject-stats', StatsTable[Pairwise]),
+             ('global-stats', StatsTable[Pairwise])],
     parameter_descriptions={
         'metadata': 'The sample metadata.',
         'time_column': 'The column within the `metadata` that the'
@@ -347,7 +348,9 @@
                       ' simulation on'
     },
     output_descriptions={
-        'stats': 'stats table for comparing if the actual PEDS values'
+        'per-subject-stats': 'stats table for comparing if the actual PEDS'
+        ' values to shuffled simulated values per subject',
+        'global-stats': 'stats table for comparing if the actual PEDS values'
         ' to Shuffled simulated values'
     },
     name='',

From 426ddbfa23e8d72f492e10096f8316056ae2f9c8 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 11 Jul 2024 14:59:29 -0700
Subject: [PATCH 21/50] more stats

---
 q2_fmt/_peds.py        | 52 +++++++++++++++++++++---------------------
 q2_fmt/plugin_setup.py |  9 ++++++--
 2 files changed, 33 insertions(+), 28 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 7f55b8b..945a3bc 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -70,10 +70,14 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
 
 
 def peds_heatmap(output_dir: str, data: pd.DataFrame,
-                 level_delimiter: str = None, stats: pd.DataFrame = None):
+                 level_delimiter: str = None, 
+                 per_subject_stats: pd.DataFrame = None,
+                 global_stats):
     _rename_features(data=data, level_delimiter=level_delimiter)
-    if stats is not None:
-        table1, stats = _make_stats(stats)
+    if  per_subject_stats is not None:
+        table1, psstats = _make_stats(per_subject_stats)
+    if  global_stats is not None:
+        table2, gstats = _make_stats(global_stats)
     J_ENV = jinja2.Environment(
         loader=jinja2.PackageLoader('q2_fmt', 'assets')
     )
@@ -111,8 +115,10 @@ def peds_heatmap(output_dir: str, data: pd.DataFrame,
 
     with open(os.path.join(output_dir, 'index.html'), 'w') as fh:
         spec_string = json.dumps(full_spec)
-        fh.write(index.render(spec=spec_string, stats=stats,
-                              table1=table1))
+        fh.write(index.render(spec=spec_string, 
+                              per_subject_stats=psstats,
+                              global_stats=gstats,
+                              table1=table1, table2=table2))
 
 
 def sample_peds(table: pd.DataFrame, metadata: qiime2.Metadata,
@@ -562,15 +568,16 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
 
     # Calculating per-subject p-values
     actual_element_wise = actual_temp.values[:, None]
-    agree_df = shuffled_donor > actual_element_wise
-    agree_series = agree_df.sum(axis=1)
+    disagree_df = shuffled_donor >= actual_element_wise
+    disagree_series = disagree_df.sum(axis=1)
+    agree_series=replicates-disagree_series
     # adding 1 here because you can mathmatically can get p-value of 0 from a
     # Monte Carlo Simulation
-    per_subject_p = (agree_series + 1)/(replicates+1)
+    per_subject_p = (disagree_series + 1)/(replicates+1)
     per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
     # Common Langauge effect size calcs in prep for stats refactor.
-    disagree = 1 - per_subject_p
-    rank_biserial = per_subject_p - disagree
+    agree = 1 - per_subject_p
+    rank_biserial = agree - per_subject_p
 
     per_sub_stats = pd.DataFrame({'A:group': "actual_values",
                                  'A:n': actual_temp.size,
@@ -582,7 +589,7 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
                                   'q-value': per_subject_q,
-                                  'disagree': disagree,
+                                  'CLES': agree,
                                   'rank-biserial-correlation': rank_biserial})
     per_sub_stats['A:group'].attrs.update({'title': 'actual_values',
                                           'description': 'PEDS values'
@@ -605,21 +612,21 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_sub_stats['B:measure'].attrs.update(measure)
     per_sub_stats['n'].attrs.update(dict(title='count', description='Number of'
                                     'comparisons'))
-    # TODO: ask about this
-    per_sub_stats['test-statistic'].attrs.update(dict(title='...',
-                                                 description='....'
-                                                             'test statistic'))
+    per_sub_stats['test-statistic'].attrs.update(
+        dict(title='Iteration',
+             description='Number of iterations that agree with the ALTERNATIVE'
+                         'hypothesis'))
     per_sub_stats['p-value'].attrs.update(dict(title='one-tail p-value',
                                           description='one-tail p-value'))
     per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
-    per_sub_stats['disagree'].attrs.update(
-        dict(title='values that agree with the null hypothesis',
-             description='values that agree with the null hypothesis'))
+    per_sub_stats['CLES'].attrs.update(
+        dict(title='Common Language Effect Size',
+             description='Common Language Effect Size'))
     per_sub_stats['rank-biserial-correlation'].attrs.update(
         dict(title='rank biserial correlation',
-             description='rank biserial correlation'))
+             description='rank biserial correlation ranging from -1 to 1'))
     return per_sub_stats
 
 
@@ -640,10 +647,3 @@ def _global_stats(p_series):
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
     return global_stats
-    # TODO: Do I need?
-    # per_sub_stats['disagree'].attrs.update(
-    #   dict(title='values that agree with the null hypothesis',
-    #         description='values that agree with the null hypothesis'))
-    # per_sub_stats['rank-biserial-correlation'].attrs.update(
-    #    dict(title='rank biserial correlation',
-    #         description='rank biserial correlation'))
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 4c3926d..a311131 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -215,7 +215,8 @@
 plugin.visualizers.register_function(
     function=q2_fmt.peds_heatmap,
     inputs={'data': Dist1D[Ordered, Matched] % Properties("peds"),
-            'stats': StatsTable[Pairwise]},
+            'per-subject-stats': StatsTable[Pairwise],
+            'global-stats': StatsTable[Pairwise]},
     parameters={'level_delimiter': Str},
     parameter_descriptions={},
     name='PEDS Heatmap',
@@ -357,7 +358,11 @@
     description='A Monte Carlo Simulation that shuffles the recipient and'
                 ' donated microbiome in order to investigate if the true'
                 ' donated microbiome is the most similar to their true FMT'
-                ' recipient',
+                ' recipient. This method is not investigating engraftment'
+                ' extent, but is instead identifying if the recipients'
+                ' engrafted unique features from their respective donated'
+                ' microbiome as opposed to features common to the gut'
+                ' microbiome',
     citations=[citations['aggarwala_precise_2021']],
     examples={
         'peds_methods': ex.simulation_peds_method

From 858ab2023a84fbe4306ff74e91e1fd0659b803e6 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Sun, 14 Jul 2024 08:22:54 -0700
Subject: [PATCH 22/50] viz messing

---
 q2_fmt/_peds.py          | 16 ++++++++--------
 q2_fmt/assets/index.html |  2 ++
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 945a3bc..c6ef41d 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -72,12 +72,12 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
 def peds_heatmap(output_dir: str, data: pd.DataFrame,
                  level_delimiter: str = None, 
                  per_subject_stats: pd.DataFrame = None,
-                 global_stats):
+                 global_stats: pd.DataFrame = None):
     _rename_features(data=data, level_delimiter=level_delimiter)
-    if  per_subject_stats is not None:
-        table1, psstats = _make_stats(per_subject_stats)
-    if  global_stats is not None:
-        table2, gstats = _make_stats(global_stats)
+    if global_stats is not None:
+        table1, gstats = _make_stats(global_stats)
+    if per_subject_stats is not None:
+        table2, psstats = _make_stats(per_subject_stats)
     J_ENV = jinja2.Environment(
         loader=jinja2.PackageLoader('q2_fmt', 'assets')
     )
@@ -115,7 +115,7 @@ def peds_heatmap(output_dir: str, data: pd.DataFrame,
 
     with open(os.path.join(output_dir, 'index.html'), 'w') as fh:
         spec_string = json.dumps(full_spec)
-        fh.write(index.render(spec=spec_string, 
+        fh.write(index.render(spec=spec_string,
                               per_subject_stats=psstats,
                               global_stats=gstats,
                               table1=table1, table2=table2))
@@ -564,13 +564,13 @@ def _shuffle_donor_associations(recipient, reference_column, donor):
     return metadata
 
 
-def _per_subject_stats(actual_temp, shuffled_donor, replicates): 
+def _per_subject_stats(actual_temp, shuffled_donor, replicates):
 
     # Calculating per-subject p-values
     actual_element_wise = actual_temp.values[:, None]
     disagree_df = shuffled_donor >= actual_element_wise
     disagree_series = disagree_df.sum(axis=1)
-    agree_series=replicates-disagree_series
+    agree_series = replicates-disagree_series
     # adding 1 here because you can mathmatically can get p-value of 0 from a
     # Monte Carlo Simulation
     per_subject_p = (disagree_series + 1)/(replicates+1)
diff --git a/q2_fmt/assets/index.html b/q2_fmt/assets/index.html
index e422cda..a42225f 100644
--- a/q2_fmt/assets/index.html
+++ b/q2_fmt/assets/index.html
@@ -70,6 +70,8 @@
     #stats tr {
       border-top: 1px solid;
     }
+
+  
   </style>
 </head>
 <body>

From c7b564de09f7618d3f703e5d8f884df546005519 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:16:40 -0700
Subject: [PATCH 23/50] unittests

---
 q2_fmt/_examples.py              |  6 +-
 q2_fmt/_peds.py                  | 21 +++----
 q2_fmt/plugin_setup.py           | 12 ++--
 q2_fmt/tests/test_engraftment.py | 98 ++++++++++++++++++++++++++------
 4 files changed, 102 insertions(+), 35 deletions(-)

diff --git a/q2_fmt/_examples.py b/q2_fmt/_examples.py
index 7ac6960..2d11bb7 100644
--- a/q2_fmt/_examples.py
+++ b/q2_fmt/_examples.py
@@ -240,7 +240,7 @@ def simulation_peds_method(use):
     md = use.init_metadata('md', peds_md_factory)
     peds_table = use.init_artifact('peds_table', peds_feature_table_factory)
 
-    peds_stats, = use.action(
+    peds_stats, global_stats = use.action(
         use.UsageAction('fmt', 'peds_simulation'),
         use.UsageInputs(
             table=peds_table,
@@ -250,9 +250,11 @@ def simulation_peds_method(use):
             subject_column='SubjectID'
         ),
         use.UsageOutputNames(
-            stats='stats'
+            per_subject_stats='per_subject_stats',
+            global_stats='global_stats'
         )
 
     )
 
     peds_stats.assert_output_type("StatsTable[Pairwise]")
+    global_stats.assert_output_type("StatsTable[Pairwise]")
diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index c6ef41d..f99ae3a 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -510,23 +510,21 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     donor = metadata_df[metadata_df.index.isin(
         _check_reference_column(metadata=metadata_df,
                                 reference_column=reference_column))]
-    try:
-        donor.index.unique() != 1
-    except AssertionError as e:
+
+    if len(donor.index.unique()) == 1:
         raise AssertionError("There is only one donated microbiome in your"
                              " data. A Monte Carlo simulation shuffles"
                              " your donated microbiome and recipient pairing"
-                             "  and needs more than one donated microbiome"
-                             " to successfully shuffle") from e
+                             " and needs more than one donated microbiome"
+                             " to successfully shuffle")
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
-    try:
-        recipient.index.unique() != 1
-    except AssertionError as e:
+
+    if len(recipient.index.unique()) == 1:
         raise AssertionError("There is only one recipient in your"
                              " data. A Monte Carlo simulation shuffles"
                              " your donated microbiome and recipient pairing"
-                             "  and needs more than one recipient"
-                             " to successfully shuffle") from e
+                             " and needs more than one recipient"
+                             " to successfully shuffle")
 
     shuffled_donor = pd.DataFrame([])
     peds = sample_peds(
@@ -539,6 +537,9 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
            drop_incomplete_timepoint=drop_incomplete_timepoint)
     actual_temp = peds["measure"]
     for i in range(replicates):
+        print(recipient)
+        print(reference_column)
+        print(donor)
         metadata = _shuffle_donor_associations(recipient, reference_column,
                                                donor)
         peds = sample_peds(
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index a311131..11bb1fe 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -215,8 +215,8 @@
 plugin.visualizers.register_function(
     function=q2_fmt.peds_heatmap,
     inputs={'data': Dist1D[Ordered, Matched] % Properties("peds"),
-            'per-subject-stats': StatsTable[Pairwise],
-            'global-stats': StatsTable[Pairwise]},
+            'per_subject_stats': StatsTable[Pairwise],
+            'global_stats': StatsTable[Pairwise]},
     parameters={'level_delimiter': Str},
     parameter_descriptions={},
     name='PEDS Heatmap',
@@ -321,8 +321,8 @@
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str],
                 'replicates': Int % Range(999, None)},
-    outputs=[('per-subject-stats', StatsTable[Pairwise]),
-             ('global-stats', StatsTable[Pairwise])],
+    outputs=[('per_subject_stats', StatsTable[Pairwise]),
+             ('global_stats', StatsTable[Pairwise])],
     parameter_descriptions={
         'metadata': 'The sample metadata.',
         'time_column': 'The column within the `metadata` that the'
@@ -349,9 +349,9 @@
                       ' simulation on'
     },
     output_descriptions={
-        'per-subject-stats': 'stats table for comparing if the actual PEDS'
+        'per_subject_stats': 'stats table for comparing if the actual PEDS'
         ' values to shuffled simulated values per subject',
-        'global-stats': 'stats table for comparing if the actual PEDS values'
+        'global_stats': 'stats table for comparing if the actual PEDS values'
         ' to Shuffled simulated values'
     },
     name='',
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 7b7f8a4..679f99f 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -20,7 +20,7 @@
                           _check_subject_column, _check_column_type,
                           _drop_incomplete_timepoints, feature_peds,
                           _check_column_missing, _rename_features,
-                          peds_simulation)
+                          peds_simulation, _shuffle_donor_associations)
 
 
 class TestBase(TestPluginBase):
@@ -1282,14 +1282,14 @@ def test_high_donor_overlap(self):
             'Feature3': [0, 0, 1, 0, 0, 1]}).set_index('id')
         metadata = Metadata(metadata_df)
 
-        stats = peds_simulation(metadata=metadata,
-                                table=table_df,
-                                time_column="group",
-                                reference_column="Ref",
-                                subject_column="subject",
-                                replicates=999)
-        real_median = stats["A:measure"].values
-        fake_median = stats["B:measure"].values
+        stats, _ = peds_simulation(metadata=metadata,
+                                   table=table_df,
+                                   time_column="group",
+                                   reference_column="Ref",
+                                   subject_column="subject",
+                                   replicates=999)
+        real_median = np.median(stats["A:measure"].values)
+        fake_median = np.median(stats["B:measure"].values)
         self.assertGreater(real_median, fake_median)
 
     def test_low_donor_overlap(self):
@@ -1314,13 +1314,77 @@ def test_low_donor_overlap(self):
             'Feature3': [0, 0, 1, 1, 1, 0]}).set_index('id')
         metadata = Metadata(metadata_df)
 
-        stats = peds_simulation(metadata=metadata,
-                                table=table_df,
-                                time_column="group",
-                                reference_column="Ref",
-                                subject_column="subject",
-                                replicates=999)
+        stats, _ = peds_simulation(metadata=metadata,
+                                   table=table_df,
+                                   time_column="group",
+                                   reference_column="Ref",
+                                   subject_column="subject",
+                                   replicates=999)
 
-        real_median = stats["A:measure"].values
-        fake_median = stats["B:measure"].values
+        real_median = np.median(stats["A:measure"].values)
+        fake_median = np.median(stats["B:measure"].values)
         self.assertGreater(fake_median, real_median)
+
+    def test_single_donor(self):
+        metadata_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1'],
+            'Ref': ['donor1', 'donor1', 'donor1', np.nan],
+            'subject': ['sub1', 'sub2', 'sub3', np.nan],
+            'group': [1, 1, 1, np.nan],
+            "Location": [np.nan, np.nan,
+                         np.nan, 'test']}).set_index('id')
+
+        table_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1'],
+            'Feature1': [1, 0, 0, 0],
+            'Feature2': [0, 1, 0, 1],
+            'Feature3': [0, 0, 1, 1]}).set_index('id')
+        metadata = Metadata(metadata_df)
+
+        with self.assertRaisesRegex(AssertionError, "There is only one"
+                                    " donated microbiome in your data. *"):
+            peds_simulation(metadata=metadata,
+                            table=table_df,
+                            time_column="group",
+                            reference_column="Ref",
+                            subject_column="subject",
+                            replicates=999)
+
+    def test_is_shuffled(self):
+        metadata_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan,
+                        np.nan],
+            'group': [1, 1, 1, np.nan, np.nan,
+                      np.nan],
+            "Location": [np.nan, np.nan,
+                         np.nan, 'test', 'test',
+                         'test']}).set_index('id')
+
+        reference = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Ref': ['donor1', 'donor2', 'donor3'],
+            'subject': ['sub1', 'sub2', 'sub3'],
+            'group': [1, 1, 1],
+            "Location": [np.nan, np.nan,
+                         np.nan]}).set_index('id')
+        reference_column = "Ref"
+        donor = pd.DataFrame({
+            'id': ['donor1', 'donor2', 'donor3'],
+            'Ref': [np.nan, np.nan, np.nan],
+            'subject': [np.nan, np.nan, np.nan],
+            'group': [np.nan, np.nan, np.nan],
+            "Location": ['test', 'test',
+                         'test']}).set_index('id')
+
+        shuffled_metadata = \
+            _shuffle_donor_associations(recipient=reference,
+                                        reference_column=reference_column,
+                                        donor=donor).to_dataframe()
+        equals = metadata_df.equals(shuffled_metadata)
+        self.assertFalse(equals)

From db0ea7f64088add4c6b6e6e9a87e8a653e6ee335 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:18:19 -0700
Subject: [PATCH 24/50] linty

---
 q2_fmt/_peds.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index f99ae3a..8e441b6 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -70,7 +70,7 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
 
 
 def peds_heatmap(output_dir: str, data: pd.DataFrame,
-                 level_delimiter: str = None, 
+                 level_delimiter: str = None,
                  per_subject_stats: pd.DataFrame = None,
                  global_stats: pd.DataFrame = None):
     _rename_features(data=data, level_delimiter=level_delimiter)

From 572f0c86fb32e9f6131ac8036f5a2a6b363d023b Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Mon, 15 Jul 2024 14:53:59 -0700
Subject: [PATCH 25/50] viz

---
 q2_fmt/_peds.py          | 34 +++++++++++++---------------
 q2_fmt/assets/index.html | 49 +++++++++++++++++++++++++++++++++-------
 2 files changed, 57 insertions(+), 26 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 8e441b6..be7cc3c 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -75,9 +75,10 @@ def peds_heatmap(output_dir: str, data: pd.DataFrame,
                  global_stats: pd.DataFrame = None):
     _rename_features(data=data, level_delimiter=level_delimiter)
     if global_stats is not None:
-        table1, gstats = _make_stats(global_stats)
+        gstats = global_stats.to_html(index=False)
+        #table2, gstats = _make_stats(global_stats)
     if per_subject_stats is not None:
-        table2, psstats = _make_stats(per_subject_stats)
+        table1, psstats = _make_stats(per_subject_stats)
     J_ENV = jinja2.Environment(
         loader=jinja2.PackageLoader('q2_fmt', 'assets')
     )
@@ -116,9 +117,9 @@ def peds_heatmap(output_dir: str, data: pd.DataFrame,
     with open(os.path.join(output_dir, 'index.html'), 'w') as fh:
         spec_string = json.dumps(full_spec)
         fh.write(index.render(spec=spec_string,
-                              per_subject_stats=psstats,
-                              global_stats=gstats,
-                              table1=table1, table2=table2))
+                              persubjectstats=psstats,
+                              globalstats=gstats,
+                              table1=table1))
 
 
 def sample_peds(table: pd.DataFrame, metadata: qiime2.Metadata,
@@ -537,9 +538,6 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
            drop_incomplete_timepoint=drop_incomplete_timepoint)
     actual_temp = peds["measure"]
     for i in range(replicates):
-        print(recipient)
-        print(reference_column)
-        print(donor)
         metadata = _shuffle_donor_associations(recipient, reference_column,
                                                donor)
         peds = sample_peds(
@@ -578,7 +576,7 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
     # Common Langauge effect size calcs in prep for stats refactor.
     agree = 1 - per_subject_p
-    rank_biserial = agree - per_subject_p
+    #rank_biserial = agree - per_subject_p
 
     per_sub_stats = pd.DataFrame({'A:group': "actual_values",
                                  'A:n': actual_temp.size,
@@ -589,9 +587,9 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
                                   'n': replicates*actual_temp.size,
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
-                                  'q-value': per_subject_q,
-                                  'CLES': agree,
-                                  'rank-biserial-correlation': rank_biserial})
+                                  'q-value': per_subject_q})
+                                  #'CLES': agree,
+                                  #'rank-biserial-correlation': rank_biserial})
     per_sub_stats['A:group'].attrs.update({'title': 'actual_values',
                                           'description': 'PEDS values'
                                            ' calculated with actual recipient'
@@ -622,12 +620,12 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
-    per_sub_stats['CLES'].attrs.update(
-        dict(title='Common Language Effect Size',
-             description='Common Language Effect Size'))
-    per_sub_stats['rank-biserial-correlation'].attrs.update(
-        dict(title='rank biserial correlation',
-             description='rank biserial correlation ranging from -1 to 1'))
+    #per_sub_stats['CLES'].attrs.update(
+     #   dict(title='Common Language Effect Size',
+      #       description='Common Language Effect Size'))
+    #per_sub_stats['rank-biserial-correlation'].attrs.update(
+    #    dict(title='rank biserial correlation',
+    #         description='rank biserial correlation ranging from -1 to 1'))
     return per_sub_stats
 
 
diff --git a/q2_fmt/assets/index.html b/q2_fmt/assets/index.html
index a42225f..fffef1b 100644
--- a/q2_fmt/assets/index.html
+++ b/q2_fmt/assets/index.html
@@ -44,7 +44,7 @@
       right: 212px !important;
       top: 85px !important;
     }
-    #stats table {
+    #persubjectstats table {
       margin-top: 50px;
       width: 100%;
       text-align: right;
@@ -53,21 +53,47 @@
       font-family: monospace;
     }
 
-    #stats thead tr {
+    #persubjectstats thead tr {
       border: none;
     }
 
-    #stats tr th {
+    #persubjectstats tr th {
       border: none;
       padding-right: 10px;
     }
 
-    #stats td {
+    #persubjectstats td {
       border: none;
       padding-right: 10px;
     }
 
-    #stats tr {
+    #persubjectstats tr {
+      border-top: 1px solid;
+    }
+    #globalstats table {
+      margin-top: 50px;
+      width: 100%;
+      text-align: right;
+      border-collapse: collapse;
+      border: 0px;
+      font-family: monospace;
+    }
+
+    #globalstats thead tr {
+      border: none;
+    }
+
+    #globalstats tr th {
+      border: none;
+      padding-right: 10px;
+    }
+
+    #globalstats td {
+      border: none;
+      padding-right: 10px;
+    }
+
+    #globalstats tr {
       border-top: 1px solid;
     }
 
@@ -80,13 +106,20 @@
 <p><strong>sampleSize:</strong> If clicked, the sampleSize box scales your heatmap boxes based on the denominator of the PEDS proportion </p>
 
 
-{% if stats %}
-<div id='stats'>
-{{ stats }}
+{% if persubjectstats %}
+<div id='persubjectstats'>
+{{ persubjectstats }}
 </div>
 <p><strong>Table 1.</strong> {{ table1 }}</p>
 {% endif %}
 
+{% if globalstats %}
+<div id='globalstats'>
+{{ globalstats }}
+</div>
+{% endif %}
+
+
 <script type="text/javascript">
   var spec = JSON.parse(document.getElementById('spec').innerHTML);
   vegaEmbed('#viz', spec, {'renderer': 'svg'}).then(function(result) {

From bed46085f1e8774350b7d9e55adb31c00040e5dc Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Mon, 15 Jul 2024 15:10:25 -0700
Subject: [PATCH 26/50] different make stats

---
 q2_fmt/_peds.py | 58 +++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 46 insertions(+), 12 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index be7cc3c..ca1a85e 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -18,7 +18,7 @@
 import json
 
 from q2_fmt._util import json_replace
-from q2_stats._visualizer import _make_stats
+from q2_stats._visualizer import _make_group_col
 
 
 def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
@@ -76,7 +76,7 @@ def peds_heatmap(output_dir: str, data: pd.DataFrame,
     _rename_features(data=data, level_delimiter=level_delimiter)
     if global_stats is not None:
         gstats = global_stats.to_html(index=False)
-        #table2, gstats = _make_stats(global_stats)
+        # table2, gstats = _make_stats(global_stats)
     if per_subject_stats is not None:
         table1, psstats = _make_stats(per_subject_stats)
     J_ENV = jinja2.Environment(
@@ -576,7 +576,7 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
     # Common Langauge effect size calcs in prep for stats refactor.
     agree = 1 - per_subject_p
-    #rank_biserial = agree - per_subject_p
+    rank_biserial = agree - per_subject_p
 
     per_sub_stats = pd.DataFrame({'A:group': "actual_values",
                                  'A:n': actual_temp.size,
@@ -587,9 +587,9 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
                                   'n': replicates*actual_temp.size,
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
-                                  'q-value': per_subject_q})
-                                  #'CLES': agree,
-                                  #'rank-biserial-correlation': rank_biserial})
+                                  'q-value': per_subject_q,
+                                  'CLES': agree,
+                                  'rank-biserial-correlation': rank_biserial})
     per_sub_stats['A:group'].attrs.update({'title': 'actual_values',
                                           'description': 'PEDS values'
                                            ' calculated with actual recipient'
@@ -620,12 +620,12 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
-    #per_sub_stats['CLES'].attrs.update(
-     #   dict(title='Common Language Effect Size',
-      #       description='Common Language Effect Size'))
-    #per_sub_stats['rank-biserial-correlation'].attrs.update(
-    #    dict(title='rank biserial correlation',
-    #         description='rank biserial correlation ranging from -1 to 1'))
+    per_sub_stats['CLES'].attrs.update(
+        dict(title='Common Language Effect Size',
+             description='Common Language Effect Size'))
+    per_sub_stats['rank-biserial-correlation'].attrs.update(
+        dict(title='rank biserial correlation',
+             description='rank biserial correlation ranging from -1 to 1'))
     return per_sub_stats
 
 
@@ -646,3 +646,37 @@ def _global_stats(p_series):
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
     return global_stats
+
+
+def _make_stats(stats):
+    method = stats['test-statistic'].attrs['title']
+    group_unit = (stats['A:group'].attrs['title']
+                  + ' vs ' + stats['B:group'].attrs['title'])
+    pval_method = stats['p-value'].attrs['title']
+    qval_method = stats['q-value'].attrs['title']
+    table1 = (f'{method} tests between groups ({group_unit}), with'
+              f' {pval_method} p-value calculations and {qval_method}'
+              f' correction for multiple comparisons (q-value).')
+    df = pd.DataFrame(index=stats.index)
+    group_a = _make_group_col('A', stats)
+    df[group_a.name] = group_a
+    group_b = _make_group_col('B', stats)
+    df[group_b.name] = group_b
+    df['A'] = stats['A:measure']
+    df['B'] = stats['B:measure']
+    df = df.merge(stats.iloc[:, 6:], left_index=True, right_index=True)
+    df.columns = pd.MultiIndex.from_tuples([
+        ('Group A', stats['A:group'].attrs['title']),
+        ('Group B', stats['B:group'].attrs['title']),
+        ('A', stats['A:measure'].attrs['title']),
+        ('B', stats['B:measure'].attrs['title']),
+        ('', 'n'),
+        ('', 'test-statistic'),
+        ('', 'p-value'),
+        ('', 'q-value'),
+        ('', 'CLES'),
+        ('', 'rank-biserial-correlation'),
+    ])
+    html = df.to_html(index=False)
+
+    return table1, html

From d3621f29107140a43660760b481a41d1969a1c5e Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 16 Jul 2024 15:37:42 -0700
Subject: [PATCH 27/50] evans stats request

---
 q2_fmt/_peds.py | 35 +++++++++++++++++++++--------------
 1 file changed, 21 insertions(+), 14 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index ca1a85e..dece1bc 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -9,6 +9,7 @@
 import qiime2
 import pandas as pd
 import numpy as np
+import random
 import warnings
 from scipy.stats import false_discovery_control, combine_pvalues
 from collections import Counter
@@ -535,7 +536,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
            subject_column=subject_column,
            filter_missing_references=filter_missing_references,
            drop_incomplete_subjects=drop_incomplete_subjects,
-           drop_incomplete_timepoint=drop_incomplete_timepoint)
+           drop_incomplete_timepoint=drop_incomplete_timepoint).set_index("id")
     actual_temp = peds["measure"]
     for i in range(replicates):
         metadata = _shuffle_donor_associations(recipient, reference_column,
@@ -547,16 +548,19 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                subject_column=subject_column,
                filter_missing_references=filter_missing_references,
                drop_incomplete_subjects=drop_incomplete_subjects,
-               drop_incomplete_timepoint=drop_incomplete_timepoint)
+               drop_incomplete_timepoint=drop_incomplete_timepoint
+               ).set_index("id")
         shuffled_donor[i] = peds["measure"]
-        per_sub_stats = _per_subject_stats(actual_temp, shuffled_donor,
-                                           replicates)
-        global_stats = _global_stats(per_sub_stats['p-value'])
+    per_sub_stats = _per_subject_stats(actual_temp, shuffled_donor,
+                                       replicates)
+    global_stats = _global_stats(per_sub_stats['p-value'])
     return per_sub_stats, global_stats
 
 
 def _shuffle_donor_associations(recipient, reference_column, donor):
-    shuffled_list = recipient[reference_column].sample(frac=1).to_list()
+    donors = recipient[reference_column].unique()
+    shuffled_list = random.choices(donors, k=len(recipient))
+    # shuffled_list = recipient[reference_column].sample(frac=1).to_list()
     recipient.loc[:, reference_column] = shuffled_list
     metadata_df = pd.concat([donor, recipient])
     metadata = qiime2.Metadata(metadata_df)
@@ -577,14 +581,15 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     # Common Langauge effect size calcs in prep for stats refactor.
     agree = 1 - per_subject_p
     rank_biserial = agree - per_subject_p
+    shuffled_names = "shuffled_" + shuffled_donor.reset_index()["id"]
 
-    per_sub_stats = pd.DataFrame({'A:group': "actual_values",
-                                 'A:n': actual_temp.size,
+    per_sub_stats = pd.DataFrame({'A:group': actual_temp.index,
+                                 'A:n': 1,
                                   'A:measure': actual_temp,
-                                  'B:group': "shuffled_values",
+                                  'B:group': shuffled_names.values,
                                   'B:n': replicates,
-                                  'B:measure': shuffled_donor.median(axis=1),
-                                  'n': replicates*actual_temp.size,
+                                  'B:measure': shuffled_donor.mean(axis=1),
+                                  'n': replicates,
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
                                   'q-value': per_subject_q,
@@ -604,10 +609,12 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_sub_stats['A:n'].attrs.update(n)
     per_sub_stats['B:n'].attrs.update(n)
     measure = {
-        'title': 'Median PEDS Value',
-        'description': 'Median PEDS Value'
+        'title': 'Mean PEDS Value',
+        'description': 'Mean PEDS Value'
     }
-    per_sub_stats['A:measure'].attrs.update(measure)
+    per_sub_stats['A:measure'].attrs.update({'title': 'PEDS Value',
+                                            'description': 'PEDS Value'
+                                             })
     per_sub_stats['B:measure'].attrs.update(measure)
     per_sub_stats['n'].attrs.update(dict(title='count', description='Number of'
                                     'comparisons'))

From d2def81649fb66febc49c9a307102607cf66d00f Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 17 Jul 2024 08:01:20 -0700
Subject: [PATCH 28/50] leftovers

---
 q2_fmt/_peds.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index dece1bc..4276550 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -559,6 +559,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
 
 def _shuffle_donor_associations(recipient, reference_column, donor):
     donors = recipient[reference_column].unique()
+    # Randomize donors using a uniform distribution
     shuffled_list = random.choices(donors, k=len(recipient))
     # shuffled_list = recipient[reference_column].sample(frac=1).to_list()
     recipient.loc[:, reference_column] = shuffled_list

From 77f0f03e463f65a956e5cfb987f35dbba5d9b967 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 18 Jul 2024 17:41:55 -0700
Subject: [PATCH 29/50] REF: recipient shuffle

---
 q2_fmt/_peds.py                  | 107 ++++++++++++++++++++-----------
 q2_fmt/tests/test_engraftment.py |  39 +----------
 2 files changed, 71 insertions(+), 75 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 4276550..00d9459 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -10,6 +10,7 @@
 import pandas as pd
 import numpy as np
 import random
+import itertools
 import warnings
 from scipy.stats import false_discovery_control, combine_pvalues
 from collections import Counter
@@ -366,7 +367,6 @@ def _check_reference_column(metadata, reference_column):
 
 def _filter_associated_reference(reference_series, metadata, time_column,
                                  filter_missing_references, reference_column):
-
     used_references = reference_series[~metadata[time_column].isna()]
 
     if used_references.isna().any():
@@ -490,14 +490,12 @@ def _create_masking(time_metadata, donor_df, recip_df, reference_column):
     donor_df = donor_df.to_numpy()
     donor_mask = donor_df[donor_index_masking]
     donor_mask = donor_mask.astype(int)
-    recipdf = recip_df.to_numpy()
-    recipdf = recipdf.astype(int)
     return donor_mask
 
 
 def _mask_recipient(donor_mask, recip_df):
-    maskedrecip = donor_mask & recip_df
-    return maskedrecip
+    recip_mask = donor_mask & recip_df
+    return recip_mask
 
 
 def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
@@ -509,11 +507,16 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                     replicates: int = 999) -> (pd.DataFrame, pd.DataFrame):
 
     metadata_df = metadata.to_dataframe()
-    donor = metadata_df[metadata_df.index.isin(
-        _check_reference_column(metadata=metadata_df,
-                                reference_column=reference_column))]
+    reference_series = _check_reference_column(metadata_df,
+                                               reference_column)
+
+    (metadata_df,
+     used_references) = _filter_associated_reference(reference_series,
+                                                     metadata_df, time_column,
+                                                     filter_missing_references,
+                                                     reference_column)
 
-    if len(donor.index.unique()) == 1:
+    if len(used_references.unique()) == 1:
         raise AssertionError("There is only one donated microbiome in your"
                              " data. A Monte Carlo simulation shuffles"
                              " your donated microbiome and recipient pairing"
@@ -528,7 +531,6 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                              " and needs more than one recipient"
                              " to successfully shuffle")
 
-    shuffled_donor = pd.DataFrame([])
     peds = sample_peds(
            table=table, metadata=metadata,
            time_column=time_column,
@@ -538,34 +540,67 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
            drop_incomplete_subjects=drop_incomplete_subjects,
            drop_incomplete_timepoint=drop_incomplete_timepoint).set_index("id")
     actual_temp = peds["measure"]
-    for i in range(replicates):
-        metadata = _shuffle_donor_associations(recipient, reference_column,
-                                               donor)
-        peds = sample_peds(
-               table=table, metadata=metadata,
-               time_column=time_column,
-               reference_column=reference_column,
-               subject_column=subject_column,
-               filter_missing_references=filter_missing_references,
-               drop_incomplete_subjects=drop_incomplete_subjects,
-               drop_incomplete_timepoint=drop_incomplete_timepoint
-               ).set_index("id")
-        shuffled_donor[i] = peds["measure"]
-    per_sub_stats = _per_subject_stats(actual_temp, shuffled_donor,
+
+    # Mismatch simulation:
+    table = table > 0
+    recip_df = _create_recipient_table(used_references, metadata_df, table)
+    donor_df = table[table.index.isin(used_references)]
+    duplicated_table, mismatch_df = \
+        _create_mismatch_pairs(recip_df,
+                               metadata_df,
+                               used_references,
+                               reference_column)
+    donor_mask = _create_sim_masking(mismatch_df, donor_df, reference_column)
+    recip_mask = donor_mask & duplicated_table
+    num_sum = np.sum(recip_mask.values, axis=1)
+    donor_sum = np.sum(donor_mask, axis=1)
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        mismatch_peds = num_sum/donor_sum
+    mismatchpairs_df = _simulate_uniform_distro(recip_df, mismatch_peds,
+                                                replicates)
+    per_sub_stats = _per_subject_stats(actual_temp,  mismatchpairs_df,
                                        replicates)
     global_stats = _global_stats(per_sub_stats['p-value'])
     return per_sub_stats, global_stats
 
 
-def _shuffle_donor_associations(recipient, reference_column, donor):
-    donors = recipient[reference_column].unique()
-    # Randomize donors using a uniform distribution
-    shuffled_list = random.choices(donors, k=len(recipient))
-    # shuffled_list = recipient[reference_column].sample(frac=1).to_list()
-    recipient.loc[:, reference_column] = shuffled_list
-    metadata_df = pd.concat([donor, recipient])
-    metadata = qiime2.Metadata(metadata_df)
-    return metadata
+# Simulation helper
+def _create_mismatch_pairs(recip_df, metadata, used_references,
+                           reference_column):
+    mismatch_pairs = []
+    for x in itertools.product(recip_df.index,
+                               metadata[reference_column].dropna()):
+        mismatch_pairs.append(x)
+    mismatched_t = list(zip(used_references.index, used_references))
+    # Removes matched donor and recipient pairs
+    filtered = [tuple for tuple in mismatch_pairs if tuple not in mismatched_t]
+    idx, values = zip(*filtered)
+    mismatched_df = pd.DataFrame({"id": idx,
+                                  reference_column: values}).set_index("id")
+    duplicated_table = \
+        mismatched_df.merge(recip_df, left_index=True, right_index=True,
+                            how='left').drop(mismatched_df.columns, axis=1)
+    return duplicated_table, mismatched_df
+
+
+def _create_sim_masking(mismatch_df, donor_df, reference_column):
+    donors = mismatch_df[reference_column]
+    donor_index_masking = []
+    for donor in donors:
+        donor_index_masking.append(donor_df.index.get_loc(donor))
+    donor_df_np = donor_df.to_numpy()
+    donor_mask = donor_df_np[donor_index_masking]
+    donor_mask = donor_mask.astype(int)
+    return donor_mask
+
+
+def _simulate_uniform_distro(recip_df, mismatch_peds, replicates):
+    mismatchpairs_df = pd.DataFrame([], columns=range(replicates))
+    for x in range(len(recip_df.index)):
+        peds_iters = random.choices(mismatch_peds, k=replicates)
+        mismatchpairs_df.loc[len(mismatchpairs_df)] = peds_iters
+    return mismatchpairs_df
 
 
 def _per_subject_stats(actual_temp, shuffled_donor, replicates):
@@ -582,12 +617,10 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     # Common Langauge effect size calcs in prep for stats refactor.
     agree = 1 - per_subject_p
     rank_biserial = agree - per_subject_p
-    shuffled_names = "shuffled_" + shuffled_donor.reset_index()["id"]
-
     per_sub_stats = pd.DataFrame({'A:group': actual_temp.index,
                                  'A:n': 1,
-                                  'A:measure': actual_temp,
-                                  'B:group': shuffled_names.values,
+                                  'A:measure': actual_temp.values,
+                                  'B:group': "shuffled_recipients",
                                   'B:n': replicates,
                                   'B:measure': shuffled_donor.mean(axis=1),
                                   'n': replicates,
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 679f99f..49cf39e 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -20,7 +20,7 @@
                           _check_subject_column, _check_column_type,
                           _drop_incomplete_timepoints, feature_peds,
                           _check_column_missing, _rename_features,
-                          peds_simulation, _shuffle_donor_associations)
+                          peds_simulation)
 
 
 class TestBase(TestPluginBase):
@@ -1351,40 +1351,3 @@ def test_single_donor(self):
                             reference_column="Ref",
                             subject_column="subject",
                             replicates=999)
-
-    def test_is_shuffled(self):
-        metadata_df = pd.DataFrame({
-            'id': ['sample1', 'sample2', 'sample3',
-                   'donor1', 'donor2', 'donor3'],
-            'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan,
-                    np.nan],
-            'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan,
-                        np.nan],
-            'group': [1, 1, 1, np.nan, np.nan,
-                      np.nan],
-            "Location": [np.nan, np.nan,
-                         np.nan, 'test', 'test',
-                         'test']}).set_index('id')
-
-        reference = pd.DataFrame({
-            'id': ['sample1', 'sample2', 'sample3'],
-            'Ref': ['donor1', 'donor2', 'donor3'],
-            'subject': ['sub1', 'sub2', 'sub3'],
-            'group': [1, 1, 1],
-            "Location": [np.nan, np.nan,
-                         np.nan]}).set_index('id')
-        reference_column = "Ref"
-        donor = pd.DataFrame({
-            'id': ['donor1', 'donor2', 'donor3'],
-            'Ref': [np.nan, np.nan, np.nan],
-            'subject': [np.nan, np.nan, np.nan],
-            'group': [np.nan, np.nan, np.nan],
-            "Location": ['test', 'test',
-                         'test']}).set_index('id')
-
-        shuffled_metadata = \
-            _shuffle_donor_associations(recipient=reference,
-                                        reference_column=reference_column,
-                                        donor=donor).to_dataframe()
-        equals = metadata_df.equals(shuffled_metadata)
-        self.assertFalse(equals)

From e889ee7083940bcfcd3d019fe7c62a38ad05249e Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Fri, 19 Jul 2024 11:42:46 -0700
Subject: [PATCH 30/50] helper-method-tests

---
 q2_fmt/tests/test_engraftment.py | 82 +++++++++++++++++++++++++++++++-
 1 file changed, 81 insertions(+), 1 deletion(-)

diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 49cf39e..53bfa75 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -20,7 +20,8 @@
                           _check_subject_column, _check_column_type,
                           _drop_incomplete_timepoints, feature_peds,
                           _check_column_missing, _rename_features,
-                          peds_simulation)
+                          peds_simulation, _create_mismatch_pairs,
+                          _simulate_uniform_distro, _create_sim_masking)
 
 
 class TestBase(TestPluginBase):
@@ -1351,3 +1352,82 @@ def test_single_donor(self):
                             reference_column="Ref",
                             subject_column="subject",
                             replicates=999)
+
+    def test_mismatch_pairs(self):
+        metadata_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan,
+                    np.nan],
+            'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan,
+                        np.nan],
+            'group': [1, 1, 1, np.nan, np.nan,
+                      np.nan],
+            "Location": [np.nan, np.nan,
+                         np.nan, 'test', 'test',
+                         'test']}).set_index('id')
+
+        recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+        used_references = pd.Series(data=['donor1', 'donor2', 'donor3'],
+                                    index=['sample1', 'sample2', 'sample3'],
+                                    name='Ref')
+        used_references.index.name = "id"
+        _, mismatched_df = _create_mismatch_pairs(recip_df, metadata_df,
+                                                  used_references,
+                                                  reference_column='Ref')
+        exp_mismatched_df = pd.DataFrame({'id': ["sample1", "sample1",
+                                                 "sample2", "sample2",
+                                                 "sample3", "sample3"],
+                                          "Ref": ["donor2", "donor3",
+                                                  "donor1", "donor3",
+                                                  "donor1", "donor2"]}
+                                         ).set_index('id')
+        pd.testing.assert_frame_equal(mismatched_df, exp_mismatched_df)
+
+    def test_uniform_distro(self):
+        recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+        mismatch_peds = [0, 0, 0, 0, 0, 0]
+
+        replicates = 999
+
+        mismatchpairs_df = _simulate_uniform_distro(recip_df, mismatch_peds,
+                                                    999)
+        self.assertEquals(mismatchpairs_df.size,
+                          (recip_df.index.size * replicates))
+
+    def test_sim_masking(self):
+
+        mismatched_df = pd.DataFrame({'id': ["sample1", "sample1",
+                                             "sample2", "sample2",
+                                             "sample3", "sample3"],
+                                      "Ref": ["donor2", "donor3",
+                                              "donor1", "donor3",
+                                              "donor1", "donor2"]}
+                                     ).set_index('id')
+
+        donor_df = pd.DataFrame({
+            'id': ['donor1', 'donor2', 'donor3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+        exp_mask = [[0, 1, 0],
+                    [0, 0, 1],
+                    [1, 0, 0],
+                    [0, 0, 1],
+                    [1, 0, 0],
+                    [0, 1, 0]]
+
+        donor_mask = _create_sim_masking(mismatched_df, donor_df,
+                                         reference_column='Ref')
+        np.testing.assert_array_equal(donor_mask, exp_mask)

From 6bdef3772645ebeec1620bad3ebdc2a301d3e08b Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Fri, 19 Jul 2024 11:58:51 -0700
Subject: [PATCH 31/50] stats-bumper

---
 q2_fmt/_peds.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 00d9459..c5b83d8 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -691,8 +691,14 @@ def _global_stats(p_series):
 
 def _make_stats(stats):
     method = stats['test-statistic'].attrs['title']
-    group_unit = (stats['A:group'].attrs['title']
-                  + ' vs ' + stats['B:group'].attrs['title'])
+    try:
+        group_unit = (stats['A:group'].attrs['title']
+                      + ' vs ' + stats['B:group'].attrs['title'])
+    except KeyError as e:
+        raise AssertionError("The wrong stats table was handed to"
+                             " per-subject-stats. Check to make sure you are"
+                             " passing in the correct stats tables") from e
+
     pval_method = stats['p-value'].attrs['title']
     qval_method = stats['q-value'].attrs['title']
     table1 = (f'{method} tests between groups ({group_unit}), with'

From 1c388ccd03af7366bd4bf933dd7c0b0857a793e4 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 23 Jul 2024 10:19:50 -0700
Subject: [PATCH 32/50] make stats

---
 q2_fmt/_peds.py        | 55 ++----------------------------------------
 q2_fmt/plugin_setup.py | 10 ++++++--
 2 files changed, 10 insertions(+), 55 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index c5b83d8..00ca4c3 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -20,7 +20,7 @@
 import json
 
 from q2_fmt._util import json_replace
-from q2_stats._visualizer import _make_group_col
+from q2_stats._visualizer import _make_stats
 
 
 def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
@@ -614,9 +614,6 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     # Monte Carlo Simulation
     per_subject_p = (disagree_series + 1)/(replicates+1)
     per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
-    # Common Langauge effect size calcs in prep for stats refactor.
-    agree = 1 - per_subject_p
-    rank_biserial = agree - per_subject_p
     per_sub_stats = pd.DataFrame({'A:group': actual_temp.index,
                                  'A:n': 1,
                                   'A:measure': actual_temp.values,
@@ -626,9 +623,7 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
                                   'n': replicates,
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
-                                  'q-value': per_subject_q,
-                                  'CLES': agree,
-                                  'rank-biserial-correlation': rank_biserial})
+                                  'q-value': per_subject_q})
     per_sub_stats['A:group'].attrs.update({'title': 'actual_values',
                                           'description': 'PEDS values'
                                            ' calculated with actual recipient'
@@ -661,12 +656,6 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
     per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
-    per_sub_stats['CLES'].attrs.update(
-        dict(title='Common Language Effect Size',
-             description='Common Language Effect Size'))
-    per_sub_stats['rank-biserial-correlation'].attrs.update(
-        dict(title='rank biserial correlation',
-             description='rank biserial correlation ranging from -1 to 1'))
     return per_sub_stats
 
 
@@ -687,43 +676,3 @@ def _global_stats(p_series):
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
              'Benjamini–Hochberg procedure'))
     return global_stats
-
-
-def _make_stats(stats):
-    method = stats['test-statistic'].attrs['title']
-    try:
-        group_unit = (stats['A:group'].attrs['title']
-                      + ' vs ' + stats['B:group'].attrs['title'])
-    except KeyError as e:
-        raise AssertionError("The wrong stats table was handed to"
-                             " per-subject-stats. Check to make sure you are"
-                             " passing in the correct stats tables") from e
-
-    pval_method = stats['p-value'].attrs['title']
-    qval_method = stats['q-value'].attrs['title']
-    table1 = (f'{method} tests between groups ({group_unit}), with'
-              f' {pval_method} p-value calculations and {qval_method}'
-              f' correction for multiple comparisons (q-value).')
-    df = pd.DataFrame(index=stats.index)
-    group_a = _make_group_col('A', stats)
-    df[group_a.name] = group_a
-    group_b = _make_group_col('B', stats)
-    df[group_b.name] = group_b
-    df['A'] = stats['A:measure']
-    df['B'] = stats['B:measure']
-    df = df.merge(stats.iloc[:, 6:], left_index=True, right_index=True)
-    df.columns = pd.MultiIndex.from_tuples([
-        ('Group A', stats['A:group'].attrs['title']),
-        ('Group B', stats['B:group'].attrs['title']),
-        ('A', stats['A:measure'].attrs['title']),
-        ('B', stats['B:measure'].attrs['title']),
-        ('', 'n'),
-        ('', 'test-statistic'),
-        ('', 'p-value'),
-        ('', 'q-value'),
-        ('', 'CLES'),
-        ('', 'rank-biserial-correlation'),
-    ])
-    html = df.to_html(index=False)
-
-    return table1, html
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 95a9777..f30064d 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -260,7 +260,13 @@
     inputs={'data': Dist1D[Ordered, Matched] % Properties("peds"),
             'per_subject_stats': StatsTable[Pairwise],
             'global_stats': StatsTable[Pairwise]},
-    input_descriptions={'data': 'PEDS output to plot'},
+    input_descriptions={'data': 'PEDS output to plot',
+                        'per_subject_stats': 'stats table for comparing if the'
+                                             ' actual PEDS values to shuffled'
+                                             ' simulated values per subject',
+                        'global_stats': 'stats table for comparing if the'
+                                        ' actual PEDS values to shuffled'
+                                        ' simulated values globally.'},
     parameters={'level_delimiter': Str},
     parameter_descriptions={
                             'level_delimiter': 'delimiter to split taxonomic'
@@ -403,7 +409,7 @@
         'per_subject_stats': 'stats table for comparing if the actual PEDS'
         ' values to shuffled simulated values per subject',
         'global_stats': 'stats table for comparing if the actual PEDS values'
-        ' to Shuffled simulated values'
+        ' to shuffled simulated values globally.'
     },
     name='',
     description='A Monte Carlo Simulation that shuffles the recipient and'

From abce83252c4f4c32923be35d6db149a572981437 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 25 Jul 2024 13:58:13 -0700
Subject: [PATCH 33/50] Apply suggestions from code review

Co-authored-by: Greg Caporaso <192372+gregcaporaso@users.noreply.github.com>
---
 q2_fmt/_peds.py        |  6 +++---
 q2_fmt/plugin_setup.py | 25 +++++++++++++------------
 2 files changed, 16 insertions(+), 15 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 00ca4c3..25c7742 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -521,7 +521,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                              " data. A Monte Carlo simulation shuffles"
                              " your donated microbiome and recipient pairing"
                              " and needs more than one donated microbiome"
-                             " to successfully shuffle")
+                             " to successfully shuffle.")
     recipient = metadata_df.loc[metadata_df[reference_column].notnull()]
 
     if len(recipient.index.unique()) == 1:
@@ -529,7 +529,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                              " data. A Monte Carlo simulation shuffles"
                              " your donated microbiome and recipient pairing"
                              " and needs more than one recipient"
-                             " to successfully shuffle")
+                             " to successfully shuffle.")
 
     peds = sample_peds(
            table=table, metadata=metadata,
@@ -565,7 +565,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     return per_sub_stats, global_stats
 
 
-# Simulation helper
+# peds_simulation helper functions
 def _create_mismatch_pairs(recip_df, metadata, used_references,
                            reference_column):
     mismatch_pairs = []
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index f30064d..e7d773b 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -372,8 +372,10 @@
     function=q2_fmt.peds_simulation,
     inputs={'table': FeatureTable[Frequency | RelativeFrequency |
                                   PresenceAbsence]},
-    parameters={'metadata': Metadata, 'time_column': Str,
-                'reference_column': Str, 'subject_column': T_subject,
+    parameters={'metadata': Metadata, 
+                'time_column': Str,
+                'reference_column': Str, 
+                'subject_column': T_subject,
                 'filter_missing_references': Bool,
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str],
@@ -382,20 +384,19 @@
              ('global_stats', StatsTable[Pairwise])],
     parameter_descriptions={
         'metadata': 'The sample metadata.',
-        'time_column': 'The column within the `metadata` that the'
+        'time_column': 'The column within `metadata` that the'
                        ' `table` should be grouped by. This column'
-                       ' should contain simple integer values.',
-        'reference_column': 'The column within the `metadata` that contains'
-                            ' the sample to use as a reference'
-                            ' for a given `table`.'
-                            ' For example, this may be the relevant donor'
-                            ' sample to compare against.',
+                       ' should contain integer values.',
+        'reference_column': 'The column in `metadata` indicating'
+                            ' the sample to use as the reference.'
+                            ' For example, this may list the relevant donor'
+                            ' sample for each recipient sample.',
         'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID to be tracked against timepoints.',
+                          ' subject ID for each sample.',
         'filter_missing_references': 'Filter out references contained within'
                                      ' the metadata that are not present'
-                                     ' in the table.'
-                                     ' Default behavior is to raise an error.',
+                                     ' in the input feature table.'
+                                     ' Default behavior is to raise an error if some reference ids are in the metadata but not the feature table.',
         'drop_incomplete_subjects': 'Filter out subjects that do not have'
                                     ' a sample at every timepoint.'
                                     ' Default behavior is to raise an error.',

From d42f10756115059e17c89ba942d4c208c4033591 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 25 Jul 2024 15:50:36 -0700
Subject: [PATCH 34/50] addresses Greg's Review

---
 q2_fmt/_peds.py                  | 88 +++++++++++++++++---------------
 q2_fmt/plugin_setup.py           | 52 +++++++++++--------
 q2_fmt/tests/test_engraftment.py | 18 +++----
 3 files changed, 85 insertions(+), 73 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 25c7742..4d805db 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -34,7 +34,6 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
 
     if peds_metric == 'sample':
         sample_peds = ctx.get_action('fmt', 'sample_peds')
-        peds_bootstrap = ctx.get_action('fmt', 'peds_bootstrap')
         peds_dist = sample_peds(
             table=table, metadata=metadata, time_column=time_column,
             subject_column=subject_column, reference_column=reference_column,
@@ -42,15 +41,6 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
             drop_incomplete_timepoint=drop_incomplete_timepoint,
             filter_missing_references=filter_missing_references)
 
-        peds_stats = peds_bootstrap(
-            table=table, metadata=metadata, time_column=time_column,
-            subject_column=subject_column,
-            reference_column=reference_column,
-            drop_incomplete_subjects=drop_incomplete_subjects,
-            drop_incomplete_timepoint=drop_incomplete_timepoint,
-            filter_missing_references=filter_missing_references)
-        stat = peds_stats[0]
-
     else:
         if drop_incomplete_subjects or drop_incomplete_timepoint:
             warnings.warn('Feature PEDS was selected as the PEDS metric, which'
@@ -64,8 +54,7 @@ def peds(ctx, table, metadata, peds_metric, time_column, reference_column,
             table=table, metadata=metadata, time_column=time_column,
             subject_column=subject_column, reference_column=reference_column,
             filter_missing_references=filter_missing_references)
-        stat = None
-    results += peds_heatmap(data=peds_dist[0], stats=stat,
+    results += peds_heatmap(data=peds_dist[0],
                             level_delimiter=level_delimiter)
 
     return tuple(results)
@@ -504,7 +493,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                     filter_missing_references: bool = False,
                     drop_incomplete_subjects: bool = False,
                     drop_incomplete_timepoint: list = None,
-                    replicates: int = 999) -> (pd.DataFrame, pd.DataFrame):
+                    iterations: int = 999) -> (pd.DataFrame, pd.DataFrame):
 
     metadata_df = metadata.to_dataframe()
     reference_series = _check_reference_column(metadata_df,
@@ -545,36 +534,45 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     table = table > 0
     recip_df = _create_recipient_table(used_references, metadata_df, table)
     donor_df = table[table.index.isin(used_references)]
-    duplicated_table, mismatch_df = \
-        _create_mismatch_pairs(recip_df,
-                               metadata_df,
-                               used_references,
-                               reference_column)
-    donor_mask = _create_sim_masking(mismatch_df, donor_df, reference_column)
+    duplicated_table, mismatched_df = \
+        _create_mismatched_pairs(recip_df,
+                                 metadata_df,
+                                 used_references,
+                                 reference_column)
+    donor_mask = _create_sim_masking(mismatched_df, donor_df, reference_column)
     recip_mask = donor_mask & duplicated_table
+    # Numerator for PEDS Calc. (Number of Donor features in the Recipient)
     num_sum = np.sum(recip_mask.values, axis=1)
+    # Denominator for PEDS Calc. (Number of unique features in the Donor)
     donor_sum = np.sum(donor_mask, axis=1)
+    # This ignores warnings that come from dividing by 0.
+    # mismatched_peds will be Nan if the denominator is 0 and thats reasonable.
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
-        mismatch_peds = num_sum/donor_sum
-    mismatchpairs_df = _simulate_uniform_distro(recip_df, mismatch_peds,
-                                                replicates)
-    per_sub_stats = _per_subject_stats(actual_temp,  mismatchpairs_df,
-                                       replicates)
+        mismatched_peds = num_sum/donor_sum
+    mismatchedpairs_df = _simulate_uniform_distro(recip_df, mismatched_peds,
+                                                  iterations)
+    per_sub_stats = _per_subject_stats(actual_temp, mismatchedpairs_df,
+                                       iterations)
     global_stats = _global_stats(per_sub_stats['p-value'])
     return per_sub_stats, global_stats
 
 
 # peds_simulation helper functions
-def _create_mismatch_pairs(recip_df, metadata, used_references,
-                           reference_column):
-    mismatch_pairs = []
+
+def _create_mismatched_pairs(recip_df, metadata, used_references,
+                             reference_column):
+    # creates mismatched_pairs of donors and recipients.
+    # Also creates a recipient feature table where there is a duplicate
+    # recipient sample for every pair.
+    # This makes numpy array math much simpler later.
+    mismatched_pairs = []
     for x in itertools.product(recip_df.index,
                                metadata[reference_column].dropna()):
-        mismatch_pairs.append(x)
+        mismatched_pairs.append(x)
     mismatched_t = list(zip(used_references.index, used_references))
     # Removes matched donor and recipient pairs
-    filtered = [tuple for tuple in mismatch_pairs if tuple not in mismatched_t]
+    filtered = [pair for pair in mismatched_pairs if pair not in mismatched_t]
     idx, values = zip(*filtered)
     mismatched_df = pd.DataFrame({"id": idx,
                                   reference_column: values}).set_index("id")
@@ -584,8 +582,11 @@ def _create_mismatch_pairs(recip_df, metadata, used_references,
     return duplicated_table, mismatched_df
 
 
-def _create_sim_masking(mismatch_df, donor_df, reference_column):
-    donors = mismatch_df[reference_column]
+def _create_sim_masking(mismatched_df, donor_df, reference_column):
+    # Creates a Donor Numpy array that duplicates donor samples to match up
+    # with duplicated_table. This will mask recipient feature that aren't
+    # in the donor.
+    donors = mismatched_df[reference_column]
     donor_index_masking = []
     for donor in donors:
         donor_index_masking.append(donor_df.index.get_loc(donor))
@@ -595,32 +596,34 @@ def _create_sim_masking(mismatch_df, donor_df, reference_column):
     return donor_mask
 
 
-def _simulate_uniform_distro(recip_df, mismatch_peds, replicates):
-    mismatchpairs_df = pd.DataFrame([], columns=range(replicates))
+def _simulate_uniform_distro(recip_df, mismatched_peds, iterations):
+    # samples a list of mismatched_peds calculations n (iterations) number
+    # of times uniformly.
+    mismatchedpairs_df = pd.DataFrame([], columns=range(iterations))
     for x in range(len(recip_df.index)):
-        peds_iters = random.choices(mismatch_peds, k=replicates)
-        mismatchpairs_df.loc[len(mismatchpairs_df)] = peds_iters
-    return mismatchpairs_df
-
+        peds_iters = random.choices(mismatched_peds, k=iterations)
+        mismatchedpairs_df.loc[len(mismatchedpairs_df)] = peds_iters
+    return mismatchedpairs_df
 
-def _per_subject_stats(actual_temp, shuffled_donor, replicates):
 
+def _per_subject_stats(actual_temp, shuffled_donor, iterations):
     # Calculating per-subject p-values
     actual_element_wise = actual_temp.values[:, None]
     disagree_df = shuffled_donor >= actual_element_wise
     disagree_series = disagree_df.sum(axis=1)
-    agree_series = replicates-disagree_series
+    agree_series = iterations-disagree_series
     # adding 1 here because you can mathmatically can get p-value of 0 from a
     # Monte Carlo Simulation
-    per_subject_p = (disagree_series + 1)/(replicates+1)
+    per_subject_p = (disagree_series + 1)/(iterations+1)
     per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
     per_sub_stats = pd.DataFrame({'A:group': actual_temp.index,
                                  'A:n': 1,
                                   'A:measure': actual_temp.values,
                                   'B:group': "shuffled_recipients",
-                                  'B:n': replicates,
+                                  'B:n': iterations,  # TODO: Change to num
+                                                      # mismatched_pairs
                                   'B:measure': shuffled_donor.mean(axis=1),
-                                  'n': replicates,
+                                  'n': iterations,
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
                                   'q-value': per_subject_q})
@@ -660,6 +663,7 @@ def _per_subject_stats(actual_temp, shuffled_donor, replicates):
 
 
 def _global_stats(p_series):
+    # Calculates global stats based off of the per-subject p-values.
     stats, p = combine_pvalues(p_series, method="stouffer")
     global_stats = pd.DataFrame([["p-values", p_series.size, stats, p,
                                   np.nan]],
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index e7d773b..91eedb7 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -372,14 +372,14 @@
     function=q2_fmt.peds_simulation,
     inputs={'table': FeatureTable[Frequency | RelativeFrequency |
                                   PresenceAbsence]},
-    parameters={'metadata': Metadata, 
+    parameters={'metadata': Metadata,
                 'time_column': Str,
-                'reference_column': Str, 
+                'reference_column': Str,
                 'subject_column': T_subject,
                 'filter_missing_references': Bool,
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str],
-                'replicates': Int % Range(999, None)},
+                'iterations': Int % Range(99, None)},
     outputs=[('per_subject_stats', StatsTable[Pairwise]),
              ('global_stats', StatsTable[Pairwise])],
     parameter_descriptions={
@@ -396,31 +396,39 @@
         'filter_missing_references': 'Filter out references contained within'
                                      ' the metadata that are not present'
                                      ' in the input feature table.'
-                                     ' Default behavior is to raise an error if some reference ids are in the metadata but not the feature table.',
+                                     ' Default behavior is to raise an error'
+                                     ' if some reference ids are in the'
+                                     ' metadata but not the feature table.',
         'drop_incomplete_subjects': 'Filter out subjects that do not have'
                                     ' a sample at every timepoint.'
-                                    ' Default behavior is to raise an error.',
-        'drop_incomplete_timepoint': 'Filter out timepoints that do not have'
-                                    ' enough samples.'
-                                    ' Default behavior is to raise an error.',
-        'replicates': 'The number of replicates to run the Monte Carlo'
+                                    ' Default behavior is to raise an error '
+                                    ' if any subject is missing a timepoint',
+        'drop_incomplete_timepoint': 'Filter out a `problematic` timepoint'
+                                     ' that is causing subject to be dropped.'
+                                     ' Default behavior is to raise an error'
+                                     ' if any subject is missing a timepoint.',
+        'iterations': 'The number of iterations to run the Monte Carlo'
                       ' simulation on'
     },
     output_descriptions={
-        'per_subject_stats': 'stats table for comparing if the actual PEDS'
-        ' values to shuffled simulated values per subject',
-        'global_stats': 'stats table for comparing if the actual PEDS values'
-        ' to shuffled simulated values globally.'
+        'per_subject_stats': 'Table describing significance of PEDS scores'
+                             ' compared to mismatched donor-recipient pairs'
+                             ' on a per-subject basis',
+        'global_stats': 'Table describing significance of PEDS scores across'
+                        ' all subjects.'
     },
-    name='',
-    description='A Monte Carlo Simulation that shuffles the recipient and'
-                ' donated microbiome in order to investigate if the true'
-                ' donated microbiome is the most similar to their true FMT'
-                ' recipient. This method is not investigating engraftment'
-                ' extent, but is instead identifying if the recipients'
-                ' engrafted unique features from their respective donated'
-                ' microbiome as opposed to features common to the gut'
-                ' microbiome',
+    name='PEDS Monte Carlo simulation',
+    description='A Monte Carlo simulation that randomizes the relationships'
+                ' between donors and recipients, to test whether the PEDS'
+                ' score between a recipient and their actual donor is'
+                ' significantly higher than PEDS scores between other'
+                ' recipients and a random donor. This is intended to only work'
+                ' in studies where there are distinct donors, and will yield'
+                ' insignificant results if there are too few donors.'
+                ' This method attempts to quantify the extent to which'
+                ' indicator features that are unique to a given donor are'
+                ' transferred to their recipients, as opposed to features that'
+                ' are not indicative of any specific donor.',
     citations=[citations['aggarwala_precise_2021']],
     examples={
         'peds_methods': ex.simulation_peds_method
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 53bfa75..acc19fd 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -20,7 +20,7 @@
                           _check_subject_column, _check_column_type,
                           _drop_incomplete_timepoints, feature_peds,
                           _check_column_missing, _rename_features,
-                          peds_simulation, _create_mismatch_pairs,
+                          peds_simulation, _create_mismatched_pairs,
                           _simulate_uniform_distro, _create_sim_masking)
 
 
@@ -1288,7 +1288,7 @@ def test_high_donor_overlap(self):
                                    time_column="group",
                                    reference_column="Ref",
                                    subject_column="subject",
-                                   replicates=999)
+                                   iterations=999)
         real_median = np.median(stats["A:measure"].values)
         fake_median = np.median(stats["B:measure"].values)
         self.assertGreater(real_median, fake_median)
@@ -1320,7 +1320,7 @@ def test_low_donor_overlap(self):
                                    time_column="group",
                                    reference_column="Ref",
                                    subject_column="subject",
-                                   replicates=999)
+                                   iterations=999)
 
         real_median = np.median(stats["A:measure"].values)
         fake_median = np.median(stats["B:measure"].values)
@@ -1351,7 +1351,7 @@ def test_single_donor(self):
                             time_column="group",
                             reference_column="Ref",
                             subject_column="subject",
-                            replicates=999)
+                            iterations=999)
 
     def test_mismatch_pairs(self):
         metadata_df = pd.DataFrame({
@@ -1377,9 +1377,9 @@ def test_mismatch_pairs(self):
                                     index=['sample1', 'sample2', 'sample3'],
                                     name='Ref')
         used_references.index.name = "id"
-        _, mismatched_df = _create_mismatch_pairs(recip_df, metadata_df,
-                                                  used_references,
-                                                  reference_column='Ref')
+        _, mismatched_df = _create_mismatched_pairs(recip_df, metadata_df,
+                                                    used_references,
+                                                    reference_column='Ref')
         exp_mismatched_df = pd.DataFrame({'id': ["sample1", "sample1",
                                                  "sample2", "sample2",
                                                  "sample3", "sample3"],
@@ -1398,12 +1398,12 @@ def test_uniform_distro(self):
 
         mismatch_peds = [0, 0, 0, 0, 0, 0]
 
-        replicates = 999
+        iterations = 999
 
         mismatchpairs_df = _simulate_uniform_distro(recip_df, mismatch_peds,
                                                     999)
         self.assertEquals(mismatchpairs_df.size,
-                          (recip_df.index.size * replicates))
+                          (recip_df.index.size * iterations))
 
     def test_sim_masking(self):
 

From c9dde880522053875d00d2a14064ec02e0088537 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 25 Jul 2024 16:23:57 -0700
Subject: [PATCH 35/50] addressing

---
 q2_fmt/_peds.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 4d805db..91ad3e1 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -540,6 +540,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                                  used_references,
                                  reference_column)
     donor_mask = _create_sim_masking(mismatched_df, donor_df, reference_column)
+    mismatched_pairs_n = len(duplicated_table.index)
     recip_mask = donor_mask & duplicated_table
     # Numerator for PEDS Calc. (Number of Donor features in the Recipient)
     num_sum = np.sum(recip_mask.values, axis=1)
@@ -553,7 +554,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     mismatchedpairs_df = _simulate_uniform_distro(recip_df, mismatched_peds,
                                                   iterations)
     per_sub_stats = _per_subject_stats(actual_temp, mismatchedpairs_df,
-                                       iterations)
+                                       iterations, mismatched_pairs_n)
     global_stats = _global_stats(per_sub_stats['p-value'])
     return per_sub_stats, global_stats
 
@@ -568,7 +569,7 @@ def _create_mismatched_pairs(recip_df, metadata, used_references,
     # This makes numpy array math much simpler later.
     mismatched_pairs = []
     for x in itertools.product(recip_df.index,
-                               metadata[reference_column].dropna()):
+                               metadata[reference_column].dropna().unique()):
         mismatched_pairs.append(x)
     mismatched_t = list(zip(used_references.index, used_references))
     # Removes matched donor and recipient pairs
@@ -606,7 +607,8 @@ def _simulate_uniform_distro(recip_df, mismatched_peds, iterations):
     return mismatchedpairs_df
 
 
-def _per_subject_stats(actual_temp, shuffled_donor, iterations):
+def _per_subject_stats(actual_temp, shuffled_donor,
+                       iterations, mismatched_pairs_n):
     # Calculating per-subject p-values
     actual_element_wise = actual_temp.values[:, None]
     disagree_df = shuffled_donor >= actual_element_wise
@@ -620,8 +622,7 @@ def _per_subject_stats(actual_temp, shuffled_donor, iterations):
                                  'A:n': 1,
                                   'A:measure': actual_temp.values,
                                   'B:group': "shuffled_recipients",
-                                  'B:n': iterations,  # TODO: Change to num
-                                                      # mismatched_pairs
+                                  'B:n': mismatched_pairs_n,
                                   'B:measure': shuffled_donor.mean(axis=1),
                                   'n': iterations,
                                   'test-statistic': agree_series,

From ee61f8a191ac3e8404418f0f9f8de4bd50f239f9 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Fri, 26 Jul 2024 08:11:36 -0700
Subject: [PATCH 36/50] ref: descriptions

---
 q2_fmt/plugin_setup.py | 261 ++++++++++++++++-------------------------
 1 file changed, 99 insertions(+), 162 deletions(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 91eedb7..f61dd6f 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -46,6 +46,50 @@
     (Str, Str % Choices("baseline", "consecutive")): Visualization
 })
 
+# Description Variables
+metadata = 'The sample `metadata`.'
+time_column = ('The column within `metadata` that the `table` should be'
+               ' grouped by. This column should contain integer values.')
+reference_column = ('The column in `metadata` indicating the sample to use as'
+                    ' the reference. For example, this may list the relevant'
+                    ' donor sample for each recipient sample.')
+subject_column = ('The column within `metadata` that contains the subject'
+                  ' ID for each sample.')
+filter_missing_references = ('Filter out references contained within the'
+                             ' metadata that are not present in the input'
+                             ' feature table. Default behavior is to raise an'
+                             ' error if some reference ids are in the metadata'
+                             ' but not the feature table.')
+drop_incomplete_subjects = ('Filter out subjects that do not have a sample at'
+                            ' every timepoint. Default behavior is to raise an'
+                            ' error if any subject is missing a timepoint.')
+drop_incomplete_timepoint = ('Filter out a `problematic` timepoint that is'
+                             ' causing subject to be dropped. Default behavior'
+                             ' is to raise an error if any subject is missing'
+                             ' a timepoint.')
+level_delimiter = 'delimiter to split taxonomic string on'
+control_column = ('The column within `metadata` that contains any relevant'
+                  ' control group IDs. Actual treatment samples should not'
+                  ' contain any value within this column.')
+distance_to = 'What reference type to calculate distance to against'
+baseline_timepoint = ('If `baseline` is selected for `distance_to`, the'
+                      ' timepoint to use as baseline reference.')
+where = ('Additional filtering for the associated `metadata` file. This can be'
+         ' used to filter by a subset of `metadata`, such as a specific value'
+         ' in one of `metadata` columns.')
+
+per_subject_stats = ('Table describing significance of PEDS scores compared to'
+                     ' mismatched donor-recipient pairs on a per-subject'
+                     ' basis.')
+global_stats = ('Table describing significance of PEDS scores across all'
+                ' subjects.')
+peds_table = 'A feature table to calculate PEDS on.'
+peds_dists = ('The distributions for the PEDS measure, grouped by the selected'
+              ' `time_column`. also contains the Numerator and Denominator for'
+              ' PEDS calulations. May also contain subject IDs, if'
+              ' `subject_column` is  provided in `metadata`.')
+
+
 plugin.pipelines.register_function(
     function=q2_fmt.engraftment,
     inputs={'diversity_measure': DistanceMatrix | SampleData[AlphaDiversity]},
@@ -67,7 +111,7 @@
                              ' microbiome',
     },
     parameter_descriptions={
-        'metadata': 'The sample `metadata`.',
+        'metadata': metadata,
         'compare': 'The type of comparison that will be used to analyze the'
                    ' input `diversity_measure`.'
                    ' The "baseline" comparison defines Group A as the'
@@ -83,30 +127,14 @@
                    ' The "all-pairwise" comparison defines Group A as all'
                    ' groups in `reference_column` and `control_column`, and'
                    ' Group B is all timepoints in `time_column`.',
-        'distance_to': 'What reference type to calculate distance to against',
-        'time_column': 'The column within the `metadata` that the'
-                       ' `diversity_measure` should be grouped by.'
-                       ' This column should contain simple integer values.',
-        'control_column': 'The column within the `metadata` that contains any'
-                          ' relevant control group IDs.'
-                          ' Actual treatment samples should not contain any'
-                          ' value within this column.',
-        'reference_column': 'The column within the `metadata` that contains'
-                            ' the sample to use as a reference for a given'
-                            ' beta `diversity_measure`. For example, this'
-                            ' may be the relevant donor sample to compare'
-                            ' against.',
-        'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID to be tracked against timepoints.',
-        'filter_missing_references': 'Filter out references contained within'
-                                     ' the metadata that are not present in'
-                                     ' the diversity measure.'
-                                     ' Default behavior is to raise an error.',
-        'baseline_timepoint': 'If `baseline` is selected for `distance_to`,'
-                              ' the timepoint to use as baseline reference.',
-        'where': 'Additional filtering for the associated `metadata` file.'
-                 ' This can be used to filter by a subset of the `metadata`,'
-                 ' such as a specific value in one of the `metadata` columns.',
+        'distance_to': distance_to,
+        'time_column': time_column,
+        'control_column': control_column,
+        'reference_column': reference_column,
+        'subject_column': subject_column,
+        'filter_missing_references': filter_missing_references,
+        'baseline_timepoint': baseline_timepoint,
+        'where': where,
         'against_group': 'Based on the selected comparison, this is the column'
                          ' that will be used to compare all samples against.',
         'alternative': 'The "two-sided" alternative hypothesis is that the'
@@ -157,34 +185,18 @@
                              ' structure'
     },
     parameter_descriptions={
-        'metadata': 'The sample metadata.',
-        'distance_to': 'What reference type to calculate distance to against',
-        'time_column': 'The column within the `metadata` that the'
-                       ' `diversity_measure` should be grouped by.'
-                       ' This column should contain simple integer values.',
-        'control_column': 'The column within the `metadata` that contains any'
-                          ' relevant control group IDs.'
-                          ' Actual treatment samples should not contain any'
-                          ' value within this column.',
-        'reference_column': 'The column within the `metadata` that contains'
-                            ' the sample to use as a reference'
-                            ' for a given beta `diversity_measure`.'
-                            ' For example, this may be the relevant donor'
-                            ' sample to compare against.',
-        'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID to be tracked against timepoints.',
+        'metadata': metadata,
+        'distance_to': distance_to,
+        'time_column': time_column,
+        'control_column': control_column,
+        'reference_column': reference_column,
+        'subject_column': subject_column,
         'group_column': 'The column within the metadata that contains'
                         ' information about groups (ex: treatment group)'
                         ' in order to compare engraftment between groups',
-        'filter_missing_references': 'Filter out references contained within'
-                                     ' the metadata that are not present'
-                                     ' in the diversity measure.'
-                                     ' Default behavior is to raise an error.',
-        'where': 'Additional filtering for the associated `metadata` file.'
-                 ' This can be used to filter by a subset of the `metadata`,'
-                 ' such as a specific value in one of the `metadata` columns.',
-        'baseline_timepoint': 'If `baseline` is selected for `distance_to`,'
-                              ' the timepoint to use as baseline reference.',
+        'filter_missing_references': filter_missing_references,
+        'where': where,
+        'baseline_timepoint': baseline_timepoint,
     },
     output_descriptions={
         'timepoint_dists': 'The distributions for the `diversity_measure`,'
@@ -222,33 +234,17 @@
                 'drop_incomplete_timepoint': List[Str],
                 'level_delimiter': Str},
     outputs=[('peds_heatmap', Visualization)],
-    input_descriptions={'table': 'A feature table to calculate PEDS on'},
+    input_descriptions={'table': peds_table},
     parameter_descriptions={
-        'metadata': 'The sample metadata.',
+        'metadata': metadata,
         'peds_metric': 'PEDS metric to run.',
-        'time_column': 'The column within the `metadata` that the'
-                       ' `table` should be grouped by. This column'
-                       ' should contain simple integer values.',
-        'reference_column': 'The column within the `metadata` that contains'
-                            ' the sample to use as a reference'
-                            ' for a given `table`.'
-                            ' For example, this may be the relevant donor'
-                            ' sample to compare against.',
-        'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID to be tracked against timepoints.',
-        'filter_missing_references': 'Filter out references contained within'
-                                     ' the metadata that are not present'
-                                     ' in the table.'
-                                     ' Default behavior is to raise an error.',
-        'drop_incomplete_subjects': 'Filter out subjects that do not have'
-                                    ' a sample at every timepoint.'
-                                    ' Default behavior is to raise an error.',
-        'drop_incomplete_timepoint': 'Filter out provided timepoint.'
-                                     ' This will be performed before'
-                                     ' drop_incomplete_subjects if the'
-                                     ' drop_incomplete_subjects parameter is'
-                                     ' passed.',
-        'level_delimiter': 'delimiter to split taxonomic string on'},
+        'time_column': time_column,
+        'reference_column': reference_column,
+        'subject_column': subject_column,
+        'filter_missing_references': filter_missing_references,
+        'drop_incomplete_subjects': drop_incomplete_subjects,
+        'drop_incomplete_timepoint': drop_incomplete_timepoint,
+        'level_delimiter': level_delimiter},
     output_descriptions={'peds_heatmap': 'PEDS heatmap visualization'},
     name='PEDS pipeline to calculate feature or sample PEDS',
     description='Runs a pipeline to calculate sample or feature PEDS,'
@@ -261,16 +257,11 @@
             'per_subject_stats': StatsTable[Pairwise],
             'global_stats': StatsTable[Pairwise]},
     input_descriptions={'data': 'PEDS output to plot',
-                        'per_subject_stats': 'stats table for comparing if the'
-                                             ' actual PEDS values to shuffled'
-                                             ' simulated values per subject',
-                        'global_stats': 'stats table for comparing if the'
-                                        ' actual PEDS values to shuffled'
-                                        ' simulated values globally.'},
+                        'per_subject_stats': per_subject_stats,
+                        'global_stats': global_stats},
     parameters={'level_delimiter': Str},
     parameter_descriptions={
-                            'level_delimiter': 'delimiter to split taxonomic'
-                                               ' string on'},
+                            'level_delimiter': level_delimiter},
     name='PEDS Heatmap',
     description='Plot heatmap for PEDS value over time'
 )
@@ -285,38 +276,18 @@
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str]},
     outputs=[('peds_dists', Dist1D[Ordered, Matched] % Properties("peds"))],
-    input_descriptions={'table': 'A feature table to calculate PEDS on'},
+    input_descriptions={'table': peds_table},
     parameter_descriptions={
-        'metadata': 'The sample metadata.',
-        'time_column': 'The column within the `metadata` that the'
-                       ' `table` should be grouped by. This column'
-                       ' should contain simple integer values.',
-        'reference_column': 'The column within the `metadata` that contains'
-                            ' the sample to use as a reference'
-                            ' for a given `table`.'
-                            ' For example, this may be the relevant donor'
-                            ' sample to compare against.',
-        'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID to be tracked against timepoints.',
-        'filter_missing_references': 'Filter out references contained within'
-                                     ' the metadata that are not present'
-                                     ' in the table.'
-                                     ' Default behavior is to raise an error.',
-        'drop_incomplete_subjects': 'Filter out subjects that do not have'
-                                    ' a sample at every timepoint.'
-                                    ' Default behavior is to raise an error.',
-        'drop_incomplete_timepoint': 'Filter out a list of provided timepoint.'
-                                     ' This will be performed before'
-                                     ' drop_incomplete_subjects if the'
-                                     ' drop_incomplete_subjects parameter is'
-                                     ' passed.'
+        'metadata': metadata,
+        'time_column': time_column,
+        'reference_column': reference_column,
+        'subject_column': subject_column,
+        'filter_missing_references': filter_missing_references,
+        'drop_incomplete_subjects': drop_incomplete_subjects,
+        'drop_incomplete_timepoint': drop_incomplete_timepoint
     },
     output_descriptions={
-        'peds_dists': 'The distributions for the PEDS measure,'
-                      ' grouped by the selected `time_column`.'
-                      ' Also contains the Numerator and Denominator for'
-                      ' PEDS calculations. May also contain subject IDs,'
-                      ' if `subject_column` is provided in the `metadata`.'
+        'peds_dists': peds_dists
     },
     name='Proportional Engraftment of Donor Strains (Features) in each'
          ' recipient sample',
@@ -335,30 +306,16 @@
                 'reference_column': Str, 'subject_column': Str,
                 'filter_missing_references': Bool},
     outputs=[('peds_dists', Dist1D[Ordered, Matched] % Properties("peds"))],
-    input_descriptions={'table': 'A feature table to calculate PEDS on.'},
+    input_descriptions={'table': peds_table},
     parameter_descriptions={
-        'metadata': 'The sample metadata.',
-        'time_column': 'The column within the `metadata` that the'
-                       ' `table` should be grouped by. This column'
-                       ' should contain simple integer values.',
-        'reference_column': 'The column within the `metadata` that contains'
-                            ' the sample to use as a reference'
-                            ' for a given `table`.'
-                            ' For example, this may be the relevant donor'
-                            ' sample to compare against.',
-        'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID to be tracked against timepoints.',
-        'filter_missing_references': 'Filter out references contained within'
-                                     ' the metadata that are not present'
-                                     ' in the table.'
-                                     ' Default behavior is to raise an error.',
+        'metadata': metadata,
+        'time_column': time_column,
+        'reference_column': reference_column,
+        'subject_column': subject_column,
+        'filter_missing_references': filter_missing_references,
     },
     output_descriptions={
-        'peds_dists': 'The distributions for the PEDS measure,'
-                      ' grouped by the selected `time_column`.'
-                      ' also contains the Numerator and Denominator for'
-                      ' PEDS calulations. May also contain subject IDs,'
-                      ' if `subject_column` is provided in the `metadata`.'
+        'peds_dists': peds_dists
     },
     name='Porportional Engraftment of Donor Strains per feature',
     description='Calculates how many recipients recieved a given'
@@ -383,39 +340,19 @@
     outputs=[('per_subject_stats', StatsTable[Pairwise]),
              ('global_stats', StatsTable[Pairwise])],
     parameter_descriptions={
-        'metadata': 'The sample metadata.',
-        'time_column': 'The column within `metadata` that the'
-                       ' `table` should be grouped by. This column'
-                       ' should contain integer values.',
-        'reference_column': 'The column in `metadata` indicating'
-                            ' the sample to use as the reference.'
-                            ' For example, this may list the relevant donor'
-                            ' sample for each recipient sample.',
-        'subject_column': 'The column within the `metadata` that contains the'
-                          ' subject ID for each sample.',
-        'filter_missing_references': 'Filter out references contained within'
-                                     ' the metadata that are not present'
-                                     ' in the input feature table.'
-                                     ' Default behavior is to raise an error'
-                                     ' if some reference ids are in the'
-                                     ' metadata but not the feature table.',
-        'drop_incomplete_subjects': 'Filter out subjects that do not have'
-                                    ' a sample at every timepoint.'
-                                    ' Default behavior is to raise an error '
-                                    ' if any subject is missing a timepoint',
-        'drop_incomplete_timepoint': 'Filter out a `problematic` timepoint'
-                                     ' that is causing subject to be dropped.'
-                                     ' Default behavior is to raise an error'
-                                     ' if any subject is missing a timepoint.',
+        'metadata': metadata,
+        'time_column': time_column,
+        'reference_column': reference_column,
+        'subject_column': subject_column,
+        'filter_missing_references': filter_missing_references,
+        'drop_incomplete_subjects': drop_incomplete_subjects,
+        'drop_incomplete_timepoint': drop_incomplete_timepoint,
         'iterations': 'The number of iterations to run the Monte Carlo'
                       ' simulation on'
     },
     output_descriptions={
-        'per_subject_stats': 'Table describing significance of PEDS scores'
-                             ' compared to mismatched donor-recipient pairs'
-                             ' on a per-subject basis',
-        'global_stats': 'Table describing significance of PEDS scores across'
-                        ' all subjects.'
+        'per_subject_stats': per_subject_stats,
+        'global_stats': global_stats
     },
     name='PEDS Monte Carlo simulation',
     description='A Monte Carlo simulation that randomizes the relationships'

From 811f661bfba646e3a8225f0b13f8b494a17f5bb2 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 30 Jul 2024 16:57:38 -0700
Subject: [PATCH 37/50] greg comments

---
 q2_fmt/_peds.py                  | 313 ++++++++++++++++++++++++++++---
 q2_fmt/tests/test_engraftment.py | 100 +++++++++-
 2 files changed, 379 insertions(+), 34 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 91ad3e1..4e01749 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -534,14 +534,15 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     table = table > 0
     recip_df = _create_recipient_table(used_references, metadata_df, table)
     donor_df = table[table.index.isin(used_references)]
-    duplicated_table, mismatched_df = \
+    mismatched_df = \
         _create_mismatched_pairs(recip_df,
                                  metadata_df,
                                  used_references,
                                  reference_column)
+    duplicated_recip_table = _create_duplicated_recip_table(mismatched_df,
+                                                            recip_df)
     donor_mask = _create_sim_masking(mismatched_df, donor_df, reference_column)
-    mismatched_pairs_n = len(duplicated_table.index)
-    recip_mask = donor_mask & duplicated_table
+    recip_mask = _mask_recipient(donor_mask, recip_df)
     # Numerator for PEDS Calc. (Number of Donor features in the Recipient)
     num_sum = np.sum(recip_mask.values, axis=1)
     # Denominator for PEDS Calc. (Number of unique features in the Donor)
@@ -551,10 +552,10 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
         mismatched_peds = num_sum/donor_sum
-    mismatchedpairs_df = _simulate_uniform_distro(recip_df, mismatched_peds,
-                                                  iterations)
-    per_sub_stats = _per_subject_stats(actual_temp, mismatchedpairs_df,
-                                       iterations, mismatched_pairs_n)
+    mismatched_pairs_n = len(duplicated_recip_table.index)
+    per_sub_stats = _per_subject_stats(mismatched_peds,
+                                       actual_temp, iterations,
+                                       mismatched_pairs_n)
     global_stats = _global_stats(per_sub_stats['p-value'])
     return per_sub_stats, global_stats
 
@@ -563,10 +564,59 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
 
 def _create_mismatched_pairs(recip_df, metadata, used_references,
                              reference_column):
-    # creates mismatched_pairs of donors and recipients.
-    # Also creates a recipient feature table where there is a duplicate
-    # recipient sample for every pair.
-    # This makes numpy array math much simpler later.
+    """Creates a Dataframe of all the incorrect Donor-Recipient pairs
+
+    Creates a list of tuples of all donor-recipient pairs then filters out
+    any real donor-recipient pairs. The result is a Dataframe that includes
+    only mismatched donor-recipient pairs.
+
+    Parameters
+    ----------
+    recip_df: pd.DataFrame
+        A feature table of FMT recipients.
+    metadata: pd.DataFrame
+        The sample metadata for the study.
+    used_references: pd.Series
+        A series with recipients as the index and the associated reference
+        as the values.
+    reference_column: Str
+        Name of the reference column in the Sample Metadata.
+
+    Returns
+    -------
+    mismatched_df: pd.DataFrame
+        A DataFrame containing all mismatched pairs of donors and recipients.
+
+    Examples
+    --------
+    recip_df = pd.DataFrame({
+        'id': ['sample1', 'sample2', 'sample3'],
+        'Feature1': [1, 0, 0],
+        'Feature2': [0, 1, 0],
+        'Feature3': [0, 0, 1]}).set_index('id')
+
+    metadata_df = pd.DataFrame({
+        'id': ['sample1', 'sample2', 'sample3','donor1', 'donor2', 'donor3'],
+        'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan, np.nan],
+        'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan, np.nan],
+        'group': [1, 1, 1, np.nan, np.nan, np.nan],
+        "Location": [np.nan, np.nan,np.nan, 'test', 'test','test']}
+        ).set_index('id')
+
+    used_references = pd.Series(data=['donor1', 'donor2', 'donor3'],
+                                    index=['sample1', 'sample2', 'sample3'],
+                                    name='Ref')
+    used_references.index.name = "id"
+
+    reference_column = "Ref"
+
+    _create_mismatched_pairs(recip_df, metadata_df, used_references,
+                             reference_column)
+    pd.DataFrame({'id': ["sample1", "sample1", "sample2", "sample2",
+                         "sample3", "sample3"],
+                  "Ref": ["donor2", "donor3", "donor1", "donor3",
+                          "donor1", "donor2"]}).set_index('id')
+    """
     mismatched_pairs = []
     for x in itertools.product(recip_df.index,
                                metadata[reference_column].dropna().unique()):
@@ -577,41 +627,216 @@ def _create_mismatched_pairs(recip_df, metadata, used_references,
     idx, values = zip(*filtered)
     mismatched_df = pd.DataFrame({"id": idx,
                                   reference_column: values}).set_index("id")
+    return mismatched_df
+
+
+def _create_duplicated_recip_table(mismatched_df, recip_df):
+    """Creates a recipient feature table that is the same dimensions as
+    the mismatched_df
+
+    Creates a recipient table that duplicates feature information so that the
+    recip_df is the same dimension as the mismatched_df. This makes numpy array
+    math easier later on.
+
+    Parameters
+    ----------
+    mismatched_df: pd.DataFrame
+        A Dataframe that contains recipient samples as the index and mismatched
+        donors as the values. A recipient sample will appear as many times as
+        there are mismatched donors to pair with.
+    recip_df: pd.DataFrame
+        A feature table of FMT recipients.
+
+    Returns
+    -------
+    duplicated_table: pd.DataFrame
+        A recipent feature table where feature information is duplicated so
+        that recipient is in the same diminsion as mismatched_df
+
+    Examples
+    --------
+    mismatched_df = pd.DataFrame({
+            'id': ["sample1", "sample1", "sample2", "sample2",
+                   "sample3", "sample3"],
+            "Ref": ["donor2", "donor3", "donor1", "donor3",
+                    "donor1", "donor2"]}
+                                ).set_index('id')
+
+    recip_df = pd.DataFrame({
+        'id': ['sample1', 'sample2', 'sample3'],
+        'Feature1': [1, 0, 0],
+        'Feature2': [0, 1, 0],
+        'Feature3': [0, 0, 1]}).set_index('id')
+
+    _create_duplicated_recip_table(mismatched_df, recip_df)
+
+    pd.DataFrame({
+            'id': ["sample1", "sample1", "sample2", "sample2",
+                   "sample3", "sample3"],
+            'Feature1': [1, 1, 0, 0, 0, 0],
+            'Feature2': [0, 0, 1, 1, 0, 0],
+            'Feature3': [0, 0, 0, 0, 1, 1]}).set_index('id')
+    """
     duplicated_table = \
         mismatched_df.merge(recip_df, left_index=True, right_index=True,
                             how='left').drop(mismatched_df.columns, axis=1)
-    return duplicated_table, mismatched_df
+    return duplicated_table
 
 
 def _create_sim_masking(mismatched_df, donor_df, reference_column):
-    # Creates a Donor Numpy array that duplicates donor samples to match up
-    # with duplicated_table. This will mask recipient feature that aren't
-    # in the donor.
+    """Create a donor mask to mask recipient feature that aren't in the donor.
+
+    Creates a Donor Numpy array that duplicates donor samples to match up
+    with duplicated_table. This will mask recipient feature that aren't
+    in the donor.
+
+    Parameters
+    ----------
+    mismatched_df: pd.DataFrame
+        A Dataframe that contains recipient samples as the index and mismatched
+        donors as the values. A recipient sample will appear as many times as
+        there are mismatched donors to pair with.
+    donor_df: pd.DataFrame
+        A feature table of FMT donors.
+    reference_column: str
+        Name of the reference column in the Sample Metadata.
+
+    Returns
+    -------
+    donor_mask: ndarray
+        a numpy array of donor feature information with the same order as
+        duplicated_recip_table. This allows for easy np array math.
+
+    Examples
+    --------
+    mismatched_df = pd.DataFrame({
+            'id': ["sample1", "sample1", "sample2", "sample2",
+                   "sample3", "sample3"],
+            "Ref": ["donor2", "donor3", "donor1", "donor3",
+                    "donor1", "donor2"]}
+                                ).set_index('id')
+    donor_df = pd.DataFrame({
+            'id': ['donor1', 'donor2', 'donor3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+    _create_sim_masking(mismatched_df, donor_df, reference_column)
+
+    ndarray[[0, 1, 0],
+            [0, 0, 1],
+            [1, 0, 0],
+            [0, 0, 1],
+            [1, 0, 0],
+            [0, 1, 0]]
+    """
     donors = mismatched_df[reference_column]
     donor_index_masking = []
     for donor in donors:
         donor_index_masking.append(donor_df.index.get_loc(donor))
     donor_df_np = donor_df.to_numpy()
     donor_mask = donor_df_np[donor_index_masking]
-    donor_mask = donor_mask.astype(int)
     return donor_mask
 
 
-def _simulate_uniform_distro(recip_df, mismatched_peds, iterations):
-    # samples a list of mismatched_peds calculations n (iterations) number
-    # of times uniformly.
-    mismatchedpairs_df = pd.DataFrame([], columns=range(iterations))
-    for x in range(len(recip_df.index)):
-        peds_iters = random.choices(mismatched_peds, k=iterations)
-        mismatchedpairs_df.loc[len(mismatchedpairs_df)] = peds_iters
-    return mismatchedpairs_df
+def _simulate_uniform_distro(mismatched_peds, iterations):
+    """Randomly samples the mismatched PEDS values.
+
+    Creates a uniform distribution of mismatched PEDS values by randomly
+    sampling N (iterations) times
+
+    Parameters
+    ----------
+    mismatched_peds: pd.DataFrame
+        A Dataframe that contains all mismatched PEDS values.
+    iterations: int
+        Number of iterations to run simulations (Number of times to randomly
+        sample mismatched_peds)
+
+    Returns
+    -------
+    peds_iters: pd.Series
+        a pd.Series with all iterations of mismatched PEDS Values. This will
+        later be compared to an actual PEDS value.
+
+    Examples
+    --------
+    mismatched_peds = [0, 0, 0,0 ]
+    iterations = 10
 
+    _simulate_uniform_distro(mismatched_peds, iterations)
 
-def _per_subject_stats(actual_temp, shuffled_donor,
+    pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
+               index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    """
+    peds_iters = random.choices(mismatched_peds, k=iterations)
+    # Tranforming to series for easy series math in _per_subject_stats()
+    peds_iters = pd.Series(peds_iters)
+    return peds_iters
+
+
+def _per_subject_stats(mismatched_peds, actual_temp,
                        iterations, mismatched_pairs_n):
-    # Calculating per-subject p-values
-    actual_element_wise = actual_temp.values[:, None]
-    disagree_df = shuffled_donor >= actual_element_wise
+    """Creates per subject PEDS stats
+
+    Creates per subject PEDS stats by sampling mismatch PEDS values 'iteration'
+    number of time. The list of mismatched PEDS values is compared to the
+    actual PEDS value. The Test Statistic for this is the iterations - the
+    number of mismatched PEDS value that are greater than the actual p-value.
+    p-values are corrected by adding 1 to each side of the faction
+    because Monte Carlo Simulation is able to produce an invalid p-value of 0.
+
+    Parameters
+    ----------
+    mismatched_peds: pd.DataFrame
+        A Dataframe that contains all mismatched PEDS values.
+    actual_temp: pd.Series
+        A Series containing Sample IDs as the index and actual PEDS values for
+        the sample as the value.
+    iterations: int
+        Number of iterations to run simulations (Number of times to randomly
+        sample mismatched_peds)
+    mismatched_peds_n: int
+        Number of mismatched donor pairs are possible.
+
+    Returns
+    -------
+    per_subject_stats: pd.DataFrame
+        a Dataframe containing test statistics, p-values, q-values and
+        sample size info. Note: this is the DataFrame that is returned
+        to the user as the per-subject-stats.qza
+
+    Examples
+    --------
+    mismatched_peds = [0, 0, 0, 0]
+    actual_temp = pd.Series([data = [1, 1, 1, 1],
+               index = [sample1, sample2, sample3, sample4]])
+    iterations = 10
+    mismatched_pairs_n = 4
+
+    _per_subject_stats(mismatched_peds, actual_temp,
+                       iterations, mismatched_pairs_n):
+
+    pd.Dataframe({"A:group": [sample1, sample2, sample3, sample4],
+                   "A:n": [1, 1, 1, 1],
+                   "A:measure" : [1, 1, 1, 1],
+                   "B:group" : [shuffled_recipients, shuffled_recipients,
+                                shuffled_recipients, shuffled_recipients],
+                   "B:n": [4, 4, 4, 4],
+                   "B:measure": [0, 0, 0, 0],
+                   "n": [10, 10, 10, 10],
+                   "test-statistic": [10, 10, 10, 10],
+                   "p-value": [0.001, 0.001, 0.001, 0.001],
+                   "q-value": [0.004, 0.002, 0.00133, 0.001]})
+    """
+    disagree_df = pd.DataFrame([])
+    peds_iter_mean = []
+    for value in actual_temp:
+        peds_iters = _simulate_uniform_distro(mismatched_peds, iterations)
+        peds_iter_mean.append(peds_iters.mean())
+        disagree_recipient = peds_iters >= value
+        # add disagree recipient as a row in the disagree df
+        disagree_df = pd.concat([disagree_df, disagree_recipient.to_frame().T])
     disagree_series = disagree_df.sum(axis=1)
     agree_series = iterations-disagree_series
     # adding 1 here because you can mathmatically can get p-value of 0 from a
@@ -623,7 +848,7 @@ def _per_subject_stats(actual_temp, shuffled_donor,
                                   'A:measure': actual_temp.values,
                                   'B:group': "shuffled_recipients",
                                   'B:n': mismatched_pairs_n,
-                                  'B:measure': shuffled_donor.mean(axis=1),
+                                  'B:measure': peds_iter_mean,
                                   'n': iterations,
                                   'test-statistic': agree_series,
                                   'p-value': per_subject_p,
@@ -664,7 +889,35 @@ def _per_subject_stats(actual_temp, shuffled_donor,
 
 
 def _global_stats(p_series):
-    # Calculates global stats based off of the per-subject p-values.
+    """Creates global PEDS stats
+
+    Uses Stouffers Method on the per-subject-values to create a global PEDS
+    statistic.
+
+    Parameters
+    ----------
+    p_series: pd.Series
+        Series of all per-subject-stats
+
+    Returns
+    -------
+    global_stats: pd.DataFrame
+        a Dataframe containing global stat statistics, p-values, and q-values.
+        Note: this is the DataFrame that is returned
+        to the user as the global-stats.qza
+
+    Examples
+    --------
+    p_series = pd.Series([data = [0.001, 0.001, 0.001, 0.001])
+
+    _global_stats(p_series)
+
+    pd.Dataframe({"measure": [p-values],
+                   "n": [4],
+                   "test-statistic": [20.93],
+                   "p-value": [0.001]
+                   "q-value": [NaN]})
+    """
     stats, p = combine_pvalues(p_series, method="stouffer")
     global_stats = pd.DataFrame([["p-values", p_series.size, stats, p,
                                   np.nan]],
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index acc19fd..0173cdf 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -21,7 +21,9 @@
                           _drop_incomplete_timepoints, feature_peds,
                           _check_column_missing, _rename_features,
                           peds_simulation, _create_mismatched_pairs,
-                          _simulate_uniform_distro, _create_sim_masking)
+                          _simulate_uniform_distro, _create_sim_masking,
+                          _mask_recipient, _create_duplicated_recip_table,
+                          _per_subject_stats, _global_stats)
 
 
 class TestBase(TestPluginBase):
@@ -1353,7 +1355,7 @@ def test_single_donor(self):
                             subject_column="subject",
                             iterations=999)
 
-    def test_mismatch_pairs(self):
+    def test_create_mismatched_pairs(self):
         metadata_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3',
                    'donor1', 'donor2', 'donor3'],
@@ -1389,7 +1391,18 @@ def test_mismatch_pairs(self):
                                          ).set_index('id')
         pd.testing.assert_frame_equal(mismatched_df, exp_mismatched_df)
 
-    def test_uniform_distro(self):
+    def test_mask_recipient(self):
+        donor_mask = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
+        recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+        recip_mask = _mask_recipient(donor_mask, recip_df)
+        exp_r_mask = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
+        np.testing.assert_array_equal(recip_mask, exp_r_mask)
+
+    def test_simulate_uniform_distro(self):
         recip_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3'],
             'Feature1': [1, 0, 0],
@@ -1405,7 +1418,7 @@ def test_uniform_distro(self):
         self.assertEquals(mismatchpairs_df.size,
                           (recip_df.index.size * iterations))
 
-    def test_sim_masking(self):
+    def test_create_sim_masking(self):
 
         mismatched_df = pd.DataFrame({'id': ["sample1", "sample1",
                                              "sample2", "sample2",
@@ -1431,3 +1444,82 @@ def test_sim_masking(self):
         donor_mask = _create_sim_masking(mismatched_df, donor_df,
                                          reference_column='Ref')
         np.testing.assert_array_equal(donor_mask, exp_mask)
+
+    def test_create_duplicated_table(self):
+        recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+        mismatched_df = pd.DataFrame({'id': ["sample1", "sample1",
+                                             "sample2", "sample2",
+                                             "sample3", "sample3"],
+                                      "Ref": ["donor2", "donor3",
+                                              "donor1", "donor3",
+                                              "donor1", "donor2"]}
+                                     ).set_index('id')
+
+        duplicated_recip_table = _create_duplicated_recip_table(mismatched_df,
+                                                                recip_df)
+
+        exp_d_r_table = pd.DataFrame({
+            'id': ['sample1', 'sample1', 'sample2', 'sample2',
+                   'sample3', 'sample3'],
+            'Feature1': [1, 1, 0, 0, 0, 0],
+            'Feature2': [0, 0, 1, 1, 0, 0],
+            'Feature3': [0, 0, 0, 0, 1, 1]}).set_index('id')
+
+        pd.testing.assert_frame_equal(duplicated_recip_table,
+                                      exp_d_r_table)
+
+    def test_per_subject_stats_labels(self):
+        mismatched_peds = [0, 0, 0, 0]
+        actual_temp = pd.Series(data=[1, 1, 1, 1],
+                                index=["sample1", "sample2", "sample3",
+                                       "sample4"])
+        iterations = 10
+        mismatched_pairs_n = 4
+
+        p_s_stats = _per_subject_stats(mismatched_peds,
+                                       actual_temp, iterations,
+                                       mismatched_pairs_n)
+
+        exp_column_names = ["A:group", "A:n", "A:measure",
+                            "B:group", "B:n", "B:measure", "n",
+                            "test-statistic", "p-value", "q-value"]
+        np.testing.assert_array_equal(p_s_stats.columns.values,
+                                      exp_column_names)
+
+    def test_per_subject_stats(self):
+        mismatched_peds = [0, 0, 0, 0]
+        actual_temp = pd.Series(data=[1, 1, 1, 1],
+                                index=["sample1", "sample2", "sample3",
+                                       "sample4"])
+        iterations = 10
+        mismatched_pairs_n = 4
+
+        p_s_stats = _per_subject_stats(mismatched_peds,
+                                       actual_temp, iterations,
+                                       mismatched_pairs_n)
+
+        exp_test_stats = pd.Series([10, 10, 10, 10])
+
+        np.testing.assert_array_equal(p_s_stats["test-statistic"].values,
+                                      exp_test_stats.values)
+
+    def test_global_stats_label(self):
+        p_series = pd.Series(data=[0.001, 0.001, 0.001, 0.001])
+
+        p = _global_stats(p_series)
+
+        exp_labels = ["Measure", "n", "test-statistic", "p-value", "q-value"]
+        np.testing.assert_array_equal(p.columns.values,
+                                      exp_labels)
+
+    def test_global_stats(self):
+        p_series = pd.Series(data=[0.001, 0.001, 0.001, 0.001])
+
+        p = _global_stats(p_series)
+
+        np.testing.assert_array_equal(p["n"].values, [4])

From 5af31bfe72760fb056dd5f6286dc36c71f3b5280 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 30 Jul 2024 17:18:26 -0700
Subject: [PATCH 38/50] fix tests

---
 q2_fmt/_peds.py                  |  2 +-
 q2_fmt/tests/test_engraftment.py | 11 +++++------
 2 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 4e01749..fca8bb8 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -542,7 +542,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     duplicated_recip_table = _create_duplicated_recip_table(mismatched_df,
                                                             recip_df)
     donor_mask = _create_sim_masking(mismatched_df, donor_df, reference_column)
-    recip_mask = _mask_recipient(donor_mask, recip_df)
+    recip_mask = _mask_recipient(donor_mask, duplicated_recip_table)
     # Numerator for PEDS Calc. (Number of Donor features in the Recipient)
     num_sum = np.sum(recip_mask.values, axis=1)
     # Denominator for PEDS Calc. (Number of unique features in the Donor)
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 0173cdf..2f9de57 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1379,9 +1379,9 @@ def test_create_mismatched_pairs(self):
                                     index=['sample1', 'sample2', 'sample3'],
                                     name='Ref')
         used_references.index.name = "id"
-        _, mismatched_df = _create_mismatched_pairs(recip_df, metadata_df,
-                                                    used_references,
-                                                    reference_column='Ref')
+        mismatched_df = _create_mismatched_pairs(recip_df, metadata_df,
+                                                 used_references,
+                                                 reference_column='Ref')
         exp_mismatched_df = pd.DataFrame({'id': ["sample1", "sample1",
                                                  "sample2", "sample2",
                                                  "sample3", "sample3"],
@@ -1413,10 +1413,9 @@ def test_simulate_uniform_distro(self):
 
         iterations = 999
 
-        mismatchpairs_df = _simulate_uniform_distro(recip_df, mismatch_peds,
+        mismatchpairs_df = _simulate_uniform_distro(mismatch_peds,
                                                     999)
-        self.assertEquals(mismatchpairs_df.size,
-                          (recip_df.index.size * iterations))
+        self.assertEquals(mismatchpairs_df.size, iterations)
 
     def test_create_sim_masking(self):
 

From a2b4cd9be4b15f9c6bae2297b6f8c4e3cf17a8a0 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 30 Jul 2024 17:30:53 -0700
Subject: [PATCH 39/50] linty

---
 q2_fmt/tests/test_engraftment.py | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 2f9de57..5596dd2 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1403,12 +1403,6 @@ def test_mask_recipient(self):
         np.testing.assert_array_equal(recip_mask, exp_r_mask)
 
     def test_simulate_uniform_distro(self):
-        recip_df = pd.DataFrame({
-            'id': ['sample1', 'sample2', 'sample3'],
-            'Feature1': [1, 0, 0],
-            'Feature2': [0, 1, 0],
-            'Feature3': [0, 0, 1]}).set_index('id')
-
         mismatch_peds = [0, 0, 0, 0, 0, 0]
 
         iterations = 999

From ae8b429630fd9aae806297dbfbb3dbbf7e61a449 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Tue, 30 Jul 2024 18:48:03 -0700
Subject: [PATCH 40/50] more tests

---
 q2_fmt/tests/test_engraftment.py | 73 +++++++++++++++++++++++++++++++-
 1 file changed, 72 insertions(+), 1 deletion(-)

diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 5596dd2..ceeafc8 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1402,13 +1402,33 @@ def test_mask_recipient(self):
         exp_r_mask = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
         np.testing.assert_array_equal(recip_mask, exp_r_mask)
 
+    def test_mask_recipient_donor_one(self):
+        donor_mask = np.array([[0, 1, 0], [0, 1, 0], [0, 1, 0]])
+        recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+        recip_mask = _mask_recipient(donor_mask, recip_df)
+        exp_r_mask = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]
+        np.testing.assert_array_equal(recip_mask, exp_r_mask)
+
     def test_simulate_uniform_distro(self):
         mismatch_peds = [0, 0, 0, 0, 0, 0]
 
         iterations = 999
 
         mismatchpairs_df = _simulate_uniform_distro(mismatch_peds,
-                                                    999)
+                                                    iterations)
+        self.assertEquals(mismatchpairs_df.size, iterations)
+
+    def test_one_iter_simulate_uniform_distro(self):
+        mismatch_peds = [0, 0, 0, 0, 0, 0]
+
+        iterations = 1
+
+        mismatchpairs_df = _simulate_uniform_distro(mismatch_peds,
+                                                    iterations)
         self.assertEquals(mismatchpairs_df.size, iterations)
 
     def test_create_sim_masking(self):
@@ -1438,6 +1458,30 @@ def test_create_sim_masking(self):
                                          reference_column='Ref')
         np.testing.assert_array_equal(donor_mask, exp_mask)
 
+    def test_create_one_donor_sim_masking(self):
+
+        mismatched_df = pd.DataFrame({'id': ["sample1",
+                                             "sample2",
+                                             "sample3"],
+                                      "Ref": ["donor2",
+                                              "donor2",
+                                              "donor2"]}
+                                     ).set_index('id')
+
+        donor_df = pd.DataFrame({
+            'id': ['donor2'],
+            'Feature1': [1],
+            'Feature2': [0],
+            'Feature3': [0]}).set_index('id')
+
+        exp_mask = [[1, 0, 0],
+                    [1, 0, 0],
+                    [1, 0, 0]]
+
+        donor_mask = _create_sim_masking(mismatched_df, donor_df,
+                                         reference_column='Ref')
+        np.testing.assert_array_equal(donor_mask, exp_mask)
+
     def test_create_duplicated_table(self):
         recip_df = pd.DataFrame({
             'id': ['sample1', 'sample2', 'sample3'],
@@ -1466,6 +1510,33 @@ def test_create_duplicated_table(self):
         pd.testing.assert_frame_equal(duplicated_recip_table,
                                       exp_d_r_table)
 
+    def test_create_no_duplicated_table(self):
+        recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+        mismatched_df = pd.DataFrame({'id': ["sample1",
+                                             "sample2",
+                                             "sample3"],
+                                      "Ref": ["donor2",
+                                              "donor2",
+                                              "donor2"]}
+                                     ).set_index('id')
+
+        duplicated_recip_table = _create_duplicated_recip_table(mismatched_df,
+                                                                recip_df)
+
+        exp_d_r_table = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+        pd.testing.assert_frame_equal(duplicated_recip_table,
+                                      exp_d_r_table)
+
     def test_per_subject_stats_labels(self):
         mismatched_peds = [0, 0, 0, 0]
         actual_temp = pd.Series(data=[1, 1, 1, 1],

From 22b61eb906abb9330c2d6b0e34ca0a6c14af71d4 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 31 Jul 2024 15:41:31 -0700
Subject: [PATCH 41/50] Apply suggestions from code review

Co-authored-by: Greg Caporaso <192372+gregcaporaso@users.noreply.github.com>
---
 q2_fmt/_peds.py        | 12 ++++++------
 q2_fmt/plugin_setup.py | 12 ++++++------
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index fca8bb8..c099722 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -684,7 +684,7 @@ def _create_duplicated_recip_table(mismatched_df, recip_df):
 
 
 def _create_sim_masking(mismatched_df, donor_df, reference_column):
-    """Create a donor mask to mask recipient feature that aren't in the donor.
+    """Create a donor mask to mask recipient features that aren't in the donor.
 
     Creates a Donor Numpy array that duplicates donor samples to match up
     with duplicated_table. This will mask recipient feature that aren't
@@ -743,7 +743,7 @@ def _simulate_uniform_distro(mismatched_peds, iterations):
     """Randomly samples the mismatched PEDS values.
 
     Creates a uniform distribution of mismatched PEDS values by randomly
-    sampling N (iterations) times
+    sampling `iterations` times with replacement.
 
     Parameters
     ----------
@@ -780,10 +780,10 @@ def _per_subject_stats(mismatched_peds, actual_temp,
     """Creates per subject PEDS stats
 
     Creates per subject PEDS stats by sampling mismatch PEDS values 'iteration'
-    number of time. The list of mismatched PEDS values is compared to the
+    number of times. The list of mismatched PEDS values is compared to the
     actual PEDS value. The Test Statistic for this is the iterations - the
     number of mismatched PEDS value that are greater than the actual p-value.
-    p-values are corrected by adding 1 to each side of the faction
+    p-values are corrected by adding 1 to each side of the fraction
     because Monte Carlo Simulation is able to produce an invalid p-value of 0.
 
     Parameters
@@ -797,7 +797,7 @@ def _per_subject_stats(mismatched_peds, actual_temp,
         Number of iterations to run simulations (Number of times to randomly
         sample mismatched_peds)
     mismatched_peds_n: int
-        Number of mismatched donor pairs are possible.
+        Number of possible mismatched donor pairs.
 
     Returns
     -------
@@ -880,7 +880,7 @@ def _per_subject_stats(mismatched_peds, actual_temp,
         dict(title='Iteration',
              description='Number of iterations that agree with the ALTERNATIVE'
                          'hypothesis'))
-    per_sub_stats['p-value'].attrs.update(dict(title='one-tail p-value',
+    per_sub_stats['p-value'].attrs.update(dict(title='one-tailed',
                                           description='one-tail p-value'))
     per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index f61dd6f..f08ce33 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -63,11 +63,11 @@
 drop_incomplete_subjects = ('Filter out subjects that do not have a sample at'
                             ' every timepoint. Default behavior is to raise an'
                             ' error if any subject is missing a timepoint.')
-drop_incomplete_timepoint = ('Filter out a `problematic` timepoint that is'
-                             ' causing subject to be dropped. Default behavior'
+drop_incomplete_timepoint = ('Filter out specified timepoint. This is useful for removing frequently missing timepoints which cause'
+                             ' many subjects to be dropped. Default behavior'
                              ' is to raise an error if any subject is missing'
                              ' a timepoint.')
-level_delimiter = 'delimiter to split taxonomic string on'
+level_delimiter = 'delimiter to split taxonomic label on'
 control_column = ('The column within `metadata` that contains any relevant'
                   ' control group IDs. Actual treatment samples should not'
                   ' contain any value within this column.')
@@ -83,9 +83,9 @@
                      ' basis.')
 global_stats = ('Table describing significance of PEDS scores across all'
                 ' subjects.')
-peds_table = 'A feature table to calculate PEDS on.'
+peds_table = 'The feature table to calculate PEDS on.'
 peds_dists = ('The distributions for the PEDS measure, grouped by the selected'
-              ' `time_column`. also contains the Numerator and Denominator for'
+              ' `time_column`. Also contains the numerator and denominator for'
               ' PEDS calulations. May also contain subject IDs, if'
               ' `subject_column` is  provided in `metadata`.')
 
@@ -359,7 +359,7 @@
                 ' between donors and recipients, to test whether the PEDS'
                 ' score between a recipient and their actual donor is'
                 ' significantly higher than PEDS scores between other'
-                ' recipients and a random donor. This is intended to only work'
+                ' recipients paired with random donors. This is intended to only work'
                 ' in studies where there are distinct donors, and will yield'
                 ' insignificant results if there are too few donors.'
                 ' This method attempts to quantify the extent to which'

From caafb9aa6e9e2238b7358af633504ad216447448 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 31 Jul 2024 18:47:04 -0700
Subject: [PATCH 42/50] review

---
 q2_fmt/_peds.py                  | 274 ++++++++++++++++++-------------
 q2_fmt/citations.bib             |  26 +++
 q2_fmt/plugin_setup.py           |  31 ++--
 q2_fmt/tests/test_engraftment.py |   9 +-
 4 files changed, 202 insertions(+), 138 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index c099722..2e0fc29 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -493,7 +493,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
                     filter_missing_references: bool = False,
                     drop_incomplete_subjects: bool = False,
                     drop_incomplete_timepoint: list = None,
-                    iterations: int = 999) -> (pd.DataFrame, pd.DataFrame):
+                    num_iterations: int = 999) -> (pd.DataFrame, pd.DataFrame):
 
     metadata_df = metadata.to_dataframe()
     reference_series = _check_reference_column(metadata_df,
@@ -528,7 +528,7 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
            filter_missing_references=filter_missing_references,
            drop_incomplete_subjects=drop_incomplete_subjects,
            drop_incomplete_timepoint=drop_incomplete_timepoint).set_index("id")
-    actual_temp = peds["measure"]
+    actual_peds = peds["measure"]
 
     # Mismatch simulation:
     table = table > 0
@@ -544,18 +544,16 @@ def peds_simulation(table: pd.DataFrame, metadata: qiime2.Metadata,
     donor_mask = _create_sim_masking(mismatched_df, donor_df, reference_column)
     recip_mask = _mask_recipient(donor_mask, duplicated_recip_table)
     # Numerator for PEDS Calc. (Number of Donor features in the Recipient)
-    num_sum = np.sum(recip_mask.values, axis=1)
+    num_engrafted_donor_features = np.sum(recip_mask.values, axis=1)
     # Denominator for PEDS Calc. (Number of unique features in the Donor)
-    donor_sum = np.sum(donor_mask, axis=1)
+    num_donor_features = np.sum(donor_mask, axis=1)
     # This ignores warnings that come from dividing by 0.
     # mismatched_peds will be Nan if the denominator is 0 and thats reasonable.
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
-        mismatched_peds = num_sum/donor_sum
-    mismatched_pairs_n = len(duplicated_recip_table.index)
+        mismatched_peds = num_engrafted_donor_features/num_donor_features
     per_sub_stats = _per_subject_stats(mismatched_peds,
-                                       actual_temp, iterations,
-                                       mismatched_pairs_n)
+                                       actual_peds, num_iterations)
     global_stats = _global_stats(per_sub_stats['p-value'])
     return per_sub_stats, global_stats
 
@@ -589,41 +587,46 @@ def _create_mismatched_pairs(recip_df, metadata, used_references,
 
     Examples
     --------
-    recip_df = pd.DataFrame({
-        'id': ['sample1', 'sample2', 'sample3'],
-        'Feature1': [1, 0, 0],
-        'Feature2': [0, 1, 0],
-        'Feature3': [0, 0, 1]}).set_index('id')
-
-    metadata_df = pd.DataFrame({
-        'id': ['sample1', 'sample2', 'sample3','donor1', 'donor2', 'donor3'],
-        'Ref': ['donor1', 'donor2', 'donor3', np.nan, np.nan, np.nan],
-        'subject': ['sub1', 'sub2', 'sub3', np.nan, np.nan, np.nan],
-        'group': [1, 1, 1, np.nan, np.nan, np.nan],
-        "Location": [np.nan, np.nan,np.nan, 'test', 'test','test']}
-        ).set_index('id')
-
-    used_references = pd.Series(data=['donor1', 'donor2', 'donor3'],
-                                    index=['sample1', 'sample2', 'sample3'],
+    >>> recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
+
+    >>> metadata_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3',
+                   'donor1', 'donor2', 'donor3'],
+            'Ref': ['donor1', 'donor2', 'donor3',
+                    np.nan, np.nan, np.nan],
+            'subject': ['sub1', 'sub2', 'sub3',
+                        np.nan, np.nan, np.nan],
+            'group': [1, 1, 1,
+                      np.nan, np.nan, np.nan],
+            'Location': [np.nan, np.nan,np.nan,
+                         'test', 'test','test']}).set_index('id')
+
+    >>> used_references = pd.Series(data=['donor1', 'donor2', 'donor3'],
+                                    index=pd.Index(['sample1', 'sample2',
+                                                    'sample3'], name='id'),
                                     name='Ref')
-    used_references.index.name = "id"
 
-    reference_column = "Ref"
+    >>> reference_column = "Ref"
 
-    _create_mismatched_pairs(recip_df, metadata_df, used_references,
-                             reference_column)
+    >>> _create_mismatched_pairs(recip_df, metadata_df, used_references,
+                                 reference_column)
     pd.DataFrame({'id': ["sample1", "sample1", "sample2", "sample2",
                          "sample3", "sample3"],
                   "Ref": ["donor2", "donor3", "donor1", "donor3",
                           "donor1", "donor2"]}).set_index('id')
     """
-    mismatched_pairs = []
-    for x in itertools.product(recip_df.index,
-                               metadata[reference_column].dropna().unique()):
-        mismatched_pairs.append(x)
-    mismatched_t = list(zip(used_references.index, used_references))
-    # Removes matched donor and recipient pairs
-    filtered = [pair for pair in mismatched_pairs if pair not in mismatched_t]
+    matched_pairs = list(zip(used_references.index, used_references))
+    donors = metadata[reference_column].dropna().unique()
+    # Generates all donor recipient pairs and then removes
+    # matched donor and recipient pairs
+    filtered = \
+        [pair
+         for pair in itertools.product(recip_df.index,
+                                       donors)if pair not in matched_pairs]
     idx, values = zip(*filtered)
     mismatched_df = pd.DataFrame({"id": idx,
                                   reference_column: values}).set_index("id")
@@ -655,20 +658,19 @@ def _create_duplicated_recip_table(mismatched_df, recip_df):
 
     Examples
     --------
-    mismatched_df = pd.DataFrame({
+    >>> mismatched_df = pd.DataFrame({
             'id': ["sample1", "sample1", "sample2", "sample2",
                    "sample3", "sample3"],
             "Ref": ["donor2", "donor3", "donor1", "donor3",
-                    "donor1", "donor2"]}
-                                ).set_index('id')
+                    "donor1", "donor2"]}).set_index('id')
 
-    recip_df = pd.DataFrame({
-        'id': ['sample1', 'sample2', 'sample3'],
-        'Feature1': [1, 0, 0],
-        'Feature2': [0, 1, 0],
-        'Feature3': [0, 0, 1]}).set_index('id')
+    >>> recip_df = pd.DataFrame({
+            'id': ['sample1', 'sample2', 'sample3'],
+            'Feature1': [1, 0, 0],
+            'Feature2': [0, 1, 0],
+            'Feature3': [0, 0, 1]}).set_index('id')
 
-    _create_duplicated_recip_table(mismatched_df, recip_df)
+    >>> _create_duplicated_recip_table(mismatched_df, recip_df)
 
     pd.DataFrame({
             'id': ["sample1", "sample1", "sample2", "sample2",
@@ -709,19 +711,19 @@ def _create_sim_masking(mismatched_df, donor_df, reference_column):
 
     Examples
     --------
-    mismatched_df = pd.DataFrame({
-            'id': ["sample1", "sample1", "sample2", "sample2",
-                   "sample3", "sample3"],
-            "Ref": ["donor2", "donor3", "donor1", "donor3",
-                    "donor1", "donor2"]}
+    >>> mismatched_df = pd.DataFrame({
+                'id': ["sample1", "sample1", "sample2", "sample2",
+                       "sample3", "sample3"],
+                "Ref": ["donor2", "donor3", "donor1", "donor3",
+                        "donor1", "donor2"]}
                                 ).set_index('id')
-    donor_df = pd.DataFrame({
-            'id': ['donor1', 'donor2', 'donor3'],
-            'Feature1': [1, 0, 0],
-            'Feature2': [0, 1, 0],
-            'Feature3': [0, 0, 1]}).set_index('id')
+    >>> donor_df = pd.DataFrame({
+              'id': ['donor1', 'donor2', 'donor3'],
+              'Feature1': [1, 0, 0],
+              'Feature2': [0, 1, 0],
+              'Feature3': [0, 0, 1]}).set_index('id')
 
-    _create_sim_masking(mismatched_df, donor_df, reference_column)
+    >>> _create_sim_masking(mismatched_df, donor_df, reference_column)
 
     ndarray[[0, 1, 0],
             [0, 0, 1],
@@ -739,65 +741,109 @@ def _create_sim_masking(mismatched_df, donor_df, reference_column):
     return donor_mask
 
 
-def _simulate_uniform_distro(mismatched_peds, iterations):
+def _simulate_uniform_distro(mismatched_peds, k):
     """Randomly samples the mismatched PEDS values.
 
     Creates a uniform distribution of mismatched PEDS values by randomly
-    sampling `iterations` times with replacement.
+    sampling `num_iterations` times with replacement.
 
     Parameters
     ----------
     mismatched_peds: pd.DataFrame
         A Dataframe that contains all mismatched PEDS values.
-    iterations: int
-        Number of iterations to run simulations (Number of times to randomly
-        sample mismatched_peds)
+    k: int
+        Number of iterations(`k`) to run simulations (Number of times to
+        randomly sample mismatched_peds)
 
     Returns
     -------
     peds_iters: pd.Series
-        a pd.Series with all iterations of mismatched PEDS Values. This will
-        later be compared to an actual PEDS value.
+        a pd.Series with all num_iterations of mismatched PEDS Values. This
+        will later be compared to an actual PEDS value.
 
     Examples
     --------
-    mismatched_peds = [0, 0, 0,0 ]
-    iterations = 10
+    >>> mismatched_peds = [0, 0, 0,0 ]
+    >>> num_iterations = 10
 
-    _simulate_uniform_distro(mismatched_peds, iterations)
+    >>> _simulate_uniform_distro(mismatched_peds, num_iterations)
 
     pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
                index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
     """
-    peds_iters = random.choices(mismatched_peds, k=iterations)
+    peds_iters = random.choices(mismatched_peds, k=k)
     # Tranforming to series for easy series math in _per_subject_stats()
     peds_iters = pd.Series(peds_iters)
     return peds_iters
 
 
-def _per_subject_stats(mismatched_peds, actual_temp,
-                       iterations, mismatched_pairs_n):
+def peds_sim_stats(value, peds_iters, num_iterations):
+    """Calculates test statistics, and p-value.
+
+    Calculates tests statistics (`count_gte` and ` count_less`) and p-value
+    for PEDS Monte Carlo Simulation
+
+    Parameters
+    ----------
+    value: float
+        A actual PEDS value to compare against
+    peds_iters: pd.Series
+        a pd.Series with all num_iterations of mismatched PEDS Values. This
+        will be compared to the actual PEDS value.
+    num_iterations: int
+        Number of iterations to run simulations (Number of times to
+        randomly sample mismatched_peds)
+
+    Returns
+    -------
+    count_gte: int
+        Count of mismatched PEDS values that were greater than the actual PEDS
+        value
+    count_less
+        Count of mismatched PEDS values that were less than the actual PEDS
+        value. This is calculated by ``num_interations - count_gte``
+    per_subject_p
+        The p-value associated with the above test stats.
+
+    Examples
+    --------
+    >>> pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                   index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    >>> num_iterations = 10
+    >>> value = 1
+
+    >>> _simulate_uniform_distro(mismatched_peds, num_iterations)
+
+    pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
+               index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    """
+    gte_series = peds_iters >= value
+    count_gte = gte_series.sum()
+    count_less = num_iterations-count_gte
+    # adding 1 here because you can mathmatically can get p-value of 0 from a
+    # Monte Carlo Simulation
+    per_subject_p = (count_gte + 1)/(num_iterations+1)
+    return count_gte, count_less, per_subject_p
+
+
+def _per_subject_stats(mismatched_peds, actual_peds,
+                       num_iterations):
     """Creates per subject PEDS stats
 
     Creates per subject PEDS stats by sampling mismatch PEDS values 'iteration'
     number of times. The list of mismatched PEDS values is compared to the
-    actual PEDS value. The Test Statistic for this is the iterations - the
-    number of mismatched PEDS value that are greater than the actual p-value.
-    p-values are corrected by adding 1 to each side of the fraction
-    because Monte Carlo Simulation is able to produce an invalid p-value of 0.
+    actual PEDS value.
 
     Parameters
     ----------
     mismatched_peds: pd.DataFrame
         A Dataframe that contains all mismatched PEDS values.
-    actual_temp: pd.Series
+    actual_peds: pd.Series
         A Series containing Sample IDs as the index and actual PEDS values for
         the sample as the value.
-    iterations: int
-        Number of iterations to run simulations (Number of times to randomly
-        sample mismatched_peds)
-    mismatched_peds_n: int
-        Number of possible mismatched donor pairs.
+    num_iterations: int
+        Number of num_iterations to run simulations (Number of times to
+        randomly sample mismatched_peds)
 
     Returns
     -------
@@ -808,20 +854,20 @@ def _per_subject_stats(mismatched_peds, actual_temp,
 
     Examples
     --------
-    mismatched_peds = [0, 0, 0, 0]
-    actual_temp = pd.Series([data = [1, 1, 1, 1],
-               index = [sample1, sample2, sample3, sample4]])
-    iterations = 10
-    mismatched_pairs_n = 4
+    >>> mismatched_peds = [0, 0, 0, 0]
+    >>> actual_peds = pd.Series([data = [1, 1, 1, 1],
+                                 index = ["sample1", "sample2",
+                                          "sample3", "sample4"]])
+    >>> num_iterations = 10
 
-    _per_subject_stats(mismatched_peds, actual_temp,
-                       iterations, mismatched_pairs_n):
+    >>> _per_subject_stats(mismatched_peds, actual_peds,
+                       num_iterations):
 
-    pd.Dataframe({"A:group": [sample1, sample2, sample3, sample4],
+    pd.Dataframe({"A:group": ["sample1", "sample2", "sample3", "sample4"],
                    "A:n": [1, 1, 1, 1],
                    "A:measure" : [1, 1, 1, 1],
-                   "B:group" : [shuffled_recipients, shuffled_recipients,
-                                shuffled_recipients, shuffled_recipients],
+                   "B:group" : ["shuffled recipients", "shuffled recipients",
+                                "shuffled recipients", "shuffled recipients"],
                    "B:n": [4, 4, 4, 4],
                    "B:measure": [0, 0, 0, 0],
                    "n": [10, 10, 10, 10],
@@ -829,29 +875,27 @@ def _per_subject_stats(mismatched_peds, actual_temp,
                    "p-value": [0.001, 0.001, 0.001, 0.001],
                    "q-value": [0.004, 0.002, 0.00133, 0.001]})
     """
-    disagree_df = pd.DataFrame([])
-    peds_iter_mean = []
-    for value in actual_temp:
-        peds_iters = _simulate_uniform_distro(mismatched_peds, iterations)
-        peds_iter_mean.append(peds_iters.mean())
-        disagree_recipient = peds_iters >= value
-        # add disagree recipient as a row in the disagree df
-        disagree_df = pd.concat([disagree_df, disagree_recipient.to_frame().T])
-    disagree_series = disagree_df.sum(axis=1)
-    agree_series = iterations-disagree_series
-    # adding 1 here because you can mathmatically can get p-value of 0 from a
-    # Monte Carlo Simulation
-    per_subject_p = (disagree_series + 1)/(iterations+1)
+    peds_iters_means = []
+    count_less_list = []
+    per_subject_p_list = []
+    for value in actual_peds:
+        peds_iters = _simulate_uniform_distro(mismatched_peds, num_iterations)
+        peds_iters_means.append(peds_iters.mean())
+        _, count_less, per_subject_p = peds_sim_stats(value, peds_iters,
+                                                      num_iterations)
+        count_less_list.append(count_less)
+        per_subject_p_list.append(per_subject_p)
+
     per_subject_q = false_discovery_control(ps=per_subject_p, method='bh')
-    per_sub_stats = pd.DataFrame({'A:group': actual_temp.index,
+    per_sub_stats = pd.DataFrame({'A:group': actual_peds.index,
                                  'A:n': 1,
-                                  'A:measure': actual_temp.values,
-                                  'B:group': "shuffled_recipients",
-                                  'B:n': mismatched_pairs_n,
-                                  'B:measure': peds_iter_mean,
-                                  'n': iterations,
-                                  'test-statistic': agree_series,
-                                  'p-value': per_subject_p,
+                                  'A:measure': actual_peds.values,
+                                  'B:group': "shuffled recipients",
+                                  'B:n': mismatched_peds.shape[0],
+                                  'B:measure': peds_iters_means,
+                                  'n': num_iterations,
+                                  'test-statistic': count_less_list,
+                                  'p-value': per_subject_p_list,
                                   'q-value': per_subject_q})
     per_sub_stats['A:group'].attrs.update({'title': 'actual_values',
                                           'description': 'PEDS values'
@@ -878,13 +922,13 @@ def _per_subject_stats(mismatched_peds, actual_temp,
                                     'comparisons'))
     per_sub_stats['test-statistic'].attrs.update(
         dict(title='Iteration',
-             description='Number of iterations that agree with the ALTERNATIVE'
-                         'hypothesis'))
+             description='Number of num_iterations that agree with the'
+                         ' ALTERNATIVE hypothesis'))
     per_sub_stats['p-value'].attrs.update(dict(title='one-tailed',
                                           description='one-tail p-value'))
     per_sub_stats['q-value'].attrs.update(
         dict(title='Benjamini–Hochberg', description='FDR corrections using'
-             'Benjamini–Hochberg procedure'))
+             ' Benjamini–Hochberg procedure'))
     return per_sub_stats
 
 
@@ -908,11 +952,11 @@ def _global_stats(p_series):
 
     Examples
     --------
-    p_series = pd.Series([data = [0.001, 0.001, 0.001, 0.001])
+    >>> p_series = pd.Series([data = [0.001, 0.001, 0.001, 0.001])
 
-    _global_stats(p_series)
+    >>> _global_stats(p_series)
 
-    pd.Dataframe({"measure": [p-values],
+    pd.Dataframe({"measure": ["p-values"],
                    "n": [4],
                    "test-statistic": [20.93],
                    "p-value": [0.001]
diff --git a/q2_fmt/citations.bib b/q2_fmt/citations.bib
index 1a6a85b..1db8063 100644
--- a/q2_fmt/citations.bib
+++ b/q2_fmt/citations.bib
@@ -15,3 +15,29 @@ @article{aggarwala_precise_2021
 	date = {2021-10},
 	langid = {english},
 }
+
+@book{stouffer_1949_american,
+  author    = {Stouffer, S. A. and Suchman, E. A. and Devinney, L. C. and Star, S. A. and Williams, R. M.},
+  title     = {The American Soldier: Adjustment during Army Life},
+  series    = {Studies in Social Psychology in World War II},
+  volume    = {I},
+  publisher = {Princeton University Press},
+  address   = {Princeton, NJ},
+  year      = {1949}
+}
+
+@article{Benjamini_fdr_1995,
+ ISSN = {00359246},
+ URL = {http://www.jstor.org/stable/2346101},
+ abstract = {The common approach to the multiplicity problem calls for controlling the familywise error rate (FWER). This approach, though, has faults, and we point out a few. A different approach to problems of multiple significance testing is presented. It calls for controlling the expected proportion of falsely rejected hypotheses-the false discovery rate. This error rate is equivalent to the FWER when all hypotheses are true but is smaller otherwise. Therefore, in problems where the control of the false discovery rate rather than that of the FWER is desired, there is potential for a gain in power. A simple sequential Bonferroni-type procedure is proved to control the false discovery rate for independent test statistics, and a simulation study shows that the gain in power is substantial. The use of the new procedure and the appropriateness of the criterion are illustrated with examples.},
+ author = {Yoav Benjamini and Yosef Hochberg},
+ journal = {Journal of the Royal Statistical Society. Series B (Methodological)},
+ number = {1},
+ pages = {289--300},
+ publisher = {[Royal Statistical Society, Oxford University Press]},
+ title = {Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing},
+ urldate = {2024-07-31},
+ volume = {57},
+ year = {1995}
+}
+
diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index f08ce33..6c3a39b 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -63,10 +63,11 @@
 drop_incomplete_subjects = ('Filter out subjects that do not have a sample at'
                             ' every timepoint. Default behavior is to raise an'
                             ' error if any subject is missing a timepoint.')
-drop_incomplete_timepoint = ('Filter out specified timepoint. This is useful for removing frequently missing timepoints which cause'
-                             ' many subjects to be dropped. Default behavior'
-                             ' is to raise an error if any subject is missing'
-                             ' a timepoint.')
+drop_incomplete_timepoint = ('Filter out specified timepoints. This is useful'
+                             ' for removing frequently missing timepoints'
+                             ' which cause many subjects to be dropped.'
+                             ' Default behavior is to raise an error if any'
+                             ' subject is missing a timepoint.')
 level_delimiter = 'delimiter to split taxonomic label on'
 control_column = ('The column within `metadata` that contains any relevant'
                   ' control group IDs. Actual treatment samples should not'
@@ -336,7 +337,7 @@
                 'filter_missing_references': Bool,
                 'drop_incomplete_subjects': Bool,
                 'drop_incomplete_timepoint': List[Str],
-                'iterations': Int % Range(99, None)},
+                'num_iterations': Int % Range(99, None)},
     outputs=[('per_subject_stats', StatsTable[Pairwise]),
              ('global_stats', StatsTable[Pairwise])],
     parameter_descriptions={
@@ -347,8 +348,8 @@
         'filter_missing_references': filter_missing_references,
         'drop_incomplete_subjects': drop_incomplete_subjects,
         'drop_incomplete_timepoint': drop_incomplete_timepoint,
-        'iterations': 'The number of iterations to run the Monte Carlo'
-                      ' simulation on'
+        'num_iterations': 'The number of iterations to run the Monte Carlo'
+                          ' simulation on'
     },
     output_descriptions={
         'per_subject_stats': per_subject_stats,
@@ -359,14 +360,14 @@
                 ' between donors and recipients, to test whether the PEDS'
                 ' score between a recipient and their actual donor is'
                 ' significantly higher than PEDS scores between other'
-                ' recipients paired with random donors. This is intended to only work'
-                ' in studies where there are distinct donors, and will yield'
-                ' insignificant results if there are too few donors.'
-                ' This method attempts to quantify the extent to which'
-                ' indicator features that are unique to a given donor are'
-                ' transferred to their recipients, as opposed to features that'
-                ' are not indicative of any specific donor.',
-    citations=[citations['aggarwala_precise_2021']],
+                ' recipients paired with random donors. This is intended to'
+                ' only work in studies where there are distinct donors,'
+                ' and will yield insignificant results if there are too'
+                ' few donors. This method attempts to quantify the extent to'
+                ' which indicator features that are unique to a given donor'
+                ' are transferred to their recipients, as opposed to features'
+                ' that are not indicative of any specific donor.',
+    citations=[citations['aggarwala_precise_2021', 'stouffer_1949_american']],
     examples={
         'peds_methods': ex.simulation_peds_method
     }
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index ceeafc8..7be72de 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1527,15 +1527,8 @@ def test_create_no_duplicated_table(self):
 
         duplicated_recip_table = _create_duplicated_recip_table(mismatched_df,
                                                                 recip_df)
-
-        exp_d_r_table = pd.DataFrame({
-            'id': ['sample1', 'sample2', 'sample3'],
-            'Feature1': [1, 0, 0],
-            'Feature2': [0, 1, 0],
-            'Feature3': [0, 0, 1]}).set_index('id')
-
         pd.testing.assert_frame_equal(duplicated_recip_table,
-                                      exp_d_r_table)
+                                      recip_df)
 
     def test_per_subject_stats_labels(self):
         mismatched_peds = [0, 0, 0, 0]

From e3d3f04b29565fa4156155a71a30d521ef2bf4f7 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 31 Jul 2024 18:53:17 -0700
Subject: [PATCH 43/50] fix doc string

---
 q2_fmt/_peds.py | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 2e0fc29..5c95659 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -777,7 +777,7 @@ def _simulate_uniform_distro(mismatched_peds, k):
     return peds_iters
 
 
-def peds_sim_stats(value, peds_iters, num_iterations):
+def _peds_sim_stats(value, peds_iters, num_iterations):
     """Calculates test statistics, and p-value.
 
     Calculates tests statistics (`count_gte` and ` count_less`) and p-value
@@ -799,23 +799,25 @@ def peds_sim_stats(value, peds_iters, num_iterations):
     count_gte: int
         Count of mismatched PEDS values that were greater than the actual PEDS
         value
-    count_less
+    count_less: int
         Count of mismatched PEDS values that were less than the actual PEDS
         value. This is calculated by ``num_interations - count_gte``
-    per_subject_p
+    per_subject_p: float
         The p-value associated with the above test stats.
 
     Examples
     --------
-    >>> pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
-                   index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
-    >>> num_iterations = 10
     >>> value = 1
+    >>> peds_iters = pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
+                                index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    >>> num_iterations = 10
 
-    >>> _simulate_uniform_distro(mismatched_peds, num_iterations)
+    >>> _peds_sim_stats()
 
-    pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
-               index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    count_gte = 0
+    count_less = 10
+    per_subject_p = (1/11) (Note: the 10 iterations is not enough to get a 
+                            significant p-value)
     """
     gte_series = peds_iters >= value
     count_gte = gte_series.sum()
@@ -881,8 +883,8 @@ def _per_subject_stats(mismatched_peds, actual_peds,
     for value in actual_peds:
         peds_iters = _simulate_uniform_distro(mismatched_peds, num_iterations)
         peds_iters_means.append(peds_iters.mean())
-        _, count_less, per_subject_p = peds_sim_stats(value, peds_iters,
-                                                      num_iterations)
+        _, count_less, per_subject_p = _peds_sim_stats(value, peds_iters,
+                                                       num_iterations)
         count_less_list.append(count_less)
         per_subject_p_list.append(per_subject_p)
 

From 43ff16775d64a461ebaf8b23330a8ca24eef9ed6 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 31 Jul 2024 19:03:42 -0700
Subject: [PATCH 44/50] fix cites

---
 q2_fmt/plugin_setup.py           | 4 +++-
 q2_fmt/tests/test_engraftment.py | 2 +-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 6c3a39b..e7365be 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -367,7 +367,9 @@
                 ' which indicator features that are unique to a given donor'
                 ' are transferred to their recipients, as opposed to features'
                 ' that are not indicative of any specific donor.',
-    citations=[citations['aggarwala_precise_2021', 'stouffer_1949_american']],
+    citations=[citations['aggarwala_precise_2021'],
+               citations['stouffer_1949_american'],
+               citations['Benjamini_fdr_1995']],
     examples={
         'peds_methods': ex.simulation_peds_method
     }
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 7be72de..b947357 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -23,7 +23,7 @@
                           peds_simulation, _create_mismatched_pairs,
                           _simulate_uniform_distro, _create_sim_masking,
                           _mask_recipient, _create_duplicated_recip_table,
-                          _per_subject_stats, _global_stats)
+                          _per_subject_stats, _global_stats, )
 
 
 class TestBase(TestPluginBase):

From 7931b0e7c8872bd72db8289458f08b78dbe00e32 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 31 Jul 2024 19:56:13 -0700
Subject: [PATCH 45/50] more tests!

---
 q2_fmt/_peds.py                  |  20 ++---
 q2_fmt/tests/test_engraftment.py | 132 ++++++++++++++++++++++++++++---
 2 files changed, 131 insertions(+), 21 deletions(-)

diff --git a/q2_fmt/_peds.py b/q2_fmt/_peds.py
index 5c95659..ccec7e0 100644
--- a/q2_fmt/_peds.py
+++ b/q2_fmt/_peds.py
@@ -749,8 +749,8 @@ def _simulate_uniform_distro(mismatched_peds, k):
 
     Parameters
     ----------
-    mismatched_peds: pd.DataFrame
-        A Dataframe that contains all mismatched PEDS values.
+    mismatched_peds: list
+        A list that contains all mismatched PEDS values.
     k: int
         Number of iterations(`k`) to run simulations (Number of times to
         randomly sample mismatched_peds)
@@ -768,8 +768,8 @@ def _simulate_uniform_distro(mismatched_peds, k):
 
     >>> _simulate_uniform_distro(mismatched_peds, num_iterations)
 
-    pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
-               index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    pd.Series(data = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+              index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
     """
     peds_iters = random.choices(mismatched_peds, k=k)
     # Tranforming to series for easy series math in _per_subject_stats()
@@ -808,15 +808,15 @@ def _peds_sim_stats(value, peds_iters, num_iterations):
     Examples
     --------
     >>> value = 1
-    >>> peds_iters = pd.Series([data = [0, 0, 0, 0, 0, 0, 0, 0, 0],
-                                index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])
+    >>> peds_iters = pd.Series(data = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                               index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
     >>> num_iterations = 10
 
     >>> _peds_sim_stats()
 
     count_gte = 0
     count_less = 10
-    per_subject_p = (1/11) (Note: the 10 iterations is not enough to get a 
+    per_subject_p = (1/11) (Note: the 10 iterations is not enough to get a
                             significant p-value)
     """
     gte_series = peds_iters >= value
@@ -838,8 +838,8 @@ def _per_subject_stats(mismatched_peds, actual_peds,
 
     Parameters
     ----------
-    mismatched_peds: pd.DataFrame
-        A Dataframe that contains all mismatched PEDS values.
+    mismatched_peds: list
+        A list that contains all mismatched PEDS values.
     actual_peds: pd.Series
         A Series containing Sample IDs as the index and actual PEDS values for
         the sample as the value.
@@ -893,7 +893,7 @@ def _per_subject_stats(mismatched_peds, actual_peds,
                                  'A:n': 1,
                                   'A:measure': actual_peds.values,
                                   'B:group': "shuffled recipients",
-                                  'B:n': mismatched_peds.shape[0],
+                                  'B:n': len(mismatched_peds),
                                   'B:measure': peds_iters_means,
                                   'n': num_iterations,
                                   'test-statistic': count_less_list,
diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index b947357..7e1bf60 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -23,7 +23,7 @@
                           peds_simulation, _create_mismatched_pairs,
                           _simulate_uniform_distro, _create_sim_masking,
                           _mask_recipient, _create_duplicated_recip_table,
-                          _per_subject_stats, _global_stats, )
+                          _per_subject_stats, _global_stats, _peds_sim_stats)
 
 
 class TestBase(TestPluginBase):
@@ -1290,7 +1290,7 @@ def test_high_donor_overlap(self):
                                    time_column="group",
                                    reference_column="Ref",
                                    subject_column="subject",
-                                   iterations=999)
+                                   num_iterations=999)
         real_median = np.median(stats["A:measure"].values)
         fake_median = np.median(stats["B:measure"].values)
         self.assertGreater(real_median, fake_median)
@@ -1322,7 +1322,7 @@ def test_low_donor_overlap(self):
                                    time_column="group",
                                    reference_column="Ref",
                                    subject_column="subject",
-                                   iterations=999)
+                                   num_iterations=999)
 
         real_median = np.median(stats["A:measure"].values)
         fake_median = np.median(stats["B:measure"].values)
@@ -1353,7 +1353,7 @@ def test_single_donor(self):
                             time_column="group",
                             reference_column="Ref",
                             subject_column="subject",
-                            iterations=999)
+                            num_iterations=999)
 
     def test_create_mismatched_pairs(self):
         metadata_df = pd.DataFrame({
@@ -1414,12 +1414,17 @@ def test_mask_recipient_donor_one(self):
         np.testing.assert_array_equal(recip_mask, exp_r_mask)
 
     def test_simulate_uniform_distro(self):
-        mismatch_peds = [0, 0, 0, 0, 0, 0]
+        # Note: This tests has a VERY small chance to have intermit failures
+        # if by random chance 1, 2, or 3 are not selected by random.choice. 
+        mismatch_peds = [1, 2, 3]
 
         iterations = 999
 
         mismatchpairs_df = _simulate_uniform_distro(mismatch_peds,
                                                     iterations)
+        self.assertIn(1, mismatchpairs_df)
+        self.assertIn(2, mismatchpairs_df)
+        self.assertIn(3, mismatchpairs_df)
         self.assertEquals(mismatchpairs_df.size, iterations)
 
     def test_one_iter_simulate_uniform_distro(self):
@@ -1536,11 +1541,9 @@ def test_per_subject_stats_labels(self):
                                 index=["sample1", "sample2", "sample3",
                                        "sample4"])
         iterations = 10
-        mismatched_pairs_n = 4
 
         p_s_stats = _per_subject_stats(mismatched_peds,
-                                       actual_temp, iterations,
-                                       mismatched_pairs_n)
+                                       actual_temp, iterations)
 
         exp_column_names = ["A:group", "A:n", "A:measure",
                             "B:group", "B:n", "B:measure", "n",
@@ -1554,11 +1557,9 @@ def test_per_subject_stats(self):
                                 index=["sample1", "sample2", "sample3",
                                        "sample4"])
         iterations = 10
-        mismatched_pairs_n = 4
 
         p_s_stats = _per_subject_stats(mismatched_peds,
-                                       actual_temp, iterations,
-                                       mismatched_pairs_n)
+                                       actual_temp, iterations)
 
         exp_test_stats = pd.Series([10, 10, 10, 10])
 
@@ -1580,3 +1581,112 @@ def test_global_stats(self):
         p = _global_stats(p_series)
 
         np.testing.assert_array_equal(p["n"].values, [4])
+
+    def test_peds_sim_stats_good_match(self):
+        value = 1
+        peds_iters = pd.Series(data=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                               index=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+        num_iterations = 10
+
+        count_gte, count_less, per_subject_p = _peds_sim_stats(value,
+                                                               peds_iters,
+                                                               num_iterations)
+        exp_count_gte = 0
+        exp_count_less = 10
+        exp_per_subject_p = (1/11)
+
+        self.assertEqual(count_gte, exp_count_gte)
+        self.assertEqual(count_less, exp_count_less)
+        self.assertEqual(per_subject_p, exp_per_subject_p)
+
+    def test_peds_sim_stats_bad_match(self):
+        value = 0
+        peds_iters = pd.Series(data=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+                               index=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+        num_iterations = 10
+
+        count_gte, count_less, per_subject_p = _peds_sim_stats(value,
+                                                               peds_iters,
+                                                               num_iterations)
+        exp_count_gte = 10
+        exp_count_less = 0
+        exp_per_subject_p = (11/11)
+
+        self.assertEqual(count_gte, exp_count_gte)
+        self.assertEqual(count_less, exp_count_less)
+        self.assertEqual(per_subject_p, exp_per_subject_p)
+
+    def test_peds_sim_stats_equal_match(self):
+        value = .5
+        peds_iters = pd.Series(data=[.5, .5, .5, .5, .5, .5, .5, .5, .5, .5],
+                               index=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+        num_iterations = 10
+
+        count_gte, count_less, per_subject_p = _peds_sim_stats(value,
+                                                               peds_iters,
+                                                               num_iterations)
+        exp_count_gte = 10
+        exp_count_less = 0
+        exp_per_subject_p = (11/11)
+
+        self.assertEqual(count_gte, exp_count_gte)
+        self.assertEqual(count_less, exp_count_less)
+        self.assertEqual(per_subject_p, exp_per_subject_p)
+
+    def test_peds_sim_stats_50_percent_bad_match(self):
+        value = .5
+        peds_iters = pd.Series(data=[.5, 1, .5, 1, .5, 1, .5, 1, .5, 1],
+                               index=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+        num_iterations = 10
+
+        count_gte, count_less, per_subject_p = _peds_sim_stats(value,
+                                                               peds_iters,
+                                                               num_iterations)
+        exp_count_gte = 10
+        exp_count_less = 0
+        exp_per_subject_p = (11/11)
+
+        self.assertEqual(count_gte, exp_count_gte)
+        self.assertEqual(count_less, exp_count_less)
+        self.assertEqual(per_subject_p, exp_per_subject_p)
+
+    def test_peds_sim_stats_50_good_match(self):
+        value = .5
+        peds_iters = pd.Series(data=[.5, 0, .5, 0, .5, 0, .5, 0, .5, 0],
+                               index=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+        num_iterations = 10
+
+        count_gte, count_less, per_subject_p = _peds_sim_stats(value,
+                                                               peds_iters,
+                                                               num_iterations)
+        exp_count_gte = 5
+        exp_count_less = 5
+        exp_per_subject_p = (6/11)
+
+        self.assertEqual(count_gte, exp_count_gte)
+        self.assertEqual(count_less, exp_count_less)
+        self.assertEqual(per_subject_p, exp_per_subject_p)
+
+    def test_peds_sim_stats_99_iters(self):
+        value = .5
+        peds_iters = pd.Series(data=[0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0, .5, 0, .5, 0, .5, 0, .5,
+                                     0, .5, 0, .5, 0], index=list(range(99)))
+        num_iterations = 99
+
+        count_gte, count_less, per_subject_p = _peds_sim_stats(value,
+                                                               peds_iters,
+                                                               num_iterations)
+        exp_count_gte = 49
+        exp_count_less = 50
+        exp_per_subject_p = (50/100)
+
+        self.assertEqual(count_gte, exp_count_gte)
+        self.assertEqual(count_less, exp_count_less)
+        self.assertEqual(per_subject_p, exp_per_subject_p)

From 28a772a5edddae3956ae41bad15ba39ae5dd20b2 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 31 Jul 2024 19:57:54 -0700
Subject: [PATCH 46/50] whitespace

---
 q2_fmt/tests/test_engraftment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/q2_fmt/tests/test_engraftment.py b/q2_fmt/tests/test_engraftment.py
index 7e1bf60..8e94332 100644
--- a/q2_fmt/tests/test_engraftment.py
+++ b/q2_fmt/tests/test_engraftment.py
@@ -1415,7 +1415,7 @@ def test_mask_recipient_donor_one(self):
 
     def test_simulate_uniform_distro(self):
         # Note: This tests has a VERY small chance to have intermit failures
-        # if by random chance 1, 2, or 3 are not selected by random.choice. 
+        # if by random chance 1, 2, or 3 are not selected by random.choice.
         mismatch_peds = [1, 2, 3]
 
         iterations = 999

From 78238f0535b2f43905ab5afe6511b6c3acaceec0 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Thu, 8 Aug 2024 11:23:24 -0700
Subject: [PATCH 47/50] lint

---
 q2_fmt/__init__.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/q2_fmt/__init__.py b/q2_fmt/__init__.py
index 3a39a4e..b2d6681 100644
--- a/q2_fmt/__init__.py
+++ b/q2_fmt/__init__.py
@@ -18,4 +18,3 @@
 
 __all__ = ['cc', 'sample_peds', 'feature_peds',
            'peds', 'peds_heatmap', 'group_timepoints', 'peds_simulation']
-

From ea88903e0fb3c0669a4ac8fc79eb90a43dddf0a5 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 14 Aug 2024 13:55:19 -0700
Subject: [PATCH 48/50] dependency note

---
 q2_fmt/plugin_setup.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 521b719..50e3601 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -367,7 +367,14 @@
                 ' few donors. This method attempts to quantify the extent to'
                 ' which indicator features that are unique to a given donor'
                 ' are transferred to their recipients, as opposed to features'
-                ' that are not indicative of any specific donor.',
+                ' that are not indicative of any specific donor. Note: '
+                ' PEDS Monte Carlo simulation may have dependency issues'
+                ' between samples and comparisons that can make the test'
+                ' overly conservative. This can be fixed by filtering down to'
+                ' a single timepoint before running this methods. Additionally'
+                ' if there are many baseline timepoints, the global test may'
+                ' be too conservative and this can be addressed by filtering'
+                ' out baseline samples prior to running this method.',
     citations=[citations['aggarwala_precise_2021'],
                citations['stouffer_1949_american'],
                citations['Benjamini_fdr_1995']],

From dad7333b1fc22ccce3cd8b323b6b3164b5b4e582 Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 14 Aug 2024 14:12:38 -0700
Subject: [PATCH 49/50] Update q2_fmt/plugin_setup.py

Co-authored-by: Greg Caporaso <192372+gregcaporaso@users.noreply.github.com>
---
 q2_fmt/plugin_setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index 50e3601..b8f7d1e 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -371,7 +371,7 @@
                 ' PEDS Monte Carlo simulation may have dependency issues'
                 ' between samples and comparisons that can make the test'
                 ' overly conservative. This can be fixed by filtering down to'
-                ' a single timepoint before running this methods. Additionally'
+                ' a single timepoint before running this method. Additionally'
                 ' if there are many baseline timepoints, the global test may'
                 ' be too conservative and this can be addressed by filtering'
                 ' out baseline samples prior to running this method.',

From 3a1eb8271b116a070a9ac597c0e9c5b19918e43e Mon Sep 17 00:00:00 2001
From: Chloe Herman <60228108+cherman2@users.noreply.github.com>
Date: Wed, 14 Aug 2024 14:12:49 -0700
Subject: [PATCH 50/50] Update q2_fmt/plugin_setup.py

Co-authored-by: Greg Caporaso <192372+gregcaporaso@users.noreply.github.com>
---
 q2_fmt/plugin_setup.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/q2_fmt/plugin_setup.py b/q2_fmt/plugin_setup.py
index b8f7d1e..fd19fdc 100644
--- a/q2_fmt/plugin_setup.py
+++ b/q2_fmt/plugin_setup.py
@@ -369,7 +369,8 @@
                 ' are transferred to their recipients, as opposed to features'
                 ' that are not indicative of any specific donor. Note: '
                 ' PEDS Monte Carlo simulation may have dependency issues'
-                ' between samples and comparisons that can make the test'
+                ' between samples and the simulated background'
+                ' distribution that can make the test'
                 ' overly conservative. This can be fixed by filtering down to'
                 ' a single timepoint before running this method. Additionally'
                 ' if there are many baseline timepoints, the global test may'