Online Statsmodels Compiler

statsmodels fits statistical models and reports them the way a statistics course expects: coefficients with standard errors, p-values and confidence intervals. It covers linear and logistic regression, ANOVA, time series models and many statistical tests, and it accepts R-style formulas such as "score ~ hours". Researchers and students use it to explain a relationship, not only to predict it. This page is an online statsmodels compiler: the code runs in your browser, so you can try it without installing anything. Run the example first, then paste any snippet below into a new cell to try it.

What the example does

The example fits a straight line to the exam scores of 20 students. sm.add_constant() adds a column of ones to the hours, which gives the model an intercept. sm.OLS() sets up an ordinary least squares model and .fit() estimates it: the line is 56.94 + 2.59 × hours. results.predict([1, hours]) takes the same two columns, a 1 for the constant and then the hours, and returns an array: 82.85 for 10 hours. results.summary() prints the full table. R-squared is 0.922, and the x1 row (the hours) holds the slope with its standard error, p-value and 95% confidence interval (2.22 to 2.97).

Fit a model with a formula

statsmodels.formula.api takes a pandas DataFrame and a formula, and adds the intercept for you. The scores rise quickly at first and slowly later, so this adds a squared term with I(hours ** 2):

import pandas as pd
import statsmodels.formula.api as smf

df = pd.DataFrame({
    "hours": range(1, 21),
    "score": [51, 55, 60, 68, 72, 75, 78, 82, 85, 88,
              90, 92, 94, 96, 97, 98, 99, 100, 101, 102],
})

linear = smf.ols("score ~ hours", data=df).fit()
curved = smf.ols("score ~ hours + I(hours ** 2)", data=df).fit()

print(curved.params.round(3))
print("R-squared:", round(linear.rsquared, 3), "->", round(curved.rsquared, 3))
print("AIC:", round(linear.aic, 1), "->", round(curved.aic, 1))

R-squared rises from 0.922 to 0.997, and AIC (lower is better) falls from 119.6 to 59.0. For a text column, C(column) fits one coefficient per category, measured against the first.

Get predictions with intervals

get_prediction() with summary_frame() gives each prediction with two 95% intervals:

import pandas as pd
import statsmodels.formula.api as smf

df = pd.DataFrame({
    "hours": range(1, 21),
    "score": [51, 55, 60, 68, 72, 75, 78, 82, 85, 88,
              90, 92, 94, 96, 97, 98, 99, 100, 101, 102],
})
model = smf.ols("score ~ hours + I(hours ** 2)", data=df).fit()

new = pd.DataFrame({"hours": [5, 10, 15]})
prediction = model.get_prediction(new).summary_frame(alpha=0.05)
print(prediction.round(1))

mean_ci bounds the average score of all students who study that long. obs_ci bounds the score of one new student, so it is wider.

Run a logistic regression

For a yes-or-no outcome, smf.logit() models the probability of a 1. Here 200 simulated students pass more often the longer they study:

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(seed=1)
hours = rng.uniform(0, 10, 200)
chance = 1 / (1 + np.exp(-(hours - 5)))          # true pass probability
df = pd.DataFrame({"hours": hours, "passed": rng.binomial(1, chance)})

model = smf.logit("passed ~ hours", data=df).fit(disp=0)
print(model.params.round(3))
print("Odds ratio per hour:", np.exp(model.params["hours"]).round(2))
print(model.predict(pd.DataFrame({"hours": [2, 5, 8]})).round(2))

The coefficients are log-odds, so np.exp() turns the slope into an odds ratio: each extra hour multiplies the odds of passing by about 3. predict() returns probabilities, and disp=0 hides the optimizer's progress message.

Forecast a time series with ARIMA

ARIMA takes a series with a date index and an order of (p, d, q): autoregressive terms, differences and moving-average terms. trend="t" adds a drift:

import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA

rng = np.random.default_rng(seed=3)
months = pd.date_range("2023-01-01", periods=36, freq="MS")
steps = rng.normal(3, 8, 36)             # random steps, 3 up on average
sales = pd.Series(200 + np.cumsum(steps), index=months, name="sales")

model = ARIMA(sales, order=(1, 1, 0), trend="t").fit()
forecast = model.get_forecast(steps=6)
table = forecast.conf_int()
table["forecast"] = forecast.predicted_mean
print(table.round(1))

The forecast grows by about 3 a month, the drift in the data. The 95% interval from conf_int() widens the further ahead it looks.

Good to know

• sm.OLS fits no intercept unless you call sm.add_constant(). Without it, the example's line is forced through zero and the slope jumps from 2.59 to 6.76. Formulas add the Intercept themselves.
• sm.OLS(y, X) takes the outcome first, the reverse of scikit-learn's fit(X, y).