# How can I assess case-level uncertainty of classification using logistic regression?

I'm hoping to fit a binary logistic regression to be used to predict the binary outcome for new cases/observations. I'm wondering if there is any way to gauge uncertainty of a prediction for individual cases - for example, if I use my model to predict the binary outcome in 500 cases (which were not used to build the model), for each case, is it possible to get an idea of how certain the predicted outcome is? Since the outcome of logistic regression is a probability between 0 and 1, with a probability < 0.5 typically corresponding to the binary outcome 0 and a probability > 0.5 corresponding to the binary outcome 1, is it logical to assume that how close the predicted probability for an individual case is to 1 or 0 corresponds to the uncertainty of the estimate for that case? If not, is there another way to achieve this, or a classification method other than logistic regression that might be more suitable?

Hope this makes sense, thanks!

• Probability < 0.5 does not correspond to a binary outcome of 1. That is a forced choice that is inconsistent with the intentions of logistic regression as a direct probability estimation model. See this for more. Just get confidence intervals for individual risk predictions. Dec 17, 2023 at 12:57

People usually distinguish two types of uncertainty: aleatoric uncertainty and epistemic uncertainty.

Aleatoric uncertainty describes the uncertainty due to pure randomness in the data-generating process, e.g., the noise term $$\varepsilon$$ in $$y=f(x)+\varepsilon$$. You will never be able to predict this noise (because it is random), thus even a perfect model that has correctly learned the true $$f(x)$$ will display random errors corresponding to $$\varepsilon$$ vs. the true outcome $$y$$. In a logistic regression context, quantifying aleatoric uncertainty is fundamentally impossible.

Epistemic uncertainty describes the uncertainty in your fitted model parameters $$\theta$$ that parametrize the model $$f(x)$$. In other words, how precisely do your data determine your fitted model? In a logistic regression context, this is what you can get from the identified covariances / confidence intervals of the parameter estimates.

Now, you might think that what you want is a prediction interval (which is distinctly different from a confidence interval) or a prediction set: an interval that contains the true outcome value for this individual data point with a high probability. To get this, one needs to combine aleatoric and epistemic uncertainty estimates. However, prediction intervals/sets are, again, rather useless for logistic regression because the aleatoric uncertainty component cannot be quantified.

Finally, the most practical answer to your question is: yes, the number between 0 and 1 that you get is a reasonable measure of model uncertainty if your model is calibrated, which is generally more likely than not for logistic regression (but you should still check it). Read up on the meaning of calibration; this is most likely what you want to look into.

Here I'm calculating standard deviation of each prediction & then take adjust prediction by positive and negative standard deviation.

import numpy as np
import statsmodels.api as sm

# Simulating some data for logistic regression
np.random.seed(42)  # for reproducibility
X = np.random.randn(100, 3)  # 100 observations and 3 features
y = np.random.binomial(1, p=1/(1 + np.exp(-np.dot(X, np.array([0.5, -0.25, 0.75])))))  # generating binary outcomes

# Adding a constant term for the intercept

# Fitting a logistic regression model
logit_model = sm.Logit(y, X).fit()

# Now, let's make predictions with the model
# Here, we'll use the same X for simplicity, but in practice, you should use new, unseen data
predictions = logit_model.predict(X)

# Compute standard errors of predictions
std_errors = np.sqrt(np.diag(np.dot(np.dot(X, logit_model.cov_params()), X.T)))

# Compute 95% confidence intervals
conf_int_lower = predictions - 1.96 * std_errors
conf_int_upper = predictions + 1.96 * std_errors

# Output the first 5 predictions and their confidence intervals
results = np.column_stack((predictions[:5], conf_int_lower[:5], conf_int_upper[:5]))
print(results)