# Standard Error of prediction for Logistic Sigmoid function

(previously: Finding the prediction interval for logistic regression)

## Update 2:

This paper describes what I am looking to implement.

## Update:

• Based upon some comments I am renaming this question. I am asking for the std error to calculate prediction bands on a Logistic sigmoid function, which is fit using regression not the regression predictions, which are binary.

• Predictions are being made using scipy.special.expit like so: prediction = [expit(x*beta + alpha) for x in a_pred]. The results look like this:

Based on Kerby Shedden's answer and this I now have to methods of calcualting the SE for the prediction interval:

$$SE = \sqrt(xSx^T)$$

and

$$SE = \sqrt(MSE + xSx^T)$$

However, is it legitimate to caluate the MSE based on residuals from binary data as shown:

## TL;DR

there are two type of SE(Standard error) - More Detail

How do I calculate the SE for a predicted indivdual (not the mean) put into a logistic regression model statsmodels.discrete.discrete_model.Logit

Like this for linear regression in statsmodels or this in R.

I can extract model parameter variance, covariance and std dev but that is it....the model does not return std error of predictions - how do I calculate these?

## Full

I am trying to fit a prediction interval for logitistic regression model. I am using statsmodels although I am happy hear answers using another package.

### My procedure so far:

Fit the model to data df:

log_mdl = statsmodels.discrete.discrete_model.Logit.from_formula ("hit ~ a",df).fit()

The models parameters returned:

                           Logit Regression Results
==============================================================================
Dep. Variable:                    hit   No. Observations:                  200
Model:                          Logit   Df Residuals:                      198
Method:                           MLE   Df Model:                            1
Date:                Tue, 10 Mar 2020   Pseudo R-squ.:                  0.6209
Time:                        14:02:57   Log-Likelihood:                -45.308
converged:                       True   LL-Null:                       -119.52
Covariance Type:            nonrobust   LLR p-value:                 3.823e-34
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -4.0571      0.532     -7.624      0.000      -5.100      -3.014
a              0.0990      0.015      6.657      0.000       0.070       0.128
==============================================================================



then I make some predictions for values a_pred so I can plot the line:

 log_pred = [expit(x* log_mdl.params[1] + log_MLE.params[0]) for x in a_pred]

I now wish to find the prediction interval for each prediction.

To do this I need the SE (Standard Error) of each prediction:

However there are two type of SE - More detail:

Ideally there would be a method like this for linear regression in statsmodels or this in R.

I have found this for 95% Confidence interval for of the true logit, which is the same as the martic operation in this. Which is apparently what R returns as the SE for a prediction.

However I am not sure if this is the SE of the predicted mean (Confidence interval) or the SE of a predicted individual (Prediction interval).

I have also found this in which there are two method of calcualting the CI, one for the funct (i assume mean) and another for an observation (I assume this means a single prediction). This reflects what is said here, that there are two SEs one for the mean and another for a single prediction which includes the variation of the signal not just the accuracy of the estimated mean.

How can I find the SE of a predicted individual?

• The response variable for Logit is binary. The prediction interval is the set {0, 1} unless either one as probability larger than the confidence level, and then it would be just one point. I don't see how prediction intervals make much sense for a binary variable. Commented Mar 10, 2020 at 18:14
• Aplogies, I am making probability predictions like so: expit(x*beta+alpha) where alpha and beta are taken from the logistic model. I have updated the question. Commented Mar 10, 2020 at 19:38
• A subject is either in one category or the other, so a prediction interval with probabilities does not make sense to me. Unless the prediction interval contains the closed interval $[0,1]$, the prediction interval is assured to exclude the correct value.
– Dave
Commented Mar 10, 2020 at 21:07
• This may be of interest. Commented Mar 10, 2020 at 22:26
• Does this answer your question? How are the standard errors computed for the fitted values from a logistic regression? Commented Mar 14, 2020 at 11:38

Something like below should work:

x = log_mdl.model.exog
c = log_mdl.cov_params()
vcov = np.dot(x, np.dot(c, x.T))
se = np.sqrt(np.diag(vcov))


If x is very large there are faster ways to do this that only compute the diagonal of vcov, i.e.

v = (x * np.dot(x, c)).sum(1)
se = np.sqrt(v)


But this is less transparent.

For terminology, I think you could refer to this as the standard error for the logit probabilities.

• I have tried this, SE = sqrt( X C Xt). it yields very close confidence bands - depending on if there is a lot of data - this seems to relate only to the confidence of the predicted mean, not the interval of 95% prediction based upon the spread of the data and the accuracy of the mean. I found an alternative definition: SE = sqrt( MSE + (X C Xt)), where MSE is the mean squared error. This seems like it would account for data variability but I how do I find the MSE for a logistic sigmoid fit? Is it legitimate to find the residuals based upon binary values? Commented Mar 11, 2020 at 9:35
• These confidence bands correctly capture the uncertainty in the probabilities, or more specifically the log odds for P(Y=1 | X=x). When the sample size is large and the number of variables is small or moderate, then these standard errors will be fairly small. This does not capture the predictive uncertainty in a future observation. There is no sensible way to talk about a predive standard deviation for logistic regression becaue what you are predicting is a category label, not a numeric value. Commented Mar 12, 2020 at 11:59

The usual way to assess predictive uncertainty for logistic regression is using cross-validation. If you want something that sort of looks like the OLS prediction interval, you can use something like below:

import numpy as np                                                                                                                                                           import statsmodels.api as sm                                                                                                                                                 import pandas as pd                                                                                                                                                          from scipy.stats.distributions import norm                                                                                                                                                                                                                                                                                                                n, p = 200, 3
# Generate data
x = np.random.normal(size=(n, p))                                                                                                                                            par = np.r_[1, 0, -1]
lpr = np.dot(x, par)
pr = 1 / (1 + np.exp(-lpr))
y = (np.random.uniform(size=n) < pr).astype(np.int)

# Fit a model
model = sm.GLM(y, x, family=sm.families.Binomial())
result = model.fit()

# Calculate the accuracy.
c = result.cov_params()
lpr_hat = np.dot(x, result.params)
pr_hat = 1 / (1 + np.exp(-lpr_hat))
sd = np.sqrt(np.diag(np.dot(x, np.dot(c, x.T))))
q = pr_hat*(1 - norm.cdf(-lpr_hat/sd)) + (1 - pr_hat)*norm.cdf(-lpr_hat/sd)


The bottom five lines are a plug-in estimate of the probability that the observed value and predicted value are the same, for each case.