Variance inflation factors

[poroc300]

I'd like to propose adding a function to calculate Variance Inflation Factors (VIFs) for assessing multicollinearity among predictor variables. While statsmodels provides this functionality, the current implementation requires unnecessary boilerplate (in my opinion) to extract VIF values across all predictors:

from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_ratios = []
for idx in range(df_num.shape[1]):
   vif_ratios.append(variance_inflation_factor(arr, idx))

In my view, the main limitations of this approach are:

1. Column name mapping - since most data analysis workflows use pandas DataFrames, additional code is needed to associate VIF values with their corresponding column names.
2. Missing data - the function does not handle missing values natively and raises a MissingDataError if any are present.
3. Categorical columns - categorical variables are not handled natively and must be excluded manually before calling the function.

I propose a wrapper function with the following behaviour:

1. Accepts both pandas DataFrames and numpy arrays as input.
2. Automatically removes rows containing missing values.
3. Discards categorical columns when the input is a pandas DataFrame.
4. Returns a pandas DataFrame with VIF values per column (when the input is a DataFrame), or a numpy array of VIF values in the original column order (when the input is an array).
5. Optionally, accepts a threshold parameter to flag column names or indices where the VIF exceeds a user-defined value.

In my private repository, I've implemented previously a function that does some of the above (the below is just for demonstration, it would require to be expanded in terms of functionality):

def vif(df):
    """
    Compute Variance Inflation Factor (VIF) for each numeric variable. Low values
    indicate low multicollinearity.

    Parameters
    ----------
    df : DataFrame
        Input DataFrame containing the variables to assess. Non-numeric columns
        are ignored automatically and rows with any NaN values and constante values
        are dropped before computation.

    Returns
    -------
    DataFrame
        Single-row DataFrame with the VIF for each numerical variable.

    """
    # get numeric columns, remove columns with constant and nan values
    df_num = df.select_dtypes("number").copy()
    df_num = df_num.dropna(axis=0, how="any")
    df_num = df_num.loc[:, df_num.nunique() > 1]
    if df_num.empty:
        raise ValueError("preprocessing made the input DataFrame empty")

    # estimate VIF for each variable
    vif_ratios = []
    for idx in range(df_num.shape[1]):
        vif_ratios.append(variance_inflation_factor(arr, idx))

    # convert vif_ratios to DataFrame
    return pd.DataFrame([vif_ratios], columns=df_num.columns, index=["VIF"])

[raphaelvallat]

Hi @poroc300 — thanks for opening the issue.

I agree about adding a vif function in Pingouin and am aligned with your proposed behaviors.

I was brainstorming this with Claude which suggested that we can bypass the statsmodels import and speed up the implementation by using the inverse correlation matrix, similar to thepcorr (partial correlation) function:

def vif(data, threshold=None, remove_na=False):
    """Variance Inflation Factor for each numeric predictor.

    VIF_i = 1 / (1 - R_i^2), where R_i^2 comes from regressing predictor i on
    all the other predictors. Equivalently, VIF_i is the i-th diagonal
    element of the inverse of the correlation matrix of the predictors.
    """
    num = data.select_dtypes("number").copy()
    if remove_na:
        num = num.dropna()
    else:
        assert not num.isna().any().any(), (
            "Data contains missing values. Use remove_na=True."
        )
    num = num.loc[:, num.nunique() > 1]  # drop constant columns
    if num.shape[1] < 2:
        raise ValueError("At least two numeric, non-constant columns are required.")

    z = (num - num.mean()) / num.std()
    corr = z.cov()  # correlation matrix, since z is standardized
    vinv = np.linalg.pinv(corr.to_numpy(), hermitian=True)
    out = pd.DataFrame({"names": corr.columns, "VIF": np.diag(vinv)})
    if threshold is not None:
        out["high_vif"] = out["VIF"] > threshold
    return out

Notes from Claude:

• The math is verified — it matches statsmodels.stats.outliers_influence.variance_inflation_factor to 3 decimal places, and it's O(1) matrix inversion instead of looping an OLS fit per column.
• This is a rough draft, not production code. Before it could land in regression.py it would need: pingouin's usual input handling (accept a plain array too, not just a DataFrame — right now select_dtypes assumes a DataFrame), the @pf.register_dataframe_method decorator to match pcorr/partial_corr, _postprocess_dataframe instead of a bare pd.DataFrame, a full numpydoc docstring with Examples/Notes/See also, and a docstring reference to the 1/(1-R²) ↔ inverse-correlation-matrix identity so it's not a mystery one-liner.
• No test coverage yet — would need a TestVIF class in tests/test_regression.py cross-checking against statsmodels (with add_constant) the way TestRegression already cross-checks against sklearn/statsmodels.

[poroc300]

Hi @raphaelvallat ,

Thank you for considering this. Would you be open to me giving it a try, following the guidelines in the contributors documentation? I've never contributed to an open source project before, but I do have experience using git. Let me know your thoughts :)