code: def multivariate_stats(
      eigenvals, r_err_sscp, r_contrast, df_resid, tolerance=1e-8
  ):
      """
      For multivariate linear model Y = X * B
      Testing hypotheses
          L*B*M = 0
      where L is contrast matrix, B is the parameters of the
      multivariate linear model and M is dependent variable transform matrix.
          T = L*inv(X'X)*L'
          H = M'B'L'*inv(T)*LBM
          E =  M'(Y'Y - B'X'XB)M
  
      Parameters
      ----------
      eigenvals : ndarray
          The eigenvalues of inv(E + H)*H
      r_err_sscp : int
          Rank of E + H
      r_contrast : int
          Rank of T matrix
      df_resid : int
          Residual degree of freedom (n_samples minus n_variables of X)
      tolerance : float
          smaller than which eigenvalue is considered 0
  
      Returns
      -------
      A DataFrame
  
      References
      ----------
      .. [*] https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_introreg_sect012.htm
      """
      v = df_resid
      p = r_err_sscp
      q = r_contrast
      s = np.min([p, q])
      ind = eigenvals > tolerance
      # n_e = ind.sum()
      eigv2 = eigenvals[ind]
      eigv1 = np.array([i / (1 - i) for i in eigv2])
      m = (np.abs(p - q) - 1) / 2
      n = (v - p - 1) / 2
  
      cols = ["Value", "Num DF", "Den DF", "F Value", "Pr > F"]
      index = [
          "Wilks' lambda",
          "Pillai's trace",
          "Hotelling-Lawley trace",
          "Roy's greatest root",
      ]
      results = pd.DataFrame(columns=cols, index=index)
  
      def fn(x):
          return np.real([x])[0]
  
      results.loc["Wilks' lambda", "Value"] = fn(np.prod(1 - eigv2))
  
      results.loc["Pillai's trace", "Value"] = fn(eigv2.sum())
  
      results.loc["Hotelling-Lawley trace", "Value"] = fn(eigv1.sum())
  
      results.loc["Roy's greatest root", "Value"] = fn(eigv1.max())
  
      r = v - (p - q + 1) / 2
      u = (p * q - 2) / 4
      df1 = p * q
      if p * p + q * q - 5 > 0:
          t = np.sqrt((p * p * q * q - 4) / (p * p + q * q - 5))
      else:
          t = 1
      df2 = r * t - 2 * u
      lmd = results.loc["Wilks' lambda", "Value"]
      lmd = np.power(lmd, 1 / t)
      F = (1 - lmd) / lmd * df2 / df1
      results.loc["Wilks' lambda", "Num DF"] = df1
      results.loc["Wilks' lambda", "Den DF"] = df2
      results.loc["Wilks' lambda", "F Value"] = F
      pval = stats.f.sf(F, df1, df2)
      results.loc["Wilks' lambda", "Pr > F"] = pval
  
      V = results.loc["Pillai's trace", "Value"]
      df1 = s * (2 * m + s + 1)
      df2 = s * (2 * n + s + 1)
      F = df2 / df1 * V / (s - V)
      results.loc["Pillai's trace", "Num DF"] = df1
      results.loc["Pillai's trace", "Den DF"] = df2
      results.loc["Pillai's trace", "F Value"] = F
      pval = stats.f.sf(F, df1, df2)
      results.loc["Pillai's trace", "Pr > F"] = pval
  
      U = results.loc["Hotelling-Lawley trace", "Value"]
      if n > 0:
          b = (p + 2 * n) * (q + 2 * n) / 2 / (2 * n + 1) / (n - 1)
          df1 = p * q
          df2 = 4 + (p * q + 2) / (b - 1)
          c = (df2 - 2) / 2 / n
          F = df2 / df1 * U / c
      else:
          df1 = s * (2 * m + s + 1)
          df2 = s * (s * n + 1)
          F = df2 / df1 / s * U
      results.loc["Hotelling-Lawley trace", "Num DF"] = df1
      results.loc["Hotelling-Lawley trace", "Den DF"] = df2
      results.loc["Hotelling-Lawley trace", "F Value"] = F
      pval = stats.f.sf(F, df1, df2)
      results.loc["Hotelling-Lawley trace", "Pr > F"] = pval
  
      sigma = results.loc["Roy's greatest root", "Value"]
      r = np.max([p, q])
      df1 = r
      df2 = v - r + q
      F = df2 / df1 * sigma
      results.loc["Roy's greatest root", "Num DF"] = df1
      results.loc["Roy's greatest root", "Den DF"] = df2
      results.loc["Roy's greatest root", "F Value"] = F
      pval = stats.f.sf(F, df1, df2)
      results.loc["Roy's greatest root", "Pr > F"] = pval
      return results
  
imports: import numpy
  from numpy.linalg import eigvals, inv, solve, matrix_rank, pinv, svd
  from scipy import stats
  import pandas
  from patsy import DesignInfo
  from statsmodels.compat.pandas import Substitution
  from statsmodels.base.model import Model, LikelihoodModelResults
  import statsmodels.base.wrapper
  from statsmodels.regression.linear_model import RegressionResultsWrapper
  from statsmodels.iolib import summary2
  from statsmodels.tools.decorators import cache_readonly
label_desc: Use image node_id 2 for calling a global function with example usage: multivariate_stats(eigenvals, r_err_sscp, r_contrast, df_resid, tolerance) and returns: results, np