Generalized Linear Models › Implementation
Machine Learning from ScratchAll chapters Generalized Linear Models › ImplementationOne section · Sections in this chapter8 / 12 ↓
Reading tools Definitions & inline links
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
In this section, we'll demonstrate how to implement Generalized Linear Models in Python using scikit-learn. We'll work through practical examples of logistic regression (for binary classification) and Poisson regression (for count data), showing how to train models, evaluate performance, and interpret results.
Logistic Regression for Binary Classification
Let's start with a logistic regression example using the breast cancer dataset. This demonstrates how to use a GLM with a binomial distribution and logit link for binary classification.
Data Preparation
We'll start by loading the breast cancer dataset and preparing it for modeling. This dataset contains features computed from digitized images of breast mass, and our goal is to predict whether a tumor is malignant or benign.
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
# Load the breast cancer dataset
data = load_breast_cancer()
X = data.data
y = data.target
# Split into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42, stratify=y
)
# Standardize features for better convergence
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)The dataset is split into 80% training and 20% testing data using stratified sampling to maintain the class distribution. Feature standardization transforms all variables to have zero mean and unit variance, helping the logistic regression optimizer converge faster and produce more stable coefficient estimates. Without standardization, features with larger scales would dominate the model.
Model Training
Now we'll train a logistic regression model using L2 regularization (ridge penalty) with the L-BFGS solver, which is efficient for this dataset size:
# Train logistic regression model
log_reg = LogisticRegression(
penalty="l2", C=1.0, solver="lbfgs", max_iter=1000, random_state=42
)
log_reg.fit(X_train_scaled, y_train)
# Make predictions
y_pred = log_reg.predict(X_test_scaled)
y_pred_proba = log_reg.predict_proba(X_test_scaled)[:, 1]
# Calculate performance metrics
accuracy = accuracy_score(y_test, y_pred)
roc_auc = roc_auc_score(y_test, y_pred_proba)Model Performance
Let's examine the model's performance on the test set:
Logistic Regression Performance ======================================== Accuracy: 0.9825 ROC AUC Score: 0.9954 Number of iterations: 19
With an accuracy above 95% and a ROC AUC score close to 1.0, the model demonstrates excellent performance in distinguishing between malignant and benign tumors. The ROC AUC score near 1.0 is particularly important.it indicates that the model can effectively rank predictions, placing malignant cases higher than benign ones across different probability thresholds. The convergence in a reasonable number of iterations confirms that the optimization process was stable and the model parameters are reliable.
Classification Report:
----------------------------------------
precision recall f1-score support
malignant 0.98 0.98 0.98 42
benign 0.99 0.99 0.99 72
accuracy 0.98 114
macro avg 0.98 0.98 0.98 114
weighted avg 0.98 0.98 0.98 114The classification report reveals balanced performance across both tumor types. Precision values above 0.95 indicate that when the model predicts a tumor is malignant (or benign), it's correct more than 95% of the time. Similarly, recall values above 0.95 show that the model successfully identifies more than 95% of actual malignant (and benign) cases. The F1-scores, which balance precision and recall, exceed 0.95 for both classes, confirming that the model doesn't favor one class over the other, an important property for medical diagnosis applications.
Confusion Matrix:
----------------------------------------
Predicted Negative Predicted Positive
Actual Negative 41 1
Actual Positive 1 71The confusion matrix provides a clear picture of where the model succeeds and where it occasionally errs. The diagonal elements (true negatives and true positives) are substantially larger than the off-diagonal elements (false positives and false negatives), indicating that most predictions are correct. In medical contexts, false negatives (predicting benign when the tumor is malignant) are particularly concerning because they could delay treatment. The very low count of false negatives demonstrates that the model is reliable for clinical screening purposes.
Coefficient Interpretation
Let's examine the most influential features:
# Get feature importance based on coefficient magnitude
feature_importance = pd.DataFrame(
{"feature": data.feature_names, "coefficient": log_reg.coef_[0]}
).sort_values("coefficient", key=abs, ascending=False)
top_features = feature_importance.head(10)Top 10 Most Influential Features: ============================================================ worst texture : -1.2551 (decreases risk) radius error : -1.0830 (decreases risk) worst concave points : -0.9537 (decreases risk) worst area : -0.9478 (decreases risk) worst radius : -0.9476 (decreases risk) worst symmetry : -0.9392 (decreases risk) area error : -0.9291 (decreases risk) worst concavity : -0.8232 (decreases risk) worst perimeter : -0.7632 (decreases risk) worst smoothness : -0.7466 (decreases risk)
Each coefficient represents how much the log-odds of malignancy change when that feature increases by one standard deviation (since we standardized the features). Positive coefficients push predictions toward malignancy, while negative coefficients push toward benign. The magnitude of each coefficient indicates its importance: features with larger absolute values exert stronger influence on the final prediction. For instance, if "worst concave points" has a large positive coefficient, it means that higher values of this feature strongly suggest malignancy.
Poisson Regression for Count Data
Now let's demonstrate Poisson regression, which is appropriate for modeling count data. We'll create a synthetic dataset representing the number of customer visits to a store based on various factors.
Data Generation
# Set random seed for reproducibility
np.random.seed(42)
# Generate synthetic data: predicting number of customer visits
n_samples = 1000
# Features: advertising spend, day of week, temperature, is_weekend
advertising = np.random.uniform(0, 100, n_samples)
day_of_week = np.random.randint(0, 7, n_samples)
temperature = np.random.normal(20, 5, n_samples)
is_weekend = (day_of_week >= 5).astype(int)
X_count = np.column_stack([advertising, temperature, is_weekend])
# Generate count data with Poisson distribution
# True relationship: log(lambda) = 2.0 + 0.02*advertising + 0.01*temperature + 0.5*is_weekend
# Equivalently: lambda = exp(2.0 + 0.02*advertising + 0.01*temperature + 0.5*is_weekend)
lambda_true = np.exp(
2.0 + 0.02 * advertising + 0.01 * temperature + 0.5 * is_weekend
)
y_count = np.random.poisson(lambda_true)
# Split data
X_count_train, X_count_test, y_count_train, y_count_test = train_test_split(
X_count, y_count, test_size=0.2, random_state=42
)This synthetic dataset simulates a realistic scenario where customer visits follow a Poisson distribution.a natural choice for count data. We've embedded a known relationship between visits and advertising spend, with additional effects from temperature and weekend status. The goal is to see if our Poisson regression model can recover these true relationships from the data.
Model Training
# Train Poisson regression model
poisson_reg = PoissonRegressor(alpha=0.1, max_iter=1000)
poisson_reg.fit(X_count_train, y_count_train)
# Make predictions
y_count_pred = poisson_reg.predict(X_count_test)
# Calculate performance metrics
mse = mean_squared_error(y_count_test, y_count_pred)
mae = mean_absolute_error(y_count_test, y_count_pred)
rmse = np.sqrt(mse)Model Performance
Poisson Regression Performance ======================================== Mean Absolute Error: 4.3530 Root Mean Squared Error: 5.7371 Mean Predicted Visits: 31.93 Mean Actual Visits: 31.67
The low MAE (Mean Absolute Error) and RMSE (Root Mean Squared Error) indicate that predictions are typically within a few visits of the actual counts.excellent performance for this type of problem. More importantly, the mean predicted visits closely match the mean actual visits, showing that the model has correctly learned both the overall scale and the relative patterns in the data. This suggests the model has successfully recovered the underlying Poisson process that generated the data.
Coefficient Analysis
# Extract coefficients
feature_names_count = ["Advertising Spend", "Temperature", "Is Weekend"]
coefficients_df = pd.DataFrame(
{
"Feature": feature_names_count,
"Coefficient": poisson_reg.coef_,
"Exp(Coefficient)": np.exp(poisson_reg.coef_),
}
)Poisson Regression Coefficients: ====================================================================== Intercept: 1.9958 (exp: 7.3578) Advertising Spend : 0.0201 (multiplier: 1.0203, +2.03% change) Temperature : 0.0102 (multiplier: 1.0103, +1.03% change) Is Weekend : 0.4968 (multiplier: 1.6435, +64.35% change)
In Poisson regression, coefficients represent the change in log-count for a one-unit increase in the predictor, holding all other variables constant. To interpret the practical effect:
- Coefficient (): The change in for a one-unit increase in the predictor
- Exponentiated coefficient (): The multiplicative effect on the expected count
- Percentage change: gives the percentage change in expected count
For example, if and , this means a one-unit increase in the predictor multiplies the expected count by 1.0202, which is a 2.02% increase. The weekend effect shows a substantial increase in expected visits, as designed in our data generation process.
Prediction Example
Let's make predictions for specific scenarios:
# Create scenarios for prediction
scenarios = pd.DataFrame(
{
"Advertising Spend": [50, 50, 80, 80],
"Temperature": [20, 20, 25, 25],
"Is Weekend": [0, 1, 0, 1],
}
)
# Make predictions
scenario_predictions = poisson_reg.predict(scenarios.values)Prediction Scenarios: ================================================================================ Scenario 1: $50 ad spend, 20°C, Weekday → Predicted visits: 24.6 Scenario 2: $50 ad spend, 20°C, Weekend → Predicted visits: 40.5 Scenario 3: $80 ad spend, 25°C, Weekday → Predicted visits: 47.3 Scenario 4: $80 ad spend, 25°C, Weekend → Predicted visits: 77.8
These scenario-based predictions illustrate how the model integrates multiple factors to produce realistic visit estimates. Comparing scenarios 1 and 2 (which differ only in weekend status) reveals the substantial weekend effect.visits increase significantly on weekends even with identical advertising and temperature. Similarly, comparing scenarios 1 and 3 shows how increased advertising spend and warmer temperature boost predicted visits. These patterns align perfectly with the coefficients we examined earlier, confirming that the model has learned interpretable and actionable relationships.
Key Parameters
Below are the main parameters that affect how GLMs work and perform in scikit-learn.
LogisticRegression Parameters
penalty: Type of regularization applied to prevent overfitting (default: 'l2'). Options include 'l1' (lasso), 'l2' (ridge), 'elasticnet', or 'none'. L2 regularization works well for most cases and is computationally efficient. L1 can perform feature selection by driving some coefficients to exactly zero.C: Inverse of regularization strength (default: 1.0). Smaller values (e.g., 0.01-0.1) apply stronger regularization and help prevent overfitting, while larger values (e.g., 10-100) reduce regularization. Start with 1.0 and adjust based on validation performance.solver: Optimization algorithm (default: 'lbfgs'). Options include 'lbfgs', 'liblinear', 'newton-cg', 'sag', and 'saga'. Use 'lbfgs' for most problems with L2 penalty, 'liblinear' for small datasets, and 'saga' for large datasets or when using L1/elasticnet penalties.max_iter: Maximum number of iterations for convergence (default: 100). Increase to 1000 or higher if you encounter convergence warnings. More iterations allow the optimizer more time to find the optimal solution.class_weight: Weights assigned to each class (default: None). Set to 'balanced' for imbalanced datasets to automatically adjust weights inversely proportional to class frequencies. This helps the model pay more attention to minority classes.random_state: Random seed for reproducibility (default: None). Set to an integer (e.g., 42) when using solvers with randomness ('sag', 'saga') to ensure consistent results across runs.
PoissonRegressor Parameters
alpha: L2 regularization strength (default: 1.0). Controls the penalty on coefficient magnitudes. Smaller values (0.01-0.1) reduce regularization for complex patterns, while larger values (10-100) increase regularization to prevent overfitting. Start with 1.0 and tune based on validation performance.max_iter: Maximum number of iterations (default: 100). Increase if the model doesn't converge within the default limit. Convergence warnings indicate the optimizer needs more iterations to find the optimal solution.tol: Convergence tolerance (default: 1e-4). The optimization stops when the change in parameters falls below this threshold. Smaller values (e.g., 1e-6) require more precise convergence but take longer to train.warm_start: Whether to reuse the previous solution as initialization (default: False). Set to True when fitting multiple times with different parameters (e.g., during hyperparameter tuning) to speed up convergence.
Key Methods
The following are the most commonly used methods for interacting with GLM models in scikit-learn.
fit(X, y): Trains the model on feature matrix X and target vector y using iteratively reweighted least squares (IRLS) to find maximum likelihood estimates. This method learns the optimal coefficients that best explain the relationship between features and the target.predict(X): Returns predicted class labels for classification (0 or 1) or predicted counts/values for regression. For logistic regression, predictions are based on whether the probability exceeds 0.5. For Poisson regression, returns the expected count for each observation.predict_proba(X): Returns probability estimates for each class (classification only). The output is a 2D array where each row sums to 1.0, with columns representing the probability of each class. Useful for setting custom decision thresholds or performing probability calibration.predict_log_proba(X): Returns natural logarithm of probability estimates (classification only). More numerically stable thanpredict_probawhen dealing with very small probabilities, which is important for avoiding underflow in downstream calculations.score(X, y): Returns the mean accuracy for classification (proportion of correct predictions) or the D² score for Poisson regression (deviance-based R² analog). Higher scores indicate better model performance on the given test data.decision_function(X): Returns the linear predictor values (η = Xβ) before applying the link function. For logistic regression, these are the log-odds. Useful for understanding model confidence and for custom threshold selection.