statsmodels cheat sheet

import numpy as np, pandas as pd
import statsmodels.api as sm              # arrays  -> sm.OLS, sm.GLM, sm.Logit
import statsmodels.formula.api as smf     # formulas -> smf.ols, smf.glm, smf.logit
import statsmodels.tsa.api as tsa         # time series

• res = smf.ols("y ~ x1 + x2", data=df).fit()
  Formula API. Lowercase names. Intercept added, NaNs dropped, factors handled.
• res = sm.OLS(y, sm.add_constant(X)).fit()
  Array API. UPPERCASE names. You do all three of those yourself.
• sm.OLS.from_formula('y ~ x', data=df)
  The same thing smf.ols is an alias for.

TipUse smf for anything with categoricals, transforms, or interactions. Use sm when the design matrix already exists (e.g. it came out of a pipeline).

• endog
  The endogenous variable — y. The thing being explained.
• exog
  The exogenous variables — X. The things doing the explaining.
• Model(endog, exog)
  y first. scikit-learn is fit(X, y). This is the reverse.
• res.model.endog / .exog
  The actual arrays that were fitted, post-patsy.
• res.model.exog_names
  Column names of the design matrix — check this when coefficients surprise you.

GotchaPassing (X, y) instead of (y, X) does not raise if the shapes happen to work. It just fits nonsense.

res = smf.ols("y ~ x1 + x2", data=df).fit()

res.params        res.bse         res.tvalues     res.pvalues
res.conf_int(alpha=0.05)
res.rsquared      res.rsquared_adj
res.aic           res.bic         res.llf
res.fvalue        res.f_pvalue
res.resid         res.fittedvalues
res.nobs          res.df_model    res.df_resid

print(res.summary())          # the human-readable table
res.summary2()                # a pandas-friendlier variant

res.params is a pandas Series indexed by term name when you use the formula API. That alone is worth the formula API.

Gotchay ~ x1 + x2**2 does not square x2 — ** is patsy's interaction-expansion operator. You want I(x2**2).

smf.ols("y ~ C(g)", data=df)                             # k-1 dummies, first level = baseline
smf.ols("y ~ C(g, Treatment(reference='b'))", data=df)   # choose the baseline
smf.ols("y ~ C(g, Sum)", data=df)                        # deviation coding (sum-to-zero)
smf.ols("y ~ C(g, Poly)", data=df)                       # ordered factors

smf.ols("y ~ bs(x, df=4) + cr(x, df=4)", data=df)        # B-spline / cyclic spline basis
smf.ols("y ~ center(x1) + standardize(x2)", data=df)     # patsy built-ins

• string / Categorical columns
  Auto-treated as categorical. C() is only needed to force it or set coding.
• integer-coded groups
  Treated as numeric unless you wrap them in C(). This is the classic silent bug.

GotchaA group column coded 1, 2, 3 will be fitted as a straight line through the group numbers unless you write C(group).

new = pd.DataFrame({"x1": [0.5], "g": ["a"]})
res.predict(new)                       # formula API: patsy re-applies the SAME transforms

# array API: you must rebuild the design matrix yourself, constant and all
res.predict(sm.add_constant(Xnew, has_constant='add'))

TipThis is the single biggest practical argument for the formula API: the transforms, the dummy coding and the baseline level all travel with the model.

GotchaArray API: if Xnew has only one row, add_constant silently does nothing unless you pass has_constant='add'.

• smf.ols("y ~ x1 + x2", data=df).fit()
  Ordinary least squares. Homoskedastic, independent errors.
• smf.wls("y ~ x", data=df, weights=w).fit()
  Weighted — known, unequal error variance. w ∝ 1/σ².
• sm.GLS(y, X, sigma=S).fit()
  Generalized — you supply the full error covariance.
• sm.GLSAR(y, X, rho=1).iterative_fit(3)
  AR(p) errors (Cochrane–Orcutt).
• smf.rlm("y ~ x", data=df, M=sm.robust.norms.HuberT()).fit()
  Robust to outliers in y. Also TukeyBiweight(), RamsayE().
• smf.quantreg("y ~ x", data=df).fit(q=0.9)
  Models a quantile, not the mean. Great for tails.
• RollingOLS(y, X, window=60).fit()
  Rolling-window betas. from statsmodels.regression.rolling import RollingOLS.
• sm.RecursiveLS(y, X).fit()
  Expanding window + CUSUM tests for parameter stability.

Heteroskedasticity of unknown form? Don't reach for WLS — just use OLS(...).fit(cov_type='HC3').

res = smf.ols("y ~ x1 + x2 + x3", data=df).fit_regularized(alpha=0.1, L1_wt=1.0)
#   L1_wt = 1.0 -> lasso     0.0 -> ridge     between -> elastic net
res.params

Gotchafit_regularized returns coefficients but no standard errors and no p-values — the usual inference is not valid after selection. If you want prediction + CV, that's scikit-learn's job.

print(res.summary())                  # the classic three-panel table
res.summary().tables[1]               # just the coefficient block
res.summary(alpha=0.10)               # 90% intervals instead of 95%

res.summary().as_latex()              # paper
res.summary().as_csv()                # spreadsheet

from statsmodels.iolib.summary2 import summary_col
summary_col([res_a, res_b, res_c], stars=True,
            info_dict={"N": lambda r: f"{int(r.nobs)}", "R2": lambda r: f"{r.rsquared:.2f}"})

Tipsummary_col is the one you want for a paper: several models side by side in one column-per-spec table, with significance stars.

• Durbin–Watson
  ≈2 is fine. <1 or >3 means autocorrelated residuals → your SEs are wrong.
• Jarque–Bera / Prob(JB)
  Residual normality. With large n this matters far less than people think.
• Cond. No. > 30
  Multicollinearity or just unscaled columns. Check VIF before you panic.
• <code>Omnibus</code>
  Another normality test. Agrees with JB most of the time.

A footnote about a strong condition number often just means one column is in units of millions. Standardize and re-look before you start deleting predictors.

res = smf.ols("y ~ x", data=df).fit(cov_type="HC3")

res = smf.ols("y ~ x", data=df).fit(cov_type="HAC",
                                    cov_kwds={"maxlags": 4})
res = smf.ols("y ~ x", data=df).fit(cov_type="cluster",
                                    cov_kwds={"groups": df.firm})
res2 = res.get_robustcov_results(cov_type="HC1")   # re-do SEs without re-fitting

GotchaPoint estimates are identical across all of these — only bse, tvalues, pvalues and conf_int() move. If a result flips from significant to not under HC3, it was never robust.

res.t_test("x1 = 2")                       # is the slope 2, not 0?
res.f_test("x1 = x2 = 0")                  # joint significance
res.f_test("x1 - x2 = 0")                  # equality of two coefficients
res.wald_test_terms()                      # one test per model term

res_small = smf.ols("y ~ x1", data=df).fit()
res_big   = smf.ols("y ~ x1 + x2 + x3", data=df).fit()
res_big.compare_f_test(res_small)          # nested models
res_big.compare_lr_test(res_small)         # likelihood ratio

Restrictions can also be given as matrices — strings are just the readable front-end.

p = res.get_prediction(new_df)
p.summary_frame(alpha=0.05)
#   mean  mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper

Gotchamean_ci is the interval for the average response — narrow. obs_ci is the interval for a new observation — much wider, and usually the one your stakeholder actually meant.

from statsmodels.stats.diagnostic import (het_breuschpagan, het_white,
     acorr_breusch_godfrey, acorr_ljungbox, linear_reset, linear_rainbow)
from statsmodels.stats.stattools import durbin_watson, jarque_bera

het_breuschpagan(res.resid, res.model.exog)    # -> lm, lm_p, f, f_p
het_white(res.resid, res.model.exog)           # also catches wrong functional form
durbin_watson(res.resid)                       # ~2 = clean
acorr_breusch_godfrey(res, nlags=4)            # higher-order autocorrelation
jarque_bera(res.resid)                         # normality
linear_reset(res, power=2, use_f=True)         # is the model linear at all?

Each returns (statistic, p-value, ...). Small p on a heteroskedasticity test is not a crisis — it is an instruction to pass cov_type='HC3'.

inf = res.get_influence()
inf.summary_frame()                    # everything, per observation

inf.cooks_distance[0]                  # Cook's D  -> rule of thumb: > 4/n
inf.hat_matrix_diag                    # leverage  -> > 2k/n is high
inf.resid_studentized_external         # |t| > 3 is a candidate outlier
inf.dffits[0]

res.outlier_test()                     # Bonferroni-corrected p per point

• high leverage, small residual
  Unusual x, but the model fits it. Mostly harmless.
• low leverage, big residual
  An outlier in y. Consider RLM.
• high leverage <b>and</b> big residual
  Influential. This single point is steering your slope.

GotchaFinding an influential point is not a licence to delete it. Report the fit with and without, or switch to RLM and let the estimator down-weight it.

from statsmodels.stats.outliers_influence import variance_inflation_factor as vif

X = sm.add_constant(df[["x1", "x2", "x3"]])          # the constant MUST be in here
[vif(X.values, i) for i in range(X.shape[1])]
#   VIF > 5  -> worth a look     VIF > 10 -> serious

GotchaRun VIF on a matrix without the constant and you get garbage. Measured on two predictors correlated at r = 0.57: with the constant, VIF = 1.5 (correct). Without it, VIF = 3480. The constant's own VIF is huge and meaningless — ignore index 0.

TipCollinearity inflates standard errors; it does not bias the coefficients. If you only care about prediction, you can often just leave it alone.

res = smf.glm("y ~ x1 + C(g)", data=df,
              family=sm.families.Binomial()).fit()            # logistic regression

res = smf.glm("count ~ x1", data=df,
              family=sm.families.Poisson(),
              exposure=df.person_years).fit()                 # a RATE model

res = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.Log())).fit()

• family=
  Gaussian, Binomial, Poisson, NegativeBinomial, Gamma, InverseGaussian, Tweedie
• link=
  Logit, Probit, CLogLog, Log, Identity, Power, Sqrt — optional; each family has a canonical default.
• exposure= / offset=
  exposure is logged for you; offset is not. Both fix the denominator of a rate.
• res.deviance / pearson_chi2
  Fit quality. pearson_chi2 / df_resid > 1 → over-dispersion.

GotchaPoisson assumes mean = variance. Real counts are almost always over-dispersed, which makes Poisson standard errors far too small. Check the dispersion, then move to NegativeBinomial.

res.params                     # on the LINK scale (log-odds, log-rate)
np.exp(res.params)             # odds ratios (Binomial) / rate ratios (Poisson)
res.fittedvalues               # on the RESPONSE scale (probabilities, means)
res.predict(new)               # response scale
res.predict(new, linear=True)  # link scale
res.resid_deviance             # the residual to plot for a GLM

Gotchaparams lives on the link scale but fittedvalues lives on the response scale. Verified: for a Poisson fit, fittedvalues[:3] = [5.16, 1.11, 3.00] while predict(linear=True) = [1.64, 0.10, 1.10]. Same rows, different scale.

res = smf.logit("y ~ x1 + C(g)", data=df).fit()      # log-odds
res = smf.probit("y ~ x1", data=df).fit()

np.exp(res.params)              # ODDS RATIOS - what you actually report
res.predict(new)                # PROBABILITIES, not 0/1 labels
res.pred_table()                # confusion matrix at the 0.5 cut
res.prsquared                   # McFadden pseudo-R2 (NOT comparable to OLS R2)
res.llr_pvalue                  # likelihood-ratio test vs the null model

sm.MNLogit(y_codes, X).fit()    # 3+ unordered categories

Gotcha.predict() returns probabilities. Verified: [0.998, 0.001, 0.789, 0.172]. Thresholding at 0.5 is your decision, and 0.5 is rarely the right one for imbalanced data.

me = res.get_margeff(at="overall", method="dydx")
me.summary()
#   'the average marginal effect of x1 on P(y=1) is +0.13'

• at=
  'overall' (average marginal effect — usually what you want) · 'mean' · 'median'
• method=
  'dydx' · 'eyex' (elasticity) · 'dyex' · 'eydx'

TipLogit coefficients are on the log-odds scale and nobody has intuition for them. Marginal effects put the answer back in probability points — report these to non-statisticians.

• smf.poisson("cnt ~ x", data=df).fit()
  Baseline count model. Assumes mean = variance.
• smf.negativebinomial("cnt ~ x", data=df).fit()
  Over-dispersed counts. The realistic default.
• sm.NegativeBinomialP(y, X, p=2).fit()
  NB1 (p=1) vs NB2 (p=2) variance form.
• sm.GeneralizedPoisson(y, X).fit()
  Handles under- and over-dispersion.
• sm.ZeroInflatedPoisson(y, X, exog_infl=Z).fit()
  Two processes: one generates the zeros, one the counts.
• sm.ZeroInflatedNegativeBinomialP(y, X).fit()
  Both problems at once.

GotchaOfficial pitfall: fitting NegativeBinomial to data that isn't over-dispersed pushes the dispersion parameter to the boundary at zero, where the likelihood is not evaluable — the optimizer stalls or returns NaN. Same for zero-inflated models on data with no excess zeros.

from statsmodels.miscmodels.ordinal_model import OrderedModel

y = pd.Categorical(df.rating, categories=[1, 2, 3, 4, 5], ordered=True)
res = OrderedModel(y, df[["x1", "x2"]], distr="logit").fit(method="bfgs")

GotchaOrderedModel takes no intercept — the cutpoints play that role. Do not add a constant, and do not use a formula with one.

res = smf.mixedlm("y ~ x1 + x2", data=df, groups=df["school"]).fit()

res = smf.mixedlm("y ~ x1", data=df, groups=df["school"],
                  re_formula="~x1").fit()          # random intercept AND random slope

res.fe_params        # fixed effects
res.random_effects   # per-group deviations (a dict)
res.cov_re           # random-effect covariance

• random intercept
  Groups differ in level. re_formula omitted.
• random slope
  Groups differ in how they respond to x. re_formula="~x1".
• <code>vc_formula=</code>
  Crossed / nested variance components.

GotchaA ConvergenceWarning from MixedLM is common and must not be ignored — often it means the random-effect variance is being estimated at ~0, i.e. you don't need the random effect. Try method=['lbfgs'], rescale your predictors, or simplify the random structure.

res = smf.gee("y ~ x1", groups="subject", data=df,
              cov_struct=sm.cov_struct.Exchangeable(),
              family=sm.families.Binomial()).fit()

• cov_struct=
  Independence() · Exchangeable() · Autoregressive() · Nested()
• MixedLM vs GEE
  MixedLM = subject-specific effects (“for a given school…”). GEE = population-average (“across all schools…”). They answer different questions and their coefficients genuinely differ for nonlinear links.

TipGEE's standard errors are robust to getting cov_struct wrong. That is the whole point of it.

from statsmodels.stats.anova import anova_lm, AnovaRM

anova_lm(smf.ols("y ~ C(g) + x1", data=df).fit(), typ=2)
anova_lm(smf.ols("y ~ C(a)*C(b)", data=df).fit(), typ=3)   # with interactions -> typ=3
anova_lm(res_small, res_big)                                # nested model comparison

AnovaRM(df, "score", "subject", within=["condition"]).fit()  # repeated measures

Gotchatyp=1 (sequential) depends on the order you wrote the terms. For unbalanced designs use typ=2, and typ=3 when interactions are present — and with type 3 you need sum-to-zero coding (C(g, Sum)) for the results to mean anything.

from statsmodels.tsa.stattools import adfuller, kpss

adfuller(y, autolag="AIC")     # H0: unit root (NON-stationary).  p < .05 -> stationary
kpss(y, regression="c")        # H0: STATIONARY.                  p < .05 -> not stationary

GotchaThe two tests have opposite null hypotheses. Reading ADF's p-value as if it were KPSS's is the most common time-series error in Python. Run both.

from statsmodels.tsa.stattools import acf, pacf, ccf, arma_order_select_ic
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

acf(y, nlags=24, alpha=0.05)              # values + confidence band
pacf(y, nlags=24, method="ywm")
plot_acf(y, lags=24); plot_pacf(y, lags=24)

arma_order_select_ic(y, max_ar=4, max_ma=4, ic=["aic", "bic"])   # grid search
ccf(x, y)                                 # cross-correlation, for lead/lag

plot_acf includes lag 0 (always 1.0). That first giant spike is not a finding.

from statsmodels.tsa.seasonal import seasonal_decompose, STL, MSTL

d = seasonal_decompose(y, model="additive", period=12)   # classic moving-average
d.trend; d.seasonal; d.resid

STL(y, period=12, robust=True).fit()      # better: handles changing seasonality + outliers
MSTL(y, periods=(24, 168)).fit()          # MULTIPLE seasonalities (hourly: daily + weekly)

from statsmodels.tsa.filters.hp_filter import hpfilter
cycle, trend = hpfilter(y, lamb=1600)     # quarterly=1600, monthly=129600, annual=6.25

Gotchaseasonal_decompose needs period= unless your index is a pandas DatetimeIndex with a set freq. And prefer STL: the classic decomposition uses a centred moving average, so it loses data at both ends.

import statsmodels.tsa.api as tsa

res = tsa.ARIMA(y, order=(1, 1, 1), trend="t").fit()

res = tsa.SARIMAX(y, exog=X,
                  order=(1, 1, 1),
                  seasonal_order=(1, 1, 1, 12),          # (P, D, Q, s)
                  enforce_stationarity=True).fit(disp=0)

res.summary(); res.aic; res.bic
res.plot_diagnostics(figsize=(10, 6))    # residual ACF + QQ + histogram in one call

• order=(p, d, q)
  AR lags, differencing, MA lags.
• seasonal_order=(P, D, Q, s)
  The same, at the seasonal period s (12 = monthly, 4 = quarterly).
• exog=
  External regressors. You must supply them for the forecast horizon too.
• trend=
  'n' none · 'c' constant · 't' linear · 'ct'

Gotchasm.tsa.ARIMA today is the state-space implementation. The old statsmodels.tsa.arima_model.ARIMA / ARMA classes were removed in 0.13 — any tutorial calling .fit(disp=-1) on them predates that.

• tsa.AutoReg(y, lags=3, seasonal=True, period=12)
  Pure AR. Fast, and select_order helps.
• tsa.ExponentialSmoothing(y, trend='add', seasonal='add', seasonal_periods=12)
  Holt–Winters. Strong baseline; hard to beat.
• ETSModel(y, error='add', trend='add', seasonal='add')
  Same family, but with likelihood, AIC and proper prediction intervals.
• sm.tsa.UnobservedComponents(y, level='local linear trend', seasonal=12)
  Structural: gives you back an interpretable trend and seasonal.
• tsa.VAR(df).fit(maxlags=4, ic='aic')
  Several series that feed back into each other.
• sm.tsa.MarkovRegression(y, k_regimes=2, switching_variance=True)
  Regime switching (recessions, volatility states).
• sm.tsa.DynamicFactor(df, k_factors=1, factor_order=1)
  Many series driven by a few latent factors.

model = tsa.VAR(df[["gdp", "cpi", "rate"]])
model.select_order(8).summary()                 # AIC / BIC / HQIC per lag
res = model.fit(maxlags=4, ic="aic")

res.irf(10).plot(orth=True)                     # impulse responses: shock one, watch the others
res.fevd(10).summary()                          # variance decomposition: who explains whom
res.test_causality("gdp", ["rate"], kind="f")   # Granger causality
res.test_whiteness(nlags=12)                    # residuals clean?
res.forecast(df.values[-4:], steps=8)

GotchaGranger causality is predictive precedence, not causation. And if the series are cointegrated, a VAR in levels is misspecified — test with coint_johansen and use VECM instead.

f = res.get_forecast(steps=12, exog=X_future)
f.predicted_mean
f.conf_int(alpha=0.05)
f.summary_frame()

res.forecast(12)                      # point forecasts only
res.predict(start=..., end=..., dynamic=False)   # IN-sample fit
res.predict(start=..., end=..., dynamic=True)    # recursive - the honest one

res2 = res.append(new_obs, refit=False)          # new data, keep the parameters

Gotchadynamic=False feeds the true lagged values back in at every step, so your “backtest” is really a one-step-ahead fit and will look far better than reality. Use dynamic=True, or hold out properly.

from statsmodels.stats.diagnostic import acorr_ljungbox

acorr_ljungbox(res.resid, lags=[10], return_df=True)   # H0: residuals are white noise
res.plot_diagnostics()

Residual ACF flat, Ljung–Box p large, QQ plot straight → you have squeezed out the structure. Anything left is noise, and no amount of extra lags will help.

from statsmodels.stats.multitest import multipletests

reject, p_adj, _, _ = multipletests(pvals, alpha=0.05, method="fdr_bh")
#   'bonferroni' 'holm' 'fdr_bh' (Benjamini-Hochberg) 'fdr_by'

TipTesting 50 hypotheses at α=0.05 gives you ~2.5 false positives by construction. fdr_bh is the pragmatic default; Bonferroni is the conservative one.

from statsmodels.stats.power import TTestIndPower, tt_ind_solve_power

TTestIndPower().solve_power(effect_size=0.5, power=0.8, alpha=0.05)   # -> n per group
tt_ind_solve_power(effect_size=0.5, nobs1=64, alpha=0.05, power=None) # -> achieved power

Leave exactly one argument as None and it solves for that one. Also FTestAnovaPower, NormalIndPower, GofChisquarePower.

• proportions_ztest([30, 40], [100, 100])
  Two-sample test of proportions.
• proportion_confint(30, 100, method='wilson')
  CI for one proportion. Wilson, not the normal approximation.
• samplesize_confint_proportion(0.3, 0.05)
  How many do I need for ±5 points?
• Table2x2(tbl).oddsratio / .riskratio
  Epidemiology two-by-two, with CIs and .summary().
• mcnemar(tbl)
  Paired binary (before/after on the same subjects).
• StratifiedTable(tables).test_null_odds()
  Mantel–Haenszel across strata.

Gotchaproportion_confint defaults to method='normal', which produces intervals outside [0,1] for small n or extreme p. Pass 'wilson'.

from statsmodels.stats.weightstats import DescrStatsW, CompareMeans, ztest

d = DescrStatsW(x, weights=w)
d.mean, d.std, d.tconfint_mean()          # WEIGHTED, with a CI

cm = CompareMeans.from_data(a, b)
cm.tconfint_diff(usevar="unequal")        # Welch CI for the difference
cm.summary()

from statsmodels.stats.multicomp import pairwise_tukeyhsd
pairwise_tukeyhsd(df.value, df.group)     # all pairs, family-wise corrected

scipy.stats gives you the p-value; statsmodels gives you the confidence interval on the difference, which is the thing worth reporting.

• sm.nonparametric.lowess(y, x, frac=0.3)
  LOESS smoother — the scatter-plot trend line.
• nonparametric.KDEUnivariate(x).fit()
  Kernel density.
• multivariate.PCA(X, standardize=True)
  PCA with scree/loadings and missing-value handling.
• MANOVA.from_formula('y1 + y2 ~ g', data=df).mv_test()
  Multivariate ANOVA.
• duration.survfunc.SurvfuncRight(time, status)
  Kaplan–Meier + survdiff log-rank.
• duration.hazard_regression.PHReg(t, X, status=e).fit()
  Cox proportional hazards.
• imputation.mice.MICEData(df)
  Multiple imputation by chained equations.
• stats.mediation.Mediation(outcome, mediator, 'x', 'm').fit()
  Causal mediation, with bootstrap CIs.

• Rank-deficient / collinear <code>exog</code>
  Models use a pseudo-inverse, so perfect collinearity does not error. Verified: duplicating a column split the true slope 2.045 into 1.022 + 1.022, rank 2 of 3 columns, no warning.
• Ill-conditioned design
  Numerically unstable results. Check res.condition_number and rescale.
• MLE “convergence”
  The criterion can fire on the objective while the parameters have not converged. Verify.
• Quasi-separation in Logit/Probit
  Perfect separation raises. Quasi-perfect separation does not — you get huge coefficients and huge SEs.
• RLM with a perfect fit
  Residuals all zero → scale = 0 → 0/0 → NaNs, no exception.
• NegBin on non-overdispersed data
  Dispersion is pushed to the boundary at 0, where the log-likelihood cannot be evaluated.
• Zero-inflated with no excess zeros
  No inflation sits on the boundary; deflation is outside the parameter space. Optimization fails.

Source: the official User Guide → Background → Pitfalls page.

# WRONG - res1 and res2 both point at the same mutated model object
mod  = sm.RLM(y, X)
res1 = mod.fit(scale_est="mad")
res2 = mod.fit(scale_est=sm.robust.scale.HuberScale())

# RIGHT - a fresh model per fit configuration
res1 = sm.RLM(y, X).fit(scale_est="mad")
res2 = sm.RLM(y, X).fit(scale_est=sm.robust.scale.HuberScale())

GotchaResults objects keep a reference back to the model. Re-fitting the same model instance with different arguments can quietly invalidate results you are still holding. Looping to collect a scalar (e.g. AIC per lag) is fine; keeping two live results is not.

• <code>sm.OLS(y, X)</code> with no constant
  No intercept is fitted. Measured: slope 1.765 vs the true 2.0, R² 0.114 vs 0.832. No error, no warning.
• R² with no intercept
  Silently becomes the uncentered R². Not comparable to any other R².
• <code>Model(X, y)</code> argument order
  Reversed vs scikit-learn. Fits nonsense if the shapes align.
• Integer-coded groups without <code>C()</code>
  Fitted as a numeric slope through the group labels.
• <code>x**2</code> in a formula
  That is patsy interaction syntax, not arithmetic. Use I(x**2).
• <code>logit.predict()</code>
  Returns probabilities, not classes. Measured: [0.998, 0.001, 0.789, 0.172].
• VIF without the constant column
  Measured on r = 0.57 predictors: VIF 3480 without the constant vs 1.5 with it.
• GLM <code>params</code> vs <code>fittedvalues</code>
  Link scale vs response scale. Do not mix them in one plot.
• <code>predict(dynamic=False)</code>
  Uses true lagged values — flatters your backtest badly.
• ADF vs KPSS null hypotheses
  They are opposite. Reading one as the other inverts your conclusion.

• Formula API
  Drops rows with NaN by default. Measured: 100 rows in, nobs = 90 out, silently.
• Array API
  Raises MissingDataError. Pass missing='drop' explicitly.
• Either way
  Check res.nobs against len(df). Every single time.

GotchaA model fitted on 90 of 100 rows is not necessarily wrong — but if you did not know it happened, you do not know what you fitted.