Formula API and Model Selection

The R-style formula API, then model selection and comparison: information criteria, nested-model tests, and cross-validation caveats for statistical models.

Formula API (R-style)

Statsmodels supports R-style formulas for intuitive model specification:

import statsmodels.formula.api as smf

# OLS with formula
results = smf.ols('y ~ x1 + x2 + x1:x2', data=df).fit()

# Categorical variables (automatic dummy coding)
results = smf.ols('y ~ x1 + C(category)', data=df).fit()

# Interactions
results = smf.ols('y ~ x1 * x2', data=df).fit()  # x1 + x2 + x1:x2

# Polynomial terms
results = smf.ols('y ~ x + I(x**2)', data=df).fit()

# Logit
results = smf.logit('y ~ x1 + x2 + C(group)', data=df).fit()

# Poisson
results = smf.poisson('count ~ x1 + x2', data=df).fit()

# ARIMA (not available via formula, use regular API)

Model Selection and Comparison

Use identical response definitions and observations with comparable likelihoods. Formula missing-value handling can silently change the sample: select/record common rows first. Compare fixed effects in mixed models under ML, not REML. GLM likelihood BIC is bic_llf; deviance BIC is a different convention. Regular LR chi-square inference excludes boundary cases such as Poisson vs NB alpha=0, separation and unidentified parameters.

Do not interpret a p-value after model/variable selection as though the model were prespecified.

Information Criteria

# Compare models using AIC/BIC
models = {
    'Model 1': model1_results,
    'Model 2': model2_results,
    'Model 3': model3_results
}

comparison = pd.DataFrame({
    'AIC': {name: res.aic for name, res in models.items()},
    'BIC': {name: res.bic for name, res in models.items()},
    'Log-Likelihood': {name: res.llf for name, res in models.items()}
})

print(comparison.sort_values('AIC'))
# Lower values rank comparable candidates; they do not establish absolute adequacy

Likelihood Ratio Test (Nested Models)

# For nested models (one is subset of the other)
from scipy import stats

lr_stat = 2 * (full_model.llf - reduced_model.llf)
df_diff = full_model.df_model - reduced_model.df_model
if df_diff <= 0:
    raise ValueError("Expected a larger full model")
p_value = stats.chi2.sf(lr_stat, df_diff)

print(f"LR statistic: {lr_stat:.4f}")
print(f"p-value: {p_value:.4f}")

if p_value < 0.05:
    print("Reject the specified restrictions under regular LR assumptions")
else:
    print("No rejection of the specified restrictions; this is not equivalence")

Cross-Validation

The example assumes independent observations and a numeric pandas design with an intercept. Use grouped or time-ordered splits for dependent data and fit preprocessing within each training fold. Scikit-learn LogisticRegression is regularized by default and is not an interchangeable statsmodels MLE wrapper.

from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error

kf = KFold(n_splits=5, shuffle=True, random_state=42)
cv_scores = []

for train_idx, val_idx in kf.split(X):
    X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
    y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]

    # Fit model
    model = sm.OLS(y_train, X_train).fit()

    # Predict
    y_pred = model.predict(X_val)

    # Score
    rmse = np.sqrt(mean_squared_error(y_val, y_pred))
    cv_scores.append(rmse)

print(f"CV RMSE: {np.mean(cv_scores):.4f} ± {np.std(cv_scores):.4f}")