Show that classification tables do not always correlate with goodness of fit for logistic regression

Question

Background

I am reading the textbook Applied Logistic Regression by David Hosmer, specifically chapter 4, which discusses logistic regression model assesment of fit. Hosmer gives an interesting example of how classification tables may not imply goodness of fit for a logistic regression model.

Hosmer introduces univariate logistic regression as follows:

Given a dichotomous random variable $Y\in\{0,1\}$, we model the expected value of $Y$ given $X$ using a logistic function.

$$E(Y|X) = \pi(x) = \frac{1}{1+e^{-(\beta_0 + \beta_1X)}}$$

and he goes on to define the logit and derive the likelihood function, but I don't think those are needed for the problem below.

Problem

Suppose that $P(Y=1) = \theta_1$ and $X \sim N(0,1)$ in the group with $Y=0$ and $X \sim N(\mu,1)$ in the group with $Y=1$.

(1) It can be shown that the slope coefficient for the logistic regression model is $\beta_1 = \mu$, and the intercept is

$$ \beta_0 = \ln\left[\frac{\theta_1}{1-\theta_1}\right] - \frac{\mu^2}{2} $$

(2) Then, the probability of misclassification (PMC), may be shown to be

$$ \mathrm{PMC} = \theta_1 \Phi\left[\frac{1}{\beta_1}\ln\left(\frac{1-\theta_1}{\theta_1}\right)-\frac{\beta_1}{2}\right] + (1-\theta_1)\Phi\left[\frac{1}{\beta_1}\ln\left(\frac{\theta_1}{1-\theta_1}\right)-\frac{\beta_1}{2}\right] $$

where $\Phi$ is the standard normal CDF. This implies that the misclassification rate depends on the slope of the model rather than the goodness of fit.

Discussion

Hosmer says that (1) is easy to show, and is shown in an old textbook named Discriminant Analysis by Lachenbruch. I was unable to find a PDF of this textbook. I think to show (1), I should use conditional probability somehow, but I'm not sure. I assume that (2) will follow using similar techniques.

How can I show (1) and (2)?

Answer 1

This notation is not perfect, but it helps keep the conditioning straight:

$$E(Y|X=x) = \frac{1}{1+e^{-(\beta_0 + \beta_1 x)}}$$

Since $X \sim N(0,1)$ and $N(\mu, 1)$,

$$f(X|Y=1) = \frac{1}{\sqrt{2 \pi} \sigma} e^{\frac{-(x-\mu)^2}{2 \sigma}} = \frac{1}{\sqrt{2 \pi}} e^{\frac{-(x-\mu)^2}{2}}$$

$$f(X|Y=0) = \frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}}$$

Recognize that

$$f(X) = P(Y=0)f(X|Y=0) + P(Y=1)f(X|Y=1)$$

$$E(Y|X=x) = 0 P(Y = 0|X=x) + 1 P(Y=1|X=x) = P(Y=1|X=x)$$

And using Bayes rule

$$ = \frac{f(X|Y=1) P(Y=1)}{f(X)} = \frac{\theta \frac{1}{\sqrt{2 \pi}} e^{\frac{-(x-\mu)^2}{2}}}{\theta \frac{1}{\sqrt{2 \pi}} e^{\frac{-(x-\mu)^2}{2}} + (1-\theta) \frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}}}$$

$$ = \frac{1}{1 + \frac{1-\theta}{\theta} e^{-x^2 / 2 + (x-\mu)^2/2}} = \frac{1}{1+e^{ln(\frac{1-\theta}{\theta}) + \mu^2/2 - \mu x}}$$

Therefore $\beta_1 = \mu$ and $\beta_0 = log(\frac{\theta}{1-\theta}) - \mu^2 / 2$

Next the probability of misclassification is when a different class would be predicted than the actual. Utilizing 0.5 as the class threshold,

$$PMC = P(Y=1) P(\frac{1}{1+e^{-(\beta_0 + \beta_1 x)}} < 0.5 | Y=1) + P(Y=0)P(\frac{1}{1+e^{-(\beta_0 + \beta_1 x)}} > 0.5 | Y = 0)$$

Note that $\frac{1}{1+e^{-(\beta_0 + \beta_1 x)}} = 1/2$ then algebra shows $x = \frac{-\beta_0}{\beta_1}$

Now we can write the PMC:

$$PMC = P(Y=1)P(X < \frac{1}{\beta_1} log(\frac{1-\theta}{\theta}) + \beta_1 / 2 < 0.5|Y=1) + P(Y=0)P(X < \frac{1}{\beta_1} log(\frac{1-\theta}{\theta}) + \beta_1 / 2 > 0.5 | Y = 0)$$

Now, change the $P(X < *)$ to be a statement for a normal CDF

$$ = \theta \Phi(\frac{1}{\beta_1} log(\frac{1-\theta}{\theta}) + \beta_1 / 2 - \mu) + (1-\theta)[1-\Phi(\frac{1}{\beta_1} log(\frac{1-\theta}{\theta}) + \beta_1 / 2)]$$

$$ = \theta \Phi(\frac{1}{\beta_1} log(\frac{1-\theta}{\theta}) - \beta_1 / 2) + (1-\theta)[1-\Phi(\frac{1}{\beta_1} log(\frac{1-\theta}{\theta}) + \beta_1 / 2)]$$

Answer 2

Many thanks to R Carnell for providing the derivation. To increase my understanding, I ran some simulations to validate the theoretical results.

First some helper functions.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import bernoulli

# expectation of Y given X -- logistic model
def EYX(x,beta0,beta1):
    denom = 1+np.exp(-(beta0 + beta1*x))
    return 1./denom

# predict class of Y given X, using logistic model
def classify(x,beta0,beta1,c):
    return np.asarray([int(EYX(xval,beta0,beta1)>c) for xval in x])

Now simulate 100 observations of the dichotomous variable $Y$ such that $P(Y=1)=\theta_1$ and $X\sim N(\mu,1)$ in the group with $Y=1$ and $X\sim N(0,1)$ in the group with $Y=0$.

# parameters
mu=1
theta=0.7

# classification cutoff
c=0.5

# generate the data using mu,theta
n=100
y=np.sort(bernoulli.rvs(theta,size=n))
x0=norm.rvs(loc= 0,scale=1,size=np.sum(y==0))
x1=norm.rvs(loc=mu,scale=1,size=np.sum(y==1))
x=np.concatenate((x0,x1))

Now we calculate the theoretical model coefficients, predict the class of each $Y_i$

# true beta0,beta1
beta0=np.log(theta/(1-theta)) - 0.5*mu**2
beta1=mu

# predict class using logistic model
yhat=classify(x,beta0,beta1,c)

Below, I plot the data, the logistic model, the predicted classes, and the undelrying distributions for $X|Y=0$ and $X|Y=1$.

plt.plot(x,y,'o',label='data');
xv = np.linspace(-3,3,100);
plt.plot(xv,(1-theta)*norm.pdf(xv),'k--',label='f(X|Y=0)');
plt.plot(xv,theta*norm.pdf(xv,loc=mu),'k--',label='f(X|Y=1)');
plt.plot(xv,EYX(xv,beta0,beta1),label='E(Y|X)');
plt.plot(xv,classify(xv,beta0,beta1,c),label='Predicted Class')
plt.legend();

We can also compute the theoretical misclassification rate.

def MCR(x,y,beta0,beta1,c):
    yhat = classify(x,beta0,beta1,c)
    num_misclass = float(np.sum(np.abs(yhat-y)>0))
    return num_misclass/float(y.size)
    
true_PMC = theta*norm.cdf(np.log((1-theta)/theta)/beta1 - beta1/2.) \
        + (1-theta)*norm.cdf(np.log(theta/(1-theta))/beta1 - beta1/2.)

obsv_PMC = MCR(x,y,beta0,beta1,c)

print(true_PMC,obsv_PMC)

which yields

0.2530043786236347 0.21

Does the observed misclassification rate match the theoretical value, on average? A quick simulation confirms.

ntimes=2000
obsv_PMC=np.zeros(ntimes)

for i in range(ntimes):
    y=np.sort(bernoulli.rvs(theta,size=n))
    x0=norm.rvs(loc= 0,scale=1,size=np.sum(y==0))
    x1=norm.rvs(loc=mu,scale=1,size=np.sum(y==1))
    x=np.concatenate((x0,x1))
    obsv_PMC[i] = MCR(x,y,beta0,beta1,c)

plt.hist(obsv_PMC,bins=100,alpha=0.7);
plt.axvline(x=true_PMC,
            linestyle='-',
            color='r',
            label='True PMC');
plt.legend();
plt.title('Misclassification Rate');

The only confusing aspect is exactly how this relates to the goodness of fit. How does this result imply that there exist cases where PMC goes to zero but goodness of fit does not increase?