1''' 2 3from pystatsmodels mailinglist 20100524 4 5Notes: 6 - unfinished, unverified, but most parts seem to work in MonteCarlo 7 - one example taken from lecture notes looks ok 8 - needs cases with non-monotonic inequality for test to see difference between 9 one-step, step-up and step-down procedures 10 - FDR does not look really better then Bonferoni in the MC examples that I tried 11update: 12 - now tested against R, stats and multtest, 13 I have all of their methods for p-value correction 14 - getting Hommel was impossible until I found reference for pvalue correction 15 - now, since I have p-values correction, some of the original tests (rej/norej) 16 implementation is not really needed anymore. I think I keep it for reference. 17 Test procedure for Hommel in development session log 18 - I have not updated other functions and classes in here. 19 - multtest has some good helper function according to docs 20 - still need to update references, the real papers 21 - fdr with estimated true hypothesis still missing 22 - multiple comparison procedures incomplete or missing 23 - I will get multiple comparison for now only for independent case, which might 24 be conservative in correlated case (?). 25 26 27some References: 28 29Gibbons, Jean Dickinson and Chakraborti Subhabrata, 2003, Nonparametric Statistical 30Inference, Fourth Edition, Marcel Dekker 31 p.363: 10.4 THE KRUSKAL-WALLIS ONE-WAY ANOVA TEST AND MULTIPLE COMPARISONS 32 p.367: multiple comparison for kruskal formula used in multicomp.kruskal 33 34Sheskin, David J., 2004, Handbook of Parametric and Nonparametric Statistical 35Procedures, 3rd ed., Chapman&Hall/CRC 36 Test 21: The Single-Factor Between-Subjects Analysis of Variance 37 Test 22: The Kruskal-Wallis One-Way Analysis of Variance by Ranks Test 38 39Zwillinger, Daniel and Stephen Kokoska, 2000, CRC standard probability and 40statistics tables and formulae, Chapman&Hall/CRC 41 14.9 WILCOXON RANKSUM (MANN WHITNEY) TEST 42 43 44S. Paul Wright, Adjusted P-Values for Simultaneous Inference, Biometrics 45 Vol. 48, No. 4 (Dec., 1992), pp. 1005-1013, International Biometric Society 46 Stable URL: http://www.jstor.org/stable/2532694 47 (p-value correction for Hommel in appendix) 48 49for multicomparison 50 51new book "multiple comparison in R" 52Hsu is a good reference but I do not have it. 53 54 55Author: Josef Pktd and example from H Raja and rewrite from Vincent Davis 56 57 58TODO 59---- 60* name of function multipletests, rename to something like pvalue_correction? 61 62 63''' 64from statsmodels.compat.python import lzip, lrange 65 66import copy 67import math 68 69import numpy as np 70from numpy.testing import assert_almost_equal, assert_equal 71from scipy import stats, interpolate 72 73from statsmodels.iolib.table import SimpleTable 74#temporary circular import 75from statsmodels.stats.multitest import multipletests, _ecdf as ecdf, fdrcorrection as fdrcorrection0, fdrcorrection_twostage 76from statsmodels.graphics import utils 77from statsmodels.tools.sm_exceptions import ValueWarning 78 79qcrit = ''' 80 2 3 4 5 6 7 8 9 10 815 3.64 5.70 4.60 6.98 5.22 7.80 5.67 8.42 6.03 8.91 6.33 9.32 6.58 9.67 6.80 9.97 6.99 10.24 826 3.46 5.24 4.34 6.33 4.90 7.03 5.30 7.56 5.63 7.97 5.90 8.32 6.12 8.61 6.32 8.87 6.49 9.10 837 3.34 4.95 4.16 5.92 4.68 6.54 5.06 7.01 5.36 7.37 5.61 7.68 5.82 7.94 6.00 8.17 6.16 8.37 848 3.26 4.75 4.04 5.64 4.53 6.20 4.89 6.62 5.17 6.96 5.40 7.24 5.60 7.47 5.77 7.68 5.92 7.86 859 3.20 4.60 3.95 5.43 4.41 5.96 4.76 6.35 5.02 6.66 5.24 6.91 5.43 7.13 5.59 7.33 5.74 7.49 8610 3.15 4.48 3.88 5.27 4.33 5.77 4.65 6.14 4.91 6.43 5.12 6.67 5.30 6.87 5.46 7.05 5.60 7.21 8711 3.11 4.39 3.82 5.15 4.26 5.62 4.57 5.97 4.82 6.25 5.03 6.48 5.20 6.67 5.35 6.84 5.49 6.99 8812 3.08 4.32 3.77 5.05 4.20 5.50 4.51 5.84 4.75 6.10 4.95 6.32 5.12 6.51 5.27 6.67 5.39 6.81 8913 3.06 4.26 3.73 4.96 4.15 5.40 4.45 5.73 4.69 5.98 4.88 6.19 5.05 6.37 5.19 6.53 5.32 6.67 9014 3.03 4.21 3.70 4.89 4.11 5.32 4.41 5.63 4.64 5.88 4.83 6.08 4.99 6.26 5.13 6.41 5.25 6.54 9115 3.01 4.17 3.67 4.84 4.08 5.25 4.37 5.56 4.59 5.80 4.78 5.99 4.94 6.16 5.08 6.31 5.20 6.44 9216 3.00 4.13 3.65 4.79 4.05 5.19 4.33 5.49 4.56 5.72 4.74 5.92 4.90 6.08 5.03 6.22 5.15 6.35 9317 2.98 4.10 3.63 4.74 4.02 5.14 4.30 5.43 4.52 5.66 4.70 5.85 4.86 6.01 4.99 6.15 5.11 6.27 9418 2.97 4.07 3.61 4.70 4.00 5.09 4.28 5.38 4.49 5.60 4.67 5.79 4.82 5.94 4.96 6.08 5.07 6.20 9519 2.96 4.05 3.59 4.67 3.98 5.05 4.25 5.33 4.47 5.55 4.65 5.73 4.79 5.89 4.92 6.02 5.04 6.14 9620 2.95 4.02 3.58 4.64 3.96 5.02 4.23 5.29 4.45 5.51 4.62 5.69 4.77 5.84 4.90 5.97 5.01 6.09 9724 2.92 3.96 3.53 4.55 3.90 4.91 4.17 5.17 4.37 5.37 4.54 5.54 4.68 5.69 4.81 5.81 4.92 5.92 9830 2.89 3.89 3.49 4.45 3.85 4.80 4.10 5.05 4.30 5.24 4.46 5.40 4.60 5.54 4.72 5.65 4.82 5.76 9940 2.86 3.82 3.44 4.37 3.79 4.70 4.04 4.93 4.23 5.11 4.39 5.26 4.52 5.39 4.63 5.50 4.73 5.60 10060 2.83 3.76 3.40 4.28 3.74 4.59 3.98 4.82 4.16 4.99 4.31 5.13 4.44 5.25 4.55 5.36 4.65 5.45 101120 2.80 3.70 3.36 4.20 3.68 4.50 3.92 4.71 4.10 4.87 4.24 5.01 4.36 5.12 4.47 5.21 4.56 5.30 102infinity 2.77 3.64 3.31 4.12 3.63 4.40 3.86 4.60 4.03 4.76 4.17 4.88 4.29 4.99 4.39 5.08 4.47 5.16 103''' 104 105res = [line.split() for line in qcrit.replace('infinity','9999').split('\n')] 106c=np.array(res[2:-1]).astype(float) 107#c[c==9999] = np.inf 108ccols = np.arange(2,11) 109crows = c[:,0] 110cv005 = c[:, 1::2] 111cv001 = c[:, 2::2] 112 113 114def get_tukeyQcrit(k, df, alpha=0.05): 115 ''' 116 return critical values for Tukey's HSD (Q) 117 118 Parameters 119 ---------- 120 k : int in {2, ..., 10} 121 number of tests 122 df : int 123 degrees of freedom of error term 124 alpha : {0.05, 0.01} 125 type 1 error, 1-confidence level 126 127 128 129 not enough error checking for limitations 130 ''' 131 if alpha == 0.05: 132 intp = interpolate.interp1d(crows, cv005[:,k-2]) 133 elif alpha == 0.01: 134 intp = interpolate.interp1d(crows, cv001[:,k-2]) 135 else: 136 raise ValueError('only implemented for alpha equal to 0.01 and 0.05') 137 return intp(df) 138 139def get_tukeyQcrit2(k, df, alpha=0.05): 140 ''' 141 return critical values for Tukey's HSD (Q) 142 143 Parameters 144 ---------- 145 k : int in {2, ..., 10} 146 number of tests 147 df : int 148 degrees of freedom of error term 149 alpha : {0.05, 0.01} 150 type 1 error, 1-confidence level 151 152 153 154 not enough error checking for limitations 155 ''' 156 from statsmodels.stats.libqsturng import qsturng 157 return qsturng(1-alpha, k, df) 158 159 160def get_tukey_pvalue(k, df, q): 161 ''' 162 return adjusted p-values for Tukey's HSD 163 164 Parameters 165 ---------- 166 k : int in {2, ..., 10} 167 number of tests 168 df : int 169 degrees of freedom of error term 170 q : scalar, array_like; q >= 0 171 quantile value of Studentized Range 172 173 ''' 174 175 from statsmodels.stats.libqsturng import psturng 176 return psturng(q, k, df) 177 178 179def Tukeythreegene(first, second, third): 180 # Performing the Tukey HSD post-hoc test for three genes 181 # qwb = xlrd.open_workbook('F:/Lab/bioinformatics/qcrittable.xls') 182 # #opening the workbook containing the q crit table 183 # qwb.sheet_names() 184 # qcrittable = qwb.sheet_by_name(u'Sheet1') 185 186 # means of the three arrays 187 firstmean = np.mean(first) 188 secondmean = np.mean(second) 189 thirdmean = np.mean(third) 190 191 # standard deviations of the threearrays 192 firststd = np.std(first) 193 secondstd = np.std(second) 194 thirdstd = np.std(third) 195 196 # standard deviation squared of the three arrays 197 firsts2 = math.pow(firststd, 2) 198 seconds2 = math.pow(secondstd, 2) 199 thirds2 = math.pow(thirdstd, 2) 200 201 # numerator for mean square error 202 mserrornum = firsts2 * 2 + seconds2 * 2 + thirds2 * 2 203 # denominator for mean square error 204 mserrorden = (len(first) + len(second) + len(third)) - 3 205 mserror = mserrornum / mserrorden # mean square error 206 207 standarderror = math.sqrt(mserror / len(first)) 208 # standard error, which is square root of mserror and 209 # the number of samples in a group 210 211 # various degrees of freedom 212 dftotal = len(first) + len(second) + len(third) - 1 213 dfgroups = 2 214 dferror = dftotal - dfgroups # noqa: F841 215 216 qcrit = 0.5 # fix arbitrary#qcrittable.cell(dftotal, 3).value 217 qcrit = get_tukeyQcrit(3, dftotal, alpha=0.05) 218 # getting the q critical value, for degrees of freedom total and 3 groups 219 220 qtest3to1 = (math.fabs(thirdmean - firstmean)) / standarderror 221 # calculating q test statistic values 222 qtest3to2 = (math.fabs(thirdmean - secondmean)) / standarderror 223 qtest2to1 = (math.fabs(secondmean - firstmean)) / standarderror 224 225 conclusion = [] 226 227 # print(qcrit 228 print(qtest3to1) 229 print(qtest3to2) 230 print(qtest2to1) 231 232 # testing all q test statistic values to q critical values 233 if qtest3to1 > qcrit: 234 conclusion.append('3to1null') 235 else: 236 conclusion.append('3to1alt') 237 if qtest3to2 > qcrit: 238 conclusion.append('3to2null') 239 else: 240 conclusion.append('3to2alt') 241 if qtest2to1 > qcrit: 242 conclusion.append('2to1null') 243 else: 244 conclusion.append('2to1alt') 245 246 return conclusion 247 248 249#rewrite by Vincent 250def Tukeythreegene2(genes): #Performing the Tukey HSD post-hoc test for three genes 251 """gend is a list, ie [first, second, third]""" 252# qwb = xlrd.open_workbook('F:/Lab/bioinformatics/qcrittable.xls') 253 #opening the workbook containing the q crit table 254# qwb.sheet_names() 255# qcrittable = qwb.sheet_by_name(u'Sheet1') 256 257 means = [] 258 stds = [] 259 for gene in genes: 260 means.append(np.mean(gene)) 261 std.append(np.std(gene)) # noqa:F821 See GH#5756 262 263 #firstmean = np.mean(first) #means of the three arrays 264 #secondmean = np.mean(second) 265 #thirdmean = np.mean(third) 266 267 #firststd = np.std(first) #standard deviations of the three arrays 268 #secondstd = np.std(second) 269 #thirdstd = np.std(third) 270 271 stds2 = [] 272 for std in stds: 273 stds2.append(math.pow(std,2)) 274 275 276 #firsts2 = math.pow(firststd,2) #standard deviation squared of the three arrays 277 #seconds2 = math.pow(secondstd,2) 278 #thirds2 = math.pow(thirdstd,2) 279 280 #mserrornum = firsts2*2+seconds2*2+thirds2*2 #numerator for mean square error 281 mserrornum = sum(stds2)*2 282 mserrorden = (len(genes[0])+len(genes[1])+len(genes[2]))-3 #denominator for mean square error 283 mserror = mserrornum/mserrorden #mean square error 284 285 286def catstack(args): 287 x = np.hstack(args) 288 labels = np.hstack([k*np.ones(len(arr)) for k,arr in enumerate(args)]) 289 return x, labels 290 291 292 293 294def maxzero(x): 295 '''find all up zero crossings and return the index of the highest 296 297 Not used anymore 298 299 300 >>> np.random.seed(12345) 301 >>> x = np.random.randn(8) 302 >>> x 303 array([-0.20470766, 0.47894334, -0.51943872, -0.5557303 , 1.96578057, 304 1.39340583, 0.09290788, 0.28174615]) 305 >>> maxzero(x) 306 (4, array([1, 4])) 307 308 309 no up-zero-crossing at end 310 311 >>> np.random.seed(0) 312 >>> x = np.random.randn(8) 313 >>> x 314 array([ 1.76405235, 0.40015721, 0.97873798, 2.2408932 , 1.86755799, 315 -0.97727788, 0.95008842, -0.15135721]) 316 >>> maxzero(x) 317 (None, array([6])) 318 ''' 319 x = np.asarray(x) 320 cond1 = x[:-1] < 0 321 cond2 = x[1:] > 0 322 #allzeros = np.nonzero(np.sign(x[:-1])*np.sign(x[1:]) <= 0)[0] + 1 323 allzeros = np.nonzero((cond1 & cond2) | (x[1:]==0))[0] + 1 324 if x[-1] >=0: 325 maxz = max(allzeros) 326 else: 327 maxz = None 328 return maxz, allzeros 329 330def maxzerodown(x): 331 '''find all up zero crossings and return the index of the highest 332 333 Not used anymore 334 335 >>> np.random.seed(12345) 336 >>> x = np.random.randn(8) 337 >>> x 338 array([-0.20470766, 0.47894334, -0.51943872, -0.5557303 , 1.96578057, 339 1.39340583, 0.09290788, 0.28174615]) 340 >>> maxzero(x) 341 (4, array([1, 4])) 342 343 344 no up-zero-crossing at end 345 346 >>> np.random.seed(0) 347 >>> x = np.random.randn(8) 348 >>> x 349 array([ 1.76405235, 0.40015721, 0.97873798, 2.2408932 , 1.86755799, 350 -0.97727788, 0.95008842, -0.15135721]) 351 >>> maxzero(x) 352 (None, array([6])) 353''' 354 x = np.asarray(x) 355 cond1 = x[:-1] > 0 356 cond2 = x[1:] < 0 357 #allzeros = np.nonzero(np.sign(x[:-1])*np.sign(x[1:]) <= 0)[0] + 1 358 allzeros = np.nonzero((cond1 & cond2) | (x[1:]==0))[0] + 1 359 if x[-1] <=0: 360 maxz = max(allzeros) 361 else: 362 maxz = None 363 return maxz, allzeros 364 365 366 367def rejectionline(n, alpha=0.5): 368 '''reference line for rejection in multiple tests 369 370 Not used anymore 371 372 from: section 3.2, page 60 373 ''' 374 t = np.arange(n)/float(n) 375 frej = t/( t * (1-alpha) + alpha) 376 return frej 377 378 379 380 381 382 383#I do not remember what I changed or why 2 versions, 384#this follows german diss ??? with rline 385#this might be useful if the null hypothesis is not "all effects are zero" 386#rename to _bak and working again on fdrcorrection0 387def fdrcorrection_bak(pvals, alpha=0.05, method='indep'): 388 '''Reject False discovery rate correction for pvalues 389 390 Old version, to be deleted 391 392 393 missing: methods that estimate fraction of true hypotheses 394 395 ''' 396 pvals = np.asarray(pvals) 397 398 399 pvals_sortind = np.argsort(pvals) 400 pvals_sorted = pvals[pvals_sortind] 401 pecdf = ecdf(pvals_sorted) 402 if method in ['i', 'indep', 'p', 'poscorr']: 403 rline = pvals_sorted / alpha 404 elif method in ['n', 'negcorr']: 405 cm = np.sum(1./np.arange(1, len(pvals))) 406 rline = pvals_sorted / alpha * cm 407 elif method in ['g', 'onegcorr']: #what's this ? german diss 408 rline = pvals_sorted / (pvals_sorted*(1-alpha) + alpha) 409 elif method in ['oth', 'o2negcorr']: # other invalid, cut-paste 410 cm = np.sum(np.arange(len(pvals))) 411 rline = pvals_sorted / alpha /cm 412 else: 413 raise ValueError('method not available') 414 415 reject = pecdf >= rline 416 if reject.any(): 417 rejectmax = max(np.nonzero(reject)[0]) 418 else: 419 rejectmax = 0 420 reject[:rejectmax] = True 421 return reject[pvals_sortind.argsort()] 422 423def mcfdr(nrepl=100, nobs=50, ntests=10, ntrue=6, mu=0.5, alpha=0.05, rho=0.): 424 '''MonteCarlo to test fdrcorrection 425 ''' 426 nfalse = ntests - ntrue 427 locs = np.array([0.]*ntrue + [mu]*(ntests - ntrue)) 428 results = [] 429 for i in range(nrepl): 430 #rvs = locs + stats.norm.rvs(size=(nobs, ntests)) 431 rvs = locs + randmvn(rho, size=(nobs, ntests)) 432 tt, tpval = stats.ttest_1samp(rvs, 0) 433 res = fdrcorrection_bak(np.abs(tpval), alpha=alpha, method='i') 434 res0 = fdrcorrection0(np.abs(tpval), alpha=alpha) 435 #res and res0 give the same results 436 results.append([np.sum(res[:ntrue]), np.sum(res[ntrue:])] + 437 [np.sum(res0[:ntrue]), np.sum(res0[ntrue:])] + 438 res.tolist() + 439 np.sort(tpval).tolist() + 440 [np.sum(tpval[:ntrue]<alpha), 441 np.sum(tpval[ntrue:]<alpha)] + 442 [np.sum(tpval[:ntrue]<alpha/ntests), 443 np.sum(tpval[ntrue:]<alpha/ntests)]) 444 return np.array(results) 445 446def randmvn(rho, size=(1, 2), standardize=False): 447 '''create random draws from equi-correlated multivariate normal distribution 448 449 Parameters 450 ---------- 451 rho : float 452 correlation coefficient 453 size : tuple of int 454 size is interpreted (nobs, nvars) where each row 455 456 Returns 457 ------- 458 rvs : ndarray 459 nobs by nvars where each row is a independent random draw of nvars- 460 dimensional correlated rvs 461 462 ''' 463 nobs, nvars = size 464 if 0 < rho and rho < 1: 465 rvs = np.random.randn(nobs, nvars+1) 466 rvs2 = rvs[:,:-1] * np.sqrt((1-rho)) + rvs[:,-1:] * np.sqrt(rho) 467 elif rho ==0: 468 rvs2 = np.random.randn(nobs, nvars) 469 elif rho < 0: 470 if rho < -1./(nvars-1): 471 raise ValueError('rho has to be larger than -1./(nvars-1)') 472 elif rho == -1./(nvars-1): 473 rho = -1./(nvars-1+1e-10) #barely positive definite 474 #use Cholesky 475 A = rho*np.ones((nvars,nvars))+(1-rho)*np.eye(nvars) 476 rvs2 = np.dot(np.random.randn(nobs, nvars), np.linalg.cholesky(A).T) 477 if standardize: 478 rvs2 = stats.zscore(rvs2) 479 return rvs2 480 481#============================ 482# 483# Part 2: Multiple comparisons and independent samples tests 484# 485#============================ 486 487def tiecorrect(xranks): 488 ''' 489 490 should be equivalent of scipy.stats.tiecorrect 491 492 ''' 493 #casting to int rounds down, but not relevant for this case 494 rankbincount = np.bincount(np.asarray(xranks,dtype=int)) 495 nties = rankbincount[rankbincount > 1] 496 ntot = float(len(xranks)) 497 tiecorrection = 1 - (nties**3 - nties).sum()/(ntot**3 - ntot) 498 return tiecorrection 499 500 501class GroupsStats(object): 502 ''' 503 statistics by groups (another version) 504 505 groupstats as a class with lazy evaluation (not yet - decorators are still 506 missing) 507 508 written this time as equivalent of scipy.stats.rankdata 509 gs = GroupsStats(X, useranks=True) 510 assert_almost_equal(gs.groupmeanfilter, stats.rankdata(X[:,0]), 15) 511 512 TODO: incomplete doc strings 513 514 ''' 515 516 def __init__(self, x, useranks=False, uni=None, intlab=None): 517 '''descriptive statistics by groups 518 519 Parameters 520 ---------- 521 x : ndarray, 2d 522 first column data, second column group labels 523 useranks : bool 524 if true, then use ranks as data corresponding to the 525 scipy.stats.rankdata definition (start at 1, ties get mean) 526 uni, intlab : arrays (optional) 527 to avoid call to unique, these can be given as inputs 528 529 530 ''' 531 self.x = np.asarray(x) 532 if intlab is None: 533 uni, intlab = np.unique(x[:,1], return_inverse=True) 534 elif uni is None: 535 uni = np.unique(x[:,1]) 536 537 self.useranks = useranks 538 539 540 self.uni = uni 541 self.intlab = intlab 542 self.groupnobs = groupnobs = np.bincount(intlab) 543 544 #temporary until separated and made all lazy 545 self.runbasic(useranks=useranks) 546 547 548 549 def runbasic_old(self, useranks=False): 550 """runbasic_old""" 551 #check: refactoring screwed up case useranks=True 552 553 #groupxsum = np.bincount(intlab, weights=X[:,0]) 554 #groupxmean = groupxsum * 1.0 / groupnobs 555 x = self.x 556 if useranks: 557 self.xx = x[:,1].argsort().argsort() + 1 #rankraw 558 else: 559 self.xx = x[:,0] 560 self.groupsum = groupranksum = np.bincount(self.intlab, weights=self.xx) 561 #print('groupranksum', groupranksum, groupranksum.shape, self.groupnobs.shape 562 # start at 1 for stats.rankdata : 563 self.groupmean = grouprankmean = groupranksum * 1.0 / self.groupnobs # + 1 564 self.groupmeanfilter = grouprankmean[self.intlab] 565 #return grouprankmean[intlab] 566 567 def runbasic(self, useranks=False): 568 """runbasic""" 569 #check: refactoring screwed up case useranks=True 570 571 #groupxsum = np.bincount(intlab, weights=X[:,0]) 572 #groupxmean = groupxsum * 1.0 / groupnobs 573 x = self.x 574 if useranks: 575 xuni, xintlab = np.unique(x[:,0], return_inverse=True) 576 ranksraw = x[:,0].argsort().argsort() + 1 #rankraw 577 self.xx = GroupsStats(np.column_stack([ranksraw, xintlab]), 578 useranks=False).groupmeanfilter 579 else: 580 self.xx = x[:,0] 581 self.groupsum = groupranksum = np.bincount(self.intlab, weights=self.xx) 582 #print('groupranksum', groupranksum, groupranksum.shape, self.groupnobs.shape 583 # start at 1 for stats.rankdata : 584 self.groupmean = grouprankmean = groupranksum * 1.0 / self.groupnobs # + 1 585 self.groupmeanfilter = grouprankmean[self.intlab] 586 #return grouprankmean[intlab] 587 588 def groupdemean(self): 589 """groupdemean""" 590 return self.xx - self.groupmeanfilter 591 592 def groupsswithin(self): 593 """groupsswithin""" 594 xtmp = self.groupdemean() 595 return np.bincount(self.intlab, weights=xtmp**2) 596 597 def groupvarwithin(self): 598 """groupvarwithin""" 599 return self.groupsswithin()/(self.groupnobs-1) #.sum() 600 601class TukeyHSDResults(object): 602 """Results from Tukey HSD test, with additional plot methods 603 604 Can also compute and plot additional post-hoc evaluations using this 605 results class. 606 607 Attributes 608 ---------- 609 reject : array of boolean, True if we reject Null for group pair 610 meandiffs : pairwise mean differences 611 confint : confidence interval for pairwise mean differences 612 std_pairs : standard deviation of pairwise mean differences 613 q_crit : critical value of studentized range statistic at given alpha 614 halfwidths : half widths of simultaneous confidence interval 615 pvalues : adjusted p-values from the HSD test 616 617 Notes 618 ----- 619 halfwidths is only available after call to `plot_simultaneous`. 620 621 Other attributes contain information about the data from the 622 MultiComparison instance: data, df_total, groups, groupsunique, variance. 623 """ 624 def __init__(self, mc_object, results_table, q_crit, reject=None, 625 meandiffs=None, std_pairs=None, confint=None, df_total=None, 626 reject2=None, variance=None, pvalues=None): 627 628 self._multicomp = mc_object 629 self._results_table = results_table 630 self.q_crit = q_crit 631 self.reject = reject 632 self.meandiffs = meandiffs 633 self.std_pairs = std_pairs 634 self.confint = confint 635 self.df_total = df_total 636 self.reject2 = reject2 637 self.variance = variance 638 self.pvalues = pvalues 639 # Taken out of _multicomp for ease of access for unknowledgeable users 640 self.data = self._multicomp.data 641 self.groups = self._multicomp.groups 642 self.groupsunique = self._multicomp.groupsunique 643 644 def __str__(self): 645 return str(self._results_table) 646 647 def summary(self): 648 '''Summary table that can be printed 649 ''' 650 return self._results_table 651 652 653 def _simultaneous_ci(self): 654 """Compute simultaneous confidence intervals for comparison of means. 655 """ 656 self.halfwidths = simultaneous_ci(self.q_crit, self.variance, 657 self._multicomp.groupstats.groupnobs, 658 self._multicomp.pairindices) 659 660 def plot_simultaneous(self, comparison_name=None, ax=None, figsize=(10,6), 661 xlabel=None, ylabel=None): 662 """Plot a universal confidence interval of each group mean 663 664 Visualize significant differences in a plot with one confidence 665 interval per group instead of all pairwise confidence intervals. 666 667 Parameters 668 ---------- 669 comparison_name : str, optional 670 if provided, plot_intervals will color code all groups that are 671 significantly different from the comparison_name red, and will 672 color code insignificant groups gray. Otherwise, all intervals will 673 just be plotted in black. 674 ax : matplotlib axis, optional 675 An axis handle on which to attach the plot. 676 figsize : tuple, optional 677 tuple for the size of the figure generated 678 xlabel : str, optional 679 Name to be displayed on x axis 680 ylabel : str, optional 681 Name to be displayed on y axis 682 683 Returns 684 ------- 685 Figure 686 handle to figure object containing interval plots 687 688 Notes 689 ----- 690 Multiple comparison tests are nice, but lack a good way to be 691 visualized. If you have, say, 6 groups, showing a graph of the means 692 between each group will require 15 confidence intervals. 693 Instead, we can visualize inter-group differences with a single 694 interval for each group mean. Hochberg et al. [1] first proposed this 695 idea and used Tukey's Q critical value to compute the interval widths. 696 Unlike plotting the differences in the means and their respective 697 confidence intervals, any two pairs can be compared for significance 698 by looking for overlap. 699 700 References 701 ---------- 702 .. [*] Hochberg, Y., and A. C. Tamhane. Multiple Comparison Procedures. 703 Hoboken, NJ: John Wiley & Sons, 1987. 704 705 Examples 706 -------- 707 >>> from statsmodels.examples.try_tukey_hsd import cylinders, cyl_labels 708 >>> from statsmodels.stats.multicomp import MultiComparison 709 >>> cardata = MultiComparison(cylinders, cyl_labels) 710 >>> results = cardata.tukeyhsd() 711 >>> results.plot_simultaneous() 712 <matplotlib.figure.Figure at 0x...> 713 714 This example shows an example plot comparing significant differences 715 in group means. Significant differences at the alpha=0.05 level can be 716 identified by intervals that do not overlap (i.e. USA vs Japan, 717 USA vs Germany). 718 719 >>> results.plot_simultaneous(comparison_name="USA") 720 <matplotlib.figure.Figure at 0x...> 721 722 Optionally provide one of the group names to color code the plot to 723 highlight group means different from comparison_name. 724 """ 725 fig, ax1 = utils.create_mpl_ax(ax) 726 if figsize is not None: 727 fig.set_size_inches(figsize) 728 if getattr(self, 'halfwidths', None) is None: 729 self._simultaneous_ci() 730 means = self._multicomp.groupstats.groupmean 731 732 733 sigidx = [] 734 nsigidx = [] 735 minrange = [means[i] - self.halfwidths[i] for i in range(len(means))] 736 maxrange = [means[i] + self.halfwidths[i] for i in range(len(means))] 737 738 if comparison_name is None: 739 ax1.errorbar(means, lrange(len(means)), xerr=self.halfwidths, 740 marker='o', linestyle='None', color='k', ecolor='k') 741 else: 742 if comparison_name not in self.groupsunique: 743 raise ValueError('comparison_name not found in group names.') 744 midx = np.where(self.groupsunique==comparison_name)[0][0] 745 for i in range(len(means)): 746 if self.groupsunique[i] == comparison_name: 747 continue 748 if (min(maxrange[i], maxrange[midx]) - 749 max(minrange[i], minrange[midx]) < 0): 750 sigidx.append(i) 751 else: 752 nsigidx.append(i) 753 #Plot the main comparison 754 ax1.errorbar(means[midx], midx, xerr=self.halfwidths[midx], 755 marker='o', linestyle='None', color='b', ecolor='b') 756 ax1.plot([minrange[midx]]*2, [-1, self._multicomp.ngroups], 757 linestyle='--', color='0.7') 758 ax1.plot([maxrange[midx]]*2, [-1, self._multicomp.ngroups], 759 linestyle='--', color='0.7') 760 #Plot those that are significantly different 761 if len(sigidx) > 0: 762 ax1.errorbar(means[sigidx], sigidx, 763 xerr=self.halfwidths[sigidx], marker='o', 764 linestyle='None', color='r', ecolor='r') 765 #Plot those that are not significantly different 766 if len(nsigidx) > 0: 767 ax1.errorbar(means[nsigidx], nsigidx, 768 xerr=self.halfwidths[nsigidx], marker='o', 769 linestyle='None', color='0.5', ecolor='0.5') 770 771 ax1.set_title('Multiple Comparisons Between All Pairs (Tukey)') 772 r = np.max(maxrange) - np.min(minrange) 773 ax1.set_ylim([-1, self._multicomp.ngroups]) 774 ax1.set_xlim([np.min(minrange) - r / 10., np.max(maxrange) + r / 10.]) 775 ylbls = [""] + self.groupsunique.astype(str).tolist() + [""] 776 ax1.set_yticks(np.arange(-1, len(means) + 1)) 777 ax1.set_yticklabels(ylbls) 778 ax1.set_xlabel(xlabel if xlabel is not None else '') 779 ax1.set_ylabel(ylabel if ylabel is not None else '') 780 return fig 781 782 783class MultiComparison(object): 784 '''Tests for multiple comparisons 785 786 Parameters 787 ---------- 788 data : ndarray 789 independent data samples 790 groups : ndarray 791 group labels corresponding to each data point 792 group_order : list[str], optional 793 the desired order for the group mean results to be reported in. If 794 not specified, results are reported in increasing order. 795 If group_order does not contain all labels that are in groups, then 796 only those observations are kept that have a label in group_order. 797 798 ''' 799 800 def __init__(self, data, groups, group_order=None): 801 802 if len(data) != len(groups): 803 raise ValueError('data has %d elements and groups has %d' % (len(data), len(groups))) 804 self.data = np.asarray(data) 805 self.groups = groups = np.asarray(groups) 806 807 # Allow for user-provided sorting of groups 808 if group_order is None: 809 self.groupsunique, self.groupintlab = np.unique(groups, 810 return_inverse=True) 811 else: 812 #check if group_order has any names not in groups 813 for grp in group_order: 814 if grp not in groups: 815 raise ValueError( 816 "group_order value '%s' not found in groups" % grp) 817 self.groupsunique = np.array(group_order) 818 self.groupintlab = np.empty(len(data), int) 819 self.groupintlab.fill(-999) # instead of a nan 820 count = 0 821 for name in self.groupsunique: 822 idx = np.where(self.groups == name)[0] 823 count += len(idx) 824 self.groupintlab[idx] = np.where(self.groupsunique == name)[0] 825 if count != self.data.shape[0]: 826 #raise ValueError('group_order does not contain all groups') 827 # warn and keep only observations with label in group_order 828 import warnings 829 warnings.warn('group_order does not contain all groups:' + 830 ' dropping observations', ValueWarning) 831 832 mask_keep = self.groupintlab != -999 833 self.groupintlab = self.groupintlab[mask_keep] 834 self.data = self.data[mask_keep] 835 self.groups = self.groups[mask_keep] 836 837 if len(self.groupsunique) < 2: 838 raise ValueError('2 or more groups required for multiple comparisons') 839 840 self.datali = [self.data[self.groups == k] for k in self.groupsunique] 841 self.pairindices = np.triu_indices(len(self.groupsunique), 1) #tuple 842 self.nobs = self.data.shape[0] 843 self.ngroups = len(self.groupsunique) 844 845 846 def getranks(self): 847 '''convert data to rankdata and attach 848 849 850 This creates rankdata as it is used for non-parametric tests, where 851 in the case of ties the average rank is assigned. 852 853 854 ''' 855 #bug: the next should use self.groupintlab instead of self.groups 856 #update: looks fixed 857 #self.ranks = GroupsStats(np.column_stack([self.data, self.groups]), 858 self.ranks = GroupsStats(np.column_stack([self.data, self.groupintlab]), 859 useranks=True) 860 self.rankdata = self.ranks.groupmeanfilter 861 862 def kruskal(self, pairs=None, multimethod='T'): 863 ''' 864 pairwise comparison for kruskal-wallis test 865 866 This is just a reimplementation of scipy.stats.kruskal and does 867 not yet use a multiple comparison correction. 868 869 ''' 870 self.getranks() 871 tot = self.nobs 872 meanranks = self.ranks.groupmean 873 groupnobs = self.ranks.groupnobs 874 875 876 # simultaneous/separate treatment of multiple tests 877 f=(tot * (tot + 1.) / 12.) / stats.tiecorrect(self.rankdata) #(xranks) 878 print('MultiComparison.kruskal') 879 for i,j in zip(*self.pairindices): 880 #pdiff = np.abs(mrs[i] - mrs[j]) 881 pdiff = np.abs(meanranks[i] - meanranks[j]) 882 se = np.sqrt(f * np.sum(1. / groupnobs[[i,j]] )) #np.array([8,8]))) #Fixme groupnobs[[i,j]] )) 883 Q = pdiff / se 884 885 # TODO : print(statments, fix 886 print(i,j, pdiff, se, pdiff / se, pdiff / se > 2.6310) 887 print(stats.norm.sf(Q) * 2) 888 return stats.norm.sf(Q) * 2 889 890 891 def allpairtest(self, testfunc, alpha=0.05, method='bonf', pvalidx=1): 892 '''run a pairwise test on all pairs with multiple test correction 893 894 The statistical test given in testfunc is calculated for all pairs 895 and the p-values are adjusted by methods in multipletests. The p-value 896 correction is generic and based only on the p-values, and does not 897 take any special structure of the hypotheses into account. 898 899 Parameters 900 ---------- 901 testfunc : function 902 A test function for two (independent) samples. It is assumed that 903 the return value on position pvalidx is the p-value. 904 alpha : float 905 familywise error rate 906 method : str 907 This specifies the method for the p-value correction. Any method 908 of multipletests is possible. 909 pvalidx : int (default: 1) 910 position of the p-value in the return of testfunc 911 912 Returns 913 ------- 914 sumtab : SimpleTable instance 915 summary table for printing 916 917 errors: TODO: check if this is still wrong, I think it's fixed. 918 results from multipletests are in different order 919 pval_corrected can be larger than 1 ??? 920 ''' 921 res = [] 922 for i,j in zip(*self.pairindices): 923 res.append(testfunc(self.datali[i], self.datali[j])) 924 res = np.array(res) 925 reject, pvals_corrected, alphacSidak, alphacBonf = \ 926 multipletests(res[:, pvalidx], alpha=alpha, method=method) 927 #print(np.column_stack([res[:,0],res[:,1], reject, pvals_corrected]) 928 929 i1, i2 = self.pairindices 930 if pvals_corrected is None: 931 resarr = np.array(lzip(self.groupsunique[i1], self.groupsunique[i2], 932 np.round(res[:,0],4), 933 np.round(res[:,1],4), 934 reject), 935 dtype=[('group1', object), 936 ('group2', object), 937 ('stat',float), 938 ('pval',float), 939 ('reject', np.bool8)]) 940 else: 941 resarr = np.array(lzip(self.groupsunique[i1], self.groupsunique[i2], 942 np.round(res[:,0],4), 943 np.round(res[:,1],4), 944 np.round(pvals_corrected,4), 945 reject), 946 dtype=[('group1', object), 947 ('group2', object), 948 ('stat',float), 949 ('pval',float), 950 ('pval_corr',float), 951 ('reject', np.bool8)]) 952 results_table = SimpleTable(resarr, headers=resarr.dtype.names) 953 results_table.title = ( 954 'Test Multiple Comparison %s \n%s%4.2f method=%s' 955 % (testfunc.__name__, 'FWER=', alpha, method) + 956 '\nalphacSidak=%4.2f, alphacBonf=%5.3f' 957 % (alphacSidak, alphacBonf)) 958 959 return results_table, (res, reject, pvals_corrected, 960 alphacSidak, alphacBonf), resarr 961 962 def tukeyhsd(self, alpha=0.05): 963 """ 964 Tukey's range test to compare means of all pairs of groups 965 966 Parameters 967 ---------- 968 alpha : float, optional 969 Value of FWER at which to calculate HSD. 970 971 Returns 972 ------- 973 results : TukeyHSDResults instance 974 A results class containing relevant data and some post-hoc 975 calculations 976 """ 977 self.groupstats = GroupsStats( 978 np.column_stack([self.data, self.groupintlab]), 979 useranks=False) 980 981 gmeans = self.groupstats.groupmean 982 gnobs = self.groupstats.groupnobs 983 # var_ = self.groupstats.groupvarwithin() 984 # #possibly an error in varcorrection in this case 985 var_ = np.var(self.groupstats.groupdemean(), ddof=len(gmeans)) 986 # res contains: 0:(idx1, idx2), 1:reject, 2:meandiffs, 3: std_pairs, 987 # 4:confint, 5:q_crit, 6:df_total, 7:reject2, 8: pvals 988 res = tukeyhsd(gmeans, gnobs, var_, df=None, alpha=alpha, q_crit=None) 989 990 resarr = np.array(lzip(self.groupsunique[res[0][0]], 991 self.groupsunique[res[0][1]], 992 np.round(res[2], 4), 993 np.round(res[8], 4), 994 np.round(res[4][:, 0], 4), 995 np.round(res[4][:, 1], 4), 996 res[1]), 997 dtype=[('group1', object), 998 ('group2', object), 999 ('meandiff', float), 1000 ('p-adj', float), 1001 ('lower', float), 1002 ('upper', float), 1003 ('reject', np.bool8)]) 1004 results_table = SimpleTable(resarr, headers=resarr.dtype.names) 1005 results_table.title = 'Multiple Comparison of Means - Tukey HSD, ' + \ 1006 'FWER=%4.2f' % alpha 1007 1008 return TukeyHSDResults(self, results_table, res[5], res[1], res[2], 1009 res[3], res[4], res[6], res[7], var_, res[8]) 1010 1011 1012def rankdata(x): 1013 '''rankdata, equivalent to scipy.stats.rankdata 1014 1015 just a different implementation, I have not yet compared speed 1016 1017 ''' 1018 uni, intlab = np.unique(x[:,0], return_inverse=True) 1019 groupnobs = np.bincount(intlab) 1020 groupxsum = np.bincount(intlab, weights=X[:,0]) 1021 groupxmean = groupxsum * 1.0 / groupnobs 1022 1023 rankraw = x[:,0].argsort().argsort() 1024 groupranksum = np.bincount(intlab, weights=rankraw) 1025 # start at 1 for stats.rankdata : 1026 grouprankmean = groupranksum * 1.0 / groupnobs + 1 1027 return grouprankmean[intlab] 1028 1029 1030#new 1031 1032def compare_ordered(vals, alpha): 1033 '''simple ordered sequential comparison of means 1034 1035 vals : array_like 1036 means or rankmeans for independent groups 1037 1038 incomplete, no return, not used yet 1039 ''' 1040 vals = np.asarray(vals) 1041 alphaf = alpha # Notation ? 1042 sortind = np.argsort(vals) 1043 pvals = vals[sortind] 1044 sortrevind = sortind.argsort() 1045 ntests = len(vals) 1046 #alphacSidak = 1 - np.power((1. - alphaf), 1./ntests) 1047 #alphacBonf = alphaf / float(ntests) 1048 v1, v2 = np.triu_indices(ntests, 1) 1049 #v1,v2 have wrong sequence 1050 for i in range(4): 1051 for j in range(4,i, -1): 1052 print(i,j) 1053 1054 1055 1056def varcorrection_unbalanced(nobs_all, srange=False): 1057 '''correction factor for variance with unequal sample sizes 1058 1059 this is just a harmonic mean 1060 1061 Parameters 1062 ---------- 1063 nobs_all : array_like 1064 The number of observations for each sample 1065 srange : bool 1066 if true, then the correction is divided by the number of samples 1067 for the variance of the studentized range statistic 1068 1069 Returns 1070 ------- 1071 correction : float 1072 Correction factor for variance. 1073 1074 1075 Notes 1076 ----- 1077 1078 variance correction factor is 1079 1080 1/k * sum_i 1/n_i 1081 1082 where k is the number of samples and summation is over i=0,...,k-1. 1083 If all n_i are the same, then the correction factor is 1. 1084 1085 This needs to be multiplied by the joint variance estimate, means square 1086 error, MSE. To obtain the correction factor for the standard deviation, 1087 square root needs to be taken. 1088 1089 ''' 1090 nobs_all = np.asarray(nobs_all) 1091 if not srange: 1092 return (1./nobs_all).sum() 1093 else: 1094 return (1./nobs_all).sum()/len(nobs_all) 1095 1096def varcorrection_pairs_unbalanced(nobs_all, srange=False): 1097 '''correction factor for variance with unequal sample sizes for all pairs 1098 1099 this is just a harmonic mean 1100 1101 Parameters 1102 ---------- 1103 nobs_all : array_like 1104 The number of observations for each sample 1105 srange : bool 1106 if true, then the correction is divided by 2 for the variance of 1107 the studentized range statistic 1108 1109 Returns 1110 ------- 1111 correction : ndarray 1112 Correction factor for variance. 1113 1114 1115 Notes 1116 ----- 1117 1118 variance correction factor is 1119 1120 1/k * sum_i 1/n_i 1121 1122 where k is the number of samples and summation is over i=0,...,k-1. 1123 If all n_i are the same, then the correction factor is 1. 1124 1125 This needs to be multiplies by the joint variance estimate, means square 1126 error, MSE. To obtain the correction factor for the standard deviation, 1127 square root needs to be taken. 1128 1129 For the studentized range statistic, the resulting factor has to be 1130 divided by 2. 1131 1132 ''' 1133 #TODO: test and replace with broadcasting 1134 n1, n2 = np.meshgrid(nobs_all, nobs_all) 1135 if not srange: 1136 return (1./n1 + 1./n2) 1137 else: 1138 return (1./n1 + 1./n2) / 2. 1139 1140def varcorrection_unequal(var_all, nobs_all, df_all): 1141 '''return joint variance from samples with unequal variances and unequal 1142 sample sizes 1143 1144 something is wrong 1145 1146 Parameters 1147 ---------- 1148 var_all : array_like 1149 The variance for each sample 1150 nobs_all : array_like 1151 The number of observations for each sample 1152 df_all : array_like 1153 degrees of freedom for each sample 1154 1155 Returns 1156 ------- 1157 varjoint : float 1158 joint variance. 1159 dfjoint : float 1160 joint Satterthwait's degrees of freedom 1161 1162 1163 Notes 1164 ----- 1165 (copy, paste not correct) 1166 variance is 1167 1168 1/k * sum_i 1/n_i 1169 1170 where k is the number of samples and summation is over i=0,...,k-1. 1171 If all n_i are the same, then the correction factor is 1/n. 1172 1173 This needs to be multiplies by the joint variance estimate, means square 1174 error, MSE. To obtain the correction factor for the standard deviation, 1175 square root needs to be taken. 1176 1177 This is for variance of mean difference not of studentized range. 1178 ''' 1179 1180 var_all = np.asarray(var_all) 1181 var_over_n = var_all *1./ nobs_all #avoid integer division 1182 varjoint = var_over_n.sum() 1183 1184 dfjoint = varjoint**2 / (var_over_n**2 * df_all).sum() 1185 1186 return varjoint, dfjoint 1187 1188def varcorrection_pairs_unequal(var_all, nobs_all, df_all): 1189 '''return joint variance from samples with unequal variances and unequal 1190 sample sizes for all pairs 1191 1192 something is wrong 1193 1194 Parameters 1195 ---------- 1196 var_all : array_like 1197 The variance for each sample 1198 nobs_all : array_like 1199 The number of observations for each sample 1200 df_all : array_like 1201 degrees of freedom for each sample 1202 1203 Returns 1204 ------- 1205 varjoint : ndarray 1206 joint variance. 1207 dfjoint : ndarray 1208 joint Satterthwait's degrees of freedom 1209 1210 1211 Notes 1212 ----- 1213 1214 (copy, paste not correct) 1215 variance is 1216 1217 1/k * sum_i 1/n_i 1218 1219 where k is the number of samples and summation is over i=0,...,k-1. 1220 If all n_i are the same, then the correction factor is 1. 1221 1222 This needs to be multiplies by the joint variance estimate, means square 1223 error, MSE. To obtain the correction factor for the standard deviation, 1224 square root needs to be taken. 1225 1226 TODO: something looks wrong with dfjoint, is formula from SPSS 1227 ''' 1228 #TODO: test and replace with broadcasting 1229 v1, v2 = np.meshgrid(var_all, var_all) 1230 n1, n2 = np.meshgrid(nobs_all, nobs_all) 1231 df1, df2 = np.meshgrid(df_all, df_all) 1232 1233 varjoint = v1/n1 + v2/n2 1234 1235 dfjoint = varjoint**2 / (df1 * (v1/n1)**2 + df2 * (v2/n2)**2) 1236 1237 return varjoint, dfjoint 1238 1239def tukeyhsd(mean_all, nobs_all, var_all, df=None, alpha=0.05, q_crit=None): 1240 '''simultaneous Tukey HSD 1241 1242 1243 check: instead of sorting, I use absolute value of pairwise differences 1244 in means. That's irrelevant for the test, but maybe reporting actual 1245 differences would be better. 1246 CHANGED: meandiffs are with sign, studentized range uses abs 1247 1248 q_crit added for testing 1249 1250 TODO: error in variance calculation when nobs_all is scalar, missing 1/n 1251 1252 ''' 1253 mean_all = np.asarray(mean_all) 1254 #check if or when other ones need to be arrays 1255 1256 n_means = len(mean_all) 1257 1258 if df is None: 1259 df = nobs_all - 1 1260 1261 if np.size(df) == 1: # assumes balanced samples with df = n - 1, n_i = n 1262 df_total = n_means * df 1263 df = np.ones(n_means) * df 1264 else: 1265 df_total = np.sum(df) 1266 1267 if (np.size(nobs_all) == 1) and (np.size(var_all) == 1): 1268 #balanced sample sizes and homogenous variance 1269 var_pairs = 1. * var_all / nobs_all * np.ones((n_means, n_means)) 1270 1271 elif np.size(var_all) == 1: 1272 #unequal sample sizes and homogenous variance 1273 var_pairs = var_all * varcorrection_pairs_unbalanced(nobs_all, 1274 srange=True) 1275 elif np.size(var_all) > 1: 1276 var_pairs, df_sum = varcorrection_pairs_unequal(nobs_all, var_all, df) 1277 var_pairs /= 2. 1278 #check division by two for studentized range 1279 1280 else: 1281 raise ValueError('not supposed to be here') 1282 1283 #meandiffs_ = mean_all[:,None] - mean_all 1284 meandiffs_ = mean_all - mean_all[:,None] #reverse sign, check with R example 1285 std_pairs_ = np.sqrt(var_pairs) 1286 1287 #select all pairs from upper triangle of matrix 1288 idx1, idx2 = np.triu_indices(n_means, 1) 1289 meandiffs = meandiffs_[idx1, idx2] 1290 std_pairs = std_pairs_[idx1, idx2] 1291 1292 st_range = np.abs(meandiffs) / std_pairs #studentized range statistic 1293 1294 df_total_ = max(df_total, 5) #TODO: smallest df in table 1295 if q_crit is None: 1296 q_crit = get_tukeyQcrit2(n_means, df_total, alpha=alpha) 1297 1298 pvalues = get_tukey_pvalue(n_means, df_total, st_range) 1299 # we need pvalues to be atleast_1d for iteration. see #6132 1300 pvalues = np.atleast_1d(pvalues) 1301 1302 reject = st_range > q_crit 1303 crit_int = std_pairs * q_crit 1304 reject2 = np.abs(meandiffs) > crit_int 1305 1306 confint = np.column_stack((meandiffs - crit_int, meandiffs + crit_int)) 1307 1308 return ((idx1, idx2), reject, meandiffs, std_pairs, confint, q_crit, 1309 df_total, reject2, pvalues) 1310 1311 1312def simultaneous_ci(q_crit, var, groupnobs, pairindices=None): 1313 """Compute simultaneous confidence intervals for comparison of means. 1314 1315 q_crit value is generated from tukey hsd test. Variance is considered 1316 across all groups. Returned halfwidths can be thought of as uncertainty 1317 intervals around each group mean. They allow for simultaneous 1318 comparison of pairwise significance among any pairs (by checking for 1319 overlap) 1320 1321 Parameters 1322 ---------- 1323 q_crit : float 1324 The Q critical value studentized range statistic from Tukey's HSD 1325 var : float 1326 The group variance 1327 groupnobs : array_like object 1328 Number of observations contained in each group. 1329 pairindices : tuple of lists, optional 1330 Indices corresponding to the upper triangle of matrix. Computed 1331 here if not supplied 1332 1333 Returns 1334 ------- 1335 halfwidths : ndarray 1336 Half the width of each confidence interval for each group given in 1337 groupnobs 1338 1339 See Also 1340 -------- 1341 MultiComparison : statistics class providing significance tests 1342 tukeyhsd : among other things, computes q_crit value 1343 1344 References 1345 ---------- 1346 .. [*] Hochberg, Y., and A. C. Tamhane. Multiple Comparison Procedures. 1347 Hoboken, NJ: John Wiley & Sons, 1987.) 1348 """ 1349 # Set initial variables 1350 ng = len(groupnobs) 1351 if pairindices is None: 1352 pairindices = np.triu_indices(ng, 1) 1353 1354 # Compute dij for all pairwise comparisons ala hochberg p. 95 1355 gvar = var / groupnobs 1356 1357 d12 = np.sqrt(gvar[pairindices[0]] + gvar[pairindices[1]]) 1358 1359 # Create the full d matrix given all known dij vals 1360 d = np.zeros((ng, ng)) 1361 d[pairindices] = d12 1362 d = d + d.conj().T 1363 1364 # Compute the two global sums from hochberg eq 3.32 1365 sum1 = np.sum(d12) 1366 sum2 = np.sum(d, axis=0) 1367 1368 if (ng > 2): 1369 w = ((ng-1.) * sum2 - sum1) / ((ng - 1.) * (ng - 2.)) 1370 else: 1371 w = sum1 * np.ones((2, 1)) / 2. 1372 1373 return (q_crit / np.sqrt(2))*w 1374 1375def distance_st_range(mean_all, nobs_all, var_all, df=None, triu=False): 1376 '''pairwise distance matrix, outsourced from tukeyhsd 1377 1378 1379 1380 CHANGED: meandiffs are with sign, studentized range uses abs 1381 1382 q_crit added for testing 1383 1384 TODO: error in variance calculation when nobs_all is scalar, missing 1/n 1385 1386 ''' 1387 mean_all = np.asarray(mean_all) 1388 #check if or when other ones need to be arrays 1389 1390 n_means = len(mean_all) 1391 1392 if df is None: 1393 df = nobs_all - 1 1394 1395 if np.size(df) == 1: # assumes balanced samples with df = n - 1, n_i = n 1396 df_total = n_means * df 1397 else: 1398 df_total = np.sum(df) 1399 1400 if (np.size(nobs_all) == 1) and (np.size(var_all) == 1): 1401 #balanced sample sizes and homogenous variance 1402 var_pairs = 1. * var_all / nobs_all * np.ones((n_means, n_means)) 1403 1404 elif np.size(var_all) == 1: 1405 #unequal sample sizes and homogenous variance 1406 var_pairs = var_all * varcorrection_pairs_unbalanced(nobs_all, 1407 srange=True) 1408 elif np.size(var_all) > 1: 1409 var_pairs, df_sum = varcorrection_pairs_unequal(nobs_all, var_all, df) 1410 var_pairs /= 2. 1411 #check division by two for studentized range 1412 1413 else: 1414 raise ValueError('not supposed to be here') 1415 1416 #meandiffs_ = mean_all[:,None] - mean_all 1417 meandiffs = mean_all - mean_all[:,None] #reverse sign, check with R example 1418 std_pairs = np.sqrt(var_pairs) 1419 1420 idx1, idx2 = np.triu_indices(n_means, 1) 1421 if triu: 1422 #select all pairs from upper triangle of matrix 1423 meandiffs = meandiffs_[idx1, idx2] # noqa: F821 See GH#5756 1424 std_pairs = std_pairs_[idx1, idx2] # noqa: F821 See GH#5756 1425 1426 st_range = np.abs(meandiffs) / std_pairs #studentized range statistic 1427 1428 return st_range, meandiffs, std_pairs, (idx1,idx2) #return square arrays 1429 1430 1431def contrast_allpairs(nm): 1432 '''contrast or restriction matrix for all pairs of nm variables 1433 1434 Parameters 1435 ---------- 1436 nm : int 1437 1438 Returns 1439 ------- 1440 contr : ndarray, 2d, (nm*(nm-1)/2, nm) 1441 contrast matrix for all pairwise comparisons 1442 1443 ''' 1444 contr = [] 1445 for i in range(nm): 1446 for j in range(i+1, nm): 1447 contr_row = np.zeros(nm) 1448 contr_row[i] = 1 1449 contr_row[j] = -1 1450 contr.append(contr_row) 1451 return np.array(contr) 1452 1453def contrast_all_one(nm): 1454 '''contrast or restriction matrix for all against first comparison 1455 1456 Parameters 1457 ---------- 1458 nm : int 1459 1460 Returns 1461 ------- 1462 contr : ndarray, 2d, (nm-1, nm) 1463 contrast matrix for all against first comparisons 1464 1465 ''' 1466 contr = np.column_stack((np.ones(nm-1), -np.eye(nm-1))) 1467 return contr 1468 1469def contrast_diff_mean(nm): 1470 '''contrast or restriction matrix for all against mean comparison 1471 1472 Parameters 1473 ---------- 1474 nm : int 1475 1476 Returns 1477 ------- 1478 contr : ndarray, 2d, (nm-1, nm) 1479 contrast matrix for all against mean comparisons 1480 1481 ''' 1482 return np.eye(nm) - np.ones((nm,nm))/nm 1483 1484def tukey_pvalues(std_range, nm, df): 1485 #corrected but very slow with warnings about integration 1486 #nm = len(std_range) 1487 contr = contrast_allpairs(nm) 1488 corr = np.dot(contr, contr.T)/2. 1489 tstat = std_range / np.sqrt(2) * np.ones(corr.shape[0]) #need len of all pairs 1490 return multicontrast_pvalues(tstat, corr, df=df) 1491 1492 1493def multicontrast_pvalues(tstat, tcorr, df=None, dist='t', alternative='two-sided'): 1494 '''pvalues for simultaneous tests 1495 1496 ''' 1497 from statsmodels.sandbox.distributions.multivariate import mvstdtprob 1498 if (df is None) and (dist == 't'): 1499 raise ValueError('df has to be specified for the t-distribution') 1500 tstat = np.asarray(tstat) 1501 ntests = len(tstat) 1502 cc = np.abs(tstat) 1503 pval_global = 1 - mvstdtprob(-cc,cc, tcorr, df) 1504 pvals = [] 1505 for ti in cc: 1506 limits = ti*np.ones(ntests) 1507 pvals.append(1 - mvstdtprob(-cc,cc, tcorr, df)) 1508 1509 return pval_global, np.asarray(pvals) 1510 1511 1512 1513 1514 1515class StepDown(object): 1516 '''a class for step down methods 1517 1518 This is currently for simple tree subset descend, similar to homogeneous_subsets, 1519 but checks all leave-one-out subsets instead of assuming an ordered set. 1520 Comment in SAS manual: 1521 SAS only uses interval subsets of the sorted list, which is sufficient for range 1522 tests (maybe also equal variance and balanced sample sizes are required). 1523 For F-test based critical distances, the restriction to intervals is not sufficient. 1524 1525 This version uses a single critical value of the studentized range distribution 1526 for all comparisons, and is therefore a step-down version of Tukey HSD. 1527 The class is written so it can be subclassed, where the get_distance_matrix and 1528 get_crit are overwritten to obtain other step-down procedures such as REGW. 1529 1530 iter_subsets can be overwritten, to get a recursion as in the many to one comparison 1531 with a control such as in Dunnet's test. 1532 1533 1534 A one-sided right tail test is not covered because the direction of the inequality 1535 is hard coded in check_set. Also Peritz's check of partitions is not possible, but 1536 I have not seen it mentioned in any more recent references. 1537 I have only partially read the step-down procedure for closed tests by Westfall. 1538 1539 One change to make it more flexible, is to separate out the decision on a subset, 1540 also because the F-based tests, FREGW in SPSS, take information from all elements of 1541 a set and not just pairwise comparisons. I have not looked at the details of 1542 the F-based tests such as Sheffe yet. It looks like running an F-test on equality 1543 of means in each subset. This would also outsource how pairwise conditions are 1544 combined, any larger or max. This would also imply that the distance matrix cannot 1545 be calculated in advance for tests like the F-based ones. 1546 1547 1548 ''' 1549 1550 def __init__(self, vals, nobs_all, var_all, df=None): 1551 self.vals = vals 1552 self.n_vals = len(vals) 1553 self.nobs_all = nobs_all 1554 self.var_all = var_all 1555 self.df = df 1556 # the following has been moved to run 1557 #self.cache_result = {} 1558 #self.crit = self.getcrit(0.5) #decide where to set alpha, moved to run 1559 #self.accepted = [] #store accepted sets, not unique 1560 1561 def get_crit(self, alpha): 1562 """ 1563 get_tukeyQcrit 1564 1565 currently tukey Q, add others 1566 """ 1567 q_crit = get_tukeyQcrit(self.n_vals, self.df, alpha=alpha) 1568 return q_crit * np.ones(self.n_vals) 1569 1570 1571 1572 def get_distance_matrix(self): 1573 '''studentized range statistic''' 1574 #make into property, decorate 1575 dres = distance_st_range(self.vals, self.nobs_all, self.var_all, df=self.df) 1576 self.distance_matrix = dres[0] 1577 1578 def iter_subsets(self, indices): 1579 """Iterate substeps""" 1580 for ii in range(len(indices)): 1581 idxsub = copy.copy(indices) 1582 idxsub.pop(ii) 1583 yield idxsub 1584 1585 1586 def check_set(self, indices): 1587 '''check whether pairwise distances of indices satisfy condition 1588 1589 ''' 1590 indtup = tuple(indices) 1591 if indtup in self.cache_result: 1592 return self.cache_result[indtup] 1593 else: 1594 set_distance_matrix = self.distance_matrix[np.asarray(indices)[:,None], indices] 1595 n_elements = len(indices) 1596 if np.any(set_distance_matrix > self.crit[n_elements-1]): 1597 res = True 1598 else: 1599 res = False 1600 self.cache_result[indtup] = res 1601 return res 1602 1603 def stepdown(self, indices): 1604 """stepdown""" 1605 print(indices) 1606 if self.check_set(indices): # larger than critical distance 1607 if (len(indices) > 2): # step down into subsets if more than 2 elements 1608 for subs in self.iter_subsets(indices): 1609 self.stepdown(subs) 1610 else: 1611 self.rejected.append(tuple(indices)) 1612 else: 1613 self.accepted.append(tuple(indices)) 1614 return indices 1615 1616 def run(self, alpha): 1617 '''main function to run the test, 1618 1619 could be done in __call__ instead 1620 this could have all the initialization code 1621 1622 ''' 1623 self.cache_result = {} 1624 self.crit = self.get_crit(alpha) #decide where to set alpha, moved to run 1625 self.accepted = [] #store accepted sets, not unique 1626 self.rejected = [] 1627 self.get_distance_matrix() 1628 self.stepdown(lrange(self.n_vals)) 1629 1630 return list(set(self.accepted)), list(set(sd.rejected)) 1631 1632 1633 1634 1635 1636 1637def homogeneous_subsets(vals, dcrit): 1638 '''recursively check all pairs of vals for minimum distance 1639 1640 step down method as in Newman-Keuls and Ryan procedures. This is not a 1641 closed procedure since not all partitions are checked. 1642 1643 Parameters 1644 ---------- 1645 vals : array_like 1646 values that are pairwise compared 1647 dcrit : array_like or float 1648 critical distance for rejecting, either float, or 2-dimensional array 1649 with distances on the upper triangle. 1650 1651 Returns 1652 ------- 1653 rejs : list of pairs 1654 list of pair-indices with (strictly) larger than critical difference 1655 nrejs : list of pairs 1656 list of pair-indices with smaller than critical difference 1657 lli : list of tuples 1658 list of subsets with smaller than critical difference 1659 res : tree 1660 result of all comparisons (for checking) 1661 1662 1663 this follows description in SPSS notes on Post-Hoc Tests 1664 1665 Because of the recursive structure, some comparisons are made several 1666 times, but only unique pairs or sets are returned. 1667 1668 Examples 1669 -------- 1670 >>> m = [0, 2, 2.5, 3, 6, 8, 9, 9.5,10 ] 1671 >>> rej, nrej, ssli, res = homogeneous_subsets(m, 2) 1672 >>> set_partition(ssli) 1673 ([(5, 6, 7, 8), (1, 2, 3), (4,)], [0]) 1674 >>> [np.array(m)[list(pp)] for pp in set_partition(ssli)[0]] 1675 [array([ 8. , 9. , 9.5, 10. ]), array([ 2. , 2.5, 3. ]), array([ 6.])] 1676 1677 1678 ''' 1679 1680 nvals = len(vals) 1681 indices_ = lrange(nvals) 1682 rejected = [] 1683 subsetsli = [] 1684 if np.size(dcrit) == 1: 1685 dcrit = dcrit*np.ones((nvals, nvals)) #example numbers for experimenting 1686 1687 def subsets(vals, indices_): 1688 '''recursive function for constructing homogeneous subset 1689 1690 registers rejected and subsetli in outer scope 1691 ''' 1692 i, j = (indices_[0], indices_[-1]) 1693 if vals[-1] - vals[0] > dcrit[i,j]: 1694 rejected.append((indices_[0], indices_[-1])) 1695 return [subsets(vals[:-1], indices_[:-1]), 1696 subsets(vals[1:], indices_[1:]), 1697 (indices_[0], indices_[-1])] 1698 else: 1699 subsetsli.append(tuple(indices_)) 1700 return indices_ 1701 res = subsets(vals, indices_) 1702 1703 all_pairs = [(i,j) for i in range(nvals) for j in range(nvals-1,i,-1)] 1704 rejs = set(rejected) 1705 not_rejected = list(set(all_pairs) - rejs) 1706 1707 return list(rejs), not_rejected, list(set(subsetsli)), res 1708 1709def set_partition(ssli): 1710 '''extract a partition from a list of tuples 1711 1712 this should be correctly called select largest disjoint sets. 1713 Begun and Gabriel 1981 do not seem to be bothered by sets of accepted 1714 hypothesis with joint elements, 1715 e.g. maximal_accepted_sets = { {1,2,3}, {2,3,4} } 1716 1717 This creates a set partition from a list of sets given as tuples. 1718 It tries to find the partition with the largest sets. That is, sets are 1719 included after being sorted by length. 1720 1721 If the list does not include the singletons, then it will be only a 1722 partial partition. Missing items are singletons (I think). 1723 1724 Examples 1725 -------- 1726 >>> li 1727 [(5, 6, 7, 8), (1, 2, 3), (4, 5), (0, 1)] 1728 >>> set_partition(li) 1729 ([(5, 6, 7, 8), (1, 2, 3)], [0, 4]) 1730 1731 ''' 1732 part = [] 1733 for s in sorted(list(set(ssli)), key=len)[::-1]: 1734 #print(s, 1735 s_ = set(s).copy() 1736 if not any(set(s_).intersection(set(t)) for t in part): 1737 #print('inside:', s 1738 part.append(s) 1739 #else: print(part 1740 1741 missing = list(set(i for ll in ssli for i in ll) 1742 - set(i for ll in part for i in ll)) 1743 return part, missing 1744 1745 1746def set_remove_subs(ssli): 1747 '''remove sets that are subsets of another set from a list of tuples 1748 1749 Parameters 1750 ---------- 1751 ssli : list of tuples 1752 each tuple is considered as a set 1753 1754 Returns 1755 ------- 1756 part : list of tuples 1757 new list with subset tuples removed, it is sorted by set-length of tuples. The 1758 list contains original tuples, duplicate elements are not removed. 1759 1760 Examples 1761 -------- 1762 >>> set_remove_subs([(0, 1), (1, 2), (1, 2, 3), (0,)]) 1763 [(1, 2, 3), (0, 1)] 1764 >>> set_remove_subs([(0, 1), (1, 2), (1,1, 1, 2, 3), (0,)]) 1765 [(1, 1, 1, 2, 3), (0, 1)] 1766 1767 ''' 1768 #TODO: maybe convert all tuples to sets immediately, but I do not need the extra efficiency 1769 part = [] 1770 for s in sorted(list(set(ssli)), key=lambda x: len(set(x)))[::-1]: 1771 #print(s, 1772 #s_ = set(s).copy() 1773 if not any(set(s).issubset(set(t)) for t in part): 1774 #print('inside:', s 1775 part.append(s) 1776 #else: print(part 1777 1778## missing = list(set(i for ll in ssli for i in ll) 1779## - set(i for ll in part for i in ll)) 1780 return part 1781 1782 1783if __name__ == '__main__': 1784 1785 examples = ['tukey', 'tukeycrit', 'fdr', 'fdrmc', 'bonf', 'randmvn', 1786 'multicompdev', 'None']#[-1] 1787 1788 if 'tukey' in examples: 1789 #Example Tukey 1790 x = np.array([[0,0,1]]).T + np.random.randn(3, 20) 1791 print(Tukeythreegene(*x)) 1792 1793 # Example FDR 1794 # ------------ 1795 if ('fdr' in examples) or ('bonf' in examples): 1796 from .ex_multicomp import example_fdr_bonferroni 1797 example_fdr_bonferroni() 1798 1799 if 'fdrmc' in examples: 1800 mcres = mcfdr(nobs=100, nrepl=1000, ntests=30, ntrue=30, mu=0.1, alpha=0.05, rho=0.3) 1801 mcmeans = np.array(mcres).mean(0) 1802 print(mcmeans) 1803 print(mcmeans[0]/6., 1-mcmeans[1]/4.) 1804 print(mcmeans[:4], mcmeans[-4:]) 1805 1806 1807 if 'randmvn' in examples: 1808 rvsmvn = randmvn(0.8, (5000,5)) 1809 print(np.corrcoef(rvsmvn, rowvar=0)) 1810 print(rvsmvn.var(0)) 1811 1812 1813 if 'tukeycrit' in examples: 1814 print(get_tukeyQcrit(8, 8, alpha=0.05), 5.60) 1815 print(get_tukeyQcrit(8, 8, alpha=0.01), 7.47) 1816 1817 1818 if 'multicompdev' in examples: 1819 #development of kruskal-wallis multiple-comparison 1820 #example from matlab file exchange 1821 1822 X = np.array([[7.68, 1], [7.69, 1], [7.70, 1], [7.70, 1], [7.72, 1], 1823 [7.73, 1], [7.73, 1], [7.76, 1], [7.71, 2], [7.73, 2], 1824 [7.74, 2], [7.74, 2], [7.78, 2], [7.78, 2], [7.80, 2], 1825 [7.81, 2], [7.74, 3], [7.75, 3], [7.77, 3], [7.78, 3], 1826 [7.80, 3], [7.81, 3], [7.84, 3], [7.71, 4], [7.71, 4], 1827 [7.74, 4], [7.79, 4], [7.81, 4], [7.85, 4], [7.87, 4], 1828 [7.91, 4]]) 1829 xli = [X[X[:,1]==k,0] for k in range(1,5)] 1830 xranks = stats.rankdata(X[:,0]) 1831 xranksli = [xranks[X[:,1]==k] for k in range(1,5)] 1832 xnobs = np.array([len(xval) for xval in xli]) 1833 meanranks = [item.mean() for item in xranksli] 1834 sumranks = [item.sum() for item in xranksli] 1835 # equivalent function 1836 #from scipy import special 1837 #-np.sqrt(2.)*special.erfcinv(2-0.5) == stats.norm.isf(0.25) 1838 stats.norm.sf(0.67448975019608171) 1839 stats.norm.isf(0.25) 1840 1841 mrs = np.sort(meanranks) 1842 v1, v2 = np.triu_indices(4,1) 1843 print('\nsorted rank differences') 1844 print(mrs[v2] - mrs[v1]) 1845 diffidx = np.argsort(mrs[v2] - mrs[v1])[::-1] 1846 mrs[v2[diffidx]] - mrs[v1[diffidx]] 1847 1848 print('\nkruskal for all pairs') 1849 for i,j in zip(v2[diffidx], v1[diffidx]): 1850 print(i,j, stats.kruskal(xli[i], xli[j])) 1851 mwu, mwupval = stats.mannwhitneyu(xli[i], xli[j], use_continuity=False) 1852 print(mwu, mwupval*2, mwupval*2<0.05/6., mwupval*2<0.1/6.) 1853 1854 1855 1856 1857 1858 uni, intlab = np.unique(X[:,0], return_inverse=True) 1859 groupnobs = np.bincount(intlab) 1860 groupxsum = np.bincount(intlab, weights=X[:,0]) 1861 groupxmean = groupxsum * 1.0 / groupnobs 1862 1863 rankraw = X[:,0].argsort().argsort() 1864 groupranksum = np.bincount(intlab, weights=rankraw) 1865 # start at 1 for stats.rankdata : 1866 grouprankmean = groupranksum * 1.0 / groupnobs + 1 1867 assert_almost_equal(grouprankmean[intlab], stats.rankdata(X[:,0]), 15) 1868 gs = GroupsStats(X, useranks=True) 1869 print('\ngroupmeanfilter and grouprankmeans') 1870 print(gs.groupmeanfilter) 1871 print(grouprankmean[intlab]) 1872 #the following has changed 1873 #assert_almost_equal(gs.groupmeanfilter, stats.rankdata(X[:,0]), 15) 1874 1875 xuni, xintlab = np.unique(X[:,0], return_inverse=True) 1876 gs2 = GroupsStats(np.column_stack([X[:,0], xintlab]), useranks=True) 1877 #assert_almost_equal(gs2.groupmeanfilter, stats.rankdata(X[:,0]), 15) 1878 1879 rankbincount = np.bincount(xranks.astype(int)) 1880 nties = rankbincount[rankbincount > 1] 1881 ntot = float(len(xranks)) 1882 tiecorrection = 1 - (nties**3 - nties).sum()/(ntot**3 - ntot) 1883 assert_almost_equal(tiecorrection, stats.tiecorrect(xranks),15) 1884 print('\ntiecorrection for data and ranks') 1885 print(tiecorrection) 1886 print(tiecorrect(xranks)) 1887 1888 tot = X.shape[0] 1889 t=500 #168 1890 f=(tot*(tot+1.)/12.)-(t/(6.*(tot-1.))) 1891 f=(tot*(tot+1.)/12.)/stats.tiecorrect(xranks) 1892 print('\npairs of mean rank differences') 1893 for i,j in zip(v2[diffidx], v1[diffidx]): 1894 #pdiff = np.abs(mrs[i] - mrs[j]) 1895 pdiff = np.abs(meanranks[i] - meanranks[j]) 1896 se = np.sqrt(f * np.sum(1./xnobs[[i,j]] )) #np.array([8,8]))) #Fixme groupnobs[[i,j]] )) 1897 print(i,j, pdiff, se, pdiff/se, pdiff/se>2.6310) 1898 1899 multicomp = MultiComparison(*X.T) 1900 multicomp.kruskal() 1901 gsr = GroupsStats(X, useranks=True) 1902 1903 print('\nexamples for kruskal multicomparison') 1904 for i in range(10): 1905 x1, x2 = (np.random.randn(30,2) + np.array([0, 0.5])).T 1906 skw = stats.kruskal(x1, x2) 1907 mc2=MultiComparison(np.r_[x1, x2], np.r_[np.zeros(len(x1)), np.ones(len(x2))]) 1908 newskw = mc2.kruskal() 1909 print(skw, np.sqrt(skw[0]), skw[1]-newskw, (newskw/skw[1]-1)*100) 1910 1911 tablett, restt, arrtt = multicomp.allpairtest(stats.ttest_ind) 1912 tablemw, resmw, arrmw = multicomp.allpairtest(stats.mannwhitneyu) 1913 print('') 1914 print(tablett) 1915 print('') 1916 print(tablemw) 1917 tablemwhs, resmw, arrmw = multicomp.allpairtest(stats.mannwhitneyu, method='hs') 1918 print('') 1919 print(tablemwhs) 1920 1921 if 'last' in examples: 1922 xli = (np.random.randn(60,4) + np.array([0, 0, 0.5, 0.5])).T 1923 #Xrvs = np.array(catstack(xli)) 1924 xrvs, xrvsgr = catstack(xli) 1925 multicompr = MultiComparison(xrvs, xrvsgr) 1926 tablett, restt, arrtt = multicompr.allpairtest(stats.ttest_ind) 1927 print(tablett) 1928 1929 1930 xli=[[8,10,9,10,9],[7,8,5,8,5],[4,8,7,5,7]] 1931 x, labels = catstack(xli) 1932 gs4 = GroupsStats(np.column_stack([x, labels])) 1933 print(gs4.groupvarwithin()) 1934 1935 1936 #test_tukeyhsd() #moved to test_multi.py 1937 1938 gmeans = np.array([ 7.71375, 7.76125, 7.78428571, 7.79875]) 1939 gnobs = np.array([8, 8, 7, 8]) 1940 sd = StepDown(gmeans, gnobs, 0.001, [27]) 1941 1942 #example from BKY 1943 pvals = [0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 1944 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000 ] 1945 1946 #same number of rejection as in BKY paper: 1947 #single step-up:4, two-stage:8, iterated two-step:9 1948 #also alpha_star is the same as theirs for TST 1949 print(fdrcorrection0(pvals, alpha=0.05, method='indep')) 1950 print(fdrcorrection_twostage(pvals, alpha=0.05, iter=False)) 1951 res_tst = fdrcorrection_twostage(pvals, alpha=0.05, iter=False) 1952 assert_almost_equal([0.047619, 0.0649], res_tst[-1][:2],3) #alpha_star for stage 2 1953 assert_equal(8, res_tst[0].sum()) 1954 print(fdrcorrection_twostage(pvals, alpha=0.05, iter=True)) 1955 print('fdr_gbs', multipletests(pvals, alpha=0.05, method='fdr_gbs')) 1956 #multicontrast_pvalues(tstat, tcorr, df) 1957 tukey_pvalues(3.649, 3, 16) 1958