"""

from pystatsmodels mailinglist 20100524

Notes:
 - unfinished, unverified, but most parts seem to work in MonteCarlo
 - one example taken from lecture notes looks ok
 - needs cases with non-monotonic inequality for test to see difference between
   one-step, step-up and step-down procedures
 - FDR does not look really better then Bonferoni in the MC examples that I tried
update:
 - now tested against R, stats and multtest,
   I have all of their methods for p-value correction
 - getting Hommel was impossible until I found reference for pvalue correction
 - now, since I have p-values correction, some of the original tests (rej/norej)
   implementation is not really needed anymore. I think I keep it for reference.
   Test procedure for Hommel in development session log
 - I have not updated other functions and classes in here.
   - multtest has some good helper function according to docs
 - still need to update references, the real papers
 - fdr with estimated true hypothesis still missing
 - multiple comparison procedures incomplete or missing
 - I will get multiple comparison for now only for independent case, which might
   be conservative in correlated case (?).


some References:

Gibbons, Jean Dickinson and Chakraborti Subhabrata, 2003, Nonparametric Statistical
Inference, Fourth Edition, Marcel Dekker
    p.363: 10.4 THE KRUSKAL-WALLIS ONE-WAY ANOVA TEST AND MULTIPLE COMPARISONS
    p.367: multiple comparison for kruskal formula used in multicomp.kruskal

Sheskin, David J., 2004, Handbook of Parametric and Nonparametric Statistical
Procedures, 3rd ed., Chapman&Hall/CRC
    Test 21: The Single-Factor Between-Subjects Analysis of Variance
    Test 22: The Kruskal-Wallis One-Way Analysis of Variance by Ranks Test

Zwillinger, Daniel and Stephen Kokoska, 2000, CRC standard probability and
statistics tables and formulae, Chapman&Hall/CRC
    14.9 WILCOXON RANKSUM (MANN WHITNEY) TEST


S. Paul Wright, Adjusted P-Values for Simultaneous Inference, Biometrics
    Vol. 48, No. 4 (Dec., 1992), pp. 1005-1013, International Biometric Society
    Stable URL: http://www.jstor.org/stable/2532694
 (p-value correction for Hommel in appendix)

for multicomparison

new book "multiple comparison in R"
Hsu is a good reference but I do not have it.


Author: Josef Pktd and example from H Raja and rewrite from Vincent Davis


TODO
----
* name of function multipletests, rename to something like pvalue_correction?


"""
from statsmodels.compat.python import lrange, lzip

import copy
import math
from typing import NamedTuple

import numpy as np
from numpy.testing import assert_almost_equal, assert_equal
import pandas as pd
from scipy import interpolate, stats

from statsmodels.graphics import utils
from statsmodels.iolib.table import SimpleTable
from statsmodels.tools.sm_exceptions import ValueWarning

try:
    # Studentized Range in SciPy 1.7+
    from scipy.stats import studentized_range
except ImportError:
    from statsmodels.stats.libqsturng import psturng, qsturng

    class StudentizedRangeTuple(NamedTuple):
        ppf: float
        sf: float

    studentized_range = StudentizedRangeTuple(ppf=qsturng, sf=psturng)


qcrit = """
  2     3     4     5     6     7     8     9     10
5   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
6   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
7   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
8   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
9   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
10  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
11  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
12  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
13  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
14  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
15  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
16  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
17  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
18  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
19  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
20  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
24  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
30  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
40  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
60  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
120   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
infinity  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
"""

res = [line.split() for line in qcrit.replace("infinity", "9999").split("\n")]
c = np.array(res[2:-1]).astype(float)
# c[c==9999] = np.inf
ccols = np.arange(2, 11)
crows = c[:, 0]
cv005 = c[:, 1::2]
cv001 = c[:, 2::2]


def get_tukeyQcrit(k, df, alpha=0.05):
    """
    return critical values for Tukey's HSD (Q)

    Parameters
    ----------
    k : int in {2, ..., 10}
        number of tests
    df : int
        degrees of freedom of error term
    alpha : {0.05, 0.01}
        type 1 error, 1-confidence level



    not enough error checking for limitations
    """
    if alpha == 0.05:
        intp = interpolate.interp1d(crows, cv005[:, k - 2])
    elif alpha == 0.01:
        intp = interpolate.interp1d(crows, cv001[:, k - 2])
    else:
        raise ValueError("only implemented for alpha equal to 0.01 and 0.05")
    return intp(df)


def get_tukeyQcrit2(k, df, alpha=0.05):
    """
    return critical values for Tukey's HSD (Q)

    Parameters
    ----------
    k : int in {2, ..., 10}
        number of tests
    df : int
        degrees of freedom of error term
    alpha : {0.05, 0.01}
        type 1 error, 1-confidence level



    not enough error checking for limitations
    """
    return studentized_range.ppf(1 - alpha, k, df)


def get_tukey_pvalue(k, df, q):
    """
    return adjusted p-values for Tukey's HSD

    Parameters
    ----------
    k : int in {2, ..., 10}
        number of tests
    df : int
        degrees of freedom of error term
    q : scalar, array_like; q >= 0
        quantile value of Studentized Range

    """
    return studentized_range.sf(q, k, df)


def Tukeythreegene(first, second, third):
    # Performing the Tukey HSD post-hoc test for three genes
    # qwb = xlrd.open_workbook('F:/Lab/bioinformatics/qcrittable.xls')
    # # opening the workbook containing the q crit table
    # qwb.sheet_names()
    # qcrittable = qwb.sheet_by_name(u'Sheet1')

    # means of the three arrays
    firstmean = np.mean(first)
    secondmean = np.mean(second)
    thirdmean = np.mean(third)

    # standard deviations of the threearrays
    firststd = np.std(first)
    secondstd = np.std(second)
    thirdstd = np.std(third)

    # standard deviation squared of the three arrays
    firsts2 = math.pow(firststd, 2)
    seconds2 = math.pow(secondstd, 2)
    thirds2 = math.pow(thirdstd, 2)

    # numerator for mean square error
    mserrornum = firsts2 * 2 + seconds2 * 2 + thirds2 * 2
    # denominator for mean square error
    mserrorden = (len(first) + len(second) + len(third)) - 3
    mserror = mserrornum / mserrorden  # mean square error

    standarderror = math.sqrt(mserror / len(first))
    # standard error, which is square root of mserror and
    # the number of samples in a group

    # various degrees of freedom
    dftotal = len(first) + len(second) + len(third) - 1
    dfgroups = 2
    dferror = dftotal - dfgroups  # noqa: F841

    qcrit = 0.5  # fix arbitrary#qcrittable.cell(dftotal, 3).value
    qcrit = get_tukeyQcrit(3, dftotal, alpha=0.05)
    # getting the q critical value, for degrees of freedom total and 3 groups

    qtest3to1 = (math.fabs(thirdmean - firstmean)) / standarderror
    # calculating q test statistic values
    qtest3to2 = (math.fabs(thirdmean - secondmean)) / standarderror
    qtest2to1 = (math.fabs(secondmean - firstmean)) / standarderror

    conclusion = []

    # print(qcrit
    print(qtest3to1)
    print(qtest3to2)
    print(qtest2to1)

    # testing all q test statistic values to q critical values
    if qtest3to1 > qcrit:
        conclusion.append("3to1null")
    else:
        conclusion.append("3to1alt")
    if qtest3to2 > qcrit:
        conclusion.append("3to2null")
    else:
        conclusion.append("3to2alt")
    if qtest2to1 > qcrit:
        conclusion.append("2to1null")
    else:
        conclusion.append("2to1alt")

    return conclusion


# rewrite by Vincent
def Tukeythreegene2(genes):  # Performing the Tukey HSD post-hoc test for three genes
    """gend is a list, ie [first, second, third]"""
    #   qwb = xlrd.open_workbook('F:/Lab/bioinformatics/qcrittable.xls')
    # opening the workbook containing the q crit table
    #   qwb.sheet_names()
    #   qcrittable = qwb.sheet_by_name(u'Sheet1')

    means = []
    stds = []
    for gene in genes:
        means.append(np.mean(gene))
        std.append(np.std(gene))  # noqa:F821  See GH#5756

    # firstmean = np.mean(first) # means of the three arrays
    # secondmean = np.mean(second)
    # thirdmean = np.mean(third)

    # firststd = np.std(first) # standard deviations of the three arrays
    # secondstd = np.std(second)
    # thirdstd = np.std(third)

    stds2 = []
    for std in stds:
        stds2.append(math.pow(std, 2))

    # firsts2 = math.pow(firststd,2) # standard deviation squared of the three arrays
    # seconds2 = math.pow(secondstd,2)
    # thirds2 = math.pow(thirdstd,2)

    # Unused result, comented out
    # mserrornum = firsts2*2+seconds2*2+thirds2*2 # numerator for mean square error
    # mserrornum = sum(stds2)*2
    # mserrorden = (len(genes[0])+len(genes[1])+len(genes[2]))-3 # denominator for mean square error
    # mserrornum/mserrorden # mean square error


def catstack(args):
    x = np.hstack(args)
    labels = np.hstack([k * np.ones(len(arr)) for k, arr in enumerate(args)])
    return x, labels


def maxzero(x):
    """find all up zero crossings and return the index of the highest

    Not used anymore


    >>> np.random.seed(12345)
    >>> x = np.random.randn(8)
    >>> x
    array([-0.20470766,  0.47894334, -0.51943872, -0.5557303 ,  1.96578057,
            1.39340583,  0.09290788,  0.28174615])
    >>> maxzero(x)
    (4, array([1, 4]))


    no up-zero-crossing at end

    >>> np.random.seed(0)
    >>> x = np.random.randn(8)
    >>> x
    array([ 1.76405235,  0.40015721,  0.97873798,  2.2408932 ,  1.86755799,
           -0.97727788,  0.95008842, -0.15135721])
    >>> maxzero(x)
    (None, array([6]))
    """
    x = np.asarray(x)
    cond1 = x[:-1] < 0
    cond2 = x[1:] > 0
    # allzeros = np.nonzero(np.sign(x[:-1])*np.sign(x[1:]) <= 0)[0] + 1
    allzeros = np.nonzero((cond1 & cond2) | (x[1:] == 0))[0] + 1
    if x[-1] >= 0:
        maxz = max(allzeros)
    else:
        maxz = None
    return maxz, allzeros


def maxzerodown(x):
    """find all up zero crossings and return the index of the highest

    Not used anymore

    >>> np.random.seed(12345)
    >>> x = np.random.randn(8)
    >>> x
    array([-0.20470766,  0.47894334, -0.51943872, -0.5557303 ,  1.96578057,
            1.39340583,  0.09290788,  0.28174615])
    >>> maxzero(x)
    (4, array([1, 4]))


    no up-zero-crossing at end

    >>> np.random.seed(0)
    >>> x = np.random.randn(8)
    >>> x
    array([ 1.76405235,  0.40015721,  0.97873798,  2.2408932 ,  1.86755799,
           -0.97727788,  0.95008842, -0.15135721])
    >>> maxzero(x)
    (None, array([6]))
    """
    x = np.asarray(x)
    cond1 = x[:-1] > 0
    cond2 = x[1:] < 0
    # allzeros = np.nonzero(np.sign(x[:-1])*np.sign(x[1:]) <= 0)[0] + 1
    allzeros = np.nonzero((cond1 & cond2) | (x[1:] == 0))[0] + 1
    if x[-1] <= 0:
        maxz = max(allzeros)
    else:
        maxz = None
    return maxz, allzeros


def rejectionline(n, alpha=0.5):
    """reference line for rejection in multiple tests

    Not used anymore

    from: section 3.2, page 60
    """
    t = np.arange(n) / float(n)
    frej = t / (t * (1 - alpha) + alpha)
    return frej


# I do not remember what I changed or why 2 versions,
# this follows german diss ???  with rline
# this might be useful if the null hypothesis is not "all effects are zero"
# rename to _bak and working again on fdrcorrection0
def fdrcorrection_bak(pvals, alpha=0.05, method="indep"):
    """Reject False discovery rate correction for pvalues

    Old version, to be deleted


    missing: methods that estimate fraction of true hypotheses

    """
    from statsmodels.stats.multitest import _ecdf as ecdf

    pvals = np.asarray(pvals)

    pvals_sortind = np.argsort(pvals)
    pvals_sorted = pvals[pvals_sortind]
    pecdf = ecdf(pvals_sorted)
    if method in ["i", "indep", "p", "poscorr"]:
        rline = pvals_sorted / alpha
    elif method in ["n", "negcorr"]:
        cm = np.sum(1.0 / np.arange(1, len(pvals)))
        rline = pvals_sorted / alpha * cm
    elif method in ["g", "onegcorr"]:  # what's this ? german diss
        rline = pvals_sorted / (pvals_sorted * (1 - alpha) + alpha)
    elif method in ["oth", "o2negcorr"]:  # other invalid, cut-paste
        cm = np.sum(np.arange(len(pvals)))
        rline = pvals_sorted / alpha / cm
    else:
        raise ValueError("method not available")

    reject = pecdf >= rline
    if reject.any():
        rejectmax = max(np.nonzero(reject)[0])
    else:
        rejectmax = 0
    reject[:rejectmax] = True
    return reject[pvals_sortind.argsort()]


def mcfdr(nrepl=100, nobs=50, ntests=10, ntrue=6, mu=0.5, alpha=0.05, rho=0.0):
    """MonteCarlo to test fdrcorrection"""
    from statsmodels.stats.multitest import fdrcorrection as fdrcorrection0

    # Unused result, commented out
    # ntests - ntrue
    locs = np.array([0.0] * ntrue + [mu] * (ntests - ntrue))
    results = []
    for _ in range(nrepl):
        # rvs = locs + stats.norm.rvs(size=(nobs, ntests))
        rvs = locs + randmvn(rho, size=(nobs, ntests))
        tt, tpval = stats.ttest_1samp(rvs, 0)
        res = fdrcorrection_bak(np.abs(tpval), alpha=alpha, method="i")
        res0 = fdrcorrection0(np.abs(tpval), alpha=alpha)
        # res and res0 give the same results
        results.append(
            [np.sum(res[:ntrue]), np.sum(res[ntrue:])]
            + [np.sum(res0[:ntrue]), np.sum(res0[ntrue:])]
            + res.tolist()
            + np.sort(tpval).tolist()
            + [np.sum(tpval[:ntrue] < alpha), np.sum(tpval[ntrue:] < alpha)]
            + [
                np.sum(tpval[:ntrue] < alpha / ntests),
                np.sum(tpval[ntrue:] < alpha / ntests),
            ]
        )
    return np.array(results)


def randmvn(rho, size=(1, 2), standardize=False):
    """create random draws from equi-correlated multivariate normal distribution

    Parameters
    ----------
    rho : float
        correlation coefficient
    size : tuple of int
        size is interpreted (nobs, nvars) where each row

    Returns
    -------
    rvs : ndarray
        nobs by nvars where each row is a independent random draw of nvars-
        dimensional correlated rvs

    """
    nobs, nvars = size
    if 0 < rho < 1:
        rvs = np.random.randn(nobs, nvars + 1)
        rvs2 = rvs[:, :-1] * np.sqrt(1 - rho) + rvs[:, -1:] * np.sqrt(rho)
    elif rho == 0:
        rvs2 = np.random.randn(nobs, nvars)
    elif rho < 0:
        if rho < -1.0 / (nvars - 1):
            raise ValueError("rho has to be larger than -1./(nvars-1)")
        elif rho == -1.0 / (nvars - 1):
            rho = -1.0 / (nvars - 1 + 1e-10)  # barely positive definite
        # use Cholesky
        A = rho * np.ones((nvars, nvars)) + (1 - rho) * np.eye(nvars)
        rvs2 = np.dot(np.random.randn(nobs, nvars), np.linalg.cholesky(A).T)
    if standardize:
        rvs2 = stats.zscore(rvs2)
    return rvs2


# ============================
#
# Part 2: Multiple comparisons and independent samples tests
#
# ============================


def tiecorrect(xranks):
    """

    should be equivalent of scipy.stats.tiecorrect

    """
    # casting to int rounds down, but not relevant for this case
    rankbincount = np.bincount(np.asarray(xranks, dtype=int))
    nties = rankbincount[rankbincount > 1]
    ntot = float(len(xranks))
    tiecorrection = 1 - (nties**3 - nties).sum() / (ntot**3 - ntot)
    return tiecorrection


class GroupsStats:
    """
    statistics by groups (another version)

    groupstats as a class with lazy evaluation (not yet - decorators are still
    missing)

    written this time as equivalent of scipy.stats.rankdata
    gs = GroupsStats(X, useranks=True)
    assert_almost_equal(gs.groupmeanfilter, stats.rankdata(X[:,0]), 15)

    TODO: incomplete doc strings

    """

    def __init__(self, x, useranks=False, uni=None, intlab=None):
        """descriptive statistics by groups

        Parameters
        ----------
        x : ndarray, 2d
            first column data, second column group labels
        useranks : bool
            if true, then use ranks as data corresponding to the
            scipy.stats.rankdata definition (start at 1, ties get mean)
        uni, intlab : arrays (optional)
            to avoid call to unique, these can be given as inputs


        """
        self.x = np.asarray(x)
        if intlab is None:
            uni, intlab = np.unique(x[:, 1], return_inverse=True)
        elif uni is None:
            uni = np.unique(x[:, 1])

        self.useranks = useranks

        self.uni = uni
        self.intlab = intlab
        self.groupnobs = np.bincount(intlab)

        # temporary until separated and made all lazy
        self.runbasic(useranks=useranks)

    def runbasic_old(self, useranks=False):
        """runbasic_old"""
        # check: refactoring screwed up case useranks=True

        # groupxsum = np.bincount(intlab, weights=X[:,0])
        # groupxmean = groupxsum * 1.0 / groupnobs
        x = self.x
        if useranks:
            self.xx = x[:, 1].argsort().argsort() + 1  # rankraw
        else:
            self.xx = x[:, 0]
        self.groupsum = groupranksum = np.bincount(self.intlab, weights=self.xx)
        # print('groupranksum', groupranksum, groupranksum.shape, self.groupnobs.shape
        # start at 1 for stats.rankdata :
        self.groupmean = grouprankmean = groupranksum * 1.0 / self.groupnobs  # + 1
        self.groupmeanfilter = grouprankmean[self.intlab]
        # return grouprankmean[intlab]

    def runbasic(self, useranks=False):
        """runbasic"""
        # check: refactoring screwed up case useranks=True

        # groupxsum = np.bincount(intlab, weights=X[:,0])
        # groupxmean = groupxsum * 1.0 / groupnobs
        x = self.x
        if useranks:
            xuni, xintlab = np.unique(x[:, 0], return_inverse=True)
            ranksraw = x[:, 0].argsort().argsort() + 1  # rankraw
            self.xx = GroupsStats(
                np.column_stack([ranksraw, xintlab]), useranks=False
            ).groupmeanfilter
        else:
            self.xx = x[:, 0]
        self.groupsum = groupranksum = np.bincount(self.intlab, weights=self.xx)
        # print('groupranksum', groupranksum, groupranksum.shape, self.groupnobs.shape
        # start at 1 for stats.rankdata :
        self.groupmean = grouprankmean = groupranksum * 1.0 / self.groupnobs  # + 1
        self.groupmeanfilter = grouprankmean[self.intlab]
        # return grouprankmean[intlab]

    def groupdemean(self):
        """groupdemean"""
        return self.xx - self.groupmeanfilter

    def groupsswithin(self):
        """groupsswithin"""
        xtmp = self.groupdemean()
        return np.bincount(self.intlab, weights=xtmp**2)

    def groupvarwithin(self):
        """groupvarwithin"""
        return self.groupsswithin() / (self.groupnobs - 1)  # .sum()


class TukeyHSDResults:
    """Results from Tukey HSD test, with additional plot methods

    Can also compute and plot additional post-hoc evaluations using this
    results class.

    Attributes
    ----------
    reject : array of boolean, True if we reject Null for group pair
    meandiffs : pairwise mean differences
    confint : confidence interval for pairwise mean differences
    std_pairs : standard deviation of pairwise mean differences
    q_crit : critical value of studentized range statistic at given alpha
    halfwidths : half widths of simultaneous confidence interval
    pvalues : adjusted p-values from the HSD test

    Notes
    -----
    halfwidths is only available after call to `plot_simultaneous`.

    Other attributes contain information about the data from the
    MultiComparison instance: data, df_total, groups, groupsunique, variance.
    """

    def __init__(
        self,
        mc_object,
        results_table,
        q_crit,
        reject=None,
        meandiffs=None,
        std_pairs=None,
        confint=None,
        df_total=None,
        reject2=None,
        variance=None,
        pvalues=None,
        alpha=None,
        group_t=None,
        group_c=None,
    ):
        self._multicomp = mc_object
        self._results_table = results_table
        self.q_crit = q_crit
        self.reject = reject
        self.meandiffs = meandiffs
        self.std_pairs = std_pairs
        self.confint = confint
        self.df_total = df_total
        self.reject2 = reject2
        self.variance = variance
        self.pvalues = pvalues
        self.alpha = alpha
        self.group_t = group_t
        self.group_c = group_c
        # Taken out of _multicomp for ease of access for unknowledgeable users
        self.data = self._multicomp.data
        self.groups = self._multicomp.groups
        self.groupsunique = self._multicomp.groupsunique

        if np.size(df_total) > 1:  # or should it be np.isscalar
            # assume we have Games-Howell, unequal var case
            self._qcrit_hsd = None
        else:
            self._qcrit_hsd = q_crit

        nobs_group = self._multicomp.groupstats.groupnobs
        self.df_total_hsd = np.sum(nobs_group - 1)

    def __str__(self):
        return str(self._results_table)

    def summary(self):
        """Summary table that can be printed"""
        return self._results_table

    def summary_frame(self):
        """Summary DataFrame

        The group columns are labeled as "group_t" and "group_c" with mean
        difference defined as treatment minus control.
        This should be less confusing than numeric labels group1 and group2.

        Returns
        -------
        pandas.DataFrame

        Notes
        -----
        The number of columns will likely increase in a future version of
        statsmodels. Do not use numeric indices for the DataFrame in order
        to be robust to the addition of columns.
        """
        frame = pd.DataFrame(
            {
                "group_t": self.group_t,
                "group_c": self.group_c,
                "meandiff": self.meandiffs,
                "p-adj": self.pvalues,
                "lower": self.confint[:, 0],
                "upper": self.confint[:, 1],
                "reject": self.reject,
            }
        )
        return frame

    def _get_q_crit(self, hsd=True, alpha=None):
        n_means = len(self.groupsunique)

        if alpha is None:
            alpha = self.alpha
            use_attr = True

        if hsd is True:
            if use_attr and self._qcrit_hsd is not None:
                q_crit = self._qcrit_hsd
            else:  # compute it
                q_crit = get_tukeyQcrit2(n_means, self.df_total_hsd, alpha=alpha)
            if use_attr:
                self._qcrit_hsd = q_crit
        else:
            raise NotImplementedError("not yet")

        return q_crit

    def _simultaneous_ci(self):
        """Compute simultaneous confidence intervals for comparison of means."""

        q_crit_hsd = self._get_q_crit(hsd=True)

        self.halfwidths = simultaneous_ci(
            q_crit_hsd,
            self.variance,
            self._multicomp.groupstats.groupnobs,
            self._multicomp.pairindices,
        )

    def plot_simultaneous(
        self, comparison_name=None, ax=None, figsize=(10, 6), xlabel=None, ylabel=None
    ):
        """Plot a universal confidence interval of each group mean

        Visualize significant differences in a plot with one confidence
        interval per group instead of all pairwise confidence intervals.

        Parameters
        ----------
        comparison_name : str, optional
            if provided, plot_intervals will color code all groups that are
            significantly different from the comparison_name red, and will
            color code insignificant groups gray. Otherwise, all intervals will
            just be plotted in black.
        ax : matplotlib axis, optional
            An axis handle on which to attach the plot.
        figsize : tuple, optional
            tuple for the size of the figure generated
        xlabel : str, optional
            Name to be displayed on x axis
        ylabel : str, optional
            Name to be displayed on y axis

        Returns
        -------
        Figure
            handle to figure object containing interval plots

        Notes
        -----
        Multiple comparison tests are nice, but lack a good way to be
        visualized. If you have, say, 6 groups, showing a graph of the means
        between each group will require 15 confidence intervals.
        Instead, we can visualize inter-group differences with a single
        interval for each group mean. Hochberg et al. [1] first proposed this
        idea and used Tukey's Q critical value to compute the interval widths.
        Unlike plotting the differences in the means and their respective
        confidence intervals, any two pairs can be compared for significance
        by looking for overlap.

        The derivation in Hochberg and Tamhane is under the equal variance
        assumption. We use the same computation in the case of unequal
        variances, however, with replacement of the common pooled variance
        by the unequal estimates of the whithin group variances.
        This provides a plot that looks more informative and plausible in the
        case where there are large differences in variances. In the equal
        sample size and equal variance case, the confidence intervals computed
        by the two methods, equal and unequal variance, are very close to
        each other in larger samples.

        References
        ----------
        .. [*] Hochberg, Y., and A. C. Tamhane. Multiple Comparison Procedures.
               Hoboken, NJ: John Wiley & Sons, 1987.

        Examples
        --------
        >>> from statsmodels.examples.try_tukey_hsd import cylinders, cyl_labels
        >>> from statsmodels.stats.multicomp import MultiComparison
        >>> cardata = MultiComparison(cylinders, cyl_labels)
        >>> results = cardata.tukeyhsd()
        >>> results.plot_simultaneous()
        <matplotlib.figure.Figure at 0x...>

        This example shows an example plot comparing significant differences
        in group means. Significant differences at the alpha=0.05 level can be
        identified by intervals that do not overlap (i.e. USA vs Japan,
        USA vs Germany).

        >>> results.plot_simultaneous(comparison_name="USA")
        <matplotlib.figure.Figure at 0x...>

        Optionally provide one of the group names to color code the plot to
        highlight group means different from comparison_name.
        """
        fig, ax1 = utils.create_mpl_ax(ax)
        if figsize is not None:
            fig.set_size_inches(figsize)
        if getattr(self, "halfwidths", None) is None:
            self._simultaneous_ci()
        means = self._multicomp.groupstats.groupmean

        sigidx = []
        nsigidx = []
        minrange = [means[i] - self.halfwidths[i] for i in range(len(means))]
        maxrange = [means[i] + self.halfwidths[i] for i in range(len(means))]

        if comparison_name is None:
            ax1.errorbar(
                means,
                lrange(len(means)),
                xerr=self.halfwidths,
                marker="o",
                linestyle="None",
                color="k",
                ecolor="k",
            )
        else:
            if comparison_name not in self.groupsunique:
                raise ValueError("comparison_name not found in group names.")
            midx = np.where(self.groupsunique == comparison_name)[0][0]
            for i in range(len(means)):
                if self.groupsunique[i] == comparison_name:
                    continue
                if (
                    min(maxrange[i], maxrange[midx]) - max(minrange[i], minrange[midx])
                    < 0
                ):
                    sigidx.append(i)
                else:
                    nsigidx.append(i)
            # Plot the main comparison
            ax1.errorbar(
                means[midx],
                midx,
                xerr=self.halfwidths[midx],
                marker="o",
                linestyle="None",
                color="b",
                ecolor="b",
            )
            ax1.plot(
                [minrange[midx]] * 2,
                [-1, self._multicomp.ngroups],
                linestyle="--",
                color="0.7",
            )
            ax1.plot(
                [maxrange[midx]] * 2,
                [-1, self._multicomp.ngroups],
                linestyle="--",
                color="0.7",
            )
            # Plot those that are significantly different
            if len(sigidx) > 0:
                ax1.errorbar(
                    means[sigidx],
                    sigidx,
                    xerr=self.halfwidths[sigidx],
                    marker="o",
                    linestyle="None",
                    color="r",
                    ecolor="r",
                )
            # Plot those that are not significantly different
            if len(nsigidx) > 0:
                ax1.errorbar(
                    means[nsigidx],
                    nsigidx,
                    xerr=self.halfwidths[nsigidx],
                    marker="o",
                    linestyle="None",
                    color="0.5",
                    ecolor="0.5",
                )

        ax1.set_title("Multiple Comparisons Between All Pairs (Tukey)")
        r = np.max(maxrange) - np.min(minrange)
        ax1.set_ylim([-1, self._multicomp.ngroups])
        ax1.set_xlim([np.min(minrange) - r / 10.0, np.max(maxrange) + r / 10.0])
        ylbls = [""] + self.groupsunique.astype(str).tolist() + [""]
        ax1.set_yticks(np.arange(-1, len(means) + 1))
        ax1.set_yticklabels(ylbls)
        ax1.set_xlabel(xlabel if xlabel is not None else "")
        ax1.set_ylabel(ylabel if ylabel is not None else "")
        return fig


class MultiComparison:
    """Tests for multiple comparisons

    Parameters
    ----------
    data : ndarray
        independent data samples
    groups : ndarray
        group labels corresponding to each data point
    group_order : list[str], optional
        the desired order for the group mean results to be reported in. If
        not specified, results are reported in increasing order.
        If group_order does not contain all labels that are in groups, then
        only those observations are kept that have a label in group_order.

    """

    def __init__(self, data, groups, group_order=None):
        if len(data) != len(groups):
            raise ValueError(
                "data has %d elements and groups has %d" % (len(data), len(groups))
            )
        self.data = np.asarray(data)
        self.groups = groups = np.asarray(groups)

