Multiple Linear Regression in Python (scikit-learn and statsmodels)

Multiple linear regression in Python predicts one number from several inputs. With scikit-learn it takes three lines: pick your feature columns, call LinearRegression().fit(X, y), and read the coefficients:

from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(homes[["sq_ft", "bedrooms", "age"]], homes["price"])
model.coef_, model.intercept_

To make the results checkable, I simulated 400 home sales where I know the true effect of each feature, then fit models with Python 3.12.5, scikit-learn 1.9.1 and statsmodels 0.15.0.

We’ll fit a model and read its coefficients, evaluate it on held-out data, add a categorical city column, get p-values with statsmodels, and plot actual against predicted prices.

How to do multiple linear regression in Python with scikit-learn

The data was built with roughly $150 per square foot, $12,000 per bedroom and minus $1,100 per year of age, plus noise. A good fit should land close to those numbers:

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
n = 400
homes = pd.DataFrame({
    "sq_ft": rng.integers(900, 3500, n),
    "bedrooms": rng.integers(1, 6, n),
    "age": rng.integers(0, 60, n),
    "city": rng.choice(["Austin", "Denver", "Phoenix"], n),
})
city_bonus = homes["city"].map({"Austin": 40_000, "Denver": 65_000, "Phoenix": 0})
homes["price"] = (60_000 + 150 * homes["sq_ft"] + 12_000 * homes["bedrooms"]
                  - 1_100 * homes["age"] + city_bonus + rng.normal(0, 25_000, n)).round(2)

from sklearn.linear_model import LinearRegression

X = homes[["sq_ft", "bedrooms", "age"]]
y = homes["price"]
model = LinearRegression().fit(X, y)

print(f"intercept: ${model.intercept_:,.2f}")
for name, coef in zip(X.columns, model.coef_):
    print(f"{name:<9} {coef:>10,.2f}")
print(f"R squared: {model.score(X, y):.3f}")

new_home = pd.DataFrame({"sq_ft": [2000], "bedrooms": [3], "age": [10]})
print(f"predicted price for 2,000 sq ft, 3 bed, 10 yrs: ${model.predict(new_home)[0]:,.2f}")

Output:

intercept: $85,827.48
sq_ft         148.28
bedrooms   14,707.49
age          -957.38
R squared: 0.905
predicted price for 2,000 sq ft, 3 bed, 10 yrs: $416,927.35

You can read each coefficient as the change in price for one extra unit of that feature, holding the others fixed. That “holding the others fixed” part is what separates multiple regression from simple one-variable regression.

How do you evaluate a multiple linear regression model?

Your model always looks good on the data it learned from, so don’t trust that score. Hold back a quarter of the rows, train on the rest, and score on the unseen part. Here I also add the city, one-hot encoded:

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
n = 400
homes = pd.DataFrame({
    "sq_ft": rng.integers(900, 3500, n),
    "bedrooms": rng.integers(1, 6, n),
    "age": rng.integers(0, 60, n),
    "city": rng.choice(["Austin", "Denver", "Phoenix"], n),
})
city_bonus = homes["city"].map({"Austin": 40_000, "Denver": 65_000, "Phoenix": 0})
homes["price"] = (60_000 + 150 * homes["sq_ft"] + 12_000 * homes["bedrooms"]
                  - 1_100 * homes["age"] + city_bonus + rng.normal(0, 25_000, n)).round(2)

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, r2_score
from sklearn.model_selection import train_test_split

X = pd.get_dummies(homes[["sq_ft", "bedrooms", "age", "city"]], columns=["city"], drop_first=True, dtype=int)
y = homes["price"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=0)

model = LinearRegression().fit(X_train, y_train)
pred = model.predict(X_test)

print("features:", list(X.columns))
print(f"test R squared: {r2_score(y_test, pred):.3f}")
print(f"test MAE:       ${mean_absolute_error(y_test, pred):,.2f}")
for name, coef in zip(X.columns, model.coef_):
    print(f"  {name:<12} {coef:>10,.2f}")

Output:

features: ['sq_ft', 'bedrooms', 'age', 'city_Denver', 'city_Phoenix']
test R squared: 0.953
test MAE:       $19,074.28
  sq_ft            148.23
  bedrooms      13,404.07
  age             -974.00
  city_Denver   23,952.60
  city_Phoenix -42,169.21

pd.get_dummies(..., drop_first=True) turns the city into 0/1 columns and drops one to avoid redundancy, so Denver and Phoenix are measured relative to Austin.

Denver comes out about $24,000 higher and Phoenix about $42,000 lower, close to the true $25,000 and $40,000. For preparing columns, see loc vs iloc in pandas.

Getting p-values with statsmodels

scikit-learn focuses on prediction. If you’d like standard errors and p-values to judge whether each feature matters, statsmodels’ OLS gives a classic regression table:

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
n = 400
homes = pd.DataFrame({
    "sq_ft": rng.integers(900, 3500, n),
    "bedrooms": rng.integers(1, 6, n),
    "age": rng.integers(0, 60, n),
    "city": rng.choice(["Austin", "Denver", "Phoenix"], n),
})
city_bonus = homes["city"].map({"Austin": 40_000, "Denver": 65_000, "Phoenix": 0})
homes["price"] = (60_000 + 150 * homes["sq_ft"] + 12_000 * homes["bedrooms"]
                  - 1_100 * homes["age"] + city_bonus + rng.normal(0, 25_000, n)).round(2)

import statsmodels.api as sm

X = sm.add_constant(homes[["sq_ft", "bedrooms", "age"]])
ols = sm.OLS(homes["price"], X).fit()
table = ols.summary2().tables[1][["Coef.", "Std.Err.", "P>|t|"]]
print(table.round(4).to_string())
print(f"\nR squared {ols.rsquared:.3f} | adjusted {ols.rsquared_adj:.3f}")

Output:

               Coef.   Std.Err.  P>|t|
const     85827.4760  7311.1069    0.0
sq_ft       148.2756     2.4518    0.0
bedrooms  14707.4876  1250.1189    0.0
age        -957.3784   104.0783    0.0

R squared 0.905 | adjusted 0.905

Remember sm.add_constant(); without it, statsmodels fits a line through zero. If you’re coming from R, this output will feel familiar; see R vs Python.

Plotting actual vs predicted values from the regression

Here’s a quick visual check you can do: if the model is good, points hug the diagonal line where predicted equals actual:

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
n = 400
homes = pd.DataFrame({
    "sq_ft": rng.integers(900, 3500, n),
    "bedrooms": rng.integers(1, 6, n),
    "age": rng.integers(0, 60, n),
    "city": rng.choice(["Austin", "Denver", "Phoenix"], n),
})
city_bonus = homes["city"].map({"Austin": 40_000, "Denver": 65_000, "Phoenix": 0})
homes["price"] = (60_000 + 150 * homes["sq_ft"] + 12_000 * homes["bedrooms"]
                  - 1_100 * homes["age"] + city_bonus + rng.normal(0, 25_000, n)).round(2)

import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

X = pd.get_dummies(homes[["sq_ft", "bedrooms", "age", "city"]], columns=["city"], drop_first=True, dtype=int)
model = LinearRegression().fit(X, homes["price"])
pred = model.predict(X)

plt.figure(figsize=(6.5, 5))
plt.scatter(homes["price"] / 1000, pred / 1000, alpha=0.5, s=14)
lims = [homes["price"].min() / 1000, homes["price"].max() / 1000]
plt.plot(lims, lims, color="black", linestyle="--", label="perfect prediction")
plt.xlabel("Actual price ($ thousands)")
plt.ylabel("Predicted price ($ thousands)")
plt.title("Multiple linear regression: actual vs predicted")
plt.legend()
plt.tight_layout()
plt.show()

Multiple linear regression checklist

• Check that relationships look roughly linear; transform skewed features like income if needed.
• Watch for highly correlated features (multicollinearity); they make individual coefficients unstable.
• One-hot encode categories, and drop one level to avoid a redundant column.
• Judge the model on test data, not the training data.
• Scale features only if you add regularization (Ridge or Lasso); plain least squares doesn’t need it.

Continue with these machine learning and data guides:

• Python libraries for machine learning
• Principal component analysis in Python
• Standard deviation in Python
• Join two DataFrames in pandas

Frequently asked questions