def ancova(dv=None, covar=None, between=None, data=None,
           export_filename=None, return_bw=False):
    """ANCOVA with one or more covariate(s).

    Parameters
    ----------
    dv : string
        Name of column containing the dependant variable.
    covar : string or list
        Name(s) of column(s) containing the covariate.
    between : string
        Name of column containing the between factor.
    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.
    return_bw : bool
        If True, return beta within parameter (used for the rm_corr function)

    Returns
    -------
    aov : DataFrame
        ANCOVA summary ::

        'Source' : Names of the factor considered
        'SS' : Sums of squares
        'DF' : Degrees of freedom
        'F' : F-values
        'p-unc' : Uncorrected p-values

    Notes
    -----
    Analysis of covariance (ANCOVA) is a general linear model which blends
    ANOVA and regression. ANCOVA evaluates whether the means of a dependent
    variable (dv) are equal across levels of a categorical independent
    variable (between) often called a treatment, while statistically
    controlling for the effects of other continuous variables that are not
    of primary interest, known as covariates or nuisance variables (covar).

    Note that in the case of one covariate, Pingouin will use a built-in
    function. However, if there are more than one covariate, Pingouin will
    use the statsmodels package to compute the ANCOVA.

    Rows with missing values are automatically removed (listwise deletion).

    See Also
    --------
    anova : One-way and two-way ANOVA

    Examples
    --------
    1. Evaluate the reading scores of students with different teaching method
    and family income as a covariate.

    >>> from pingouin import ancova, read_dataset
    >>> df = read_dataset('ancova')
    >>> ancova(data=df, dv='Scores', covar='Income', between='Method')
         Source           SS  DF          F     p-unc
    0    Method   571.030045   3   3.336482  0.031940
    1    Income  1678.352687   1  29.419438  0.000006
    2  Residual  1768.522365  31        NaN       NaN

    2. Evaluate the reading scores of students with different teaching method
    and family income + BMI as a covariate.

    >>> ancova(data=df, dv='Scores', covar=['Income', 'BMI'], between='Method')
         Source        SS  DF       F     p-unc
    0    Method   552.284   3   3.233  0.036113
    1    Income  1573.952   1  27.637  0.000011
    2       BMI    60.014   1   1.054  0.312842
    3  Residual  1708.509  30     NaN       NaN
    """
    # Safety checks
    assert isinstance(data, pd.DataFrame)
    assert dv in data, '%s is not in data.' % dv
    assert between in data, '%s is not in data.' % between
    assert isinstance(covar, (str, list)), 'covar must be a str or a list.'

    # Drop missing values
    data = data[_flatten_list([dv, between, covar])].dropna()

    # Check the number of covariates
    if isinstance(covar, list):
        if len(covar) > 1:
            return ancovan(dv=dv, covar=covar, between=between, data=data,
                           export_filename=export_filename)
        else:
            covar = covar[0]

    # Assert that covariate is numeric
    assert data[covar].dtype.kind in 'fi'

    def linreg(x, y):
        return np.corrcoef(x, y)[0, 1] * np.std(y, ddof=1) / np.std(x, ddof=1)

    # Compute slopes
    groups = data[between].unique()
    slopes = np.zeros(shape=groups.shape)
    ss = np.zeros(shape=groups.shape)
    for i, b in enumerate(groups):
        dt_covar = data[data[between] == b][covar].values
        ss[i] = ((dt_covar - dt_covar.mean())**2).sum()
        slopes[i] = linreg(dt_covar, data[data[between] == b][dv].values)
    ss_slopes = ss * slopes
    bw = ss_slopes.sum() / ss.sum()
    bt = linreg(data[covar], data[dv])

    # Run the ANOVA
    aov_dv = anova(data=data, dv=dv, between=between, detailed=True)
    aov_covar = anova(data=data, dv=covar, between=between, detailed=True)

    # Create full ANCOVA
    ss_t_dv = aov_dv.loc[0, 'SS'] + aov_dv.loc[1, 'SS']
    ss_t_covar = aov_covar.loc[0, 'SS'] + aov_covar.loc[1, 'SS']
    # Sums of squares
    ss_t = ss_t_dv - bt**2 * ss_t_covar
    ss_w = aov_dv.loc[1, 'SS'] - bw**2 * aov_covar.loc[1, 'SS']
    ss_b = ss_t - ss_w
    ss_c = ss_slopes.sum() * bw
    # DOF
    df_c = 1
    df_b = aov_dv.loc[0, 'DF']
    df_w = aov_dv.loc[1, 'DF'] - 1
    # Mean squares
    ms_c = ss_c / df_c
    ms_b = ss_b / df_b
    ms_w = ss_w / df_w
    # F-values
    f_c = ms_c / ms_w
    f_b = ms_b / ms_w
    # P-values
    p_c = f(df_c, df_w).sf(f_c)
    p_b = f(df_b, df_w).sf(f_b)

    # Create dataframe
    aov = pd.DataFrame({'Source': [between, covar, 'Residual'],
                        'SS': [ss_b, ss_c, ss_w],
                        'DF': [df_b, df_c, df_w],
                        'F': [f_b, f_c, np.nan],
                        'p-unc': [p_b, p_c, np.nan],
                        })

    # Export to .csv
    if export_filename is not None:
        _export_table(aov, export_filename)

    if return_bw:
        return aov, bw
    else:
        return aov
ANCOVA with one or more covariate(s).

    Parameters
    ----------
    dv : string
        Name of column containing the dependant variable.
    covar : string or list
        Name(s) of column(s) containing the covariate.
    between : string
        Name of column containing the between factor.
    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.
    return_bw : bool
        If True, return beta within parameter (used for the rm_corr function)

    Returns
    -------
    aov : DataFrame
        ANCOVA summary ::

        'Source' : Names of the factor considered
        'SS' : Sums of squares
        'DF' : Degrees of freedom
        'F' : F-values
        'p-unc' : Uncorrected p-values

    Notes
    -----
    Analysis of covariance (ANCOVA) is a general linear model which blends
    ANOVA and regression. ANCOVA evaluates whether the means of a dependent
    variable (dv) are equal across levels of a categorical independent
    variable (between) often called a treatment, while statistically
    controlling for the effects of other continuous variables that are not
    of primary interest, known as covariates or nuisance variables (covar).

    Note that in the case of one covariate, Pingouin will use a built-in
    function. However, if there are more than one covariate, Pingouin will
    use the statsmodels package to compute the ANCOVA.

    Rows with missing values are automatically removed (listwise deletion).

    See Also
    --------
    anova : One-way and two-way ANOVA

    Examples
    --------
    1. Evaluate the reading scores of students with different teaching method
    and family income as a covariate.

    >>> from pingouin import ancova, read_dataset
    >>> df = read_dataset('ancova')
    >>> ancova(data=df, dv='Scores', covar='Income', between='Method')
         Source           SS  DF          F     p-unc
    0    Method   571.030045   3   3.336482  0.031940
    1    Income  1678.352687   1  29.419438  0.000006
    2  Residual  1768.522365  31        NaN       NaN

    2. Evaluate the reading scores of students with different teaching method
    and family income + BMI as a covariate.

    >>> ancova(data=df, dv='Scores', covar=['Income', 'BMI'], between='Method')
         Source        SS  DF       F     p-unc
    0    Method   552.284   3   3.233  0.036113
    1    Income  1573.952   1  27.637  0.000011
    2       BMI    60.014   1   1.054  0.312842
    3  Residual  1708.509  30     NaN       NaN