def anova2(dv=None, between=None, data=None, export_filename=None):
    """Two-way ANOVA.

    This is an internal function. The main call to this function should be done
    by the :py:func:`pingouin.anova` function.

    Parameters
    ----------
    dv : string
        Name of column containing the dependant variable.
    between : list of string
        Name of column containing the two between factors. Must contain exactly
        two values (e.g. ['factor1', 'factor2'])
    data : pandas DataFrame
        DataFrame
    export_filename : string
        Filename (without extension) for the output file.
        If None, do not export the table.
        By default, the file will be created in the current python console
        directory. To change that, specify the filename with full path.

    Returns
    -------
    aov : DataFrame
        ANOVA summary ::

        'Source' : Factor names
        'SS' : Sums of squares
        'DF' : Degrees of freedom
        'MS' : Mean squares
        'F' : F-values
        'p-unc' : uncorrected p-values
        'np2' : Partial eta-square effect sizes
    """
    # Validate the dataframe
    _check_dataframe(dv=dv, between=between, data=data, effects='between')

    fac1, fac2 = between

    # Drop missing values
    data = data[[dv, fac1, fac2]].dropna()

    # Reset index (avoid duplicate axis error)
    data = data.reset_index(drop=True)

    grp_both = data.groupby(between)[dv]
    ng1, ng2 = data[fac1].nunique(), data[fac2].nunique()

    if grp_both.count().nunique() == 1:
        # BALANCED DESIGN
        aov_fac1 = anova(data=data, dv=dv, between=fac1, detailed=True)
        aov_fac2 = anova(data=data, dv=dv, between=fac2, detailed=True)
        # Sums of squares
        ss_fac1 = aov_fac1.loc[0, 'SS']
        ss_fac2 = aov_fac2.loc[0, 'SS']
        ss_tot = ((data[dv] - data[dv].mean())**2).sum()
        ss_resid = np.sum(grp_both.apply(lambda x: (x - x.mean())**2))
        ss_inter = ss_tot - (ss_resid + ss_fac1 + ss_fac2)
        # Degrees of freedom
        df_fac1 = aov_fac1.loc[0, 'DF']
        df_fac2 = aov_fac2.loc[0, 'DF']
        df_inter = (ng1 - 1) * (ng2 - 1)
        df_resid = data[dv].size - (ng1 * ng2)
    else:
        # UNBALANCED DESIGN
        import statsmodels.api as sm
        from statsmodels.formula.api import ols
        warnings.warn("Groups are unbalanced. Type II ANOVA will be computed "
                      "using statsmodels")
        formula = "%s ~ C(%s) * C(%s)" % (dv, fac1, fac2)
        lm = ols(formula, data=data).fit()
        sts = sm.stats.anova_lm(lm, typ=2)
        # Sums of squares
        ss_fac1 = sts.iloc[0, 0]
        ss_fac2 = sts.iloc[1, 0]
        ss_inter = sts.iloc[2, 0]
        ss_resid = sts.iloc[3, 0]
        # Degrees of freedom
        df_fac1 = sts.iloc[0, 1]
        df_fac2 = sts.iloc[1, 1]
        df_inter = sts.iloc[2, 1]
        df_resid = sts.iloc[3, 1]

    # Mean squares
    ms_fac1 = ss_fac1 / df_fac1
    ms_fac2 = ss_fac2 / df_fac2
    ms_inter = ss_inter / df_inter
    ms_resid = ss_resid / df_resid

    # F-values
    fval_fac1 = ms_fac1 / ms_resid
    fval_fac2 = ms_fac2 / ms_resid
    fval_inter = ms_inter / ms_resid

    # P-values
    pval_fac1 = f(df_fac1, df_resid).sf(fval_fac1)
    pval_fac2 = f(df_fac2, df_resid).sf(fval_fac2)
    pval_inter = f(df_inter, df_resid).sf(fval_inter)

    # Partial eta-square
    np2_fac1 = ss_fac1 / (ss_fac1 + ss_resid)
    np2_fac2 = ss_fac2 / (ss_fac2 + ss_resid)
    np2_inter = ss_inter / (ss_inter + ss_resid)

    # Create output dataframe
    aov = pd.DataFrame({'Source': [fac1, fac2, fac1 + ' * ' + fac2,
                                   'residual'],
                        'SS': np.round([ss_fac1, ss_fac2, ss_inter,
                                        ss_resid], 3),
                        'DF': [df_fac1, df_fac2, df_inter, df_resid],
                        'MS': np.round([ms_fac1, ms_fac2, ms_inter,
                                        ms_resid], 3),
                        'F': [fval_fac1, fval_fac2, fval_inter, np.nan],
                        'p-unc': [pval_fac1, pval_fac2, pval_inter, np.nan],
                        'np2': [np2_fac1, np2_fac2, np2_inter, np.nan]
                        })
    col_order = ['Source', 'SS', 'DF', 'MS', 'F', 'p-unc', 'np2']

    aov = aov.reindex(columns=col_order)
    aov.dropna(how='all', axis=1, inplace=True)
    # Export to .csv
    if export_filename is not None:
        _export_table(aov, export_filename)
    return aov