Source code for twosample

"""
.. py:currentmodule:: twosample

.. module:: twosample
   :platform: Unix, Windows
   :synopsis: Module to calculate two sample tests: Brunner-Munzel, Wilcoxon-Mann-Whitney

.. moduleauthor:: Martin Happ <martin.happ@aon.at>
"""

import numpy as np
import math
import scipy
import scipy.stats
import scipy.special
from collections import namedtuple
import PyNonpar.pseudorank
import functools
from functools import lru_cache


[docs]def brunner_munzel_test(x: list, y: list, alternative="two.sided", quantile="t") -> list: """ Function to calculate the Brunner-Munzel test. It is recommended to use the t-approximation for small samples. Args: x (list(float)): data from first group \n y (list(float)): data from second group \n alternative (str): either 'two.sided', 'less' or 'greater' \n quantile (str): either 't' for the small sample approximation with a t-distribution or 'normal' for the large sample procedure. \n Returns: namedtuple('BrunnerMunzelResult', ('alternative', 'statistic', 'df', 'pvalue')): \n chosen alternative (str) \n test statistic (float)\n estimated degrees of freedom for the t-approximation, only when quantile = 't' (float) \n p-value (float) """ # Check inputs if not isinstance(x, list): raise TypeError("x must be a list") if not isinstance(y, list): raise TypeError("y must be a list") if (not isinstance(alternative, str)) or (alternative not in ['two.sided', 'less', 'greater']): raise TypeError('alternative must be either two.sided, less or greater') if (not isinstance(quantile, str)) or (quantile not in ['t', 'normal']): raise TypeError('quantile must be either t or normal') n = (len(x), len(y)) N = n[0] + n[1] data = x + y grp = [1 for i in range(len(data))] ranks = PyNonpar.pseudorank.psrank(data, grp) x = ranks[0:n[0]] y = ranks[n[0]:N] xi = PyNonpar.pseudorank.psrank(x, grp[0:n[0]]) yi = PyNonpar.pseudorank.psrank(y, grp[0:n[1]]) R2 = np.mean(y) R1 = np.mean(x) # calculate variance estimator S1_squared = 0.0 for i in range(n[0]): S1_squared += (x[i] - xi[i] - R1 + (n[0] + 1) * 1 / 2) ** 2 S1_squared = S1_squared * 1 / (n[0] - 1) S2_squared = 0.0 for i in range(n[1]): S2_squared += (y[i] - yi[i] - R2 + (n[1] + 1) * 1 / 2) ** 2 S2_squared = S2_squared * 1 / (n[1] - 1) sigma_BF_squared = N * S1_squared * 1 / (N - n[0]) + N * S2_squared * 1 / (N - n[1]) # calculate test statistic W_N_BF = (R2 - R1) * 1 / math.sqrt(sigma_BF_squared) * math.sqrt(n[0] * n[1] * 1 / N) p_value = 0.0 if quantile == "normal": if alternative == "greater": p_value = scipy.stats.norm.cdf(W_N_BF) elif alternative == "less": p_value = 1 - scipy.stats.norm.cdf(W_N_BF) elif alternative == "two.sided": p_value = 2 * min(scipy.stats.norm.cdf(W_N_BF), 1 - scipy.stats.norm.cdf(W_N_BF)) result = namedtuple('BrunnerMunzelResult', ('alternative', 'statistic', 'pvalue')) output = result(alternative, W_N_BF, p_value) else: f_hat = ((1 / N * sigma_BF_squared) ** 2) * 1 / ( (S1_squared * 1 / (N - n[0])) ** 2 * 1 / (n[0] - 1) + (S2_squared * 1 / (N - n[1])) ** 2 * 1 / ( n[1] - 1)) if alternative == "greater": p_value = scipy.stats.t.cdf(W_N_BF, f_hat) elif alternative == "less": p_value = 1 - scipy.stats.t.cdf(W_N_BF, f_hat) elif alternative == "two.sided": p_value = 2 * min(scipy.stats.t.cdf(W_N_BF, f_hat), 1 - scipy.stats.t.cdf(W_N_BF, f_hat)) result = namedtuple('BrunnerMunzelResult', ('alternative', 'statistic', 'df', 'pvalue')) output = result(alternative, W_N_BF, f_hat, p_value) return output
[docs]def hodges_lehmann(x: list, y: list, alpha = 0.05): """ Function to calculate the Hodges-Lehmann estimator.\n Let us assume a location shift effect, i.e., F_1(x) = F(x - mu_1), i = 1,2 and theta = mu_2 - mu_1. \n The Hodges-Lehmann estimator theta_hat is an asymptotically unbiased estimator for theta. \n Note that the calculated confidence interval is only valid if there are no ties in the data. Args: x (list(float)): data from first group \n y (list(float)): data from second group \n Returns: namedtuple('HodgesLehmannEstimator', ('estimate', 'lowerCI', 'upperCI')):\n Hodges-Lehmann estimator theta_hat (float) \n lower bound for the asymptotic 1-alpha/2 confidence interval for theta (float)\n upper bound for the 1-alpha/2 confidence interval for theta (float) \n """ n = len(x) m = len(y) diff = [0 for i in range(n * m)] for i in range(n): for j in range(m): diff[i * m + j] = y[j] - x[i] hl = np.median(diff) # asymptotic 1-alpha/2 confidence interval L = 1 + n*m/2 - scipy.stats.norm.ppf(1-alpha/2)*math.sqrt(n*m*(n + m + 1)*1/12) U = n*m/2 + scipy.stats.norm.ppf(1-alpha/2)*math.sqrt(n*m*(n + m + 1)*1/12) diff.sort() lower = diff[int(round(L) - 1)] upper = diff[int(round(U) - 1)] result = namedtuple('HodgesLehmannEstimator', ('estimate', 'lowerCI', 'upperCI')) output = result(hl, lower, upper) return output
@functools.lru_cache(maxsize=1000) def _h(s, m, N, ranks): r = 0.0 partial_sum_ranks = 0.0 for i in range(m): partial_sum_ranks += ranks[i] if s < 0: r = 0 elif m == N: if s == partial_sum_ranks: r = 1 else: r = 0 elif m == 0: if s == 0: r = 1 else: r = 0 else: r = _h(s,m,N-1,ranks) + _h(s-ranks[N-1],m-1,N-1,ranks) return r
[docs]def wilcoxon_mann_whitney_test(x: list, y: list, alternative = "two.sided", alpha = 0.05, method = "asymptotic"): """ Function to calculate the Wilcoxon-Mann-Whitney test. Args: x (list(float)): data from first group \n y (list(float)): data from second group \n alternative (str): either 'two.sided', 'less' (x less y) or 'greater' (x greater y) \n alpha (float): 1-alpha confidence interval (only valid if there are no ties) \n method (str): use 'asymptotic' for the asymptotic test (normal distribution) or the exact test with 'exact'. Returns: namedtuple('WilcoxonMannWhitneyResult', ('alternative', 'statistic', 'HodgesLehmann', 'lowerCI', 'upperCI', 'pvalue')): \n chosen alternative (str) \n test statistic (float)\n hodges lehmann estimate for a location shift effect y - x (float) \n lower bound for 1-alpha CI for location shift effect (float) \n upper bound for 1-alpha CI for location shift effect (float) \n p-value """ n = len(x) m = len(y) data = x + y N = len(data) hl = hodges_lehmann(x, y, alpha) ranks = [0 for i in range(N)] ranks = PyNonpar.pseudorank.psrank(data, ranks) x = ranks[:n] y = ranks[n:] R1 = np.mean(x) R2 = np.mean(y) # case asymptotic test if method == "asymptotic": # variance estimator sigma_R_squared = 0.0 for i in range(N): sigma_R_squared += (ranks[i] - (N + 1)*1/2 ) ** 2 sigma_R_squared = sigma_R_squared*1/(N - 1) W_N = math.sqrt(n*m*1/N)*(R2 - R1)*1/math.sqrt(sigma_R_squared) if alternative == "greater": p_value = scipy.stats.norm.cdf(W_N) elif alternative == "less": p_value = 1 - scipy.stats.norm.cdf(W_N) elif alternative == "two.sided": p_value = 2 * min(scipy.stats.norm.cdf(W_N), 1 - scipy.stats.norm.cdf(W_N)) result = namedtuple('WilcoxonMannWhitneyResult', ('alternative', 'statistic', 'HodgesLehmann', 'lowerCI', 'upperCI', 'pvalue')) output = result(alternative, W_N, hl[0], hl[1], hl[2], p_value) elif method == "exact": R_W = np.sum(y) combinations = 0.0 total_combinations = scipy.special.binom(N, m) ranks.sort() ranks2 = [2*ranks[i] for i in range(N)] ranks2 = tuple(ranks2) for i in range(int(2*R_W)): combinations += _h(i, m, N, ranks2) if alternative == "two.sided": p_value = 2*min(combinations/total_combinations + _h(2*R_W,m,N,ranks2)/total_combinations, 1 - combinations/total_combinations) elif alternative == "greater": p_value = combinations / total_combinations + _h(2*R_W,m,N,ranks2)/total_combinations elif alternative == "less": p_value = 1 - combinations / total_combinations result = namedtuple('WilcoxonMannWhitneyResult', ('alternative', 'statistic', 'HodgesLehmann', 'lowerCI', 'upperCI', 'pvalue')) output = result(alternative, R_W, hl[0], hl[1], hl[2], p_value) return output
[docs]def wilcoxon_mann_whitney_ssp(x: list, y: list, power = 0.8, alpha = 0.05, t = 1/2): """ Function to do sample size planning for the WMW test. Args: x (list(float)): prior information from first group \n y (list(float)): prior information from second group \n power (float): probablity to detect an effect based on prior information \n alpha (float): type-I error probability \n t (float): ratio of subjects assigned to first group. Returns: namedtuple('WilcoxonMannWhitneyResult', ('alpha', 'power', 'relEffect', 'N', 't', 'n1', 'n2', 'Nrounded', 'n1rouned', 'n2rouned')): \n alpha (float) \n power(float)\n relEffect (float): calculated relative Effect \n N (float): minimal sample size required \n n1 (float): sample size for group 1 required; N = t*n1 \n n2 (float): sample size for group 2 required; N = (1-t)*n1 """ m1 = len(x) m2 = len(y) data = x + y N = len(data) ranks = [0 for i in range(N)] ranks = PyNonpar.pseudorank.psrank(data, ranks) R1 = ranks[:m1] R2 = ranks[m1:] # ranks within samples: R11 = PyNonpar.pseudorank.psrank(x, [0 for i in range(m1)]) R22 = PyNonpar.pseudorank.psrank(y, [0 for i in range(m2)]) # placements P1 = [0 for i in range(m1)] P2 = [0 for i in range(m2)] for i in range(m1): P1[i] = R1[i] - R11[i] for i in range(m2): P2[i] = R2[i] - R22[i] # effect size: pStar = (np.mean(R2) - np.mean(R1)) / (m1 + m2) + 0.5 # variances: sigmaStar = 0.0 for i in range(m1+m2): sigmaStar += (ranks[i] - (m1 + m2) / 2) ** 2 sigmaStar = np.sqrt( sigmaStar * 1/(m1 + m2) ** 3 ) sigma1Star = 0.0 for i in range(m1): sigma1Star += (P1[i] - np.mean(P1)) ** 2 sigma1Star = np.sqrt( sigma1Star * 1/(m1 * m2 * m2) ) sigma2Star = 0.0 for i in range(m2): sigma2Star += (P2[i] - np.mean(P2)) ** 2 sigma2Star = np.sqrt( sigma2Star * 1/(m1 * m1 * m2) ) # estimated sample size: N = (sigmaStar * scipy.stats.norm.ppf(1 - alpha / 2) + scipy.stats.norm.ppf(power) * np.sqrt(t * sigma2Star ** 2 + (1 - t) * sigma1Star ** 2)) ** 2 * 1/ (t * (1 - t) * (pStar - 0.5) ** 2) n1 = N * t n2 = N * (1 - t) result = namedtuple('WilcoxonMannWhitneySampleSize', ('alpha', 'power', 'relEffect', 'N', 't', 'n1', 'n2', 'Nrounded', 'n1rouned', 'n2rouned')) output = result(alpha, power, pStar, N, t, n1, n2, math.ceil(N), math.ceil(n1), math.ceil(n2)) return output