Pointwise standard errors for a logistic regression fit with statsmodels

Question

Source

An Introduction to Statistical Learning with Applications in R, found here: https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf

Task

I'm trying to replicate the example of a polynomial logistic regression on the "Wage" dataset on page 267/8.

Theory outline

According to the book, once predictions have been made, confidence intervals can be calculated like so. For a model of the form $$\hat{f}(x_0)=\hat{\beta_0}+\hat{\beta_1}x_0+\hat{\beta_2}x_0^2+\hat{\beta_3}x_0^3+\hat{\beta_4}x_0^4,$$ with a $5\times 5$ covariance matrix $C$ and vector $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$, the pointwise standard error is the square root of $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. So for every $x_0$ in our dataset we have a plot of predictions $\hat{f}(x_0)$ and a plot of upper and lower confidence intervals $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.

For a logistic regression, the same principal can be applied, but the confidence is around the conditional probability logit function, as opposed to the predictions that come straight from the formula above.

Data and approach / reusable code

First of all, this is the code for generating the logistic regression model and plotting the results. This bit is fine and I've successfully reproduced what is in the book:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from statsmodels.datasets import get_rdataset
from statsmodels.discrete import discrete_model

from sklearn.preprocessing import PolynomialFeatures
polynomial_feat = PolynomialFeatures(degree=4)

# Get dataset from the R package
data = get_rdataset("Wage", package="ISLR") 
df = data.data.reset_index() 

# Split data into wage (response, y) and age (predictor, X_orig)
y = df.wage
X_orig = df.filter(['age'], axis=1)
# Get the polynomial features from the predictor variable
X = polynomial_feat.fit_transform(X_orig) 

# Set up the test ages for a smooth results plot
X_test = np.linspace(18, 80, 1000)
X_test = X_test[:,np.newaxis] 
X_test_poly = polynomial_feat.fit_transform(X_test) 

# Create a dummy response variable, 1 if wage > 250k and 0 otherwise
y_dummy = pd.DataFrame({'wage': y[:]}) 
y_dummy['wage_split'] = np.where(y_dummy['wage'] > 250, 1, 0) 
y_dummy = y_dummy.drop(['wage'], axis=1)

# Fit a logistic regression model with statsmodels
logit_model = discrete_model.Logit(y_dummy, X).fit() 
# Get predictions, i.e. Pr(Wage > 250 | Age)
y_preds = logit_model.predict(X_test_poly)

# Plot the results
plt.figure(figsize=(8, 8)) 
plt.plot(X_test, y_preds, 'b-') 
plt.ylim(top=0.2) 
plt.xlabel("Age")
plt.ylabel("P(Wage > 250 | Age)")
plt.title("Probability of Earning > 250k with Logistic Regression")

So now I attempt to plot the confidence intervals. I don't think there is a method to do this directly in statsmodels (please correct me if I'm wrong).

My issue

My issue here is in the calculation of the pointwise standard errors and the confidence intervals. We know that the response values for the logistic regression model must be $y\in [0, 1]$, since it is a conditional probability.

The problem is that for every $x_0$, the value of $$\sqrt{l_0^TCl_0}$$ is going to be relatively large. I can demonstrate this by using the first age value, $x_0=18$:

# Get the covariance matrix from the model class
C = logit_model.normalized_cov_params
x = 18.
L_T = np.array([1, x, x**2, x**3, x**4])

# Compute the pointwise standard error, as outlined above
L_T = np.matrix(L_T)
L = np.transpose(L_T)
C = np.matrix(C)

var_f = np.matmul(np.matmul(L_T, C), L)
var_f = np.asarray(var_f)[0][0]
pointwise_se = np.sqrt(var_f) 
print(pointwise_se)

The output of this is pointwise_se = 6.14.

From the plot above, I can see that the prediction of $\text{Pr}(\text{Wage} > 250 | x=18)$ is close to zero, and from the example provided in the book I can see that the confidence interval around this value is not wide, and definitely doesn't go negative or greater than 1.

If I was to get a confidence interval from a pointwise standard error of $6.14$, the plot would be silly, and not a replication of that in the book.

My question

What am I doing wrong in my calculation of the pointwise standard error?

Answer 1

Since you are doing logistic regression and not simple linear regression, the equation $\hat f(x_0)=\hat\beta_0+\hat\beta_1x_0+\hat\beta_2x_0^2+\hat\beta_3x_0^3+\hat\beta_4x_0^4$ does not refer to the probability of earning >250K, but to the logit of that probability. This is the same as saying that logistic regression is a linear model that uses logit as a link function.

So, you have to define functions to convert between probabilities and logits (maybe they are already implemented in Numpy or something, but they are simple enough to type):

def logit(p):
    return np.log(p/(1-p))

def invlogit(x):
    # inverse function of logit
    return 1/(1+np.exp(-x))

Now, we have to apply the pointwise SE you calculated to the logit of the point estimates, and then convert back to probabilities:

upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)

Where std_err is an array with the standard errors of $\hat f(x)$ that you correctly calculated. Then, upper_limit and lower_limit will give an interval around the estimated probability.