The Toolkit

Three features where the third is almost a mix of the other two. No pair is correlated above 0.73, so a correlation heatmap looks unremarkable. Yet the variance inflation factors are 48, 50 and 103 — far past the usual worry level of 10. Multicollinearity can involve several features at once, which pairwise correlations cannot show.

• VIF rules of thumb (5, 10) are conventions, not tests; what matters is whether the coefficients you interpret are stable (O’Brien 2007).
• Collinearity hurts the interpretation of coefficients far more than predictions.
• Dropping one of two correlated features changes the meaning of the other’s coefficient.

rng = np.random.default_rng(6)
n = 1_000
a = rng.normal(size=n); b = rng.normal(size=n)
c = 0.7 * a + 0.7 * b + rng.normal(0, 0.1, n)      # almost a combination of a and b
X = np.c_[a, b, c]

def vif(X, j):
    others = np.c_[np.delete(X, j, axis=1), np.ones(len(X))]
    fit = others @ np.linalg.lstsq(others, X[:, j], rcond=None)[0]
    r2 = 1 - ((X[:, j] - fit) ** 2).sum() / ((X[:, j] - X[:, j].mean()) ** 2).sum()
    return 1 / (1 - r2)

vifs = [round(vif(X, j), 1) for j in range(3)]
max_pairwise = np.abs(np.corrcoef(X.T)[np.triu_indices(3, 1)]).max()

Assumes import numpy as np and import pandas as pd (and SciPy where imported). Each snippet simulates or makes up its own data, as the comments say.

# Python — correlation & VIF
import pandas as pd
from statsmodels.stats.outliers_influence import variance_inflation_factor

# Correlation heatmap
corr = df.corr()
sns.heatmap(corr, annot=True, cmap='RdBu_r', center=0)

# VIF for each feature
vif = pd.DataFrame()
vif['Feature'] = X.columns
vif['VIF'] = [variance_inflation_factor(X.values, i)
            for i in range(X.shape[1])]