code: def rm_anova2(
      data=None, dv=None, within=None, subject=None, effsize="ng2"
  ):
      """Two-way repeated measures ANOVA.
  
      This is an internal function. The main call to this function should be done
      by the :py:func:`pingouin.rm_anova` function.
      """
      a, b = within
  
      # Validate the dataframe
      data = _check_dataframe(
          dv=dv,
          within=within,
          data=data,
          subject=subject,
          effects="within",
      )
  
      assert not data[a].isnull().any(), (
          "Cannot have missing values in %s" % a
      )
      assert not data[b].isnull().any(), (
          "Cannot have missing values in %s" % b
      )
      assert not data[subject].isnull().any(), (
          "Cannot have missing values in %s" % subject
      )
  
      # Pivot and melt the table. This has several effects:
      # 1) Force missing values to be explicit (a NaN cell is created)
      # 2) Automatic collapsing to the mean if multiple within factors are present
      # 3) If using dropna, remove rows with missing values (listwise deletion).
      # The latter is the same behavior as JASP (= strict complete-case analysis).
      data_piv = data.pivot_table(
          index=subject, columns=within, values=dv, observed=True
      )
      data_piv = data_piv.dropna()
      data = data_piv.melt(
          ignore_index=False, value_name=dv
      ).reset_index()
  
      # Group sizes and grandmean
      n_a = data[a].nunique()
      n_b = data[b].nunique()
      n_s = data[subject].nunique()
      mu = data[dv].mean(numeric_only=True)
  
      # Groupby means
      # I think that observed=True is actually not needed here since we have already used
      # `observed=True` in pivot_table.
      grp_s = data.groupby(subject, observed=True)[dv].mean(
          numeric_only=True
      )
      grp_a = data.groupby([a], observed=True)[dv].mean(
          numeric_only=True
      )
      grp_b = data.groupby([b], observed=True)[dv].mean(
          numeric_only=True
      )
      grp_ab = data.groupby([a, b], observed=True)[dv].mean(
          numeric_only=True
      )
      grp_as = data.groupby([a, subject], observed=True)[dv].mean(
          numeric_only=True
      )
      grp_bs = data.groupby([b, subject], observed=True)[dv].mean(
          numeric_only=True
      )
  
      # Sums of squares
      ss_tot = np.sum((data[dv] - mu) ** 2)
      ss_s = (n_a * n_b) * np.sum((grp_s - mu) ** 2)
      ss_a = (n_b * n_s) * np.sum((grp_a - mu) ** 2)
      ss_b = (n_a * n_s) * np.sum((grp_b - mu) ** 2)
      ss_ab_er = n_s * np.sum((grp_ab - mu) ** 2)
      ss_ab = ss_ab_er - ss_a - ss_b
      ss_as_er = n_b * np.sum((grp_as - mu) ** 2)
      ss_as = ss_as_er - ss_s - ss_a
      ss_bs_er = n_a * np.sum((grp_bs - mu) ** 2)
      ss_bs = ss_bs_er - ss_s - ss_b
      ss_abs = ss_tot - ss_a - ss_b - ss_s - ss_ab - ss_as - ss_bs
  
      # DOF
      df_a = n_a - 1
      df_b = n_b - 1
      df_s = n_s - 1
      df_ab_er = n_a * n_b - 1
      df_ab = df_ab_er - df_a - df_b
      df_as_er = n_a * n_s - 1
      df_as = df_as_er - df_s - df_a
      df_bs_er = n_b * n_s - 1
      df_bs = df_bs_er - df_s - df_b
      df_tot = n_a * n_b * n_s - 1
      df_abs = df_tot - df_a - df_b - df_s - df_ab - df_as - df_bs
  
      # Mean squares
      ms_a = ss_a / df_a
      ms_b = ss_b / df_b
      ms_ab = ss_ab / df_ab
      ms_as = ss_as / df_as
      ms_bs = ss_bs / df_bs
      ms_abs = ss_abs / df_abs
  
      # F-values
      f_a = ms_a / ms_as
      f_b = ms_b / ms_bs
      f_ab = ms_ab / ms_abs
  
      # P-values
      p_a = f(df_a, df_as).sf(f_a)
      p_b = f(df_b, df_bs).sf(f_b)
      p_ab = f(df_ab, df_abs).sf(f_ab)
  
      # Effect sizes
      if effsize == "n2":
          # ..Eta-squared
          ef_a = ss_a / ss_tot
          ef_b = ss_b / ss_tot
          ef_ab = ss_ab / ss_tot
      elif effsize == "ng2":
          # .. Generalized eta-squared (from Bakeman 2005 Table 1) -- default
          ef_a = ss_a / (ss_a + ss_s + ss_as + ss_bs + ss_abs)
          ef_b = ss_b / (ss_b + ss_s + ss_as + ss_bs + ss_abs)
          ef_ab = ss_ab / (ss_ab + ss_s + ss_as + ss_bs + ss_abs)
      else:
          # .. Partial eta squared
          ef_a = (f_a * df_a) / (f_a * df_a + df_as)
          ef_b = (f_b * df_b) / (f_b * df_b + df_bs)
          ef_ab = (f_ab * df_ab) / (f_ab * df_ab + df_abs)
  
      # Epsilon
      piv_a = data.pivot_table(
          index=subject, columns=a, values=dv, observed=True
      )
      piv_b = data.pivot_table(
          index=subject, columns=b, values=dv, observed=True
      )
      eps_a = epsilon(piv_a, correction="gg")
      eps_b = epsilon(piv_b, correction="gg")
      # Note that the GG epsilon of the interaction slightly differs between
      # R and Pingouin. An alternative is to use the lower bound, which is
      # very conservative (same behavior as described on real-statistics.com).
      eps_ab = epsilon(data_piv, correction="gg")
  
      # Greenhouse-Geisser correction
      df_a_c, df_as_c = (
          np.maximum(d * eps_a, 1.0) for d in (df_a, df_as)
      )
      df_b_c, df_bs_c = (
          np.maximum(d * eps_b, 1.0) for d in (df_b, df_bs)
      )
      df_ab_c, df_abs_c = (
          np.maximum(d * eps_ab, 1.0) for d in (df_ab, df_abs)
      )
      p_a_corr = f(df_a_c, df_as_c).sf(f_a)
      p_b_corr = f(df_b_c, df_bs_c).sf(f_b)
      p_ab_corr = f(df_ab_c, df_abs_c).sf(f_ab)
  
      # Create dataframe
      aov = pd.DataFrame(
          {
              "Source": [a, b, a + " * " + b],
              "SS": [ss_a, ss_b, ss_ab],
              "ddof1": [df_a, df_b, df_ab],
              "ddof2": [df_as, df_bs, df_abs],
              "MS": [ms_a, ms_b, ms_ab],
              "F": [f_a, f_b, f_ab],
              "p-unc": [p_a, p_b, p_ab],
              "p-GG-corr": [p_a_corr, p_b_corr, p_ab_corr],
              effsize: [ef_a, ef_b, ef_ab],
              "eps": [eps_a, eps_b, eps_ab],
          }
      )
      return _postprocess_dataframe(aov)
  
file_path: pingouin/pingouin/parametric.py
imports: import warnings
  from collections.abc import Iterable
  import numpy
  import pandas
  from scipy.stats import f
  import pandas_flavor
  from pingouin import _check_dataframe, remove_na, _flatten_list, bayesfactor_ttest, epsilon, sphericity, _postprocess_dataframe
label_desc: Use image node_id 3 for calling a global function with example usage: rm_anova2(data, dv, within, subject, effsize) and returns: _postprocess_dataframe