code: def tukeyhsd(
      mean_all, nobs_all, var_all, df=None, alpha=0.05, q_crit=None
  ):
      """simultaneous Tukey HSD
  
  
      check: instead of sorting, I use absolute value of pairwise differences
      in means. That's irrelevant for the test, but maybe reporting actual
      differences would be better.
      CHANGED: meandiffs are with sign, studentized range uses abs
  
      q_crit added for testing
  
      TODO: error in variance calculation when nobs_all is scalar, missing 1/n
  
      """
      mean_all = np.asarray(mean_all)
      # check if or when other ones need to be arrays
  
      n_means = len(mean_all)
  
      if df is None:
          df = nobs_all - 1
  
      if (
          np.size(df) == 1
      ):  # assumes balanced samples with df = n - 1, n_i = n
          df_total = n_means * df
          df = np.ones(n_means) * df
      else:
          df_total = np.sum(df)
  
      if (np.size(nobs_all) == 1) and (np.size(var_all) == 1):
          # balanced sample sizes and homogenous variance
          var_pairs = (
              1.0 * var_all / nobs_all * np.ones((n_means, n_means))
          )
  
      elif np.size(var_all) == 1:
          # unequal sample sizes and homogenous variance
          var_pairs = var_all * varcorrection_pairs_unbalanced(
              nobs_all, srange=True
          )
      elif np.size(var_all) > 1:
          var_pairs, df_sum = varcorrection_pairs_unequal(
              nobs_all, var_all, df
          )
          var_pairs /= 2.0
          # check division by two for studentized range
  
      else:
          raise ValueError("not supposed to be here")
  
      # meandiffs_ = mean_all[:,None] - mean_all
      meandiffs_ = (
          mean_all - mean_all[:, None]
      )  # reverse sign, check with R example
      std_pairs_ = np.sqrt(var_pairs)
  
      # select all pairs from upper triangle of matrix
      idx1, idx2 = np.triu_indices(n_means, 1)
      meandiffs = meandiffs_[idx1, idx2]
      std_pairs = std_pairs_[idx1, idx2]
  
      st_range = (
          np.abs(meandiffs) / std_pairs
      )  # studentized range statistic
  
      df_total_ = max(df_total, 5)  # TODO: smallest df in table
      if q_crit is None:
          q_crit = get_tukeyQcrit2(n_means, df_total, alpha=alpha)
  
      pvalues = get_tukey_pvalue(n_means, df_total, st_range)
      # we need pvalues to be atleast_1d for iteration. see #6132
      pvalues = np.atleast_1d(pvalues)
  
      reject = st_range > q_crit
      crit_int = std_pairs * q_crit
      reject2 = np.abs(meandiffs) > crit_int
  
      confint = np.column_stack(
          (meandiffs - crit_int, meandiffs + crit_int)
      )
  
      return (
          (idx1, idx2),
          reject,
          meandiffs,
          std_pairs,
          confint,
          q_crit,
          df_total,
          reject2,
          pvalues,
      )
  
imports: from collections import namedtuple
  from statsmodels.compat.python import lzip, lrange
  import copy
  import math
  import numpy
  from numpy.testing import assert_almost_equal, assert_equal
  from scipy import stats, interpolate
  from statsmodels.iolib.table import SimpleTable
  from statsmodels.stats.multitest import multipletests, _ecdf, fdrcorrection, fdrcorrection_twostage
  from statsmodels.graphics import utils
  from statsmodels.tools.sm_exceptions import ValueWarning
label_desc: Use image node_id 20 for calling a global function with example usage: tukeyhsd(mean_all, nobs_all, var_all, df, alpha, q_crit) and returns: idx1,  idx2,  reject,  meandiffs,  std_pairs,  confint,  q_crit,  df_total,  reject2,  pvalues