<commit_before>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 microbiome data. It also provides support to subsample data to faciliate
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, u such that :math:`\mu_{1} \neq \mu_{2} \neq ...
\neq \mu_{n}`. For each of the populations, a sample, S can be drawn, with a
parameter, x where :math:`x \approx \mu` and for the samples, we can use a
test, f, to show that :math:`x_{1} \neq x_{2} \neq ... \neq x_{n}`.

Since we known 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 comitted 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.

We can then use the rate of false negatives and use the statsmodels.power
package to solve for an effect size. This can be used to extrapolate a power
curve for the data.

The general format for functions in this module is to define a statistical
test function which will take a list of ids, and return a p value. The test is
then evaluated over a series of subsample sizes.

With microbiome data, there are three ways we can approach selecting our
sample. We may choose to simply draw n observations at random from the two
underlying samples. Alternatively, we can draw subsamples which are
significantly different. Finally, we can try to match samples based on a set
of control categories.

Functions
---------

.. autosummary::
    :toctree: generated/

    get_subsampled_power
    confidence_bound
    bootstrap_power_curve
    get_significant_subsample
    get_paired_subsamples

"""

# -----------------------------------------------------------------------------
# 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.
# -----------------------------------------------------------------------------

from __future__ import division
from future.utils import viewitems
from numpy import (array, zeros, ones, round as nround, hstack, isnan,
                   sqrt, arange)
from numpy.random import choice
from scipy.stats import t, nanstd


def get_subsampled_power(mode, test, meta=None, cat=None, control_cats=None,
                         order=None, strict=True, samples=None, sub_size=None,
                         scaling=5, alpha_pwr=0.05, min_counts=20,
                         max_counts=50, counts_interval=10, counts_start=None,
                         num_iter=500, num_runs=10):
    r"""Subsamples data to iterative calculate power

      Parameters
    ----------
    test : function
        the statistical test which accepts an array-like of sample ids
        (list of lists) and returns a p-value.
    meta : {None, dataframe}
        the metadata associated with the samples. Required for "PAIRED" mode.
    cat : {None, str}
        the metadata categories for comparison. Required for "PAIRED" mode.
    control_cats : {None, list}
        the metadata categories to be used as controls. For example, if you
        wanted to control age (`cat` = "AGE"), you might want to control for
        gender and health status (i.e. `control_cats` = ["SEX", "HEALTHY"]).
        Required for "PAIRED" mode.
    order : {None, list}, optional
        Default is None. 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: bool, optional
        Default is True. If a missing value (nan) is encountered, the group
        will be skipped when 'strict' is True.
    samples : {None, 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. Required for "ALL" and "SIG"
        mode.
    sub_size : {None, int}, optional
        the maximum number of samples to select from a group. If no value is
        provided, this will be the same as the size of the smallest group.
        Otherwise, this will be compared to the size of the smallest group, and
        which ever is lower will be used.
    alpha_pwr : float, optional
        Default is 0.05. The alpha value used to calculate the power.
    min_counts : unsigned int, optional
        Default is 20. The minimum number of paired samples which must exist
        for a category and set of control categories to be able to subsample
        and make power calculations.
    max_counts : unsigned int, optional
        Default is 50. The maximum number of samples per group to draw for
        effect size calculation.
    counts_interval : unsigned int, optional
        Default is 10.
    counts_start : {None, unsigned int}, optional
        Defualt is None. How many samples should be drawn for the smallest
        subsample. If this is None, the `counts_interval` will be used.
    num_iter : unsigned int
        Default is 1000. The number of p-values to generate for each point
        on the curve.
    num_runs : unsigned int
        Default is 10. The number of times to calculate each curve.

    Returns
    -------
    power : array
        power calculated for each subsample at each count
    sample_counts : array
        the number of samples drawn at each power calculation

    Raises
    ------
    ValueError
        if mode is PAIRED and meta, cat or control_cats is None
    ValueError
        if mode is ALL or SIG and samples is None
    RuntimeError
        if there are fewer samples than the minimum count
    RuntimeError
        if the `counts_interval` is greater than the difference between the
        sample start and the max value.

    Examples
    --------
    Suppose we wanted to look at the power curve for two variables, `ind` and
    `dep`, using completely random subsampling. To control for the pseudo
    random number generation, we will use a seed.

    >>> import numpy as np
    >>> np.random.seed(20)
    >>> ind = np.random.randint(0, 20, 15)
    >>> print ind
    [ 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)
    >>> print dep
    [ 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 the `pearsonr` function from scipy, which
    returns the pearson correlation coefficient and the probability value
    that the data is not correlated. The function takes two vectors as its
    input.

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

    Now, let's use random sampling to estimate the power of our test on
    the first distribution. Since our test picks pairs, the "samples" vector
    will just be a list of positions in each array.

    >>> samples = [np.arange(0, 15, 1)]
    >>> print f(samples)
    3.64594525966e-08

    Since we know the two samples are correlated overall, let's try using
    completely random subsampling. This is recommended when sample populations
    are of simillar size, giving each sampled contained in the population a
    simillar probability of being drawn. If one sample is much larger than the
    other, signifigant subsampling can help decrease some of the noise.

    >>> from skbio.stats.power import get_subsampled_power
    >>> pwr_ests, counts = get_subsampled_power(mode="ALL",
    ...                                         test=f,
    ...                                         samples=samples,
    ...                                         min_counts=3,
    ...                                         max_counts=10,
    ...                                         counts_start=3,
    ...                                         counts_interval=1)
    >>> print counts
    [3 4 5 6 7 8 9]
    >>> print pwr_ests
    [[ 0.22   0.652  0.89   0.958  0.992  1.     1.   ]
     [ 0.234  0.642  0.876  0.96   0.99   1.     1.   ]
     [ 0.242  0.654  0.848  0.946  0.998  1.     1.   ]
     [ 0.244  0.664  0.884  0.946  0.988  1.     1.   ]
     [ 0.248  0.666  0.866  0.948  0.986  1.     1.   ]
     [ 0.242  0.658  0.9    0.94   0.99   1.     1.   ]
     [ 0.242  0.638  0.874  0.952  0.992  1.     1.   ]
     [ 0.24   0.66   0.904  0.95   0.988  1.     1.   ]
     [ 0.232  0.64   0.912  0.972  0.988  1.     1.   ]
     [ 0.256  0.646  0.854  0.952  0.992  1.     1.   ]]

    The power_est can then be used to fit an effect_size using the power module
    if the statsmodel library, or can be average and plotted.

    """

    # Checks the mode arguments
    if mode == "PAIRED":
        if meta is None or cat is None or control_cats is None:
            raise ValueError("PAIRED mode requires a meta dataframe, a "
                             "cat to vary and a set of control_cats.")
        else:
            sub_ids = get_paired_subsamples(meta, cat, control_cats, order,
                                            strict)
    elif mode == 'SIG':
        if samples is None:
            raise ValueError("SIG mode requires samples be defined.")
        else:
            sub_ids = get_significant_subsample([test], samples, sub_size,
                                                alpha_pwr, num_iter, scaling)
    elif mode == "ALL":
        if samples is None:
            raise ValueError("ALL mode requires samples be defined.")
        else:
            sub_ids = samples
    else:
        raise ValueError('%s is not a supported mode. Modes are "ALL", "SIG", '
                         'and "PAIRED".' % mode)

    # Determines the minium number of ids in a category
    num_ids = array([len(id_) for id_ in sub_ids]).min()

    # Checks there are enough samples to subsample
    if num_ids <= min_counts:
        raise RuntimeError('There are not enough samples for subsampling.')

    # Calculates the effect size vector
    if counts_start is None:
        counts_start = counts_interval

    if (max_counts - counts_start) < counts_interval:
        raise RuntimeError("No subsamples of the specified size can be drawn.")

    sample_counts = arange(counts_start,
                           min(max_counts, num_ids),
                           counts_interval)

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

    # Calculates the first power curve instance
    power[0, :] = _calculate_power_curve(test, sub_ids, sample_counts,
                                         num_iter=num_iter, alpha=alpha_pwr)

    # Calculates the power instances
    for id1 in arange(1, num_runs):
        # Gets the subsample
        if mode == "PAIRED":
            sub_ids = get_paired_subsamples(meta, cat, control_cats, order,
                                            strict)
        elif mode == 'SIGNIFICANT':
            sub_ids = get_significant_subsample([test], samples, sub_size,
                                                alpha_pwr, num_iter, scaling)
        else:
            sub_ids = samples
        # Calculates the power curve
        power[id1, :] = _calculate_power_curve(test, sub_ids, sample_counts,
                                               num_iter=num_iter,
                                               alpha=alpha_pwr)
<commit_msg>fixes line error<commit_after>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 microbiome data. It also provides support to subsample data to faciliate
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, u such that :math:`\mu_{1} \neq \mu_{2} \neq ...
\neq \mu_{n}`. For each of the populations, a sample, S can be drawn, with a
parameter, x where :math:`x \approx \mu` and for the samples, we can use a
test, f, to show that :math:`x_{1} \neq x_{2} \neq ... \neq x_{n}`.

Since we known 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 comitted 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.

We can then use the rate of false negatives and use the `statsmodels.power`
package to solve for an effect size. This can be used to extrapolate a power
curve for the data.

The general format for functions in this module is to define a statistical
test function which will take a list of ids, and return a p value. The test is
then evaluated over a series of subsample sizes.

With microbiome data, there are three ways we can approach selecting our
sample. We may choose to simply draw n observations at random from the two
underlying samples. Alternatively, we can draw subsamples which are
significantly different. Finally, we can try to match samples based on a set
of control categories.

Functions
---------

.. autosummary::
    :toctree: generated/

    get_subsampled_power
    confidence_bound
    bootstrap_power_curve
    get_significant_subsample
    get_paired_subsamples

"""

# -----------------------------------------------------------------------------
# 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.
# -----------------------------------------------------------------------------

from __future__ import division
from future.utils import viewitems
from numpy import (array, zeros, ones, round as nround, hstack, isnan,
                   sqrt, arange)
from numpy.random import choice
from scipy.stats import t, nanstd


def get_subsampled_power(mode, test, meta=None, cat=None, control_cats=None,
                         order=None, strict=True, samples=None, sub_size=None,
                         scaling=5, alpha_pwr=0.05, min_counts=20,
                         max_counts=50, counts_interval=10, counts_start=None,
                         num_iter=500, num_runs=10):
    r"""Subsamples data to iterative calculate power

      Parameters
    ----------
    test : function
        the statistical test which accepts an array-like of sample ids
        (list of lists) and returns a p-value.
    meta : {None, dataframe}
        the metadata associated with the samples. Required for "PAIRED" mode.
    cat : {None, str}
        the metadata categories for comparison. Required for "PAIRED" mode.
    control_cats : {None, list}
        the metadata categories to be used as controls. For example, if you
        wanted to control age (`cat` = "AGE"), you might want to control for
        gender and health status (i.e. `control_cats` = ["SEX", "HEALTHY"]).
        Required for "PAIRED" mode.
    order : {None, list}, optional
        Default is None. 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: bool, optional
        Default is True. If a missing value (nan) is encountered, the group
        will be skipped when 'strict' is True.
    samples : {None, 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. Required for "ALL" and "SIG"
        mode.
    sub_size : {None, int}, optional
        the maximum number of samples to select from a group. If no value is
        provided, this will be the same as the size of the smallest group.
        Otherwise, this will be compared to the size of the smallest group, and
        which ever is lower will be used.
    alpha_pwr : float, optional
        Default is 0.05. The alpha value used to calculate the power.
    min_counts : unsigned int, optional
        Default is 20. The minimum number of paired samples which must exist
        for a category and set of control categories to be able to subsample
        and make power calculations.
    max_counts : unsigned int, optional
        Default is 50. The maximum number of samples per group to draw for
        effect size calculation.
    counts_interval : unsigned int, optional
        Default is 10.
    counts_start : {None, unsigned int}, optional
        Defualt is None. How many samples should be drawn for the smallest
        subsample. If this is None, the `counts_interval` will be used.
    num_iter : unsigned int
        Default is 1000. The number of p-values to generate for each point
        on the curve.
    num_runs : unsigned int
        Default is 10. The number of times to calculate each curve.

    Returns
    -------
    power : array
        power calculated for each subsample at each count
    sample_counts : array
        the number of samples drawn at each power calculation

    Raises
    ------
    ValueError
        if mode is PAIRED and meta, cat or control_cats is None
    ValueError
        if mode is ALL or SIG and samples is None
    RuntimeError
        if there are fewer samples than the minimum count
    RuntimeError
        if the `counts_interval` is greater than the difference between the
        sample start and the max value.

    Examples
    --------
    Suppose we wanted to look at the power curve for two variables, `ind` and
    `dep`, using completely random subsampling. To control for the pseudo
    random number generation, we will use a seed.

    >>> import numpy as np
    >>> np.random.seed(20)
    >>> ind = np.random.randint(0, 20, 15)
    >>> print ind
    [ 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)
    >>> print dep
    [ 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 the `pearsonr` function from scipy, which
    returns the pearson correlation coefficient and the probability value
    that the data is not correlated. The function takes two vectors as its
    input.

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

    Now, let's use random sampling to estimate the power of our test on
    the first distribution. Since our test picks pairs, the "samples" vector
    will just be a list of positions in each array.

    >>> samples = [np.arange(0, 15, 1)]
    >>> print f(samples)
    3.64594525966e-08

    Since we know the two samples are correlated overall, let's try using
    completely random subsampling. This is recommended when sample populations
    are of simillar size, giving each sampled contained in the population a
    simillar probability of being drawn. If one sample is much larger than the
    other, signifigant subsampling can help decrease some of the noise.

    >>> from skbio.stats.power import get_subsampled_power
    >>> pwr_ests, counts = get_subsampled_power(mode="ALL",
    ...                                         test=f,
    ...                                         samples=samples,
    ...                                         min_counts=3,
    ...                                         max_counts=10,
    ...                                         counts_start=3,
    ...                                         counts_interval=1)
    >>> print counts
    [3 4 5 6 7 8 9]
    >>> print pwr_ests
    [[ 0.22   0.652  0.89   0.958  0.992  1.     1.   ]
     [ 0.234  0.642  0.876  0.96   0.99   1.     1.   ]
     [ 0.242  0.654  0.848  0.946  0.998  1.     1.   ]
     [ 0.244  0.664  0.884  0.946  0.988  1.     1.   ]
     [ 0.248  0.666  0.866  0.948  0.986  1.     1.   ]
     [ 0.242  0.658  0.9    0.94   0.99   1.     1.   ]
     [ 0.242  0.638  0.874  0.952  0.992  1.     1.   ]
     [ 0.24   0.66   0.904  0.95   0.988  1.     1.   ]
     [ 0.232  0.64   0.912  0.972  0.988  1.     1.   ]
     [ 0.256  0.646  0.854  0.952  0.992  1.     1.   ]]

    The power_est can then be used to fit an effect_size using the power module
    if the statsmodel library, or can be average and plotted.

    """

    # Checks the mode arguments
    if mode == "PAIRED":
        if meta is None or cat is None or control_cats is None:
            raise ValueError("PAIRED mode requires a meta dataframe, a "
                             "cat to vary and a set of control_cats.")
        else:
            sub_ids = get_paired_subsamples(meta, cat, control_cats, order,
                                            strict)
    elif mode == 'SIG':
        if samples is None:
            raise ValueError("SIG mode requires samples be defined.")
        else:
            sub_ids = get_significant_subsample([test], samples, sub_size,
                                                alpha_pwr, num_iter, scaling)
    elif mode == "ALL":
        if samples is None:
            raise ValueError("ALL mode requires samples be defined.")
        else:
            sub_ids = samples
    else:
        raise ValueError('%s is not a supported mode. Modes are "ALL", "SIG", '
                         'and "PAIRED".' % mode)

    # Determines the minium number of ids in a category
    num_ids = array([len(id_) for id_ in sub_ids]).min()

    # Checks there are enough samples to subsample
    if num_ids <= min_counts:
        raise RuntimeError('There are not enough samples for subsampling.')

    # Calculates the effect size vector
    if counts_start is None:
        counts_start = counts_interval

    if (max_counts - counts_start) < counts_interval:
        raise RuntimeError("No subsamples of the specified size can be drawn.")

    sample_counts = arange(counts_start,
                           min(max_counts, num_ids),
                           counts_interval)

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

    # Calculates the first power curve instance
    power[0, :] = _calculate_power_curve(test, sub_ids, sample_counts,
                                         num_iter=num_iter, alpha=alpha_pwr)

    # Calculates the power instances
    for id1 in arange(1, num_runs):
        # Gets the subsample
        if mode == "PAIRED":
            sub_ids = get_paired_subsamples(meta, cat, control_cats, order,
                                            strict)
        elif mode == 'SIG':
            sub_ids = get_significant_subsample([test], samples, sub_size,
                                                alpha_pwr, num_iter, scaling)
        else:
            sub_ids = samples
        # Calculates the power curve
        power[id1, :] = _calculate_power_curve(test, sub_ids, sample_counts,
                                               num_iter=num_iter,
                                               alpha=alpha_pwr)

{'fraction_non_alphanumeric': 0.08949944474149632, 'fraction_numerical': 0.032246123473039116, 'mean_word_length': 2.045722159589127, 'pattern_counts': {'":': 6, '<': 7, '<?xml version=': 0, '>': 87, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 52, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2808220', 'n_tokens_mistral': 7115, 'n_tokens_neox': 6546, 'n_words': 2902}