How do I run a linear regression in Python and read the output?

The short version

Use statsmodels with a formula like "sales ~ temp" and call .fit(). Read the coefficient (the change in the outcome for a one-unit rise in the predictor), its p-value (evidence against the coefficient being zero) and its confidence interval (the plausible range). Prefer statsmodels over scikit-learn when you want to interpret the model rather than predict with it.

Fitting a model

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

df = pd.DataFrame({
    "temp": [10, 15, 20, 25, 30],
    "sales": [20, 25, 33, 38, 46],
})

model = smf.ols("sales ~ temp", data=df).fit()
print(model.summary())

The formula reads "sales explained by temp". An intercept is included automatically.

Reading the key numbers

model.params["temp"]        # 1.30 , the slope
model.params["Intercept"]   # 6.40
model.rsquared              # 0.9937
model.pvalues["temp"]       # 0.000215
model.conf_int().loc["temp"]  # [1.11, 1.49]

The slope of 1.30 means that each extra degree is associated with about 1.3 more sales, on average, in this data. The intercept is the predicted value when temp is zero, which is often meaningless in itself (there may be no such thing as zero temperature in your context), so do not over-interpret it.

What the p-value does and does not say

A small p-value (here 0.0002) means the data would be unlikely if the true coefficient were exactly zero. It does not tell you the effect is large, important, or causal. A tiny effect can have a tiny p-value in a big sample. Always look at the size of the coefficient and its confidence interval, which say how big the effect plausibly is, alongside the p-value, which only speaks to whether it is distinguishable from zero.

R squared is not a verdict

rsquared is the share of variation in the outcome explained by the model. A high value does not mean the model is correct: a badly misspecified model can still fit closely, and a high R squared on a curved relationship fitted with a straight line is misleading. Always plot the data and the residuals.

Adding a categorical predictor

Wrap a categorical column in C() so it is encoded as indicator variables rather than treated as a number.

smf.ols("sales ~ temp + C(region)", data=df).fit()

Without the C(), a numeric-looking category code would be treated as a quantity, which is a silent and serious error.

The assumptions that matter

Ordinary least squares assumes a roughly linear relationship, residuals with constant spread, and independent observations. Check by plotting residuals against fitted values: a curve means the relationship is not linear, and a fan shape means the spread is not constant. Time-series data usually violates independence, which needs a different model.

Correlation, still not causation

A coefficient describes an association in your data. It becomes a causal claim only with a design that supports one (an experiment, or careful adjustment for confounders). A regression by itself cannot tell you that changing the predictor would change the outcome, and presenting it as if it can is the most common misuse of the technique.