Supervised Machine Learning with scikit-learn

Supervised learning is the branch of machine learning where every training example carries a label. You show the model inputs paired with known outputs, it learns the mapping between them, and then you ask it to predict outputs for inputs it has never seen. This article is a practical guide to the complete supervised learning workflow in scikit-learn: fitting models, evaluating them honestly, preventing overfitting with regularization, searching for the best hyperparameters, and packaging everything into a reusable pipeline.

Please see a 10 code along for practice here: https://datalad.co.uk/supervised-machine-learning-with-scikit-learn-10-code-along-examples/.

A free cheatsheet here: https://datalad.co.uk/supervised-machine-learning-with-scikit-learn-cheatsheet/

The Workflow

Every supervised learning project follows the same spine:

1. Prepare features (X) and target (y), encode categorical variables, handle missing values, and scale.
2. Split into training and test sets before touching any model.
3. Fit a model on the training set.
4. Evaluate on the held-out test set.
5. Improve with cross-validation, regularization, and hyperparameter search.
6. Package the full process in a pipeline.

K-Nearest Neighbours

K-Nearest Neighbors (KNN) classifies a new point by finding the k most similar training examples and returning a majority vote. There is no mathematical training step: the algorithm stores the training data and, at prediction time, computes distances to find neighbors.

Fit and predict

import numpy as np
from sklearn.neighbors import KNeighborsClassifier
y = customers_df["churned"].values
X = customers_df[["tenure", "support_calls"]].values
knn = KNeighborsClassifier(n_neighbors=6)
knn.fit(X, y)
X_new = np.array([[30.0, 3],
                  [107.0, 8],
                  [213.0, 2]])
y_pred = knn.predict(X_new)
print("Predictions: {}".format(y_pred))

y = ...values extracts the target as a 1D NumPy array; X = ...values extracts features as a 2D array. Scikit-learn requires both formats. .fit() stores the training data for distance lookups later; .predict() finds the 6 nearest training examples for each row of X_new and returns the majority class.

Train-test split

Evaluating a model on the same data used to train it is like grading a student on questions they already memorized the answers to. Split the data first, and only ever touch the test set once, at the end.

from sklearn.model_selection import train_test_split
X = customers_df.drop("churned", axis=1).values
y = customers_df["churned"].values
X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.20,
    random_state=42,
    stratify=y        # preserves class proportions in both splits
)
knn = KNeighborsClassifier(n_neighbors=5)
knn.fit(X_train, y_train)
print(knn.score(X_test, y_test))

stratify=y ensures that if the original dataset is 70% class 0 and 30% class 1, both the training and test splits carry the same proportions. Without it, a random split could accidentally concentrate one class in training and leave very little of it in the test set.

Finding the right k: the model complexity curve

A very small k makes the model hypersensitive to individual training points, fitting noise rather than signal. A very large k smooths the decision boundary too aggressively, missing real patterns. The right value sits between those extremes.

neighbors = np.arange(1, 13)
train_accuracies = {}
test_accuracies = {}
for neighbor in neighbors:
    knn = KNeighborsClassifier(n_neighbors=neighbor)
    knn.fit(X_train, y_train)
    train_accuracies[neighbor] = knn.score(X_train, y_train)
    test_accuracies[neighbor]  = knn.score(X_test, y_test)
plt.title("KNN: Varying Number of Neighbors")
plt.plot(neighbors, train_accuracies.values(), label="Training Accuracy")
plt.plot(neighbors, test_accuracies.values(),  label="Testing Accuracy")
plt.legend()
plt.xlabel("Number of Neighbors")
plt.ylabel("Accuracy")
plt.show()

The loop trains a fresh KNN for each k from 1 to 12 and records both training and test accuracy. On the resulting plot, the best k is where test accuracy peaks and the gap between the two curves is smallest: high test accuracy indicates a good fit, and a small gap means the model generalizes rather than memorizing.

Linear Regression

Simple regression (one feature)

from sklearn.linear_model import LinearRegression
X = ad_spend_df['radio_budget'].values.reshape(-1, 1)
y = ad_spend_df['revenue'].values
reg = LinearRegression()
reg.fit(X, y)
predictions = reg.predict(X)
plt.scatter(X, y, color="blue")
plt.plot(X, predictions, color="red")
plt.xlabel("Radio Budget ($)")
plt.ylabel("Revenue ($)")
plt.show()

.reshape(-1, 1) converts a 1D array of shape (n,) into a 2D column of shape (n, 1). Scikit-learn always expects features to be 2D, even when there is only one. .fit() finds the slope and intercept that minimize the sum of squared distances between the data points and the line. The scatter shows the raw data in blue; the fitted line in red gives a visual check of how well the model captures the trend.

Multiple regression

X = ad_spend_df.drop("revenue", axis=1).values
y = ad_spend_df["revenue"].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
reg = LinearRegression()
reg.fit(X_train, y_train)
y_pred = reg.predict(X_test)
print("Predictions: {}, Actual Values: {}".format(y_pred[:2], y_test[:2]))

The API is identical to simple regression. With multiple features, the model fits a hyperplane through higher-dimensional space, estimating a separate coefficient for each input.

Regression metrics

from sklearn.metrics import root_mean_squared_error
r_squared = reg.score(X_test, y_test)
rmse      = root_mean_squared_error(y_test, y_pred)
print("R²: {}".format(r_squared))
print("RMSE: {}".format(rmse))

R² is a proportion: 0.85 means the model explains 85% of the variance in revenue. A value of 1.0 is a perfect fit; 0.0 means the model does no better than always predicting the mean.

RMSE expresses the average prediction error in the original units of the target. If revenue is measured in thousands of dollars and RMSE is 4.2, the model is off by roughly $4,200 on average. Squaring the residuals before averaging means RMSE penalizes large errors more heavily than small ones.

Cross-Validation

A single train-test split produces a single score. A different random seed could give a noticeably different result, making it hard to know whether you are measuring the model’s quality or the luck of the split. Cross-validation runs the evaluation on multiple different splits and reports the distribution of results.

from sklearn.model_selection import cross_val_score, KFold
kf = KFold(n_splits=6, shuffle=True, random_state=5)
reg = LinearRegression()
cv_scores = cross_val_score(reg, X, y, cv=kf)
print(cv_scores)
print(np.mean(cv_scores))                       # average performance
print(np.std(cv_scores))                        # how consistent the scores are
print(np.quantile(cv_scores, [0.025, 0.975]))   # 95% interval

KFold(n_splits=6) divides the data into 6 equal chunks. cross_val_score trains on 5 and tests on the remaining 1, rotates through all 6 combinations, and returns an array of 6 scores. shuffle=True randomizes the data order before splitting so that any natural ordering does not bias the folds.

A model with high mean and low standard deviation across folds is consistently good. A model with high mean but high variance is only reliable on some slices of the data.

Regularization: Ridge and Lasso

Standard linear regression minimizes prediction error with no constraints on the size of coefficients. When the model is complex and data is noisy, it can memorize training idiosyncrasies and generalize poorly. Regularization adds a penalty for large coefficients, discouraging overfitting. The penalty strength is controlled by alpha: higher alpha means heavier penalization.

Ridge: searching over alpha values

from sklearn.linear_model import Ridge
alphas = [0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0]
ridge_scores = []
for alpha in alphas:
    ridge = Ridge(alpha=alpha)
    ridge.fit(X_train, y_train)
    score = ridge.score(X_test, y_test)
    ridge_scores.append(score)
print(ridge_scores)

At low alpha, Ridge behaves like ordinary linear regression. At high alpha, it forces all coefficients close to zero and the model becomes more conservative. This loop tests six values and collects R² for each, making it easy to spot where performance peaks before the penalty becomes too aggressive.

Lasso: built-in feature selection

from sklearn.linear_model import Lasso
lasso = Lasso(alpha=0.3)
lasso.fit(X, y)
lasso_coef = lasso.coef_
plt.bar(feature_names, lasso_coef)
plt.xticks(rotation=45)
plt.show()

Lasso’s penalty is strong enough to push individual coefficients all the way to zero, removing those features from the model entirely. The bar chart makes this visible: tall bars are the features that drive the prediction; any bar at zero was eliminated. This makes Lasso useful not just for prediction but for understanding which inputs actually matter.

Classification Metrics

Accuracy counts the fraction of correct predictions, but it can be misleading when classes are imbalanced. A model that always predicts “no churn” on a dataset where 95% of customers stay would report 95% accuracy while being completely useless for the minority class. The confusion matrix and derived metrics give a more complete picture.

from sklearn.metrics import confusion_matrix, classification_report
knn = KNeighborsClassifier(n_neighbors=6)
knn.fit(X_train, y_train)
y_pred = knn.predict(X_test)
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))

The confusion matrix layout:

                 Predicted 0      Predicted 1
Actual 0       True Negative    False Positive
Actual 1       False Negative   True Positive

From those four numbers, three core metrics follow:

Precision and recall pull in opposite directions. Lowering the classification threshold catches more true positives (higher recall) but also more false positives (lower precision). Which error type is more costly depends on the application: in fraud detection, missing a fraud is worse than a false alarm; in spam filtering, the opposite may be true.

classification_report prints all three metrics for each class in a single formatted table.

Logistic Regression and ROC Curves

Logistic regression fits an S-shaped curve that outputs a probability between 0 and 1 rather than an unbounded number. Those probabilities are then compared against a threshold to produce a class prediction.

Getting probability estimates

from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
logreg.fit(X_train, y_train)
# predict_proba returns [[P(class=0), P(class=1)], ...]
# [:, 1] selects P(class=1) for every test sample
y_pred_probs = logreg.predict_proba(X_test)[:, 1]

predict_proba returns a two-column array. Column 0 holds P(class=0) and column 1 holds P(class=1) for each sample. [:, 1]slices out the positive-class probability across all rows. These raw scores are what the ROC curve operates on.

The ROC curve

from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y_test, y_pred_probs)
plt.plot([0, 1], [0, 1], 'k--')   # diagonal = random guessing
plt.plot(fpr, tpr)
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve for Customer Churn')
plt.show()

The ROC curve sweeps through every possible classification threshold and records the resulting false positive rate and true positive rate at each point. The dashed diagonal is the baseline: a model that guesses randomly at every threshold. A curve that bows toward the top-left corner indicates a model that catches most true positives before accumulating many false positives.

AUC

from sklearn.metrics import roc_auc_score
print(roc_auc_score(y_test, y_pred_probs))

AUC (area under the curve) collapses the ROC curve to a single number. The intuitive interpretation: AUC is the probability that the model assigns a higher score to a randomly chosen positive example than to a randomly chosen negative one.

Hyperparameter Tuning

GridSearchCV: exhaustive search

GridSearchCV tests every combination in a defined parameter grid. It is the right tool when the search space is small enough that exhaustive enumeration is affordable.

from sklearn.model_selection import GridSearchCV
param_grid = {"alpha": np.linspace(0.00001, 1, 20)}
lasso_cv = GridSearchCV(lasso, param_grid, cv=kf)
lasso_cv.fit(X_train, y_train)
print("Tuned parameters: {}".format(lasso_cv.best_params_))
print("Tuned score: {}".format(lasso_cv.best_score_))

np.linspace(0.00001, 1, 20) generates 20 evenly spaced alpha values. For each value, GridSearchCV runs full cross-validation and records the average score. With 20 values and a 6-fold CV setup, this fits 120 models. best_params_ and best_score_report the winner.

RandomizedSearchCV: sampling large spaces

When the parameter space is large, exhaustive enumeration becomes expensive. RandomizedSearchCV samples a random subset of combinations instead.

from sklearn.model_selection import RandomizedSearchCV
params = {
    "penalty":      ["l1", "l2"],
    "tol":          np.linspace(0.0001, 1.0, 50),
    "C":            np.linspace(0.1, 1.0, 50),
    "class_weight": ["balanced", {0: 0.8, 1: 0.2}]
}
logreg_cv = RandomizedSearchCV(logreg, params, cv=kf)
logreg_cv.fit(X_train, y_train)
print("Best Parameters: {}".format(logreg_cv.best_params_))
print("Best Accuracy: {}".format(logreg_cv.best_score_))

The full grid here contains 2 × 50 × 50 × 2 = 10,000 combinations. Testing all of them would be prohibitive. By default, RandomizedSearchCV samples 10 combinations and returns the best. The trade-off: a small chance of missing the absolute optimum in exchange for a large reduction in compute time. In practice, the benefit of exhaustive search over a well-sampled random search is rarely worth the cost.

Preprocessing

Encoding categorical variables

tracks_encoded = pd.get_dummies(tracks_df, drop_first=True)

Models require numerical inputs. pd.get_dummies converts each categorical column into binary indicator columns: a stylecolumn with three values (Pop, Rock, Jazz) becomes two binary columns, style_Rock and style_Jazz. If both are 0, the model infers the track is Pop. drop_first=True removes one column per categorical variable, because knowing two out of three categories implies the third. Keeping all three would introduce perfect multicollinearity, which can destabilize linear models.

Handling missing values

# Option A: drop rows with missing values
tracks_df = tracks_df.dropna(subset=["style", "popularity", "loudness"])
# Option B: impute (typically used inside a Pipeline)
from sklearn.impute import SimpleImputer
imputer = SimpleImputer()   # default strategy='mean'

Dropping rows is simple and makes no assumptions, but it wastes data. Imputing fills each gap with the column mean, preserving every row at the cost of introducing artificial values. The right choice depends on how many rows you can afford to lose and whether the mean is a meaningful fill-in for the feature in question.

Scaling features

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)  # learn mean and std from training data, then apply
X_test_scaled  = scaler.transform(X_test)        # apply the same training statistics to test data

StandardScaler converts each feature to mean 0 and standard deviation 1. Without scaling, features measured in large units (annual revenue in dollars) dominate distance-based models over features in small units (number of calls) purely because of scale, not importance.

The critical rule: call fit_transform on training data only, then transform on the test set. Calling fit on test data would use information from the test set to determine the scaling parameters, which is data leakage.

Pipelines

A pipeline chains preprocessing and modeling into a single object. .fit() runs every step in order; .predict() on new data applies those same steps automatically, using the statistics learned from training.

Imputer + classifier

from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
steps = [
    ("imputer", SimpleImputer()),
    ("knn",     KNeighborsClassifier(n_neighbors=3))
]
pipeline = Pipeline(steps)
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)
print(confusion_matrix(y_test, y_pred))

Each step is a (name, object) tuple. The name is used later to reference parameters in grid search. .fit() runs imputer.fit_transform(X_train) and feeds the result into knn.fit(). .predict() imputes X_test with the same column means learned from training, then runs KNN. The correct order and correct statistics are applied automatically on every call.

Scaler + regressor

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Lasso
steps = [
    ("scaler", StandardScaler()),
    ("lasso",  Lasso(alpha=0.5))
]
pipeline = Pipeline(steps)
pipeline.fit(X_train, y_train)
print(pipeline.score(X_test, y_test))

Pipeline with GridSearchCV

To tune a parameter inside a pipeline step, use double-underscore notation: stepname__parameter.

steps = [
    ("scaler", StandardScaler()),
    ("logreg", LogisticRegression())
]
pipeline = Pipeline(steps)
parameters = {"logreg__C": np.linspace(0.001, 1.0, 20)}
cv = GridSearchCV(pipeline, param_grid=parameters)
cv.fit(X_train, y_train)
print(cv.best_score_, cv.best_params_)

"logreg__C" means “go to the step named logreg and adjust its C parameter.” For each value of C, GridSearchCV runs the full pipeline with cross-validation, scaling the data before fitting logistic regression each time.

Full pipeline: imputer + scaler + model + grid search

steps = [
    ("imp_mean", SimpleImputer()),
    ("scaler",   StandardScaler()),
    ("logreg",   LogisticRegression())
]
pipeline = Pipeline(steps)
params = {
    "logreg__solver": ["newton-cg", "saga", "lbfgs"],
    "logreg__C":      np.linspace(0.001, 1.0, 10)
}
tuning = GridSearchCV(pipeline, param_grid=params)
tuning.fit(X_train, y_train)
y_pred = tuning.predict(X_test)
print("Best Parameters: {}, Accuracy: {}".format(
    tuning.best_params_,
    tuning.score(X_test, y_test)
))

This is the complete supervised ML workflow packaged as a single reusable object. Pass raw, potentially missing training data to .fit() and it handles imputation, scaling, and hyperparameter search automatically. When you call .predict() on new data, the same imputation and scaling learned from training are applied, with no risk of forgetting a step or leaking test statistics.

Comparing Models

Regression: cross-validation comparison

models = {
    "Linear Regression": LinearRegression(),
    "Ridge":             Ridge(alpha=0.1),
    "Lasso":             Lasso(alpha=0.1)
}
results = []
for model in models.values():
    kf = KFold(n_splits=6, random_state=42, shuffle=True)
    cv_scores = cross_val_score(model, X_train, y_train, cv=kf)
    results.append(cv_scores)
plt.boxplot(results, labels=models.keys())
plt.show()

A boxplot of six cross-validation scores captures more than a single mean. A model with a slightly lower median but a much smaller interquartile range may be the safer production choice, because it performs consistently rather than relying on a fortunate split.