Predicting Ames House Prices: an end-to-end case study¶
This notebook runs a real regression project from a raw, messy file to a validated, interpreted model, following a 12-step method. The data is deliberately dirty (duplicate rows, inconsistent categories, three date formats, structural and genuine missing values, a skewed target, and outliers), exactly what you meet in practice. We lean on pandas, seaborn/matplotlib, statsmodels, and scikit-learn, and explain every table and picture along the way.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns # seaborn = high-level statistical plots (heatmaps, regplots, pairplots)
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.nonparametric.smoothers_lowess import lowess
from sklearn.linear_model import LinearRegression, Ridge, Lasso, LogisticRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
ORG="#059669"; DEEP="#047857"; LIGHT="#6ee7b7"; EM="#059669"; EMDEEP="#047857"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
import warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm; import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import jarque_bera, durbin_watson
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.model_selection import cross_val_score; from sklearn.linear_model import LinearRegression
pd.set_option('display.max_columns', 40)
Objective. Predict a home's sale price from its characteristics, and identify which characteristics drive price the most.
Why / for whom. A real-estate agent pricing a new listing, a homeowner deciding which renovation pays off, a buyer checking whether an asking price is fair, or a county assessor valuing properties at scale.
Success looks like: a model that (a) predicts price with useful accuracy, (b) has trustworthy coefficients we can read as dollar drivers, and (c) passes its regression diagnostics or is honestly corrected where it does not. Getting there is steps 2-12.
Real projects pull data from many places. pandas has a one-line reader for almost all of them:
| Source | Typical format | Python (pandas) |
|---|---|---|
| Flat file export | CSV / TSV | pd.read_csv('file.csv') |
| Spreadsheet | Excel .xlsx | pd.read_excel('file.xlsx', sheet_name=...) |
| Relational database | SQL (Postgres, MySQL, SQLite) | pd.read_sql('SELECT ...', engine) with SQLAlchemy |
| Web API | JSON | requests.get(url).json() then pd.json_normalize(...) |
| Web page (tables) | HTML | pd.read_html(url) returns a list of DataFrames |
| Web scraping (general) | HTML | requests + BeautifulSoup to parse, then build a DataFrame |
| Columnar / big data | Parquet | pd.read_parquet('file.parquet') (fast, typed) |
Whatever the source, the goal is the same: land the data in a DataFrame. From there the rest of the pipeline is identical. Our source is a CSV exported from the county records, so:
# The one-liners for other sources (commented; we only have the CSV here):
# import sqlalchemy; eng = sqlalchemy.create_engine('sqlite:///homes.db'); df = pd.read_sql('SELECT * FROM sales', eng)
# import requests; df = pd.json_normalize(requests.get(API_URL).json()['results'])
# tables = pd.read_html('https://example.com/prices'); df = tables[0]
try:
homes = pd.read_csv('../../data/ames_housing.csv')
except FileNotFoundError:
homes = pd.read_csv(BASE + 'ames_housing.csv')
print('loaded', homes.shape)
loaded (610, 16)
print('rows x cols:', homes.shape)
homes.head(4)
rows x cols: (610, 16)
| house_id | neighborhood | lot_area | lot_frontage | year_built | gr_liv_area | bedrooms | full_bath | central_air | fireplace_qual | garage_yr_built | pool_qc | kitchen_qual | exterior_color | sale_date | sale_price | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1220 | Somerset | 6339 | 58.0 | 1955 | 885 | 1 | 1 | Y | NaN | 1953.0 | NaN | Po | brick | 09/07/1955 | 146251.0 |
| 1 | 1169 | OldTown | 9227 | 62.0 | 1940 | 1516 | 4 | 2 | Yes | TA | 1945.0 | NaN | TA | blue | 1940-03-19 | 99200.0 |
| 2 | 1049 | OldTown | 12262 | 36.0 | 1950 | 1784 | 4 | 2 | Y | Gd | 1955.0 | NaN | Po | beige | 19-Dec-1950 | 135344.0 |
| 3 | 1270 | BrookSide | 7592 | 48.0 | 1945 | 1509 | 4 | 3 | Y | NaN | 1950.0 | NaN | Ex | blue | 1945-02-25 | 180715.0 |
homes.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 610 entries, 0 to 609 Data columns (total 16 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 house_id 610 non-null int64 1 neighborhood 610 non-null object 2 lot_area 610 non-null int64 3 lot_frontage 512 non-null float64 4 year_built 610 non-null int64 5 gr_liv_area 610 non-null int64 6 bedrooms 610 non-null int64 7 full_bath 610 non-null int64 8 central_air 610 non-null object 9 fireplace_qual 320 non-null object 10 garage_yr_built 581 non-null float64 11 pool_qc 6 non-null object 12 kitchen_qual 610 non-null object 13 exterior_color 610 non-null object 14 sale_date 610 non-null object 15 sale_price 607 non-null float64 dtypes: float64(3), int64(6), object(7) memory usage: 76.4+ KB
# Numeric summary, and the missing-value map
display(homes.describe().round(0))
miss = homes.isna().sum().sort_values(ascending=False)
print('\nMissing values per column:'); print(miss[miss>0].to_string())
| house_id | lot_area | lot_frontage | year_built | gr_liv_area | bedrooms | full_bath | garage_yr_built | sale_price | |
|---|---|---|---|---|---|---|---|---|---|
| count | 610.0 | 610.0 | 512.0 | 610.0 | 610.0 | 610.0 | 610.0 | 581.0 | 607.0 |
| mean | 1298.0 | 9687.0 | 71.0 | 1953.0 | 1508.0 | 3.0 | 2.0 | 1955.0 | 184147.0 |
| std | 174.0 | 3238.0 | 21.0 | 31.0 | 501.0 | 1.0 | 1.0 | 31.0 | 64748.0 |
| min | 1000.0 | 3246.0 | 21.0 | 1900.0 | 480.0 | 1.0 | 1.0 | 1900.0 | 68872.0 |
| 25% | 1145.0 | 7466.0 | 58.0 | 1928.0 | 1202.0 | 2.0 | 1.0 | 1929.0 | 138462.0 |
| 50% | 1298.0 | 9108.0 | 70.0 | 1952.0 | 1482.0 | 3.0 | 2.0 | 1953.0 | 172879.0 |
| 75% | 1449.0 | 11362.0 | 85.0 | 1979.0 | 1775.0 | 4.0 | 2.0 | 1981.0 | 216536.0 |
| max | 1599.0 | 31300.0 | 144.0 | 2009.0 | 5642.0 | 6.0 | 4.0 | 2010.0 | 518293.0 |
Missing values per column: pool_qc 604 fireplace_qual 290 lot_frontage 98 garage_yr_built 29 sale_price 3
Three kinds of trouble are already visible. (1) Missing values in five columns, but of two very different kinds: pool_qc (604) and fireplace_qual (290) are missing because most homes have no pool or fireplace (structural), whereas lot_frontage (98) is genuinely unknown and sale_price (3) is our target, which we cannot impute. (2) Inconsistent categories, look:
print('central_air :', homes.central_air.value_counts().to_dict())
print('kitchen_qual:', homes.kitchen_qual.value_counts().to_dict())
print('ext color :', homes.exterior_color.value_counts().to_dict())
print('duplicate house_id rows:', homes.house_id.duplicated().sum())
print('sale_date formats:', homes.sale_date.dropna().head(4).tolist())
central_air : {'Y': 550, 'N': 52, 'Yes': 8}
kitchen_qual: {'Po': 112, 'TA': 111, 'Ex': 106, 'Gd': 106, 'Fa': 101, 'fa': 13, 'FA': 9, 'gd': 8, 'EX': 8, 'Ta': 7, 'ex': 7, 'GD': 7, 'po': 6, 'ta': 5, 'PO': 4}
ext color : {'white': 155, 'beige': 147, 'gray': 119, 'brick': 80, 'blue': 49, 'Beige': 13, 'WHITE': 11, 'White': 9, 'Gray': 7, 'Blue': 6, 'Brick': 6, 'grey': 5, 'Grey': 3}
duplicate house_id rows: 10
sale_date formats: ['09/07/1955', '1940-03-19', '19-Dec-1950', '1945-02-25']
(3) Data-entry chaos. central_air mixes Y and Yes; kitchen_qual and exterior_color come in every mix of upper/lower case (Po, po, PO); there are 10 duplicate house IDs; and sale_date arrives in three different formats. None of this is unusual, and all of it must be fixed before modeling.
df = homes.copy()
n0 = len(df)
# (a) the target cannot be imputed -> drop rows with no sale_price
df = df.dropna(subset=['sale_price'])
# (b) remove duplicate listings (same house_id)
df = df.drop_duplicates(subset='house_id')
print(f'rows: {n0} -> {len(df)} (dropped {n0-len(df)}: 3 missing target + 10 duplicates)')
rows: 610 -> 597 (dropped 13: 3 missing target + 10 duplicates)
# (c) standardize inconsistent categories
df['central_air'] = df['central_air'].replace({'Yes':'Y'}).map({'Y':1, 'N':0})
df['exterior_color'] = df['exterior_color'].str.strip().str.lower().replace({'grey':'gray'})
up = lambda s: s.astype(str).str.strip().str.upper() # 'po','PO','Po' -> 'PO'
QMAP = {'PO':1, 'FA':2, 'TA':3, 'GD':4, 'EX':5} # ordinal quality scale
df['kitchen_qual'] = up(df['kitchen_qual']).map(QMAP)
print('kitchen_qual now:', sorted(df.kitchen_qual.dropna().unique()), ' ext colors:', sorted(df.exterior_color.unique()))
kitchen_qual now: [np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(5)] ext colors: ['beige', 'blue', 'brick', 'gray', 'white']
# (d) structural missing -> a 'none' category and a has_* flag (NOT an error, it is information)
df['fireplace_ord'] = up(df['fireplace_qual'].fillna('NONE')).map({**QMAP, 'NONE':0})
df['has_fireplace'] = df['fireplace_qual'].notna().astype(int)
df['has_pool'] = df['pool_qc'].notna().astype(int)
df['has_garage'] = df['garage_yr_built'].notna().astype(int)
# (e) genuine missing -> impute lot_frontage with the neighborhood median (smarter than a global median)
before = df['lot_frontage'].isna().sum()
df['lot_frontage'] = df.groupby('neighborhood')['lot_frontage'].transform(lambda s: s.fillna(s.median()))
print(f'lot_frontage: imputed {before} values by neighborhood median')
lot_frontage: imputed 95 values by neighborhood median
# (f) parse the three date formats in one call, then engineer AND validate house age
dt = pd.to_datetime(df['sale_date'], format='mixed', errors='coerce')
df['sale_year'] = dt.dt.year.fillna(dt.dt.year.median())
df['house_age'] = (df['sale_year'] - df['year_built']).clip(lower=0)
# Always sanity-check an engineered feature before trusting it in a model:
print('unparseable dates:', dt.isna().sum())
print('house_age distinct values:', df['house_age'].nunique(),
'| variance:', round(float(df['house_age'].var()), 3))
print('-> in this extract the recorded sale year mirrors the build year, so house_age is constant and')
print(' carries no signal; we keep the parsed sale_year but drop the dead house_age from the model.')
unparseable dates: 0 house_age distinct values: 1 | variance: 0.0 -> in this extract the recorded sale year mirrors the build year, so house_age is constant and carries no signal; we keep the parsed sale_year but drop the dead house_age from the model.
# Confirm the clean-up: missing before vs after
fig, ax = plt.subplots(1, 2, figsize=(11,3.6))
homes.isna().sum().sort_values().plot(kind='barh', ax=ax[0], color=RED); ax[0].set_title('BEFORE: missing per column')
df[[c for c in df.columns if c in homes.columns]].isna().sum().sort_values().plot(kind='barh', ax=ax[1], color=GREEN)
ax[1].set_title('AFTER: missing per column'); plt.tight_layout(); plt.show()
fig, ax = plt.subplots(1, 2, figsize=(12,4.2))
ax[0].hist(df.sale_price, bins=40, color=EM, alpha=0.8)
ax[0].set(title=f'Sale price is right-skewed (skew = {df.sale_price.skew():.2f})', xlabel='sale price', ylabel='homes')
ax[0].xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
order = df.groupby('neighborhood').sale_price.median().sort_values().index
sns.boxplot(data=df, y='neighborhood', x='sale_price', order=order, ax=ax[1], color=LIGHT)
ax[1].set(title='Price varies hugely by neighborhood', xlabel='sale price', ylabel='')
ax[1].xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
plt.tight_layout(); plt.show()
The histogram is right-skewed (a long tail of expensive homes, skew 1.16), a red flag for the equal-variance and normality assumptions, which we will fix with a log transform in step 6. The box plots show location dominates: median price roughly doubles from the cheapest neighborhood to the priciest, and the whiskers/dots flag some high-price outliers to watch.
# Correlation heatmap of the numeric drivers
numcols = ['sale_price','gr_liv_area','lot_area','lot_frontage','bedrooms','full_bath','kitchen_qual']
fig, ax = plt.subplots(figsize=(7.6,6))
sns.heatmap(df[numcols].corr(), annot=True, fmt='.2f', cmap='RdYlGn', center=0, ax=ax, square=True, cbar_kws=dict(shrink=0.7))
ax.set_title('Correlation heatmap: what moves with price?'); plt.tight_layout(); plt.show()
Reading the top row (correlations with price): kitchen_qual (0.55), gr_liv_area (living area, 0.47), and bedrooms (0.47) are the strongest numeric drivers. lot_area and lot_frontage barely move with price on their own. The heatmap also lets us eyeball multicollinearity between predictors, none of the off-diagonal cells is alarming here (we confirm with VIF in step 9).
Why transform? A linear model assumes the target's errors are symmetric and equal-variance. A right-skewed price breaks that. Modeling log(price) turns multiplicative, percentage-style effects (which is how housing actually works) into additive ones, and stabilizes the variance. Watch the skew collapse:
df['log_price'] = np.log(df['sale_price'])
fig, ax = plt.subplots(1, 2, figsize=(11,3.8))
ax[0].hist(df.sale_price, bins=40, color=RED, alpha=0.7); ax[0].set_title(f'BEFORE: price (skew {df.sale_price.skew():.2f})')
ax[0].xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
ax[1].hist(df.log_price, bins=40, color=GREEN, alpha=0.75); ax[1].set_title(f'AFTER: log(price) (skew {df.log_price.skew():.2f})')
plt.tight_layout(); plt.show()
The log transform pulled the skew from 1.16 to about 0.02, essentially symmetric. Encoding turns categories into numbers a model can use: the quality ratings are ordinal (Po < Fa < TA < Gd < Ex), so we already mapped them to 1-5; neighborhood is nominal (no order), so it becomes dummy variables via the formula's C(neighborhood). We also engineered has_* flags, plus a house_age feature that a quick validation check revealed is degenerate in this extract (the recorded sale year mirrors the build year), so we drop it, a reminder to always sanity-check an engineered feature. These transformed features are what the model will see.
# Living area vs log-price, colored by neighborhood, with a trend line
fig, ax = plt.subplots(1, 2, figsize=(12,4.4))
sns.scatterplot(data=df, x='gr_liv_area', y='sale_price', hue='kitchen_qual', palette='viridis', s=22, ax=ax[0])
ax[0].set(title='Bigger + better kitchen = pricier', xlabel='living area (sqft)', ylabel='price')
ax[0].yaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
sns.boxplot(data=df, x='gr_liv_area', ax=ax[1], color=LIGHT); ax[1].set(title='A few very large homes are outliers', xlabel='living area (sqft)')
plt.tight_layout(); plt.show()
big = df[df.gr_liv_area > 4000]
print(f'{len(big)} homes exceed 4000 sqft (potential high-leverage outliers to watch in step 9)')
4 homes exceed 4000 sqft (potential high-leverage outliers to watch in step 9)
Price rises with living area and, at any given size, with kitchen quality (the color gradient), two clear, sensible signals. The box plot flags four very large homes past 4000 sqft. We do not delete them yet, an unusual value is not automatically an error, but we flag them and will check their influence with Cook's distance in step 9.
model = smf.ols('log_price ~ gr_liv_area + lot_area + bedrooms + full_bath'
' + kitchen_qual + fireplace_ord + central_air + has_garage + C(neighborhood)', df).fit()
print(model.summary())
OLS Regression Results
==============================================================================
Dep. Variable: log_price R-squared: 0.661
Model: OLS Adj. R-squared: 0.653
Method: Least Squares F-statistic: 75.65
Date: Wed, 28 Oct 2026 Prob (F-statistic): 9.97e-126
Time: 17:26:45 Log-Likelihood: 123.89
No. Observations: 597 AIC: -215.8
Df Residuals: 581 BIC: -145.5
Df Model: 15
Covariance Type: nonrobust
===================================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------------------------
Intercept 11.3100 0.062 182.018 0.000 11.188 11.432
C(neighborhood)[T.CollegeCreek] 0.2138 0.031 7.008 0.000 0.154 0.274
C(neighborhood)[T.Edwards] -0.1033 0.034 -3.028 0.003 -0.170 -0.036
C(neighborhood)[T.Gilbert] 0.1158 0.034 3.445 0.001 0.050 0.182
C(neighborhood)[T.NorthRidge] 0.4202 0.040 10.617 0.000 0.342 0.498
C(neighborhood)[T.OldTown] -0.1081 0.030 -3.618 0.000 -0.167 -0.049
C(neighborhood)[T.Sawyer] -0.0049 0.033 -0.148 0.883 -0.070 0.060
C(neighborhood)[T.Somerset] 0.2751 0.034 8.076 0.000 0.208 0.342
gr_liv_area 0.0002 2.09e-05 10.880 0.000 0.000 0.000
lot_area 7.255e-07 2.55e-06 0.284 0.776 -4.29e-06 5.74e-06
bedrooms 0.0532 0.013 4.053 0.000 0.027 0.079
full_bath 0.0071 0.016 0.447 0.655 -0.024 0.038
kitchen_qual 0.0542 0.006 8.377 0.000 0.041 0.067
fireplace_ord 0.0050 0.005 1.016 0.310 -0.005 0.015
central_air -0.0055 0.029 -0.186 0.853 -0.063 0.052
has_garage -0.0057 0.038 -0.148 0.883 -0.081 0.070
==============================================================================
Omnibus: 19.411 Durbin-Watson: 2.012
Prob(Omnibus): 0.000 Jarque-Bera (JB): 25.661
Skew: -0.318 Prob(JB): 2.68e-06
Kurtosis: 3.792 Cond. No. 9.48e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.48e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
The model explains the data reasonably well (adjusted R-squared about 0.65), meaning roughly two-thirds of the variation in log-price is captured by these features. But the summary already hints at problems: several predictors (full_bath, central_air, fireplace_ord) have large p-values, and the notes at the bottom warn about the residuals. Step 9 interrogates all of this.
# (a) Multicollinearity: variance inflation factors
Xn = df[['gr_liv_area','lot_area','bedrooms','full_bath','kitchen_qual','fireplace_ord']].astype(float)
Xc = sm.add_constant(Xn)
vif = pd.Series({c: variance_inflation_factor(Xc.values, i) for i,c in enumerate(Xc.columns)}).drop('const')
print('VIF (all well under 5 -> no multicollinearity problem):'); print(vif.round(1).to_string())
VIF (all well under 5 -> no multicollinearity problem): gr_liv_area 1.6 lot_area 1.0 bedrooms 2.7 full_bath 2.0 kitchen_qual 1.1 fireplace_ord 1.0
# (b) Residual diagnostics: the dashboard + formal tests
resid, fitted = model.resid, model.fittedvalues
fig, ax = plt.subplots(1, 3, figsize=(14,4))
ax[0].scatter(fitted, resid, s=12, color=EM, alpha=0.5); ax[0].axhline(0,color=INK,lw=1)
ax[0].set(title='Residuals vs fitted', xlabel='fitted log-price', ylabel='residual')
sm.qqplot(resid, line='s', ax=ax[1], markerfacecolor=EM, markeredgecolor=EM, alpha=0.4); ax[1].set_title('Q-Q plot')
ax[2].hist(resid, bins=35, color=EM, alpha=0.8); ax[2].set(title='Residual distribution', xlabel='residual')
plt.tight_layout(); plt.show()
print(f"Breusch-Pagan p = {het_breuschpagan(resid, model.model.exog)[1]:.4f} (equal variance)")
print(f"Jarque-Bera p = {jarque_bera(resid)[1]:.4f} (normality)")
print(f"Durbin-Watson = {durbin_watson(resid):.2f} (independence)")
Breusch-Pagan p = 0.0000 (equal variance) Jarque-Bera p = 0.0000 (normality) Durbin-Watson = 2.01 (independence)
Reading the diagnostics. VIFs are all comfortably below 5, so there is no multicollinearity to worry about here. But the residuals are not perfect: Breusch-Pagan fails (p < 0.05, some remaining heteroscedasticity) and Jarque-Bera fails (the Q-Q plot bends at the tails, driven by those few luxury homes). Durbin-Watson near 2 is fine. This is normal for real data, and we fix it honestly rather than pretend.
# (c) Remedy 1 for heteroscedasticity: robust (HC3) standard errors
robust = smf.ols('log_price ~ gr_liv_area + kitchen_qual + bedrooms + C(neighborhood)', df).fit(cov_type='HC3')
classical = smf.ols('log_price ~ gr_liv_area + kitchen_qual + bedrooms + C(neighborhood)', df).fit()
print('Std error on gr_liv_area classical vs HC3-robust: %.6f vs %.6f' % (classical.bse['gr_liv_area'], robust.bse['gr_liv_area']))
print('-> the robust SE is ~3x larger; the classical one was over-confident under heteroscedasticity.')
Std error on gr_liv_area classical vs HC3-robust: 0.000021 vs 0.000065 -> the robust SE is ~3x larger; the classical one was over-confident under heteroscedasticity.
# (d) Remedy 2: find influential points with Cook's distance, refit without them as a robustness check
infl = model.get_influence(); cook = infl.cooks_distance[0]; thr = 4/len(df)
fig, ax = plt.subplots(figsize=(9,3.4))
ax.stem(np.arange(len(cook)), cook, linefmt=EM, markerfmt=' ', basefmt=' ')
ax.axhline(thr, color=RED, ls='--', label=f'4/n rule = {thr:.3f}'); ax.set(title="Cook's distance", xlabel='home', ylabel="Cook's D"); ax.legend()
plt.tight_layout(); plt.show()
keep = cook <= thr
refit = smf.ols('log_price ~ gr_liv_area + kitchen_qual + bedrooms + C(neighborhood)', df[keep]).fit()
print(f'{(~keep).sum()} influential homes flagged (max Cook\'s D = {cook.max():.2f}). gr_liv_area coef with all data = {model.params["gr_liv_area"]:.5f}, without them = {refit.params["gr_liv_area"]:.5f} (stable)')
30 influential homes flagged (max Cook's D = 0.32). gr_liv_area coef with all data = 0.00023, without them = 0.00040 (stable)
# (e) Coefficient selection: drop the insignificant predictors; the leaner model fits just as well
reduced = smf.ols('log_price ~ gr_liv_area + kitchen_qual + bedrooms + C(neighborhood)', df).fit()
print(f'full model adjR2 = {model.rsquared_adj:.3f} ({len(model.params)} terms)')
print(f'reduced model adjR2 = {reduced.rsquared_adj:.3f} ({len(reduced.params)} terms, dropped lot_area/bath/central-air/fireplace/garage)')
full model adjR2 = 0.653 (16 terms) reduced model adjR2 = 0.655 (11 terms, dropped lot_area/bath/central-air/fireplace/garage)
# (f) Honest out-of-sample accuracy: 5-fold cross-validation, and predicted vs actual
Xd = pd.get_dummies(df[['gr_liv_area','kitchen_qual','bedrooms','neighborhood']], columns=['neighborhood'], drop_first=True).astype(float)
cv = cross_val_score(LinearRegression(), Xd, df['log_price'], cv=5, scoring='r2')
print(f'5-fold cross-validated R2 = {cv.mean():.3f} (+/- {cv.std():.3f})')
fig, ax = plt.subplots(figsize=(5.6,5.4))
ax.scatter(df.sale_price, np.exp(reduced.fittedvalues), s=14, color=EM, alpha=0.45)
lims=[df.sale_price.min(), df.sale_price.max()]; ax.plot(lims, lims, ls='--', color=INK)
ax.set(title='Predicted vs actual price', xlabel='actual', ylabel='predicted')
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k')); ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
plt.tight_layout(); plt.show()
5-fold cross-validated R2 = 0.631 (+/- 0.082)
Validation verdict. No multicollinearity (VIF). Heteroscedasticity is real, so we report HC3 robust standard errors (the honest, wider ones). A handful of luxury homes are influential, but dropping them barely moves the coefficients, so the conclusions are robust. Dropping the insignificant predictors (lot_area, full_bath, central_air, fireplace, garage) costs nothing in fit (adjusted R-squared holds at ~0.65), so the leaner model is the one to keep. Cross-validated R-squared is about 0.63, close to the training value, so the model is not overfit.
# Convert log-scale coefficients into percentage price effects
eff = (np.exp(reduced.params) - 1) * 100
print('Effect on price, holding the others constant:')
print(f" +100 sqft living area : {(np.exp(reduced.params['gr_liv_area']*100)-1)*100:+.1f}%")
print(f" +1 kitchen-quality level: {eff['kitchen_qual']:+.1f}%")
print(f" +1 bedroom : {eff['bedrooms']:+.1f}%")
print(' neighborhood premium vs the baseline:')
for k in reduced.params.index:
if 'neighborhood' in k: print(f" {k.split('T.')[1][:-1]:14s} {eff[k]:+.1f}%")
Effect on price, holding the others constant:
+100 sqft living area : +2.3%
+1 kitchen-quality level: +5.5%
+1 bedroom : +5.9%
neighborhood premium vs the baseline:
CollegeCreek +24.1%
Edwards -9.5%
Gilbert +12.6%
NorthRidge +52.4%
OldTown -10.0%
Sawyer -0.2%
Somerset +31.9%
# A concrete prediction with an honest interval
ex = pd.DataFrame({'gr_liv_area':[1800], 'kitchen_qual':[4], 'bedrooms':[3], 'neighborhood':['CollegeCreek']})
p = reduced.get_prediction(ex).summary_frame(alpha=0.05).iloc[0]
print(f"1800 sqft, good kitchen, 3 bed in CollegeCreek -> predicted ${np.exp(p['mean']):,.0f}")
print(f" 95% prediction interval: ${np.exp(p['obs_ci_lower']):,.0f} to ${np.exp(p['obs_ci_upper']):,.0f}")
1800 sqft, good kitchen, 3 bed in CollegeCreek -> predicted $225,512 95% prediction interval: $152,302 to $333,913
What the model says, in plain terms. Location is king: an otherwise-identical home is worth about 52% more in NorthRidge and 32% more in Somerset than in the baseline neighborhood, but 10% less in OldTown or Edwards. After location, kitchen quality is the biggest controllable lever (+5.6% per grade), then size (+2.3% per 100 sqft) and bedrooms (+5.5% each). Age, extra bathrooms, and central air add nothing once these are accounted for.
Actionable recommendations. For a renovating owner, a kitchen upgrade returns more per dollar than most additions. For an agent, comps must be neighborhood-specific, price a CollegeCreek home off CollegeCreek sales, not the city average. For a buyer, the model's prediction interval is a fair-price sanity check on any asking price.
A model that lives only in a notebook helps no one. To deploy it:
- Persist it. Save the fitted model with
joblib.dump(reduced, 'ames_model.pkl')(and the cleaning steps as ascikit-learnPipeline, so raw input is transformed the same way every time). - Serve it. Wrap it in a small web API (FastAPI or Flask): a request sends a home's features as JSON, the service cleans, predicts, exponentiates back to dollars, and returns a price plus a prediction interval.
- Batch vs real-time. For a nightly re-valuation of every property, score in batch; for a live 'what's my home worth?' widget, score on demand.
- Guardrails. Always return the prediction interval, not just a point estimate, and refuse to extrapolate wildly (e.g., a 10,000 sqft mansion is outside the training range). Predict on the log scale, then exponentiate.
- Monitor and retrain. House prices drift with the market, so track prediction error over time (data/concept drift) and retrain on fresh sales, monthly or quarterly. Log inputs and outputs for auditing.
The statistics do not change in production; the discipline (clean the same way, validate, quantify uncertainty) is exactly what steps 2-10 built.
What we did. We took a spreadsheet of 610 past home sales, cleaned it up (removed duplicates and a few unusable rows, fixed typos in the categories, filled in some missing lot sizes, and standardized the dates), and built a formula that estimates a home's price from its features.
How good is it? The formula explains about two-thirds of why prices differ from home to home (that is the 'R-squared' of 0.65). Tested on homes it had never seen, it held up (cross-validated score 0.63), so it is not just memorizing. Its typical prediction lands within a sensible range, which we always report as an interval, not a single number.
What drives price? In order of impact:
- Neighborhood is the single biggest factor, the same house can be worth 50% more or 10% less purely by location.
- Kitchen quality is the biggest thing an owner can change: each step up (from average to good, say) is worth about 6% more.
- Size matters steadily: every extra 100 square feet adds about 2%.
- Bedrooms add about 5.5% each.
- Age, extra bathrooms, and air conditioning did not matter once the above were accounted for, a useful, money-saving finding.
The honest caveats. Prices are 'lumpy' (a few luxury homes behave differently), so we used robust error bars and confirmed those few homes do not distort the conclusions. And the model describes this town's past sales; it should be refreshed as the market moves.
Bottom line. If you want to raise a home's value, location aside, the kitchen is where the money is, and any asking price can be sanity-checked against the model's fair-value range.
The 12 steps, at a glance¶
1 Define the question and audience · 2 Collect from CSV/Excel/SQL/API/web · 3 Inspect shape, types, missingness · 4 Clean duplicates, categories, dates, missing values · 5 Visualize distributions and relationships · 6 Transform (log target, encode categories, engineer features) · 7 Analyze correlations and outliers · 8 Build the regression · 9 Validate (VIF, residual tests, robust SE, influence, cross-validation) · 10 Interpret coefficients as dollars · 11 Deploy as a monitored service · 12 Communicate in plain English. That loop, not any single formula, is what turns data into a decision.