r"""
Empirical Power Estimation (:mod:`skbio.stats.power`)
=====================================================

.. currentmodule:: skbio.stats.power

The purpose of this module is to provide empirical, post-hoc power estimation
of normally and non-normally distributed data. It also provides support to
subsample data to facilitate this analysis.

The underlying principle is based on subsampling and Monte Carlo simulation.
Assume that there is some set of populations, :math:`K_{1}, K_{2}, ... K_{n}`
which have some property, :math:`\mu` such that :math:`\mu_{1} \neq \mu_{2}
\neq ... \neq \mu_{n}`. For each of the populations, a sample, :math:`S` can be
drawn, with a parameter, :math:`x` where :math:`x \approx \mu` and for the
samples, we can use a test, :math:`f`, to show that :math:`x_{1} \neq x_{2}
\neq ... \neq x_{n}`.

Since we know that :math:`\mu_{1} \neq \mu_{2} \neq ... \neq \mu_{n}`,
we know we should reject the null hypothesis. If we fail to reject the null
hypothesis, we have committed a Type II error and our result is a false
negative. We can estimate the frequency of Type II errors at various sampling
depths by repeatedly subsampling the populations and observing how often we
see a false negative. If we repeat this several times for each subsampling
depth, and vary the depths we use, we can start to approximate a relationship
between the number of samples we use and the rate of false negatives, also
called the statistical power of the test.

To generate complete power curves from data which appears underpowered, the
`statsmodels.stats.power` package can be used to solve for an effect size. The
effect size can be used to extrapolate a power curve for the data.

Most functions in this module accept a statistical test function which takes a
list of samples and returns a p value. The test is then evaluated over a series
of subsamples.

Sampling may be handled in two ways. For any set of samples, we may simply
choose to draw :math:`n` observations at random for each sample. Alternatively,
if metadata is available, samples can be matched based on a set of control
categories so that paired samples are drawn at random from the set of available
matches.

Functions
---------

.. autosummary::
    :toctree:

    subsample_power
    subsample_paired_power
    confidence_bound
    paired_subsamples

Examples
--------
Suppose we wanted to test that there's a relationship between two random
variables, `ind` and `dep`. Let's use random subsampling to estimate the
statistical power of our test with an alpha of 0.1, 0.01, and 0.001.

To control for the pseudo-random number generation, we will use a seed.
When using these functions with your own data, you don't need to include the
step.

>>> import numpy as np
>>> np.random.seed(20)
>>> ind = np.random.randint(0, 20, 15)
>>> ind
array([ 3, 15,  9, 11,  7,  2,  0,  8, 19, 16,  6,  6, 16,  9,  5])
>>> dep = (3 * ind + 5 + np.random.randn(15) * 5).round(3)
>>> dep
array([ 15.617,  47.533,  28.04 ,  33.788,  19.602,  12.229,   4.779,
        36.838,  67.256,  55.032,  22.157,   7.051,  58.601,  38.664,
        18.783])

Let's define a test that will draw a list of sample pairs and determine
if they're correlated. We'll use `scipy.stats.pearsonr` which takes two arrays
and returns a correlation coefficient and a p-value representing the
probability the two distributions are correlated.

>>> from scipy.stats import pearsonr
>>> f = lambda x: pearsonr(x[0], x[1])[1]

Now, let's use random sampling to estimate the power of our test on
the first distribution.

>>> samples = [ind, dep]
>>> print("%.3e" % f(samples))
3.646e-08

In `subsample_power`, we can maintain a paired relationship between samples
by setting `draw_mode` to "matched". We can also set our critical value, so
that we estimate power for a critical value of :math:`\alpha = 0.05`, an
estimate for the critical value of 0.01, and a critical value of 0.001.

>>> from skbio.stats.power import subsample_power
>>> pwr_100, counts_100 = subsample_power(test=f,
...                                       samples=samples,
...                                       max_counts=10,
...                                       min_counts=3,
...                                       counts_interval=1,
...                                       draw_mode="matched",
...                                       alpha_pwr=0.1,
...                                       num_iter=25)
>>> pwr_010, counts_010 = subsample_power(test=f,
...                                       samples=samples,
...                                       max_counts=10,
...                                       min_counts=3,
...                                       counts_interval=1,
...                                       draw_mode="matched",
...                                       alpha_pwr=0.01,
...                                       num_iter=25)
>>> pwr_001, counts_001 = subsample_power(test=f,
...                                       samples=samples,
...                                       max_counts=10,
...                                       min_counts=3,
...                                       counts_interval=1,
...                                       draw_mode="matched",
...                                       alpha_pwr=0.001,
...                                       num_iter=25)
>>> counts_100
array([3, 4, 5, 6, 7, 8, 9])
>>> pwr_100.mean(0)
array([ 0.484,  0.844,  0.932,  0.984,  1.   ,  1.   ,  1.   ])
>>> pwr_010.mean(0)
array([ 0.044,  0.224,  0.572,  0.836,  0.928,  0.996,  1.   ])
>>> pwr_001.mean(0)
array([ 0.   ,  0.016,  0.108,  0.332,  0.572,  0.848,  0.956])

Based on this power estimate, as we increase our confidence that we have not
committed a type I error and identified a false positive, the number of samples
we need to be confident that we have not committed a type II error increases.

"""

# ----------------------------------------------------------------------------
# Copyright (c) 2013--, scikit-bio development team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file COPYING.txt, distributed with this software.
# ----------------------------------------------------------------------------

import collections
import copy

import numpy as np
import scipy.stats

from skbio.util._decorator import experimental


@experimental(as_of="0.4.0")
def subsample_power(test, samples, draw_mode='ind', alpha_pwr=0.05, ratio=None,
                    max_counts=50, counts_interval=10, min_counts=None,
                    num_iter=500, num_runs=10):
    r"""Subsamples data to iteratively calculate power

    Parameters
    ----------
    test : function
        The statistical test which accepts a list of arrays of values
        (sample ids or numeric values) and returns a p value or one-dimensional
        array of p values.
    samples : array_like
        `samples` can be a list of lists or a list of arrays where each
        sublist or row in the array corresponds to a sampled group.
    draw_mode : {"ind", "matched"}, optional
        "matched" samples should be used when observations in
        samples have corresponding observations in other groups. For instance,
        this may be useful when working with regression data where
        :math:`x_{1}, x_{2}, ..., x_{n}` maps to
        :math:`y_{1}, y_{2}, ..., y_{n}`. Sample vectors must be the same
        length in "matched" mode.
        If there is no reciprocal relationship between samples, then
        "ind" mode should be used.
    alpha_pwr : float, optional
        The critical value used to calculate the power.
    ratio : 1-D array, optional
        The fraction of the sample counts which should be
        assigned to each group. If this is a 1-D array, it must be the same
        length as `samples`. If no value is supplied (`ratio` is None),
        then an equal number of observations will be drawn for each sample. In
        `matched` mode, this will be set to one.
    max_counts : positive int, optional
        The maximum number of samples per group to draw for effect size
        calculation.
    counts_interval : positive int, optional
        The difference between each subsampling count.
    min_counts : positive int, optional
        How many samples should be drawn for the smallest
        subsample. If this is None, the `counts_interval` will be used.
    num_iter : positive int, optional
        The number of p-values to generate for each point
        on the curve.
    num_runs : positive int, optional
        The number of times to calculate each curve.

    Returns
    -------
    power : array
        The power calculated for each subsample at each count. The array has
        `num_runs` rows, a length with the same number of elements as
        `sample_counts` and a depth equal to the number of p values returned by
        `test`. If `test` returns a float, the returned array will be
        two-dimensional instead of three.
    sample_counts : array
        The number of samples drawn at each power calculation.

    Raises
    ------
    ValueError
        If the `mode` is "matched", an error will occur if the arrays in
        `samples` are not the same length.
    ValueError
        There is a ValueError if there are fewer samples than the minimum
        count.
    ValueError
        If the `counts_interval` is greater than the difference between the
        sample start and the max value, the function raises a ValueError.
    ValueError
        There are not an equal number of groups in `samples` and in `ratios`.
    TypeError
        `test` does not return a float or a 1-dimensional numpy array.


    Examples
    --------
    Let's say we wanted to look at the relationship between the presence of a
    specific bacteria, *Gardnerella vaginalis* in the vaginal community, and
    the probability of a pre or post menopausal woman experiencing a urinary
    tract infection (UTI). Healthy women were enrolled in the study either
    before or after menopause, and followed for eight weeks. Participants
    submitted fecal samples at the beginning of the study, and were then
    followed for clinical symptoms of a UTI. A confirmed UTI was an endpoint
    in the study.

    Using available literature and 16S sequencing, a set of candidate taxa were
    identified as correlated with UTIs, including *G. vaginalis*. In the 100
    women (50 premenopausal and 50 postmenopausal samples) who had UTIs, the
    presence or absence of *G. vaginalis* was confirmed with quantitative PCR.

    We can model the probability that detectable *G. vaginalis* was found in
    these samples using a binomial model. (*Note that this is a simulation.*)

    >>> import numpy as np
    >>> np.random.seed(25)
    >>> pre_rate = np.random.binomial(1, 0.85, size=(50,))
    >>> pre_rate.sum()
    45
    >>> pos_rate = np.random.binomial(1, 0.40, size=(50,))
    >>> pos_rate.sum()
    21

    Let's set up a test function, so we can test the probability of
    finding a difference in frequency between the two groups. We'll use
    `scipy.stats.chisquare` to look for the difference in frequency between
    groups.

    >>> from scipy.stats import chisquare
    >>> test = lambda x: chisquare(np.array([x[i].sum() for i in
    ...     range(len(x))]))[1]

    Let's make sure that our two distributions are different.

    >>> print(round(test([pre_rate, pos_rate]), 3))
    0.003

    Since there are an even number of samples, and we don't have enough
    information to try controlling the data, we'll use
    `skbio.stats.power.subsample_power` to compare the two groups. If we had
    metadata about other risk factors, like a reproductive history, BMI,
    tobacco use, we might want to use
    `skbio.stats.power.subsample_paired_power`.
    We'll also use "ind" `draw_mode`, since there is no linkage between the
    two groups of samples.

    >>> from skbio.stats.power import subsample_power
    >>> pwr_est, counts = subsample_power(test=test,
    ...                                   samples=[pre_rate, pos_rate],
    ...                                   num_iter=100,
    ...                                   num_runs=5,
    ...                                   counts_interval=5)
    >>> counts
    array([ 5, 10, 15, 20, 25, 30, 35, 40, 45])
    >>> np.nanmean(pwr_est, axis=0) # doctest: +NORMALIZE_WHITESPACE
    array([ 0.056,  0.074,  0.226,  0.46 ,  0.61 ,  0.806,  0.952,  1.   ,
            1.   ])
    >>> counts[np.nanmean(pwr_est, axis=0) > 0.8].min()
    30

    So, we can estimate that we will see a significant difference in the
    presence of *G. vaginalis* in the stool of pre and post women with UTIs if
    we have at least 30 samples per group.

    If we wanted to test the relationship of a second candidate taxa which is
    more rare in the population, but may have a similar effect, based on
    available literature, we might also start by trying to identify 30
    samples per group where the second candidate taxa is present.

    Suppose, now, that we want to test that a secondary metabolite seen only in
    the presence of *G vaginalis* to see if it is also correlated with UTIs. We
    can model the abundance of the metabolite as a normal distribution.

    >>> met_pos = (np.random.randn(pre_rate.sum() + pos_rate.sum()) * 2000 +
    ...     2500)
    >>> met_pos[met_pos < 0] = 0
    >>> met_neg = met_neg = (np.random.randn(100 - (pre_rate.sum() +
    ...     pos_rate.sum())) * 2000 + 500)
    >>> met_neg[met_neg < 0] = 0

    Let's compare the populations with a kruskal-wallis test. Physically, there
    cannot be a negative concentration of a chemical, so we've set the lower
    bound at 0. This means that we can no longer assume our distribution is
    normal.

    >>> from scipy.stats import kruskal
    >>> def metabolite_test(x):
    ...     return kruskal(x[0], x[1])[1]
    >>> print(round(metabolite_test([met_pos, met_neg]), 3))
    0.005

    When we go to perform the statistical test on all the data, you might
    notice that there are twice as many samples from women with *G. vaginalis*
    than those without. It might make sense to account for this difference when
    we're testing power. So, we're going to set the `ratio` parameter, which
    lets us draw twice as many samples from women with *G. vaginalis*.

    >>> pwr_est2, counts2 = subsample_power(test=metabolite_test,
    ...                                     samples=[met_pos, met_neg],
    ...                                     counts_interval=5,
    ...                                     num_iter=100,
    ...                                     num_runs=5,
    ...                                     ratio=[2, 1])
    >>> counts2
    array([  5.,  10.,  15.,  20.,  25.,  30.])
    >>> np.nanmean(pwr_est2, axis=0)
    array([ 0.14 ,  0.272,  0.426,  0.646,  0.824,  0.996])
    >>> counts2[np.nanmean(pwr_est2, axis=0) > 0.8].min()
    25.0

    When we consider the number of samples per group needed in the power
    analysis, we need to look at the ratio. The analysis says that we need 25
    samples in the smallest group, in this case, the group of women without
    *G. vaginalis* and 50 samples from women with *G. vaginalis* to see a
    significant difference in the abundance of our secondary metabolite at 80%
    power.

    """

    # Checks the inputs
    ratio, num_p, sample_counts = \
        _check_subsample_power_inputs(test=test,
                                      samples=samples,
                                      draw_mode=draw_mode,
                                      ratio=ratio,
                                      min_counts=min_counts,
                                      max_counts=max_counts,
                                      counts_interval=counts_interval)

    # Prealocates the power array
    power = np.zeros((num_runs, len(sample_counts), num_p))

    # Calculates the power instances
    for id2, c in enumerate(sample_counts):
        count = np.round(c * ratio, 0).astype(int)
        for id1 in range(num_runs):
            ps = _compare_distributions(test=test,
                                        samples=samples,
                                        num_p=num_p,
                                        counts=count,
                                        num_iter=num_iter,
                                        mode=draw_mode)
            power[id1, id2, :] = _calculate_power(ps, alpha_pwr)

    power = power.squeeze()

    return power, sample_counts


@experimental(as_of="0.4.0")
def subsample_paired_power(test, meta, cat, control_cats, order=None,
                           strict_match=True, alpha_pwr=0.05,
                           max_counts=50, counts_interval=10, min_counts=None,
                           num_iter=500, num_runs=10):
    r"""Estimates power iteratively using samples with matching metadata

    Parameters
    ----------
    test : function
        The statistical test which accepts a list of arrays sample ids and
        returns a p value.
    meta : pandas.DataFrame
        The metadata associated with the samples.
    cat : str
        The metadata category being varied between samples.
    control_cats : list
        The metadata categories to be used as controls. For example, if
        you wanted to vary age (`cat` = "AGE"), you might want to control
        for gender and health status (i.e. `control_cats` = ["SEX",
        "HEALTHY"]).
    order : list, optional
        The order of groups in the category. This can be used
        to limit the groups selected. For example, if there's a category with
        groups 'A', 'B' and 'C', and you only want to look at A vs B, `order`
        would be set to ['A', 'B'].
    strict_match : bool, optional
        This determines how data is grouped using
        `control_cats`. If a sample within `meta` has an undefined value (NaN)
        for any of the columns in `control_cats`, the sample will not be
        considered as having a match and will be ignored when `strict_match`
        is True. If `strict_match` is False, missing values (NaN) in the
        `control_cats` can be considered matches.
    alpha_pwr : float, optional
        The critical value used to calculate the power.
    max_counts : positive int, optional
        The maximum number of observations per sample to draw
        for effect size calculation.
    counts_interval : positive int, optional
        The difference between each subsampling count.
    min_counts : positive int, optional
        How many samples should be drawn for the smallest
        subsample. If this is None, the `counts_interval` will be used.
    num_iter : positive int, optional
        The number of p-values to generate for each point on the curve.
    num_runs : positive int, optional
        The number of times to calculate each curve.

    Returns
    -------
    power : array
        The power calculated for each subsample at each count. The array is
        `num_runs` rows, a length with the same number of elements as
        `sample_counts` and a depth equal to the number of p values returned by
        `test`. If `test` returns a float, the returned array will be
        two-dimensional instead of three.
    sample_counts : array
        The number of samples drawn at each power calculation.

    Raises
    ------
    ValueError
        There is a ValueError if there are fewer samples than the minimum
        count.
    ValueError
        If the `counts_interval` is greater than the difference between the
        sample start and the max value, the function raises a ValueError.
    TypeError
        `test` does not return a float or a 1-dimensional numpy array.


    Examples
    --------
    Assume you are interested in the role of a specific cytokine of protein
    translocation in myeloid-lineage cells. You are able to culture two
    macrophage lineages (bone marrow derived phagocytes and
    peritoneally-derived macrophages). Due to unfortunate circumstances, your
    growth media must be acquired from multiple sources (lab, company A,
    company B). Also unfortunate, you must use labor-intensive low throughput
    assays. You have some preliminary measurements, and you'd like to
    predict how many (more) cells you need to analyze for 80% power.

    You have information about 60 cells, which we'll simulate below. Note
    that we are setting a random seed value for consistency.

    >>> import numpy as np
    >>> import pandas as pd
    >>> np.random.seed(25)
    >>> data = pd.DataFrame.from_dict({
    ...     'CELL_LINE': np.random.binomial(1, 0.5, size=(60,)),
    ...     'SOURCE': np.random.binomial(2, 0.33, size=(60,)),
    ...     'TREATMENT': np.hstack((np.zeros((30)), np.ones((30)))),
    ...     'INCUBATOR': np.random.binomial(1, 0.2, size=(60,))})
    >>> data['OUTCOME'] = (0.25 + data.TREATMENT * 0.25) + \
    ...     np.random.randn(60) * (0.1 + data.SOURCE/10 + data.CELL_LINE/5)
    >>> data.loc[data.OUTCOME < 0, 'OUTCOME'] = 0
    >>> data.loc[data.OUTCOME > 1, 'OUTCOME'] = 1

    We will approach this by assuming that the distribution of our outcome is
    not normally distributed, and apply a kruskal-wallis test to compare
    between the cytokine treated and untreated cells.

    >>> from scipy.stats import kruskal
    >>> f = lambda x: kruskal(*[data.loc[i, 'OUTCOME'] for i in x])[1]

    Let's check that cytokine treatment has a significant effect across all
    the cells.

    >>> treatment_stat = [g for g in data.groupby('TREATMENT').groups.values()]
    >>> f(treatment_stat)
    0.0019386336266250209

    Now, let's pick the control categories. It seems reasonable to assume there
    may be an effect of cell line on the treatment outcome, which may be
    attributed to differences in receptor expression. It may also be possible
    that there are differences due cytokine source. Incubators were maintained
    under the same conditions throughout the experiment, within one degree of
    temperature difference at any given time, and the same level of CO2.
    So, at least initially, let's ignore differences due to the incubator.

    It's recommended that as a first pass analysis, control variables be
    selected based on an idea of what may be biologically relevant to the
    system, although further iteration might encourage the consideration of
    variable with effect sizes similar, or larger than the variable of
    interest.

    >>> control_cats = ['SOURCE', 'CELL_LINE']
    >>> from skbio.stats.power import subsample_paired_power
    >>> pwr, cnt = subsample_paired_power(test=f,
    ...                                   meta=data,
    ...                                   cat='TREATMENT',
    ...                                   control_cats=control_cats,
    ...                                   counts_interval=5,
    ...                                   num_iter=25,
    ...                                   num_runs=5)
    >>> cnt
    array([  5.,  10.,  15.,  20.])
    >>> pwr.mean(0)
    array([ 0.24 ,  0.528,  0.68 ,  0.88 ])
    >>> pwr.std(0).round(3)
    array([ 0.088,  0.127,  0.168,  0.08 ])

    Estimating off the power curve, it looks like 20 cells per group may
    provide adequate power for this experiment, although the large variance
    in power might suggest extending the curves or increasing the number of
    samples per group.

    """

    # Handles the order argument
    if order is None:
        order = sorted(meta.groupby(cat).groups.keys())
    order = np.array(order)

    # Checks for the number of sampling pairs available
    meta_pairs, index = _identify_sample_groups(meta, cat, control_cats, order,
                                                strict_match)
    min_obs = min([_get_min_size(meta, cat, control_cats, order, strict_match),
                  np.floor(len(index)*0.9)])
    sub_ids = _draw_paired_samples(meta_pairs, index, min_obs)

    ratio, num_p, sample_counts = \
        _check_subsample_power_inputs(test=test,
                                      samples=sub_ids,
                                      draw_mode='matched',
                                      min_counts=min_counts,
                                      max_counts=max_counts,
                                      counts_interval=counts_interval)

    # Prealocates the power array
    power = np.zeros((num_runs, len(sample_counts), num_p))

    # Calculates power instances
    for id2, c in enumerate(sample_counts):
        for id1 in range(num_runs):
            ps = np.zeros((num_p, num_iter))
            for id3 in range(num_iter):
                subs = _draw_paired_samples(meta_pairs, index, c)
                ps[:, id3] = test(subs)
            power[id1, id2, :] = _calculate_power(ps, alpha_pwr)

    power = power.squeeze()

    return power, sample_counts


@experimental(as_of="0.4.0")
def confidence_bound(vec, alpha=0.05, df=None, axis=None):
    r"""Calculates a confidence bound assuming a normal distribution

    Parameters
    ----------
    vec : array_like
        The array of values to use in the bound calculation.
    alpha : float, optional
        The critical value, used for the confidence bound calculation.
    df : float, optional
        The degrees of freedom associated with the
        distribution. If None is given, df is assumed to be the number of
        elements in specified axis.
    axis : positive int, optional
        The axis over which to take the deviation. When axis
        is None, a single value will be calculated for the whole matrix.

    Returns
    -------
    bound : float
        The confidence bound around the mean. The confidence interval is
        [mean - bound, mean + bound].

    """

    # Determines the number of non-nan counts
    vec = np.asarray(vec)
    vec_shape = vec.shape
    if axis is None and len(vec_shape) == 1:
        num_counts = vec_shape[0] - np.isnan(vec).sum()
    elif axis is None:
        num_counts = vec_shape[0] * vec_shape[1] - np.isnan(vec).sum()
    else:
        num_counts = vec_shape[axis] - np.isnan(vec).sum() / \
            (vec_shape[0] * vec_shape[1])

    # Gets the df if not supplied
    if df is None:
        df = num_counts - 1

    # Calculates the bound
    # In the conversion from scipy.stats.nanstd -> np.nanstd `ddof=1` had to be
    # added to match the scipy default of `bias=False`.
    bound = np.nanstd(vec, axis=axis, ddof=1) / np.sqrt(num_counts - 1) * \
        scipy.stats.t.ppf(1 - alpha / 2, df)

    return bound


@experimental(as_of="0.4.0")
def paired_subsamples(meta, cat, control_cats, order=None, strict_match=True):
    r"""Draws a list of samples varied by `cat` and matched for `control_cats`

    This function is designed to provide controlled samples, based on a
    metadata category. For example, one could control for age, sex, education
    level, and diet type while measuring exercise frequency.

    Parameters
    ----------
    meta : pandas.DataFrame
        The metadata associated with the samples.
    cat : str, list
        The metadata category (or a list of categories) for comparison.
    control_cats : list
        The metadata categories to be used as controls. For example, if you
        wanted to vary age (`cat` = "AGE"), you might want to control for
        gender and health status (i.e. `control_cats` = ["SEX", "HEALTHY"])
    order : list, optional
        The order of groups in the category. This can be used
        to limit the groups selected. For example, if there's a category with
        groups 'A', 'B' and 'C', and you only want to look at A vs B, `order`
        would be set to ['A', 'B'].
    strict_match: bool, optional
        This determines how data is grouped using `control_cats`. If a sample
        within `meta` has an undefined value (`NaN`) for any of the columns in
        `control_cats`, the sample will not be considered as having a match and
        will be ignored when `strict_match` is True. If `strict_match` is
        False, missing values (NaN) in the `control_cats` can be considered
        matches.


    Returns
    -------
    ids : array
        a set of ids which satisfy the criteria. These are not grouped by
        `cat`. An empty array indicates there are no sample ids which satisfy
        the requirements.

    Examples
    --------
    If we have a mapping file for a set of random individuals looking at
    housing, sex, age and antibiotic use.

    >>> import pandas as pd
    >>> import numpy as np
    >>> meta = {'SW': {'HOUSING': '2', 'SEX': 'M', 'AGE': np.nan, 'ABX': 'Y'},
    ...         'TS': {'HOUSING': '2', 'SEX': 'M', 'AGE': '40s', 'ABX': 'Y'},
    ...         'CB': {'HOUSING': '3', 'SEX': 'M', 'AGE': '40s', 'ABX': 'Y'},
    ...         'BB': {'HOUSING': '1', 'SEX': 'M', 'AGE': '40s', 'ABX': 'Y'}}
    >>> meta = pd.DataFrame.from_dict(meta, orient="index")
    >>> meta #doctest: +SKIP
       ABX HOUSING  AGE SEX
    BB   Y       1  40s   M
    CB   Y       3  40s   M
    SW   Y       2  NaN   M
    TS   Y       2  40s   M

    We may want to vary an individual's housing situation, while holding
    constant their age, sex and antibiotic use so we can estimate the effect
    size for housing, and later compare it to the effects of other variables.

    >>> from skbio.stats.power import paired_subsamples
    >>> ids = paired_subsamples(meta, 'HOUSING', ['SEX', 'AGE', 'ABX'])
    >>> np.hstack(ids) #doctest: +ELLIPSIS
    array(['BB', 'TS', 'CB']...)

    So, for this set of data, we can match TS, CB, and BB based on their age,
    sex, and antibiotic use. SW cannot be matched in either group because
    `strict_match` was true, and there is missing AGE data for this sample.

    """

    # Handles the order argument
    if order is None:
        order = sorted(meta.groupby(cat).groups.keys())
    order = np.array(order)

    # Checks the groups in the category
    min_obs = _get_min_size(meta, cat, control_cats, order, strict_match)

    # Identifies all possible subsamples
    meta_pairs, index = _identify_sample_groups(meta=meta,
                                                cat=cat,
                                                control_cats=control_cats,
                                                order=order,
                                                strict_match=strict_match)

    # Draws paired ids
    ids = _draw_paired_samples(meta_pairs=meta_pairs,
                               index=index,
                               num_samps=min_obs)

    return ids


def _get_min_size(meta, cat, control_cats, order, strict_match):
    """Determines the smallest group represented"""

    if strict_match:
        all_cats = copy.deepcopy(control_cats)
        all_cats.append(cat)
        meta = meta[all_cats].dropna()

    return meta.groupby(cat).count().loc[order, control_cats[0]].min()


def _check_nans(x, switch=False):
    r"""Returns False if x is a nan and True is x is a string or number
    """
    if isinstance(x, str):
        return True
    elif isinstance(x, (float, int)):
        return not np.isnan(x)
    elif switch and isinstance(x, (list, tuple)) and np.nan in x:
        return False
    elif switch and isinstance(x, (list, tuple)):
        return True
    else:
        raise TypeError('input must be a string, float or a nan')


def _calculate_power(p_values, alpha=0.05):
    r"""Calculates statistical power empirically

    Parameters
    ----------
    p_values : 1-D array
        A 1-D numpy array with the test results.

    alpha : float
        The critical value for the power calculation.

    Returns
    -------
    power : float
        The empirical power, or the fraction of observed p values below the
        critical value.

    """

    p_values = np.atleast_2d(p_values)

    w = (p_values < alpha).sum(axis=1)/p_values.shape[1]

    return w


def _compare_distributions(test, samples, num_p, counts=5, mode="ind",
                           num_iter=100):
    r"""Compares two distribution arrays iteratively

    Parameters
    ----------
    test : function
        The statistical test which accepts an array_like of sample ids
        (list of lists) and returns a p-value. This can be a one-dimensional
        array, or a float.
    samples : list of arrays
        A list where each 1-d array represents a sample. If `mode` is
        "matched", there must be an equal number of observations in each
        sample.
    num_p : positive int, optional
        The number of p-values returned by the test.
    counts : positive int or 1-D array, optional
        The number of samples to draw from each distribution.
        If this is a 1-D array, the length must correspond to the number of
        samples. The function will not draw more observations than are in a
        sample. In "matched" `mode`, the same number of observations will be
        drawn from each group.
    mode : {"ind", "matched", "paired"}, optional
        "matched" samples should be used when observations in
        samples have corresponding observations in other groups. For instance,
        this may be useful when working with regression data where
        :math:`x_{1}, x_{2}, ..., x_{n}` maps to :math:`y_{1}, y_{2}, ... ,
        y_{n}`.
    num_iter : positive int, optional
        Default 1000. The number of p-values to generate for each point on the
        curve.

    Returns
    -------
    p_values : array
        The p-values for the subsampled tests. If `test` returned a single
        p value, p_values is a one-dimensional array. If `test` returned an
        array, `p_values` has dimensions `num_iter` x `num_p`

    Raises
    ------
    ValueError
        If mode is not "ind" or "matched".
    ValueError
        If the arrays in samples are not the same length in "matched" mode.
    ValueError
        If counts is a 1-D array and counts and samples are different lengths.

    """

    # Prealocates the pvalue matrix
    p_values = np.zeros((num_p, num_iter))

    # Determines the number of samples per group
    num_groups = len(samples)
    samp_lens = [len(sample) for sample in samples]

    if isinstance(counts, int):
        counts = np.array([counts] * num_groups)

    for idx in range(num_iter):
        if mode == "matched":
            pos = np.random.choice(np.arange(0, samp_lens[0]), counts[0],
                                   replace=False)
            subs = [sample[pos] for sample in samples]
        else:
            subs = [np.random.choice(np.array(pop), counts[i], replace=False)
                    for i, pop in enumerate(samples)]

        p_values[:, idx] = test(subs)

    if num_p == 1:
        p_values = p_values.squeeze()

    return p_values


def _check_subsample_power_inputs(test, samples, draw_mode='ind', ratio=None,
                                  max_counts=50, counts_interval=10,
                                  min_counts=None):
    r"""Makes sure that everything is sane before power calculations

    Parameters
    ----------
    test : function
        The statistical test which accepts a list of arrays of values
        (sample ids or numeric values) and returns a p value or one-dimensional
        array of p values.
    samples : array_like
        `samples` can be a list of lists or a list of arrays where each
        sublist or row in the array corresponds to a sampled group.
    draw_mode : {"ind", "matched"}, optional
        "matched" samples should be used when observations in
        samples have corresponding observations in other groups. For instance,
        this may be useful when working with regression data where
        :math:`x_{1}, x_{2}, ..., x_{n}` maps to
        :math:`y_{1}, y_{2}, ..., y_{n}`. Sample vectors must be the same
        length in "matched" mode.
        If there is no reciprocal relationship between samples, then
        "ind" mode should be used.
    ratio : 1-D array, optional
        The fraction of the sample counts which should be
        assigned to each group. If this is a 1-D array, it must be the same
        length as `samples`. If no value is supplied (`ratio` is None),
        then an equal number of observations will be drawn for each sample. In
        `matched` mode, this will be set to one.
    max_counts : positive int, optional
        The maximum number of samples per group to draw for effect size
        calculation.
    counts_interval : positive int, optional
        The difference between each subsampling count.
    min_counts : positive int, optional
        How many samples should be drawn for the smallest
        subsample. If this is None, the `counts_interval` will be used.

    Returns
    -------
    ratio : 1-D array
        The fraction of the sample counts which should be assigned to each
        group.
    num_p : positive integer
        The number of p values returned by `test`.
    sample_counts : array
        The number of samples drawn at each power calculation.

    Raises
    ------
    ValueError
        If the `mode` is "matched", an error will occur if the arrays in
        `samples` are not the same length.
    ValueError
        There is a ValueError if there are fewer samples than the minimum
        count.
    ValueError
        If the `counts_interval` is greater than the difference between the
        sample start and the max value, the function raises a ValueError.
    ValueError
        There are not an equal number of groups in `samples` and in `ratios`.
    TypeError
        `test` does not return a float or a 1-dimensional numpy array.

    """

    if draw_mode not in {'ind', 'matched'}:
        raise ValueError('mode must be "matched" or "ind".')

    # Determines the minimum number of ids in a category
    id_counts = np.array([len(id_) for id_ in samples])
    num_ids = id_counts.min()
    # Determines the number of groups
    num_groups = len(samples)

    # Checks that "matched" mode is handled appropriately
    if draw_mode == "matched":
        for id_ in samples:
            if not len(id_) == num_ids:
                raise ValueError('Each vector in samples must be the same '
                                 'length in "matched" draw_mode.')

    # Checks the number of counts is appropriate
    if min_counts is None:
        min_counts = counts_interval
    if (max_counts - min_counts) < counts_interval:
        raise ValueError("No subsamples of the specified size can be drawn.")

    # Checks the ratio argument is sane
    if ratio is None or draw_mode == 'matched':
        ratio = np.ones((num_groups))
    else:
        ratio = np.asarray(ratio)
    if not ratio.shape == (num_groups,):
        raise ValueError('There must be a ratio for each group.')

    ratio_counts = np.array([id_counts[i] / ratio[i]
                             for i in range(num_groups)])
    largest = ratio_counts.min()

    # Determines the number of p values returned by the test
    p_return = test(samples)
    if isinstance(p_return, float):
        num_p = 1
    elif isinstance(p_return, np.ndarray) and len(p_return.shape) == 1:
        num_p = p_return.shape[0]
    else:
        raise TypeError('test must return a float or one-dimensional array.')

    # Calculates the same counts
    sample_counts = np.arange(min_counts,
                              min(max_counts, largest),
                              counts_interval)

    return ratio, num_p, sample_counts


def _identify_sample_groups(meta, cat, control_cats, order, strict_match):
    """Aggregates samples matches for `control_cats` that vary by `cat`

    Parameters
    ----------
    meta : pandas.DataFrame
        The metadata associated with the samples.
    cat : str, list
        The metadata category (or a list of categories) for comparison.
    control_cats : list
        The metadata categories to be used as controls. For example, if you
        wanted to vary age (`cat` = "AGE"), you might want to control for
        gender and health status (i.e. `control_cats` = ["SEX", "HEALTHY"])
    order : list
        The order of groups in the category. This can be used
        to limit the groups selected. For example, if there's a category with
        groups 'A', 'B' and 'C', and you only want to look at A vs B, `order`
        would be set to ['A', 'B'].
    ctrl_pos : int
        The location of the smallest group in `order`.
    strict_match: bool, optional
        This determines how data is grouped using `control_cats`. If a sample
        within `meta` has an undefined value (`NaN`) for any of the columns in
        `control_cats`, the sample will not be considered as having a match and
        will be ignored when `strict_match` is True. If `strict_match` is
        False, missing values (NaN) in the `control_cats` can be considered
        matches.

    Returns
    -------
    meta_pairs : dict
        Describes the categories matched for metadata. The
        `control_cat`-grouped samples are numbered, corresponding to the
        second list in `index`. The group is keyed to the list of sample arrays
        with the same length of `order`.
    index : list
        A list of numpy arrays describing the positions of samples to be drawn.
        The first array is an index array. The second gives an integer
        corresponding to the `control_cat`-group, and the third lists the
        position of the reference group sample in the list of samples.

    """

    # Sets up variables to be filled
    meta_pairs = {}
    index = []
    i1 = 0

    # Groups the data by the control groups
    ctrl_groups = meta.groupby(control_cats).groups
    # Identifies the samples that satisfy the control pairs. Keys are iterated
    # in sorted order so that results don't change with different dictionary
    # ordering.
    for g in sorted(ctrl_groups, key=lambda k: str(k)):
        ids = ctrl_groups[g]
        # If strict_match, Skips over data that has nans
        if not _check_nans(g, switch=True) and strict_match:
            continue
        # Draws the samples that are matched for control cats
        m_ids = meta.loc[ids].groupby(cat).groups
        # Checks if samples from the cat groups are represented in those
        # Samples
        id_vecs = [sorted(m_ids[o]) for o in order if o in
                   m_ids]
        # If all groups are represented, the index and results are retained
        if len(id_vecs) == len(order):
            min_vec = np.array([len(v) for v in id_vecs])
            loc_vec = np.arange(0, min_vec.min())
            meta_pairs[i1] = id_vecs
            index.append(np.zeros(loc_vec.shape) + i1)
            i1 = i1 + 1
        # If the groups are not represented, an empty array gets passed
        else:
            index.append(np.array([]))

    # Converts index to a 1d array
    index = np.hstack(index)

    # If index is empty, sets up meta_paris with a no key.
    if not meta_pairs:
        meta_pairs['no'] = order

    return meta_pairs, index


def _draw_paired_samples(meta_pairs, index, num_samps):
    """Draws a random set of ids from a matched list

    Parameters
    ----------
    meta_pairs : dict
        Describes the categories matched for metadata. The
        `control_cat`-grouped samples are numbered, corresponding to the
        second list in `index`. The group is keyed to the list of sample arrays
        with the same length of `order`.
    index : list
        A list of numpy arrays describing the positions of samples to be drawn.
        The first array is an index array. The second gives an integer
        corresponding to the `control_cat`-group, and the third lists the
        position of the reference group sample in the list of samples.

    Returns
    -------
    ids : list
        A set of randomly selected ids groups from each group.
    """

    # Handles an empty paired vector
    if 'no' in meta_pairs:
        return [np.array([]) for o in meta_pairs['no']]

    # Identifies the absolute positions of the control group being drawn
    set_pos = np.random.choice(index, int(num_samps),
                               replace=False).astype(int)

    subs = []

    # Draws the other groups. Get a collection.Counter object for simplicity
    counter = collections.Counter(set_pos)

    # counter.items() order isn't guaranteed in python 3.6 and then the random
    # choice isn't reproducible between python version, even specifying seed;
    # so we access such values through sets.
    set_list = set(set_pos)

    # then, as stated by @RNAer, since we can't assure that items in sets are
    # ordered, we choose to order set_list before accessing values
    set_list = sorted(set_list)

    # now set_list is ordered and we can iterate over it to get counter obj
    for set_ in set_list:
        num_ = counter[set_]
        r2 = [np.random.choice(col, num_, replace=False) for col in
              meta_pairs[set_]]
        subs.append(r2)

    ids = [np.hstack(ids) for ids in zip(*subs)]

    return ids


def _calculate_power_curve(test, samples, sample_counts, ratio=None,
                           mode='ind', num_iter=1000, alpha=0.05):
    r"""Generates an empirical power curve for the samples.

    Parameters
    ----------
    test : function
        The statistical test which accepts a list of arrays of values and
        returns a p value.
    samples : array_like
        `samples` can be a list of lists or an array where each sublist or row
        in the array corresponds to a sampled group.
    sample_counts : 1-D array
        A vector of the number of samples which should be sampled in each
        curve.
    mode : {"ind", "matched"}, optional
        "matched" samples should be used when observations in
        samples have corresponding observations in other groups. For instance,
        this may be useful when working with regression data where
        :math:`x_{1}, x_{2}, ..., x_{n}` maps to :math:`y_{1}, y_{2}, ... ,
        y_{n}`.
    ratio : 1-D array, optional
        The fraction of the sample counts which should be
        assigned to each group. If this is a 1-D array, it must be the same
        length as `samples`. If no value is supplied (`ratio` is None),
        then an equal number of observations will be drawn for each sample.
    num_iter : int
        The default is 1000. The number of p-values to generate for each point
        on the curve.

    Returns
    -------
    p_values : array
        The p-values associated with the input sample counts.
    Raises
    ------
    ValueError
        If ratio is an array and ratio is not the same length as samples
    """
    # Casts array-likes to arrays
    sample_counts = np.asarray(sample_counts)
    # Determines the number of groups
    num_groups = len(samples)
    num_samps = len(sample_counts)
    if isinstance(alpha, float):
        vec = True
        pwr = np.zeros((num_samps))
        alpha = np.array([alpha])
    else:
        vec = False
        num_crit = alpha.shape[0]
        pwr = np.zeros((num_crit, num_samps))
    # Checks the ratio argument
    if ratio is None:
        ratio = np.ones((num_groups))
    ratio = np.asarray(ratio)
    if not ratio.shape == (num_groups,):
        raise ValueError('There must be a ratio for each group.')

    # Loops through the sample sizes
    for id2, s in enumerate(sample_counts):
        count = np.round(s * ratio, 0).astype(int)
        for id1, a in enumerate(alpha):
            ps = _compare_distributions(test=test,
                                        samples=samples,
                                        counts=count,
                                        num_p=1,
                                        num_iter=num_iter,
                                        mode=mode)
            if vec:
                pwr[id2] = _calculate_power(ps, a)
            else:
                pwr[id1, id2] = _calculate_power(ps, a)

    return pwr