Logistic Regression – Predicting Prescriptions Using Logistic Regression

This guide introduces logistic regression through a real-world healthcare analytics example: predicting which prescribers are most likely to issue prescriptions for Ozempic (generic name Semaglutide). You’ll follow a full modeling workflow — from cleaning and preparing the data to interpreting the final model’s coefficients, evaluating predictive accuracy, and making new predictions.

By walking through this case step by step, you’ll learn how logistic regression can reveal the relationships between features like prescriber specialty, region, and cost patterns, and how these factors influence the likelihood of prescribing Semaglutide.

This project focuses on using logistic regression to uncover the patterns behind Semaglutide prescriptions in U.S. Medicare data. You’ll explore which prescriber features most strongly predict whether a prescription involves Semaglutide.

We’ll use a subset of the Centers for Medicare & Medicaid Services (CMS) 2023 Medicare Part D Prescriber Public Use File. Each row represents a prescriber and a specific diabetes drug.

This dataset contains real-world prescribing behavior. Using logistic regression, we’ll explore questions such as:

Class Distributions

Before we can analyze or model anything, we need to clean the data and add few new features (do some feature engineering). One of the first things you will do is create a binary target variable to indicate whether a record involves Semaglutide. This kind of variable only takes on two values: 1 if the drug is Semaglutide, and 0 if it is not. In Python, you can create a binary variable using a Boolean expression and then convert the result to integers using .astype(int).

For example:

This creates a column where any row with Semaglutide is marked as 1, and all others are marked as 0. This will eventually serve as the outcome we try to predict.

Feature Engineering - Creating New Columns

Another important task we will perform is called feature engineering. Feature engineering is the process of transforming raw data into meaningful features that better represent the underlying problem to the predictive model. This includes creating, modifying, or grouping existing variables to enhance model performance or interpretability.

One simple but powerful example is creating a new column that provides more insight than the raw data alone. For instance, to understand how costly each prescription claim is, we can divide Tot_Drug_Cst (total drug cost) by Tot_Clms (total number of claims) to compute a new column called Cost_per_claim:

This gives us a per-claim cost that may be more meaningful for comparison across prescribers or regions.

Another useful technique is mapping existing values into broader categories. For example, individual state abbreviations in the Prscrbr_State_Abrvtn column can be mapped to broader U.S. Census regions. This can help us spot regional trends more easily than looking at each state individually:

By grouping states into regions, we reduce granularity while retaining useful geographic patterns, a valuable preprocessing step before modeling.

However, not all states in the dataset may appear in the dictionary. In those cases, .map() returns a missing value (NaN). To handle this, you can fill in those missing values with a default label like "Missing" using .fillna("Missing").

Value Counts

Finally, it is often useful to explore the distribution of values in a column. For example, how many rows fall into each region or how many are Semaglutide prescriptions. You can do this using .value_counts() for counts and .value_counts(normalize=True) for proportions.

These basic data preparation steps: binary indicators, feature engineering, mapping, and counting are critical in setting up a dataset for modeling. Now, you are ready for exploratory analysis.

Before building the model, we’ll create our target variable, add informative features, and prepare categorical groupings.

Let’s begin by loading the dataset and exploring its first few rows.

We’ll now define our target variable: Semaglutide_drug, which equals 1 when the generic name is Semaglutide and 0 otherwise.

Next, we’ll engineer new columns that reveal more about prescriber behavior. One useful metric is cost per claim:

We can also group prescribers by U.S. region for easier interpretation.

At this stage, we’ve created a structured dataset suitable for exploratory data analysis.

Before we build any models, we need to take a step back and get to know our data. This is the purpose of exploratory data analysis, (EDA). You look for patterns, inconsistencies, and clues that might be helpful to understand how different variables relate to the outcome you’re trying to predict. In our case, that outcome is Semaglutide_drug, a column that tells us whether or not a prescription record involves the diabetes drug Semaglutide.

Each of the questions below helps uncover different aspects of the data structure and variability that may influence our model’s performance or interpretability.

Missing Values

We will start by checking for missing values in the numeric columns. This might seem minor, but it is actually really important. If a feature is missing data for just a few rows, we might be able to fill it in or move on. But if half the values are missing, it could be a sign that the variable is not reliable for modeling.

To calculate the percentage of missing values, we can use the .isna() method to find missing entries, then apply .mean() to calculate the proportion of missing values in each column. Multiplying this by 100 gives us the percentage. You can also use .sum() if you want the raw count.

Group Statistics: Mean and Standard Deviation

Next, we will split the data into two groups: one where Semaglutide was prescribed and one where it was not. Then we will calculate the mean (average) and standard deviation (a measure of spread) for key numeric features in each group. This gives us insight into whether, for example, Semaglutide is associated with higher costs or larger patient counts. These statistics can help us start forming hypotheses, like whether prescribers who use Semaglutide tend to prescribe more 30-day fills or serve more patients.

To do this, we can use .groupby("Semaglutide_drug") followed by .agg(["mean", "std"]) on the numeric columns. This will create a side-by-side comparison table showing the average and variability for each group.

Correlation Between Numeric Variables

We will also explore how the numeric features relate to one another using correlation. Correlation is a measure that ranges from -1 to 1 and tells us how strongly two variables move together. A value close to 1 means they move together in the same direction, while -1 indicates that one tends to decrease as the other increases. A value near 0 suggests no clear relationship.

Heatmap Visualization

To visualize these relationships, we will use a heatmap. A heatmap is a color-coded grid where darker or brighter colors represent stronger relationships. This allows us to quickly see which variables are closely linked and potentially redundant.

To create the heatmap, we first use .corr() on the numeric columns to compute all pairwise correlations. Then we pass that matrix into sns.heatmap(), a Seaborn function that creates the visualization. By setting annot=True, we can print the correlation values directly on the plot, which makes it easier to interpret.

Geographic Patterns

Lastly, we will examine regional prescribing patterns. We want to know: do certain prescribers in certain regions prescribe Semaglutide more often? A good way to check this is with a count plot, which shows how many records come from each region — and whether Semaglutide was prescribed.

Using sns.countplot(), we can plot the number of prescribers in each Prscrbr_State_Region and split the bars by Semaglutide_drug using the hue parameter. This lets us compare across regions in one chart and spot any geographic trends in prescription behavior.

To explore relationships between numeric variables (like cost and total fills), we can use .corr() to compute pairwise correlations and sns.heatmap() to visualize them as a color-coded matrix.

Before modeling, we’ll explore patterns, missing values, and relationships between variables.

Check missing numeric data:

Next, compare group averages to spot differences between prescribers who do and don’t prescribe Semaglutide.

Let’s visualize correlation between numeric variables using a heatmap.

Finally, explore geographic patterns:

In predictive modeling, one of the first steps is to distinguish between predictors (also known as features or independent variables) and response (or target). The predictors are the pieces of information the model will use to make its decisions, while the response is the variable we wish to predict. In this context, we are interested in predicting whether a prescriber issued a prescription for Semaglutide which is a binary outcome that will form the basis of our classification model.

Splitting the Data

Models are not trained on entire datasets. Instead, we partition the data into multiple subsets to serve distinct roles in the model development process. The most common partitioning scheme involves three subsets:

Understanding the Subsets

In supervised learning, our dataset is split into predictors (X) and a target variable (y). We further divide these into training, validation, and test subsets to properly evaluate model performance and prevent overfitting.

Here is what each of these variables means:

These splits are crucial to simulate how the model will perform in real-world settings and ensure that we’re not simply memorizing the data.

Stratified Sampling

One subtle but essential consideration is that we must maintain the distribution of the response variable, particularly in classification settings with imbalanced classes. To achieve this, we use stratified sampling, which ensures that the proportion of cases (e.g., Semaglutide = 1 vs. 0) remains consistent across the training, validation, and test sets. This avoids the model performing poorly simply because the subsets are not represented in the data.

Finally, it is good to inspect each of the resulting subsets. How many observations are in each split? Is the class balance preserved? These simple diagnostics are foundational checks that ensure the integrity of downstream modeling efforts which you will perform in the questions below.

To evaluate the model fairly, we’ll split data into train, validation, and test sets.

Check that class proportions remain balanced across subsets:

Before we can fit our logistic regression model, we need to make sure our dataset is clean and formatted correctly. This stage, called preprocessing, ensures that our features are in a numerical format, have no missing values, are properly scaled, and are aligned across all datasets. Logistic regression, like many models, assumes that the data has been prepared in a certain way. If we skip these steps or do them incorrectly, our model may perform poorly or fail to train altogether.

These question will walk you through five key preprocessing steps, some of which have partially completed code to help guide you.

Handling Missing Values in Categorical Variables

Missing values can cause errors during modeling and interfere with scaling or encoding. For categorical columns like Prscrbr_State_Region and Prscrbr_Type_Grouped, we’ll fill in missing values with the string "Missing". This way, even rows with unknown data are retained and can be captured as their own category during encoding.

For numeric columns like Tot_30day_Fills, Tot_Day_Suply, Cost_per_claim, and Total_Patients, we’ll fill missing values using the median from the training set. This is preferred over the mean because the median is less sensitive to outliers. You’ll use .fillna() to perform this replacement.

For one-hot encoded (binary) columns, we’ll fill missing values with 0. These columns represent the presence or absence of a category, so 0 safely indicates that the feature was not activated for that row.

One-Hot Encoding Categorical Variables

Machine learning models can’t interpret text categories directly. We convert them into numeric form using one-hot encoding, which creates a separate binary column for each unique category. You may hear them as dummy variables, too. For example, the column Prscrbr_State_Region might be transformed into:

We use pd.get_dummies() to apply one-hot encoding. The option drop_first=True tells pandas to omit the first category — this prevents duplicate, which is especially important in models like logistic regression.

Why We Encode Train, Validation, and Test Separately

We always apply one-hot encoding to X_train first. That’s because we want the model to learn from the structure of the training data, including which categories exist. We then apply the same process to X_val and X_test — but here’s the tricky part:

To resolve this, we:

This guarantees that the model sees inputs with the same structure at all stages (training, validation, testing), even if the underlying data varies.

Standardizing Numeric Features

Features that are on very different numeric scales can cause issues for models like logistic regression. For example, Tot_Day_Suply might be in the hundreds while Cost_per_claim could be in the thousands. If we don’t scale them, the model might assign disproportionate importance to the larger features.

To address this, we use StandardScaler() from sklearn.preprocessing. This function subtracts the mean and divides by the standard deviation, resulting in a column with mean 0 and standard deviation 1. We fit the scaler on X_train[numeric_cols], and apply the transformation to X_train, X_val, and X_test.

Converting Boolean Columns

Some features may be stored as True/False. Most models, including logistic regression, expect numeric input. We use .astype(int) to convert all Boolean columns into 1/0 format, which the model can interpret as binary indicators.

Final Structure Check

After all these steps, it’s important to verify that X_train, X_val, and X_test all have the same number and order of columns. This ensures the model receives a consistent input structure during training and evaluation.

Now we’ll clean and standardize features for modeling.

Standardize numeric variables:

Finally, confirm column consistency:

X_train shape: (29091, 26)

X_val shape: (9697, 26)

X_test shape: (9698, 26)

Logistic Regression and the Sigmoid Function

In binary classification problems, our goal is to predict the probability of a binary outcome: such as success/failure or 1/0. Unlike linear regression, which can produce any real number, logistic regression bounds the output between 0 and 1 by applying the sigmoid function. This lets us model probabilities directly using the equation:

$p = 1 / (1 + e^{(-(\beta_0 + \beta_1 \times X)})$

where

The result is an S-shaped curve that flattens near 0 and 1, making it ideal for modeling probabilities.

Log Odds (Logit) Transformation

Modeling probability with a linear equation (like in linear regression) does not work because probabilities must stay between 0 and 1. To make logistic regression behave like linear regression, we apply a transformation to the probability using log-odds, or the logit function:

Three Equivalent Forms of the Logistic Model

Each form is mathematically equivalent, and which one you use depends on the context:

Key Features of the Logistic Curve

How to Interpret Coefficients

In a logistic regression model, each coefficient (β) represents the change in the log-odds of the outcome for a one-unit increase in the predictor, holding all else constant.

Checking Multicollinearity with VIF

Before fitting our model, we use Variance Inflation Factor (VIF) to check for multicollinearity:

VIF(Xᵢ) = 1 / (1 – R²ᵢ)

where ${R_i}^2$ is the $R^2$ from a regression of $X_i$ onto all of other predictors. You can easily see that having ${R_i}^2$ close to one refer to collinearity and so the VIF will be large.

A VIF above 10 suggests the variable is highly collinear and may need to be removed.

Feature Selection with AIC and Forward Selection

To reduce the number of features, we use forward selection guided by Akaike Information Criterion (AIC):

AIC = 2·k – 2·log(L),

where

The model with the lowest AIC fits the data by striking a balance between fit and the number of parameters (features) used. If we pick the model with the smallest AIC, we are choosing the model with a low k (fewer features) while still ensuring it has a high likelihood log(L).

Forward selection begins with no predictors and adds them one at a time, at each step choosing the variable that leads to the greatest reduction in AIC.

Interpreting Model Coefficients with Odds Ratios

Once the model is fit, we convert coefficients into odds ratios to interpret them:

Odds Ratio = exp(β)

Evaluating Model Performance

Confusion Matrix

A confusion matrix compares the model’s predicted classes with the actual outcomes. It is used to calculate accuracy, precision, recall, and more.

From the confusion matrix, we can compute key metrics:

ROC Curve and AUC

A Receiver Operating Characteristic (ROC) curve plots the tradeoff between:

It shows how the model performs across all classification thresholds.

The Area Under the Curve (AUC) is a summary metric:

You should compute and compare AUC scores for:

This helps check for overfitting, which occurs when a model learns the noise or specific quirks of the training data rather than the underlying patterns. An overfitted model may perform very well on the training set but poorly on new, unseen data (test and validation dataset!). By evaluating performance on validation and test sets, we can ensure the model generalizes well to other data.

Ready to Model

Now that you’ve reviewed the key concepts, proceed with training your logistic regression model and interpreting the results using this knowledge!

We’ll now fit a logistic regression model, interpret coefficients, and assess performance.

Now that you’ve trained your final model, let’s use it to predict how likely new prescribers are to prescribe Semaglutide.

The Sigmoid Function and Likelihood in Logistic Regression (Semaglutide Example)

Let’s say we have a prescriber with the following characteristics:

Using the sigmoid function, the model calculates p = 0.83, meaning there’s an 83% chance this provider prescribes Semaglutide.

So:

Maximum Likelihood Estimation (MLE)

To train the model, we use maximum likelihood estimation.

This means we:

Example:

If the model predicts:

Then the total likelihood is:

Likelihood = 0.75 × 0.80 × 0.55

We want to find the coefficients that maximize this product across all rows in the dataset.

At this point, our model is trained and evaluated. Let’s see how it performs on new, unseen prescribers.

We’ll sample 20 new records from the test set and compute their predicted probability of prescribing Semaglutide.

You can now interpret the highest probabilities. Typically, top prescribers share similar characteristics — such as being in Southern or Western regions, having higher cost per claim, or belonging to Primary Care specialties.

This consistent pattern suggests that both geographic and specialty factors influence Semaglutide prescribing behavior.

You’ve completed a full logistic regression pipeline:

This workflow represents a robust and interpretable approach to understanding binary outcomes — a method widely used in both healthcare and business analytics.

By mastering this process, you can now apply logistic regression confidently to new datasets where understanding the probability of an event is key.