How to use the pymer4.utils._welch_ingredients function in pymer4

def welch_dof(x, y):
    """
    Compute adjusted dof via Welch-Satterthwaite equation
    
    Args:
        x (np.ndarray): 1d numpy array
        y (np.ndarray): 1d numpy array

    Returns:
        float: degrees of freedom
    """

    if isinstance(x, np.ndarray) and isinstance(y, np.ndarray):

        x_numerator, x_denominator = _welch_ingredients(x)
        y_numerator, y_denominator = _welch_ingredients(y)

        return np.power(x_numerator + y_numerator, 2) / (x_denominator + y_denominator)
    else:
        raise TypeError("Both x and y must be 1d numpy arrays")

def welch_dof(x, y):
    """
    Compute adjusted dof via Welch-Satterthwaite equation
    
    Args:
        x (np.ndarray): 1d numpy array
        y (np.ndarray): 1d numpy array

    Returns:
        float: degrees of freedom
    """

    if isinstance(x, np.ndarray) and isinstance(y, np.ndarray):

        x_numerator, x_denominator = _welch_ingredients(x)
        y_numerator, y_denominator = _welch_ingredients(y)

        return np.power(x_numerator + y_numerator, 2) / (x_denominator + y_denominator)
    else:
        raise TypeError("Both x and y must be 1d numpy arrays")

if cluster is not None:
            # Cluster corrected dof (num clusters - num coef)
            # Differs from stats and statsmodels which do num cluster - 1
            # Ref: http://cameron.econ.ucdavis.edu/research/Cameron_Miller_JHR_2015_February.pdf
            df = cluster.nunique() - x.shape[1]
        else:
            df = x.shape[0] - x.shape[1]
            if isinstance(weights, str) and wls_dof_correction:
                if weight_groups.ngroups != 2:
                    w = "Welch-Satterthwait DOF correction only supported for 2 groups in the data"
                    warnings.warn(w)
                    self.warnings.append(w)
                else:
                    welch_ingredients = np.array(
                        self.data.groupby(weights)[dv]
                        .apply(_welch_ingredients)
                        .values.tolist()
                    )
                    df = (
                        np.power(welch_ingredients[:, 0].sum(), 2)
                        / welch_ingredients[:, 1].sum()
                    )

        p = 2 * (1 - t_dist.cdf(np.abs(t), df))
        df = np.array([df] * len(t))
        sig = np.array([_sig_stars(elem) for elem in p])

        if conf_int == "boot":

            # Parallelize bootstrap computation for CIs
            par_for = Parallel(n_jobs=n_jobs, backend="multiprocessing")