# scipy-yli: Helpful SciPy utilities and recipes # Copyright © 2022 Lee Yingtong Li (RunasSudo) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Affero General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Affero General Public License for more details. # # You should have received a copy of the GNU Affero General Public License # along with this program. If not, see . import numpy as np import pandas as pd from scipy import stats import statsmodels.api as sm import functools import warnings from .config import config from .utils import Estimate, Interval, PValueStyle, as_2groups, as_numeric, check_nan, convert_pandas_nullable, fmt_p # ---------------- # Student's t test class TTestResult: """ Result of a Student's *t* test See :func:`yli.ttest_ind`. """ def __init__(self, *, statistic, dof, pvalue, dep, ind, group1, group2, mu1, mu2, sd1, sd2, delta, delta_direction): #: *t* statistic (*float*) self.statistic = statistic #: Degrees of freedom of the *t* distribution (*int*) self.dof = dof #: *p* value for the *t* statistic (*float*) self.pvalue = pvalue #: Name of the dependent variable (*str*) self.dep = dep #: Name of the independent variable (*str*) self.ind = ind #: Name of the first group (*str*) self.group1 = group1 #: Name of the second group (*str*) self.group2 = group2 #: Mean of the first group (*float*) self.mu1 = mu1 #: Mean of the second group (*float*) self.mu2 = mu2 #: Standard deviation of the first group (*float*) self.sd1 = sd1 #: Standard deviation of the second group (*float*) self.sd2 = sd2 #: Absolute value of the mean difference (:class:`yli.utils.Estimate`) self.delta = delta #: Description of the direction of the effect (*str*) self.delta_direction = delta_direction def _comparison_table(self, html): """Return a table showing the means/SDs for each group""" table_data = [[ '{:.2f} ({:.2f})'.format(self.mu1, self.sd1), '{:.2f} ({:.2f})'.format(self.mu2, self.sd2), ]] if html: table = pd.DataFrame(table_data, index=pd.Index(['\ue000 (SD)'], name=self.dep), columns=pd.Index([self.group1, self.group2], name=self.ind)) # U+E000 is in Private Use Area, mark μ symbol table_str = table._repr_html_() return table_str.replace('\ue000', 'μ') else: table = pd.DataFrame(table_data, index=pd.Index(['μ (SD)'], name=self.dep), columns=pd.Index([self.group1, self.group2], name=self.ind)) return str(table) def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): return '{}
t({:.0f}) = {:.2f}; p {}
Δμ ({:g}% CI) = {}, {}'.format(self._comparison_table(True), self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML), (1-config.alpha)*100, self.delta.summary(), self.delta_direction) def summary(self): """ Return a stringified summary of the *t* test :rtype: str """ return '{}\n\nt({:.0f}) = {:.2f}; p {}\nΔμ ({:g}% CI) = {}, {}'.format(self._comparison_table(False), self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION), (1-config.alpha)*100, self.delta.summary(), self.delta_direction) def summary_short(self, html): """ Return a stringified summary of the *t* test (*t* statistic only) :rtype: str """ if html: return 't({:.0f}) = {:.2f}'.format(self.dof, self.statistic) else: return 't({:.0f}) = {:.2f}'.format(self.dof, self.statistic) def ttest_ind(df, dep, ind, *, nan_policy='warn'): """ Perform an independent 2-sample Student's *t* test :param df: Data to perform the test on :type df: DataFrame :param dep: Column in *df* for the dependent variable (numeric) :type dep: str :param ind: Column in *df* for the independent variable (dichotomous) :type ind: str :param nan_policy: How to handle *nan* values (see :ref:`nan-handling`) :type nan_policy: str :rtype: :class:`yli.sig_tests.TTestResult` **Example:** .. code-block:: df = pd.DataFrame({ 'Type': ['Fresh'] * 10 + ['Stored'] * 10, 'Potency': [10.2, 10.5, 10.3, ...] }) yli.ttest_ind(df, 'Potency', 'Type') .. code-block:: text Type Fresh Stored Potency μ (SD) 10.37 (0.32) 9.83 (0.24) t(18) = 4.24; p < 0.001* Δμ (95% CI) = 0.54 (0.27–0.81), Fresh > Stored The output states that the value of the *t* statistic is 4.24, the *t* distribution has 18 degrees of freedom, and the test is significant with *p* value < 0.001. The mean difference is 0.54 in favour of the *Fresh* group, with 95% confidence interval 0.27–0.81. """ # Check for/clean NaNs df = check_nan(df[[ind, dep]], nan_policy) # Ensure 2 groups for ind group1, data1, group2, data2 = as_2groups(df, dep, ind) # Do t test # Use statsmodels rather than SciPy because this provides the mean difference automatically d1 = sm.stats.DescrStatsW(data1, ddof=1) d2 = sm.stats.DescrStatsW(data2, ddof=1) cm = sm.stats.CompareMeans(d1, d2) statistic, pvalue, dof = cm.ttest_ind() delta = d1.mean - d2.mean ci0, ci1 = cm.tconfint_diff(config.alpha) # t test is symmetric so take absolute value if d2.mean > d1.mean: delta, ci0, ci1 = -delta, -ci1, -ci0 d1, d2 = d2, d1 group1, group2 = group2, group1 # Now group1 > group2 return TTestResult( statistic=abs(statistic), dof=dof, pvalue=pvalue, dep=dep, ind=ind, group1=group1, group2=group2, mu1=d1.mean, mu2=d2.mean, sd1=d1.std, sd2=d2.std, delta=Estimate(delta, ci0, ci1), delta_direction='{} > {}'.format(group1, group2)) # ------------- # One-way ANOVA class FTestResult: """ Result of an *F* test for ANOVA/regression See :func:`yli.anova_oneway` and :meth:`yli.regress.RegressionResult.ftest`. """ def __init__(self, statistic, dof1, dof2, pvalue): #: *F* statistic (*float*) self.statistic = statistic #: Degrees of freedom in the *F* distribution numerator (*int*) self.dof1 = dof1 #: Degrees of freedom in the *F* distribution denominator (*int*) self.dof2 = dof2 #: *p* value for the *F* statistic (*float*) self.pvalue = pvalue def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): return 'F({:.0f}, {:.0f}) = {:.2f}; p {}'.format(self.dof1, self.dof2, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML)) def summary(self): """ Return a stringified summary of the *F* test :rtype: str """ return 'F({:.0f}, {:.0f}) = {:.2f}; p {}'.format(self.dof1, self.dof2, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION)) def anova_oneway(df, dep, ind, *, nan_policy='warn'): """ Perform one-way ANOVA :param df: Data to perform the test on :type df: DataFrame :param dep: Column in *df* for the dependent variable (numeric) :type dep: str :param ind: Column in *df* for the independent variable (categorical) :type ind: str :param nan_policy: How to handle *nan* values (see :ref:`nan-handling`) :type nan_policy: str :rtype: :class:`yli.sig_tests.FTestResult` **Example:** .. code-block:: df = pd.DataFrame({ 'Method': [1]*8 + [2]*7 + [3]*9, 'Score': [96, 79, 91, ...] }) yli.anova_oneway(df, 'Score', 'Method') .. code-block:: text F(2, 21) = 29.57; p < 0.001* The output states that the value of the *F* statistic is 29.57, the *F* distribution has 2 degrees of freedom in the numerator and 21 in the denominator, and the test is significant with *p* value < 0.001. """ # Check for/clean NaNs df = check_nan(df[[ind, dep]], nan_policy) # Group by independent variable groups = df.groupby(ind)[dep] # Perform one-way ANOVA result = stats.f_oneway(*[groups.get_group(k) for k in groups.groups]) # See stats.f_oneway implementation dfbn = len(groups.groups) - 1 dfwn = len(df) - len(groups.groups) return FTestResult(result.statistic, dfbn, dfwn, result.pvalue) # ----------------- # Mann-Whitney test class MannWhitneyResult: """ Result of a Mann–Whitney *U* test See :func:`yli.mannwhitney`. """ def __init__(self, *, statistic, pvalue, dep, ind, group1, group2, med1, med2, iqr1, iqr2, range1, range2, rank_biserial, direction, brunnermunzel=None): #: Lesser of the two Mann–Whitney *U* statistics (*float*) self.statistic = statistic #: *p* value for the *U* statistic (*float*) self.pvalue = pvalue #: Name of the dependent variable (*str*) self.dep = dep #: Name of the independent variable (*str*) self.ind = ind #: Name of the first group (*str*) self.group1 = group1 #: Name of the second group (*str*) self.group2 = group2 #: Median of the first group (*float*) self.med1 = med1 #: Median of the second group (*float*) self.med2 = med2 #: Interquartile range of the first group (:class:`yli.utils.Interval`) self.iqr1 = iqr1 #: Interquartile range of the second group (:class:`yli.utils.Interval`) self.iqr2 = iqr2 #: Range of the first group (:class:`yli.utils.Interval`) self.range1 = range2 #: Range of the second group (:class:`yli.utils.Interval`) self.range2 = range2 #: Absolute value of the rank-biserial correlation (*float*) self.rank_biserial = rank_biserial #: Description of the direction of the effect (*str*) self.direction = direction #: :class:`BrunnerMunzelResult` on the same data, or *None* if N/A self.brunnermunzel = brunnermunzel def _comparison_table(self, html): """Return a table showing the medians/IQRs/ranges for each group""" table_data = [ ['{:.2f} ({})'.format(self.med1, self.iqr1.summary()), '{:.2f} ({})'.format(self.med2, self.iqr2.summary())], ['{:.2f} ({})'.format(self.med1, self.range1.summary()), '{:.2f} ({})'.format(self.med2, self.range2.summary())], ] table = pd.DataFrame(table_data, index=pd.Index(['Median (IQR)', 'Median (range)'], name=self.dep), columns=pd.Index([self.group1, self.group2], name=self.ind)) if html: return table._repr_html_() else: return str(table) def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): line1 = '{}
U = {:.1f}; p {}
r = {:.2f}, {}'.format(self._comparison_table(True), self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML), self.rank_biserial, self.direction) if self.brunnermunzel: return line1 + '
' + self.brunnermunzel._repr_html_() else: return line1 def summary(self): """ Return a stringified summary of the Mann–Whitney test :rtype: str """ line1 = '{}\n\nU = {:.1f}; p {}\nr = {:.2f}, {}'.format(self._comparison_table(False), self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION), self.rank_biserial, self.direction) if self.brunnermunzel: return line1 + '\n' + self.brunnermunzel.summary() else: return line1 def summary_short(self, html): """ Return a stringified summary of the Mann–Whitney test (*U* statistic only) :rtype: str """ if html: return 'U = {:.1f}'.format(self.statistic) else: return 'U = {:.1f}'.format(self.statistic) class BrunnerMunzelResult: """ Result of a Brunner–Munzel test See :func:`yli.mannwhitney`. This library calls the Brunner–Munzel test statistic *W*. """ """Result of a Brunner-Munzel test""" def __init__(self, statistic, pvalue): #: *W* statistic (*float*) self.statistic = statistic #: *p* value for the *W* statistic (*float*) self.pvalue = pvalue def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): return 'W = {:.1f}; p {}'.format(self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML)) def summary(self): """ Return a stringified summary of the Brunner–Munzel test :rtype: str """ return 'W = {:.1f}; p {}'.format(self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION)) def mannwhitney(df, dep, ind, *, nan_policy='warn', brunnermunzel=True, use_continuity=False, alternative='two-sided', method='auto'): """ Perform a Mann–Whitney *U* test By default, this function performs a Brunner–Munzel test if the Mann–Whitney test is significant. If the Mann–Whitney test is significant but the Brunner–Munzel test is not, a warning is raised. The Brunner–Munzel test is returned only if non-significant. :param df: Data to perform the test on :type df: DataFrame :param dep: Column in *df* for the dependent variable (numeric) :type dep: str :param ind: Column in *df* for the independent variable (dichotomous) :type ind: str :param nan_policy: How to handle *nan* values (see :ref:`nan-handling`) :type nan_policy: str :param brunnermunzel: Whether to compute the Brunner–Munzel test if the Mann–Whitney test is significant :type brunnermunzel: bool :param use_continuity: See *scipy.stats.mannwhitneyu* :param alternative: See *scipy.stats.mannwhitneyu* :param method: See *scipy.stats.mannwhitneyu* :return: The result of the Mann–Whitney test. The result of a Brunner–Munzel test is included in the result object if and only if *brunnermunzel* is *True*, *and* the Mann–Whitney test is significant, *and* the Brunner–Munzel test is non-significant. :rtype: :class:`yli.sig_tests.MannWhitneyResult` **Example:** .. code-block:: df = pd.DataFrame({ 'Sample': ['Before'] * 12 + ['After'] * 12, 'Oxygen': [11.0, 11.2, 11.2, ...] }) yli.mannwhitney(df, 'Oxygen', 'Sample', method='asymptotic', alternative='less') .. code-block:: text Sample After Before Oxygen Median (IQR) 10.75 (10.55–10.95) 11.55 (11.20–11.83) Median (range) 10.75 (11.00–12.10) 11.55 (11.00–12.10) U = 6.0; p < 0.001* r = 0.92, Before > After The output states that the value of the Mann–Whitney *U* statistic is 6.0, and the one-sided test is significant with asymptotic *p* value < 0.001. The rank-biserial correlation is 0.92 in favour of the *Before* group. """ # Check for/clean NaNs df = check_nan(df[[ind, dep]], nan_policy) # Convert pandas nullable types for independent variables as this breaks mannwhitneyu df = convert_pandas_nullable(df) # Ensure 2 groups for ind group1, data1, group2, data2 = as_2groups(df, dep, ind) # Do Mann-Whitney test # Stata does not perform continuity correction result = stats.mannwhitneyu(data1, data2, use_continuity=use_continuity, alternative=alternative, method=method) u1 = result.statistic u2 = len(data1) * len(data2) - u1 r = abs(2*u1 / (len(data1) * len(data2)) - 1) # rank-biserial # If significant, perform a Brunner-Munzel test for our interest if result.pvalue < 0.05 and brunnermunzel: result_bm = stats.brunnermunzel(data1, data2) if result_bm.pvalue >= 0.05: warnings.warn('Mann-Whitney test is significant but Brunner-Munzel test is not. This could be due to a difference in shape, rather than location.') return MannWhitneyResult( statistic=min(u1, u2), pvalue=result.pvalue, #med1=data1.median(), med2=data2.median(), rank_biserial=r, direction=('{1} > {0}' if u1 < u2 else '{0} > {1}').format(group1, group2), brunnermunzel=BrunnerMunzelResult(statistic=result_bm.statistic, pvalue=result_bm.pvalue)) return MannWhitneyResult( statistic=min(u1, u2), pvalue=result.pvalue, dep=dep, ind=ind, group1=group1, group2=group2, med1=data1.median(), med2=data2.median(), iqr1=Interval(data1.quantile(0.25), data1.quantile(0.75)), iqr2=Interval(data2.quantile(0.25), data2.quantile(0.75)), range1=Interval(data1.min(), data1.max()), range2=Interval(data2.min(), data2.max()), rank_biserial=r, direction=('{1} > {0}' if u1 < u2 else '{0} > {1}').format(group1, group2)) # ------------------------ # Pearson chi-squared test class ChiSquaredResult: # TODO: Documentation def __init__(self, statistic, dof, pvalue): #: *χ*:sup:`2` statistic (*float*) self.statistic = statistic #: Degrees of freedom for the *χ*:sup:`2` distribution (*int*) self.dof = dof #: *p* value for the *χ*:sup:`2` test (*float*) self.pvalue = pvalue def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): return 'χ2({}) = {:.2f}; p {}'.format(self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML)) def summary(self): """ Return a stringified summary of the *χ*:sup:`2` test :rtype: str """ return 'χ²({}) = {:.2f}; p {}'.format(self.ct, self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION)) class PearsonChiSquaredResult(ChiSquaredResult): """ Result of a Pearson *χ*:sup:`2` test See :func:`yli.chi2`. """ def __init__(self, ct, statistic, dof, pvalue, oddsratio=None, riskratio=None): super().__init__(statistic, dof, pvalue) #: Contingency table for the observations (*DataFrame*) self.ct = ct #: Odds ratio (*float*; *None* if not a 2×2 table) self.oddsratio = oddsratio #: Risk ratio (*float*; *None* if not a 2×2 table) self.riskratio = riskratio def _repr_html_(self): if self.oddsratio is not None: return '{0}
χ2({1}) = {2:.2f}; p {3}
OR ({4:g}% CI) = {5}
RR ({4:g}% CI) = {6}'.format( self.ct._repr_html_(), self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML), (1-config.alpha)*100, self.oddsratio.summary(), self.riskratio.summary()) else: return '{}
χ2({}) = {:.2f}; p {}'.format( self.ct._repr_html_(), self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML)) def summary(self): """ Return a stringified summary of the *χ*:sup:`2` test :rtype: str """ if self.oddsratio is not None: return '{0}\n\nχ²({1}) = {2:.2f}; p {3}\nOR ({4:g}% CI) = {5}\nRR ({4:g}% CI) = {6}'.format( self.ct, self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION), (1-config.alpha)*100, self.oddsratio.summary(), self.riskratio.summary()) else: return '{}\n\nχ²({}) = {:.2f}; p {}'.format( self.ct, self.dof, self.statistic, fmt_p(self.pvalue, PValueStyle.RELATION)) def summary_short(self, html): """ Return a stringified summary of the *χ*:sup:`2` test (*χ*:sup:`2` statistic only) :rtype: str """ if html: return 'χ2({}) = {:.2f}'.format(self.dof, self.statistic) else: return 'χ²({}) = {:.2f}'.format(self.dof, self.statistic) def chi2(df, dep, ind, *, nan_policy='warn'): """ Perform a Pearson *χ*:sup:`2` test If a 2×2 contingency table is obtained (i.e. if both variables are dichotomous), an odds ratio and risk ratio are calculated. The ratios are calculated for the higher-valued value in each variable (i.e. *True* compared with *False* for a boolean). The risk ratio is calculated relative to the independent variable (rows of the contingency table). :param df: Data to perform the test on :type df: DataFrame :param dep: Column in *df* for the dependent variable (categorical) :type dep: str :param ind: Column in *df* for the independent variable (categorical) :type ind: str :param nan_policy: How to handle *nan* values (see :ref:`nan-handling`) :type nan_policy: str :rtype: :class:`yli.sig_tests.PearsonChiSquaredResult` **Example:** .. code-block:: df = pd.DataFrame({ 'Response': np.repeat([False, True, False, True], [250, 750, 400, 1600]), 'Stress': np.repeat([False, False, True, True], [250, 750, 400, 1600]) }) yli.chi2(df, 'Stress', 'Response') .. code-block:: text Stress False True Response False 250 400 True 750 1600 χ²(1) = 9.82; p = 0.002* OR (95% CI) = 1.33 (1.11–1.60) RR (95% CI) = 1.11 (1.03–1.18) The output shows the contingency table, and states that the value of the Pearson *χ*:sup:`2` statistic is 9.82, the *χ*:sup:`2` distribution has 1 degree of freedom, and the test is significant with *p* value 0.002. The odds of *Stress* in the *Response* = *True* group are 1.33 times that in the *Response* = *False* group, with 95% confidence interval 1.11–1.60. The risk of *Stress* in the *Response* = *True* group is 1.11 that in the *Response* = *False* group, with 95% confidence interval 1.03–1.18. """ # Check for/clean NaNs df = check_nan(df[[ind, dep]], nan_policy) # Compute contingency table ct = pd.crosstab(df[ind], df[dep]) # Get expected counts expected = stats.contingency.expected_freq(ct) # Warn on low expected counts if (expected < 5).sum() / expected.size > 0.2: warnings.warn('{} of {} cells ({:.0f}%) have expected count < 5'.format((expected < 5).sum(), expected.size, (expected < 5).sum() / expected.size * 100)) if (expected < 1).any(): warnings.warn('{} cells have expected count < 1'.format((expected < 1).sum())) if ct.shape == (2,2): # 2x2 table # Use statsmodels to get OR and RR smct = sm.stats.Table2x2(np.flip(ct.to_numpy()), shift_zeros=False) result = smct.test_nominal_association() ORci = smct.oddsratio_confint(config.alpha) RRci = smct.riskratio_confint(config.alpha) return PearsonChiSquaredResult( ct=ct, statistic=result.statistic, dof=result.df, pvalue=result.pvalue, oddsratio=Estimate(smct.oddsratio, ORci[0], ORci[1]), riskratio=Estimate(smct.riskratio, RRci[0], RRci[1])) else: # rxc table # Just use SciPy result = stats.chi2_contingency(ct, correction=False) return PearsonChiSquaredResult(ct=ct, statistic=result[0], dof=result[2], pvalue=result[1]) # ------------------- # Pearson correlation class PearsonRResult: """ Result of Pearson correlation See :func:`yli.pearsonr`. """ def __init__(self, statistic, pvalue): #: Pearson *r* correlation statistic (*float*) self.statistic = statistic #: *p* value for the *r* statistic (*float*) self.pvalue = pvalue def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): return 'r ({:g}% CI) = {}; p {}'.format((1-config.alpha)*100, self.statistic.summary(), fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML)) def summary(self): """ Return a stringified summary of the Pearson correlation :rtype: str """ return 'r ({:g}% CI) = {}; p {}'.format((1-config.alpha)*100, self.statistic.summary(), fmt_p(self.pvalue, PValueStyle.RELATION)) def pearsonr(df, dep, ind, *, nan_policy='warn'): """ Compute the Pearson correlation coefficient (Pearson's *r*) :param df: Data to perform the test on :type df: DataFrame :param dep: Column in *df* for the dependent variable (numerical) :type dep: str :param ind: Column in *df* for the independent variable (numerical) :type ind: str :param nan_policy: How to handle *nan* values (see :ref:`nan-handling`) :type nan_policy: str :rtype: :class:`yli.sig_tests.PearsonRResult` **Example:** .. code-block:: df = pd.DataFrame({ 'y': [41, 39, 47, 51, 43, 40, 57, 46, 50, 59, 61, 52], 'x': [24, 30, 33, 35, 36, 36, 37, 37, 38, 40, 43, 49] }) yli.pearsonr(df, 'y', 'x') .. code-block:: text r (95% CI) = 0.65 (0.11–0.89); p = 0.02* The output states that the value of the Pearson correlation coefficient is 0.65, with 95% confidence interval 0.11–0.89, and the test is significant with *p* value 0.02. """ # Check for/clean NaNs df = check_nan(df[[ind, dep]], nan_policy) # Compute Pearson's r result = stats.pearsonr(df[ind], df[dep]) ci = result.confidence_interval() return PearsonRResult(statistic=Estimate(result.statistic, ci.low, ci.high), pvalue=result.pvalue) # -------------------- # Spearman correlation class SpearmanResult: def __init__(self, statistic, pvalue): self.statistic = statistic self.pvalue = pvalue def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): return 'ρ ({:g}% CI) = {}; p {}'.format((1-config.alpha)*100, self.statistic.summary(), fmt_p(self.pvalue, PValueStyle.RELATION | PValueStyle.HTML)) def summary(self): """ Return a stringified summary of the Spearman correlation :rtype: str """ return 'ρ ({:g}% CI) = {}; p {}'.format((1-config.alpha)*100, self.statistic.summary(), fmt_p(self.pvalue, PValueStyle.RELATION)) def spearman(df, dep, ind, *, nan_policy='warn'): # Check for/clean NaNs df = check_nan(df[[ind, dep]], nan_policy) # Ensure numeric, factorising categorical variables as required ind, _ = as_numeric(df[ind]) dep, _ = as_numeric(df[dep]) # Compute Spearman's rho result = stats.spearmanr(ind, dep) # Compute confidence interval ci = stats._stats_py._pearsonr_fisher_ci(result.correlation, len(dep), 1 - config.alpha, 'two-sided') return SpearmanResult(statistic=Estimate(result.correlation, ci.low, ci.high), pvalue=result.pvalue) # ---------------------------- # Automatic selection of tests class AutoBinaryResult: """ Result of automatically computed univariable tests of association for a dichotomous dependent variable See :func:`yli.auto_univariable`. Results data stored within instances of this class is not intended to be directly accessed. """ def __init__(self, *, dep, group1, group2, result_data, result_labels): #: Name of the dependent variable (*str*) self.dep = dep #: Name of the first group (*str*) self.group1 = group1 #: Name of the second group (*str*) self.group2 = group2 # List of tuples (first group summary, second group summary, test result) self._result_data = result_data # List of tuples (plaintext label, HTML label) self._result_labels = result_labels def __repr__(self): if config.repr_is_summary: return self.summary() return super().__repr__() def _repr_html_(self): result = ''.format(self.dep, self.group1, self.group2) for data, label in zip(self._result_data, self._result_labels): result += ''.format(label[1], data[0], data[1], data[2].summary_short(True), fmt_p(data[2].pvalue, PValueStyle.TABULAR | PValueStyle.HTML)) result += '
{}{}{}p
{}{}{}{}{}
' return result def summary(self): """ Return a stringified summary of the tests of association :rtype: str """ # Format data for output result_data_fmt = [ [r[0], r[1], r[2].summary_short(False), fmt_p(r[2].pvalue, PValueStyle.TABULAR)] for r in self._result_data ] result_labels_fmt = [r[0] for r in self._result_labels] table = pd.DataFrame(result_data_fmt, index=result_labels_fmt, columns=pd.Index([self.group1, self.group2, '', 'p'], name=self.dep)) return str(table) def auto_univariable(df, dep, inds, *, ordinal=[], nan_policy='warn'): """ Automatically compute univariable tests of association for a dichotomous dependent variable The tests performed are: * For a dichotomous independent variable – :func:`yli.chi2` * For a continuous independent variable – :func:`yli.ttest_ind` * For an ordinal independent variable – :func:`yli.mannwhitney` If *nan_policy* is *warn* or *omit*, rows with *nan* values are omitted only from the individual tests of association for the missing variables. :param df: Data to perform the test on :type df: DataFrame :param dep: Column in *df* for the dependent variable (dichotomous) :type dep: str :param inds: Columns in *df* for the independent variables :type inds: List[str] :param ordinal: Columns in *df* to treat as ordinal rather than continuous :type ordinal: List[str] :param nan_policy: How to handle *nan* values (see :ref:`nan-handling`) :type nan_policy: str :rtype: :class:`yli.sig_tests.AutoBinaryResult` """ # Check for/clean NaNs in dependent variable df = check_nan(df[inds + [dep]], nan_policy, cols=[dep]) # Ensure 2 groups for dep # TODO: Work for non-binary dependent variables? group1, data1, group2, data2 = as_2groups(df, inds, dep) result_data = [] result_labels = [] for ind in inds: # Check for/clean NaNs in independent variable # Following this, we pass nan_policy='raise' to assert no NaNs remaining df_cleaned = check_nan(df, nan_policy, cols=[ind]) if df_cleaned[ind].dtype in ('bool', 'boolean', 'category', 'object'): # Pearson chi-squared test result = chi2(df_cleaned, dep, ind, nan_policy='raise') values = sorted(df_cleaned[ind].unique()) # Value counts result_labels.append(( '{}, {}'.format(ind, ':'.join(str(v) for v in values)), '{}, {}'.format(ind, ':'.join(str(v) for v in values)), )) result_data.append(( ':'.join(str((data1[ind] == v).sum()) for v in values), ':'.join(str((data2[ind] == v).sum()) for v in values), result )) elif df_cleaned[ind].dtype in ('float64', 'int64', 'Float64', 'Int64'): if ind in ordinal: # Mann-Whitney test result = mannwhitney(df_cleaned, ind, dep, nan_policy='raise') result_labels.append(( '{}, median (IQR)'.format(ind), '{}, median (IQR)'.format(ind), )) result_data.append(( '{:.2f} ({})'.format(result.med1, result.iqr1.summary()), '{:.2f} ({})'.format(result.med2, result.iqr2.summary()), result )) else: # t test result = ttest_ind(df_cleaned, ind, dep, nan_policy='raise') result_labels.append(( '{}, μ (SD)'.format(ind), '{}, μ (SD)'.format(ind), )) result_data.append(( '{:.2f} ({:.2f})'.format(result.mu1, result.sd1), '{:.2f} ({:.2f})'.format(result.mu2, result.sd2), result )) else: raise Exception('Unsupported independent dtype for auto_univariable, {}'.format(df[ind].dtype)) return AutoBinaryResult(dep=dep, group1=group1, group2=group2, result_data=result_data, result_labels=result_labels)