Decoding Logistic Regression Using Maximum Likelihood Estimation
Logistic regression is one of the most widely used statistical models in machine learning, particularly for binary classification problems. Despite its popularity, many practitioners treat logistic regression as a black box – they understand how to use it but not how it actually works under the hood.
In this post, we‘ll pull back the curtain and take a deep dive into the mathematics of logistic regression, focusing on how the model coefficients are estimated using the statistical principle of maximum likelihood estimation (MLE). By the end, you‘ll have a solid understanding of the theory and assumptions behind logistic regression and be better equipped to use it effectively in your machine learning work.
The Logistic Function
At the core of logistic regression is the logistic function, sometimes called the sigmoid function due to its characteristic S-shape. The logistic function maps any real-valued input to a value between 0 and 1, making it useful for modeling probabilities. It has the following form:
$$ f(z) = \frac{1}{1+e^{-z}} $$
Where $z$ is a linear combination of the feature variables and coefficients:
$$ z = \beta_0 + \beta_1x_1 + \beta_2x_2 + \ldots + \beta_px_p $$
Here, $x_i$ represents the $i$-th feature variable and $\beta_i$ is the corresponding coefficient that we need to estimate from the data.
The logistic function has some nice properties that make it suitable for modeling probabilities:
- It always outputs a value between 0 and 1 for any input $z$
- It approaches 0 as $z$ approaches $-\infty$ and approaches 1 as $z$ approaches $+\infty$
- It equals 0.5 when $z=0$
Interestingly, the logistic function first arose not in statistics but in ecology as a model of population growth [1]. It can be derived from the following differential equation:
$$ \frac{dP}{dt} = rP\left(1-\frac{P}{K}\right) $$
Where $P$ is the population size, $r$ is the growth rate, and $K$ is the carrying capacity. The solution to this equation is the logistic function, where $z$ is a function of time.
In the context of logistic regression, we‘re using the logistic function to model the probability that an observation belongs to the "success" or "positive" class, denoted mathematically as $P(Y=1|X)$. By taking the log of both sides of the logistic function, we can express the log-odds of success as a linear function of the features:
$$ \log\left(\frac{P(Y=1|X)}{1-P(Y=1|X)}\right) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \ldots + \beta_px_p $$
This is where the coefficients get their interpretation as the change in log-odds of success per unit change in the corresponding feature variable, holding all other features constant.
Why Not Linear Regression?
A natural question is why we need logistic regression at all – why not just use linear regression to model probabilities? The key issue is that linear regression assumes the output is unbounded and can take on any real value, positive or negative. This doesn‘t work for probabilities, which must be between 0 and 1.
If we used linear regression to model a binary outcome, we might get impossible predicted probabilities like -0.3 or 1.5. The logistic function squashes the output of linear regression to the [0, 1] range, ensuring we get valid probabilities.
Logistic regression also doesn‘t assume that the errors (residuals) are normally distributed, which is a key assumption of linear regression [2]. The errors in logistic regression follow a binomial distribution, which is more appropriate for binary data.
Maximum Likelihood Estimation
Now that we understand the form of logistic regression, how do we actually estimate the coefficients $\beta_i$ from data? The most common approach is maximum likelihood estimation (MLE).
MLE is a general statistical principle for estimating parameters of a model given observed data. The basic idea is to choose the parameter values that maximize the likelihood function, which quantifies how likely the observed data are given the model parameters.
Mathematically, for data $\mathbf{X} = (X_1, \ldots, X_n)$ and parameters $\theta$, the likelihood function is defined as:
$$ L(\theta|\mathbf{X}) = P(\mathbf{X}|\theta) = \prod_{i=1}^n P(X_i | \theta) $$
The second equality holds if we assume the observations are independent. We often work with the log-likelihood $\ell(\theta|\mathbf{X}) = \log L(\theta|\mathbf{X})$ instead of the likelihood itself since it‘s easier to maximize (products become sums).
For logistic regression, the log-likelihood takes the following form:
$$
\begin{aligned}
\ell(\beta) &= \sum_{i=1}^n \left[ y_i \log(p(x_i)) + (1-y_i)\log(1-p(xi)) \right] \
&= \sum{i=1}^n \left[ y_i (\beta_0 + \beta1x{i1} + \ldots + \betapx{ip}) – \log(1+e^{\beta_0 + \beta1x{i1} + \ldots + \betapx{ip}}) \right]
\end{aligned}
$$
Where $p(x_i) = P(Y=1|X=x_i)$ is the probability of success for the $i$-th observation, $yi$ is the observed class label (0 or 1), and $x{ij}$ is the value of the $j$-th feature for the $i$-th observation.
To find the maximum likelihood estimates of the coefficients, we need to maximize this log-likelihood function with respect to $\beta$. This requires setting the gradient vector of first partial derivatives to zero:
$$
\begin{aligned}
\frac{\partial \ell}{\partial \beta0} &= \sum{i=1}^n (y_i – p(x_i)) = 0 \
\frac{\partial \ell}{\partial \betaj} &= \sum{i=1}^n x_{ij}(y_i – p(x_i)) = 0, \quad j=1,\ldots,p
\end{aligned}
$$
These equations are nonlinear in $\beta$ and don‘t have a closed-form solution. So we have to use an iterative numerical optimization algorithm like Newton-Raphson or Fisher scoring to solve them.
Most statistical software packages have reliable implementations of these algorithms for fitting logistic regression models. Here‘s how we would fit a model in R using the glm() function:
model <- glm(y ~ x1 + x2 + x3, family = binomial(), data = mydata)The family = binomial() argument specifies that we want to fit a logistic regression model, and the data argument passes our data frame of observations. The estimated coefficients are stored in model$coefficients.
Example: Titanic Survival
To illustrate the MLE logistic regression workflow with a real example, let‘s analyze the classic Titanic survival data set. This data set contains information on passengers aboard the Titanic, including whether they survived, their age, sex, passenger class, and more.
We‘ll model the probability of survival based on two features: sex (male or female) and passenger class (1st, 2nd, or 3rd). Here‘s a preview of the data:
| Survived | Sex | Pclass |
|---|---|---|
| 0 | male | 3 |
| 1 | female | 1 |
| 1 | female | 3 |
| 1 | female | 1 |
| 0 | male | 3 |
We can fit the logistic regression model in R as follows:
titanic_model <- glm(Survived ~ Sex + Pclass, family = binomial(), data = titanic_data)Here are the estimated coefficients:
| Estimate | Std. Error | z value | P(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 3.5210 | 0.3537 | 9.950 | < 0.0001 |
| Sexmale | -2.5137 | 0.1602 | -15.695 | < 0.0001 |
| Pclass2 | -1.0225 | 0.2140 | -4.779 | < 0.0001 |
| Pclass3 | -2.2617 | 0.2121 | -10.665 | < 0.0001 |
To interpret these coefficients, we need to think in terms of log-odds. The intercept of 3.5210 means that for a female passenger in 1st class (the reference categories), the log-odds of survival are 3.5210. Exponentiating this gives an odds ratio of $\exp(3.5210) = 33.8$, or equivalently, a probability of survival of $33.8/(1+33.8) = 0.97$.
For male passengers, the log-odds of survival are lower by 2.5137 compared to female passengers, holding passenger class constant. In terms of an odds ratio, this is a multiplicative effect of $\exp(-2.5137) = 0.08$. So the odds of survival for males are 0.08 times the odds for females, i.e. much lower.
Similarly, the coefficients for Pclass2 and Pclass3 show that the log-odds of survival decrease as we move from 1st to 2nd to 3rd class, compared to the reference level of 1st class.
We can assess the statistical significance of the coefficients using p-values from the Wald test, shown in the last column. The small p-values indicate that all the coefficients are highly significant.
Another way to assess the overall fit of the model is the likelihood ratio test, which compares the full model to a reduced model with no predictors. In R:
test_model <- glm(Survived ~ 1, family = binomial(), data = titanic_data)
anova(test_model, titanic_model, test="Chisq")This gives a chi-square statistic of 346.21 with 3 degrees of freedom and a p-value < 0.0001, indicating that our full model with Sex and Pclass fits significantly better than a null model predicting the overall mean survival rate.
Limitations and Extensions
While logistic regression is a powerful and widely used model, it‘s important to be aware of its assumptions and limitations:
-
It assumes a linear relationship between the log-odds of the response and the predictor variables. If this assumption is violated, the model may fit poorly. Polynomial terms or splines can help model nonlinear relationships.
-
It assumes no multicollinearity (high correlation) among the predictors. Multicollinearity can cause unstable estimates and large standard errors.
-
It assumes independence of observations. If observations are correlated (e.g. repeated measures data), we need to use a method that accounts for this, such as generalized estimating equations (GEE).
-
It can suffer from complete or quasi-complete separation, where the outcome can be perfectly or nearly perfectly predicted by a predictor. This causes the MLE to blow up or fail to converge. Penalization methods like ridge regression can help mitigate this.
-
It assumes that binary outcomes come from a binomial distribution. If the data are overdispersed (more variable than expected), a quasi-likelihood approach may be needed.
Several extensions to logistic regression aim to address these limitations:
-
Firth‘s penalized likelihood method shrinks the MLE coefficients to avoid separation issues [3].
-
Lasso and ridge regression add a penalty term to the likelihood function to perform variable selection and shrink the coefficients [4].
-
Quasi-binomial regression uses a dispersion parameter to account for overdispersion in binary data [5].
Understanding the math behind logistic regression and the assumptions involved helps you recognize when these extensions might be needed and interpret their results correctly.
Conclusion
In this post, we took a deep dive into the theory and mathematics behind logistic regression, with a focus on how the coefficients are estimated using maximum likelihood. We derived the logistic regression log-likelihood function and showed how Newton‘s method can be used to find the coefficients that maximize it. We also worked through a real data example predicting Titanic passenger survival.
While modern software makes it easy to fit logistic regression models without worrying about the details, understanding the underlying statistical theory is still valuable. It allows you to diagnose when the model‘s assumptions are violated, recognize when extensions may be needed, and interpret the model‘s coefficients and output more insightfully.
Logistic regression is just one of many models that rely on maximum likelihood estimation. Similar principles apply to Poisson regression, Cox proportional hazards models, mixed effects models, and more. Once you understand MLE, you‘ll start to see it everywhere.
I hope this post helped illuminate what‘s happening under the hood of logistic regression. The coefficients may seem mysterious at first, but by deriving them using maximum likelihood, we see that they arise from a principled statistical approach to learning from data. Far from a black box, logistic regression is built on elegant mathematics and a powerful statistical theory of estimation.
