docs: def posthoc_durbin(a, y_col=None, block_col=None, group_col=None, melted=False, sort=False, p_adjust=None):
  
      '''Pairwise post hoc test for multiple comparisons of rank sums according to
      Durbin and Conover for a two-way balanced incomplete block design (BIBD). See
      references for additional information [1]_, [2]_.
  
      Parameters
      ----------
      a : array_like or pandas DataFrame object
          An array, any object exposing the array interface or a pandas
          DataFrame.
  
          If `melted` is set to False (default), `a` is a typical matrix of block design,
          i.e. rows are blocks, and columns are groups. In this case you do
          not need to specify col arguments.
  
          If `a` is an array and `melted` is set to True,
          y_col, block_col and group_col must specify the indices of columns
          containing elements of correspondary type.
  
          If `a` is a Pandas DataFrame and `melted` is set to True,
          y_col, block_col and group_col must specify columns names (string).
  
      y_col : str or int
          Must be specified if `a` is a pandas DataFrame object.
          Name of the column that contains y data.
  
      block_col : str or int
          Must be specified if `a` is a pandas DataFrame object.
          Name of the column that contains blocking factor values.
  
      group_col : str or int
          Must be specified if `a` is a pandas DataFrame object.
          Name of the column that contains treatment (group) factor values.
  
      melted : bool, optional
          Specifies if data are given as melted columns "y", "blocks", and
          "groups".
  
      sort : bool, optional
          If True, sort data by block and group columns.
  
      p_adjust : str, optional
          Method for adjusting p values. See statsmodels.sandbox.stats.multicomp
          for details. Available methods are:
          'bonferroni' : one-step correction
          'sidak' : one-step correction
          'holm-sidak' : step-down method using Sidak adjustments
          'holm' : step-down method using Bonferroni adjustments
          'simes-hochberg' : step-up method  (independent)
          'hommel' : closed method based on Simes tests (non-negative)
          'fdr_bh' : Benjamini/Hochberg  (non-negative)
          'fdr_by' : Benjamini/Yekutieli (negative)
          'fdr_tsbh' : two stage fdr correction (non-negative)
          'fdr_tsbky' : two stage fdr correction (non-negative)
  
      Returns
      -------
      Pandas DataFrame containing p values.
  
      References
      ----------
      .. [1] W. J. Conover and R. L. Iman (1979), On multiple-comparisons procedures,
          Tech. Rep. LA-7677-MS, Los Alamos Scientific Laboratory.
      .. [2] W. J. Conover (1999), Practical nonparametric Statistics,
          3rd. edition, Wiley.
  
      Examples
      --------
      >>> x = np.array([[31,27,24],[31,28,31],[45,29,46],[21,18,48],[42,36,46],[32,17,40]])
      >>> sp.posthoc_durbin(x)
  
      '''
  
      if melted and not all([block_col, group_col, y_col]):
          raise ValueError('block_col, group_col, y_col should be explicitly specified if using melted data')
  
      def compare_stats(i, j):
          dif = np.abs(Rj[groups[i]] - Rj[groups[j]])
          tval = dif / denom
          pval = 2. * ss.t.sf(np.abs(tval), df = df)
          return pval
  
      x, y_col, group_col, block_col = __convert_to_block_df(a, y_col, group_col, block_col, melted)
  
      if not sort:
          x[group_col] = Categorical(x[group_col], categories=x[group_col].unique(), ordered=True)
          x[block_col] = Categorical(x[block_col], categories=x[block_col].unique(), ordered=True)
      x.sort_values(by=[block_col, group_col], ascending=True, inplace=True)
      x.dropna(inplace=True)
  
      groups = x[group_col].unique()
      t = len(groups)
      b = x[block_col].unique().size
      r = b
      k = t
      x['y_ranked'] = x.groupby(block_col)[y_col].rank()
      Rj = x.groupby(group_col)['y_ranked'].sum()
      A = (x['y_ranked'] ** 2).sum()
      C = (b * k * (k + 1) ** 2) / 4.
      D = (Rj ** 2).sum() - r * C
      T1 = (t - 1) / (A - C) * D
      denom = np.sqrt(((A - C) * 2 * r) / (b * k - b - t + 1) * (1 - T1 / (b * (k -1))))
      df = b * k - b - t + 1
  
      vs = np.zeros((t, t), dtype=np.float)
      combs = it.combinations(range(t), 2)
  
      tri_upper = np.triu_indices(vs.shape[0], 1)
      tri_lower = np.tril_indices(vs.shape[0], -1)
      vs[:,:] = 0
  
      for i, j in combs:
          vs[i, j] = compare_stats(i, j)
  
      if p_adjust:
          vs[tri_upper] = multipletests(vs[tri_upper], method = p_adjust)[1]
  
      vs[tri_lower] = vs.T[tri_lower]
      np.fill_diagonal(vs, -1)
      return DataFrame(vs, index=groups, columns=groups)