Making linear to logistic regression with sigmoid function - why is a transformation of predicted y needed?

Question

I noticed that one can run a linear regression for binary outcomes and get the same predictions as from a logistic regression after using a sigmoid function. That is what I awaited. But the surprising part for me is that I needed to transform the predicted values in advance. Here is what I did:

1. Run a logistic regression for y ~ x, where y is binary and x a continuous random variable.
2. Made predictions of logits (y_logit) over the range of x with the model I received in step 1.
3. Applied inverse logit to transform the y_logit values from step 2 to probabilities (y_prob): y_prob <- 1 / (1 + exp(-y_logit)).

Afterwards I started with the OLS:

4. Run a OLS linear regression for y ~ x.
5. Made predictions of y over the range of x with the model coefficients I received in step 4: y_pred <- b0 + b1* x.
6. Run an OLS where I regressed the predicted logits y_logit from step 2 on the predicted y (y_pred) from step 5: y_logit ~ y_pred. This was a perfect linear fit.
7. Transformed predicted y_pred values from step 5 with the coefficients from model of step 6: y_transformed <- b0' + b1'*y_pred.
8. Applied a sigmoid function on y_transformed from step 7: y_prob_lm <- 1/ (1 + e^(-y_transformed)).
9. Compared whether probabilities y_prob from step 3 received by logistic regression are same as the y_prob_lm probabilities recieved by OLS from step 8: Yes, they are the same. y_prob = y_prob_lm.

I thought I could jump from step 5 to step 8, i.e. I thought I could apply the sigmoid function directly on the predicted y values of the OLS (y_pred). But I had to transform them in advance (steps 6 and 7). Why do we need to transform predicted y values before applying the sigmoid function when predicting probabilities with OLS?

---

Whole story with R. Before getting to steps 1 to 9 we need to generate data:

library(magrittr)
# Generate data
# Make it reproducible
set.seed(1)
# Number of observations.
n <- 100000
# Random x variable.
x <- rnorm(n, 50, 10)
# Coefficients.
b0 <- log(.001)
b_x <- log(1.2)
b_noise <- log(.9)
# Formula.
z <- b0 + b_x* x + b_noise* rnorm(n, 0, 1)
# y.
y <- {1/ (1 + exp(-z))} %>%
  rbinom(n= n, size= 1, prob= .) %>%
  as.logical()

# Make it a dataframe.
df <- data.frame(x, y, y_num= as.numeric(y) - 1) 

Here are the steps 1 to 3 with the logistic regression:

# Fit model (STEP 1).
mod_log <- glm(formula= y ~ x, family= binomial(link="logit"), data= df)

# What values to predict for.
df_pred <- data.frame(x= seq(0, max(x), .1))

# Predict.
# logit predictions (STEP 2)
df_pred$log_logit <- predict(mod_log, newdata= df_pred)
# Apply inverse logit to transform to probabilities (STEP 3)
df_pred$log_prob <- 1 / (1 + exp(-df_pred$log_logit))

plot(df_pred$x, df_pred$log_prob)
# This results in a plausible looking plot.

Steps 4 to 9:

# Run linear model (STEP 4)
mod_lm <- lm(formula= y_num ~ x, data= df)
# Predict y with linear model fit (STEP 5)
df_pred$lm_y <- coef(mod_lm)[1] + coef(mod_lm)[2] *df_pred$x
# Transform y.
df_pred$lm_y_transformed <- data.frame(lm_y= df_pred$lm_y, log_logit=  df_pred$log_logit) %>%
  lm(log_logit ~ lm_y, .) %>% # (STEP 6)
  coef() %>%
  {.[1] + .[2] *df_pred$lm_y} # (STEP 7)
# Apply sigmoid function on y (STEP 8)
df_pred$lm_prob <- 1/(1 + exp(-df_pred$lm_y_transformed))
# Test whether the probabilities are the same (STEP 9).
all.equal(df_pred$log_prob, df_pred$lm_prob)
# TRUE!

Answer 1

Spending time doing ad hoc approaches when logistic regression works just fine is a bit hard to understand. More importantly, your approach breaks down when an estimated probability is 0 or 1 and you attempt to compute its logit for later steps, as you’ll get infinities.