Getting started
Install
pip install checkmodel
Check a model
Fit a model as usual and hand it to check_model():
import pandas as pd
import statsmodels.formula.api as smf
from checkmodel import check_model
mtcars = pd.read_csv(
"https://raw.githubusercontent.com/pablo-legeren/checkmodel/main/tests/data/mtcars.csv",
index_col=0,
)
fit = smf.ols("mpg ~ wt + hp", data=mtcars).fit()
report = check_model(fit)
print(report)
# Model checks: OLS (gaussian, identity link), n = 32
Posterior predictive check: OK: Model-predicted data resemble the observed data (mean p = 0.400, sd p = 0.460).
Normality of residuals: Warning: Non-normality of residuals detected (p = 0.021).
Homogeneity of variance: OK: Error variance appears to be homoscedastic (p = 0.402).
Collinearity: OK: Low collinearity (all VIF < 5; max 1.77).
Influential observations: OK: No outliers detected (Cook's distance > 0.808).
Autocorrelation is not in the list: the cars in mtcars have no natural
order, and a Durbin-Watson test on unordered rows only reacts to how they
happen to be sorted. check_model() adds it when the data are indexed by time
(a pandas DatetimeIndex or PeriodIndex). Request it explicitly otherwise:
check_model(fit, checks=["normality", "autocorrelation"])
report.plot() draws all the diagnostics in one figure. Each panel says what
a good model looks like in its subtitle.
Formal tests and plots
A small p-value on 32 observations and a clean-looking Q-Q plot can disagree. Read the plots first: the tests are a prompt to look, not a verdict.
Read individual results
The report behaves like a dictionary of check results:
normality = report["normality"]
normality.passed # False
normality.p_value # 0.0214
normality.method # 'Shapiro-Wilk'
normality.plot() # Q-Q plot and histogram
report.summary() # one row per check, as a DataFrame
report.skipped # checks that do not apply, with the reason
Every check is also available on its own:
from checkmodel import check_collinearity, check_outliers
check_collinearity(fit).to_dataframe() # VIF table, one row per term
check_outliers(fit, threshold="4/n") # use the 4/n rule instead of R's default
A model matrix with exactly collinear columns, such as a dummy for every category plus an intercept, fails the collinearity check whatever its VIFs:
Warning: Perfect collinearity: the model matrix is rank deficient (aliased terms: cyl_8).
result.rank_deficient and result.aliased hold the same information; the
aliased columns are the ones R reports as NA coefficients.
scikit-learn models
scikit-learn estimators do not store their training data, so pass it:
from sklearn.linear_model import LinearRegression
X, y = mtcars[["wt", "hp"]], mtcars["mpg"]
model = LinearRegression().fit(X, y)
check_model(model, X=X, y=y)
The numbers are the same as for the statsmodels fit above. Penalized models
(Ridge, Lasso, ElasticNet, penalized LogisticRegression,
PoissonRegressor and GammaRegressor with alpha > 0) are supported too.
They have no likelihood-based criteria, so AIC and BIC are left out of
their metric tables. Their leverage, Cook's distance and standardized
residuals are computed as if the model were unpenalized, so the checks that
use them come with a warning, and the report says so in its header:
# Model checks: Ridge (gaussian, identity link), n = 32
Note: `Ridge` is a penalized fit; leverage, Cook's distance and standardized residuals are computed as for the unpenalized model and are approximate.
Logistic regression
import statsmodels.api as sm
logit = smf.glm("vs ~ mpg + wt", data=mtcars, family=sm.families.Binomial()).fit()
check_model(logit)
Normality and homoscedasticity are not expected for a binary outcome, so they
are not run. Binned residuals take their place. If every fitted probability is
0 or 1 (complete separation), the coefficients are not identified and
check_model() warns that the diagnostics are unreliable.
Count data
Poisson regressions are checked for overdispersion; zero-inflation has its own check:
from checkmodel import check_zeroinflation
bio = pd.read_csv(
"https://raw.githubusercontent.com/pablo-legeren/checkmodel/main/tests/data/biochemists.csv"
)
pois = smf.glm(
"art ~ fem + mar + kid5 + phd + ment", data=bio, family=sm.families.Poisson()
).fit()
check_model(pois)
check_zeroinflation(pois)
# Model checks: GLM (poisson, log link), n = 915
Posterior predictive check: Warning: Observed data fall outside the model-predicted range for: sd (p < .001), proportion of 0s (p < .001).
Collinearity: OK: Low collinearity (all VIF < 5; max 1.29).
Influential observations: OK: No outliers detected (Cook's distance > 0.84).
Overdispersion: Warning: Overdispersion detected (dispersion ratio = 1.829, p < .001).
Warning: Model is underfitting zeros (probable zero-inflation): 275 observed, 191 predicted, ratio 0.69.
Mixed models
statsmodels MixedLM results get the residual checks plus checks on the random
effects, singularity and convergence:
from checkmodel import icc
sleep = pd.read_csv(
"https://raw.githubusercontent.com/pablo-legeren/checkmodel/main/tests/data/sleepstudy.csv"
)
mixed = smf.mixedlm("Reaction ~ Days", sleep, groups="Subject", re_formula="~Days").fit()
check_model(mixed)
# Model checks: MixedLM (gaussian, identity link), n = 180
Posterior predictive check: OK: Model-predicted data resemble the observed data (mean p = 0.460, sd p = 0.520).
Normality of residuals: Warning: Non-normality of residuals detected (p < .001).
Normality of random effects: OK: Random effects 'Subject: (Intercept)' appear as normally distributed (p = 0.428).
OK: Random effects 'Subject: Days' appear as normally distributed (p = 0.898).
Homogeneity of variance: Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
Influential observations: Warning: 4 outliers detected: cases 7, 56, 59, 130 (Cook's distance > 0.696).
Singularity: OK: The random-effects structure is not singular.
Convergence: OK: Model converged.
Skipped:
collinearity: Not enough model terms to check for multicollinearity (need at least two).
icc() and r2() decompose the variance as R's performance does:
icc(mixed) # Adjusted ICC: 0.722, Unadjusted ICC: 0.521
r2(mixed) # Conditional R2: 0.799, Marginal R2: 0.279
Bayesian models
Install the extra (Python 3.12+), fit with PyMC or bambi, and wrap the result
with bayesian():
pip install "checkmodel[bayes]"
import bambi as bmb
from checkmodel import bayesian, check_model, model_performance, r2_bayes
model = bmb.Model("Reaction ~ Days + (1|Subject)", sleep)
fit = bayesian(model.fit(), model)
check_model(fit) # posterior predictive check and convergence (R-hat, ESS)
r2_bayes(fit) # conditional and marginal Bayesian R² with HDI
model_performance(fit) # ELPD, LOOIC, WAIC, R2, RMSE, Sigma
With bambi, bayesian() computes the expected values, posterior predictive
draws and pointwise log-likelihood itself. With PyMC, add a
pm.Deterministic("mu", ...) for the expected value of each observation and
pass the model so the rest can be computed:
import pymc as pm
with pm.Model() as model:
a = pm.Normal("a", 20, 10)
b = pm.Normal("b", 0, 10)
sigma = pm.HalfNormal("sigma", 5)
mu = pm.Deterministic("mu", a + b * mtcars["wt"].to_numpy())
pm.Normal("mpg", mu, sigma, observed=mtcars["mpg"].to_numpy())
idata = pm.sample()
fit = bayesian(idata, model)
Singular fits
If a random-effect variance belongs at zero, statsmodels' default
optimizer can stop just short of it. check_convergence() reports this;
refit with .fit(method="powell") and check_singularity() will then
flag the singular fit, as lme4 does.
Performance metrics
from checkmodel import model_performance, compare_performance, r2
r2(fit)
model_performance(fit)
AIC AICc BIC R2 R2_adjusted RMSE Sigma
0 156.652 158.134 162.515 0.827 0.815 2.469 2.593
Compare several models; keyword names label the rows:
small = smf.ols("mpg ~ wt", data=mtcars).fit()
large = smf.ols("mpg ~ wt + hp + qsec", data=mtcars).fit()
compare_performance(small=small, medium=fit, large=large, rank=True)
With rank=True the table gains a Performance_Score (the mean of every
index rescaled to 0–1, higher is better) and is sorted by it.