Why can you not estimate Beta with least squares in logistic regression? [duplicate]

Question

As the title says, why can you not estimate Beta (coefficients) in logistic regression with least-squares? I read in a book that it is not possible and that we use maximum likelihood instead, but I have found no motivation for this anywhere.

Answer 1

Linear regression is, as the name suggests, a linear model. This enables us to use linear algebra to find its parameters, this is called ordinary least squares. We cannot use OLS for generalized linear models like logistic regression, because they are non-linear. GLMs are defined in terms of a linear predictor

$$ \eta = \boldsymbol{X} \beta $$

that is passed through the link function $g$ to obtain the prediction

$$ E(Y\,|\,\boldsymbol{X} ) = \mu = g^{-1}(\eta) $$

Because the link function is non-linear, we cannot use linear algebra to find the parameters, but we need an optimization algorithm. Since generalized linear models are defined in terms of conditional distributions, we can fit them using maximum likelihood. On another hand, you can fit non-linear curves to the data by minimizing squared error, but because of the non-linearity, to do it you would use an optimization algorithm as well.

Answer 2

There are two possibilities:

1. Take the usual $\hat\beta=(X^TX)^{-1}X^Ty$ as your coefficient estimates.

2. Minimize $\sum (y_i-\hat p_i)^2$ square loss instead of $-\sum \bigg[y_i\log(\hat p_i)+(1-y_i)\log(1-\hat p_i) \bigg]$ negative log likelihood loss.

($y_i$ is the $i$-th observed binary $0/1$ outcome; $\hat p_i$ is the $i$-th estimated probability of category $1$.)

That these are not equivalent might not be obvious, but they’re not.

For the former, logistic regression also involves the nonlinear link function, so there is no reason to believe that the resulting coefficients are close to the right ones on the nonlinear logistic regression. I would encourage you to simulate some data; calculate coefficient estimates this way and using traditional logistic regression (maximum likelihood); and compare estimates. Contrary to what the book says, however, nothing will stop you from calculating such a $\hat\beta$ and asserting it as your estimate of the regression coefficients. (At the same time, a vector of all $1$s is an estimate, even if a silly estimate.)

For the latter, that involves the minimization of one of the good scoring-rules for a logistic regression, known as Brier score in this setting. That this is not the preferred method (due to the efficiency of the likelihood-based method) has not stopped people from using square loss. On Yann LeCun’s MNIST page, there is an example that uses square loss.

2-layer NN, 300 hidden units, mean square error

EDIT

I will do an R example demonstrating what happens when you estimate the logistic regression $\beta$ via $\hat\beta=(X^TX)^{-1}X^Ty$. This simulation is based on another Cross Validated answer for how to simulate a logistic regression.

set.seed(2022)
N <- 800
x <- runif(N, -3, 3)
pr <- 1/(1 + exp(-x))
y <- rbinom(N, 1, pr)
L <- glm(y ~ x, family = binomial)

# Now calculate the OLS-style estimate
#
X <- cbind(1, x)
beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y

# Now plot the compared probability predictions
#
plot(
    pr, 
    1/(1 + exp(-predict(L))),
    xlim = c(0, 1),
    ylim = c(0, 1),
    xlab = "True Probability",
    ylab = "Estimated Probabilit",
    col = "black"
)
points(
    pr,
    1/(1 + exp(-(beta_hat[1] + beta_hat[2] * x))),
    col = 'red'
)
abline(0, 1)

When you plot this, you see that the standard logistic regression gives close estimates of the true probability values, while the OLS-style estimates give awful estimates of the true probability values. Consequently, it isn't so much that you can't use the OLS-style $\hat\beta=(X^TX)^{-1}X^Ty$, but that it is a terrible estimation method.

If you prefer ggplot2, you can use this code.

d1 <- data.frame(
    true_probability = pr,
    estimated_probability = 1/(1 + exp(-predict(L))), 
    method = "Standard Logistic Regression"
)
d2 <- data.frame(
    true_probability = pr,
    estimated_probability = 1/(1 + exp(-(beta_hat[1] + beta_hat[2] * x))),
    method = "OLS-style"

)
d3 <- data.frame(
    true_probability = pr,
    estimated_probability = pr,
    method = "Ideal"
)
d <- rbind(d1, d2, d3)
ggplot2::ggplot(d, ggplot2::aes(x = true_probability, y = estimated_probability, col = method)) +
    ggplot2::geom_line() +
    ggplot2::theme_bw()