Mastering Logistic Regression: Interpretation and Fitting

The Logistic Function

The magic of this method starts with the shape of the function itself, with a specific parameterization¹ for regression purposes. Since this function is the core of the method, we’ll spend some time developing some intuition for how this function behaves.

Now with Probability!

Having a function that maps values between 0 and 1 sets the stage perfectly for introducing probability, which will also take values between 0 and 1! To make this into a concrete example, suppose I’m interested in modeling how many fruit snacks one can eat (explanatory variable, $x$) before getting a stomach ache (response variable, $y$). It’s reasonable to think about the probability of such an event since I could probably have 100 fruit snacks one day and get a stomach ache and do the same another day and not get a stomach ache. It is a little ambitious to think that there are no other variables that might affect this probability, and that the entire population will have the same tolerance for fruit snacks. For example, if kids can tolerate more fruit snacks than adults, we’d likely need to model more covariates, but we’ll just keep it simple for now.

Introducing the mathematical description, we’re modeling a single data point as a dichotomous (Bernoulli) random variable $Y$ with parameter $p$. The value of $Y$ can either be $Y=1$ for stomach ache, and $Y=0$ for not stomach ache.² $p$ is the probability that we observe a stomach ache. In this case, $p$ will depend on the number of fruit snacks that someone eats, so we should instead write $p(x)$, where $x$ is the number of fruit snacks. The probabilities are related to $x$ by the logistic function as $p(x)$ is defined above. If we’re making multiple observations of $Y_i$ given some values of $x_i$, we’ll need to introduce a subscript $i$ that indexes all the observations and write $p(x_i)$ since the probability is different for each observation. Finally, the coefficients for the logistic regression are $\beta_0$ and $\beta_1$, that (for now) we assume are given to us. We can include the dependence explicitly in $p(x_i; \beta_0, \beta_1)$, using a ”;” to separate our observations from our parameters. In symbols, this is all summarized by the statement

If you’re curious, I read this, “Y sub i conditional on x sub i is distributed as a Bernoulli random variable with parameter p of x sub i given beta 0 and beta 1”. Often we’ll just remember that $p$ depends on our parameters and data, so we can shorten the notation to just $p_i = p(x_i; \beta_0, \beta_1)$. I’ll shorten this to $Y_i \sim \text{Bernoulli}(p_i)$ since we’re thinking of $x_i$ as fixed values as well, sacrificing some precision for clarity.

This is our data model. We have now attached a probabilistic nature to the data events that we observe. Our data model provides a description of the probability of observing each event, known as a probability mass function (pmf). For the bernoulli distribution, $\text{Bernoulli}(p_i)$, the pmf is given by:

This is a notationally heavy way to say the probability of getting a stomach ache ($Y_i = 1$) is $p_i$ and not getting one is 1 minus that. The precision of such a statement takes some getting used to, but trace the function through to make sure you can map the english back to the mathematical statement.

• $Y_i = y_i$ is a common point of confusion. Just think that $Y_i$ is before we observe the data, so we can talk about the probability of such a random variable. $y_i$ is the actual value, and has no probability associated with it, you either had the stomach ache or not. For interpretation, you can always think about “plugging” in the value of the lower case variable, and we can read $Pr(Y_i = 1)$ as “probability of getting a stomach ache”.
• $Pr(Y_i = y_i ; p_i)$ again separates data from parameters, and sometimes the parameters of the pmf are left off. We need to be given our parameter $p_i$ (which in turn is a function of number of fruit snacks and $\beta_0, \beta_1$), in order to state the probability of observing an event $Y_i$.

Interpretation of Coefficients

With $p(x)$ understood as the probability of an event, we can determine the meaning of $\beta_0, \beta_1$ on the probability of an event. We need to invert the logistic function to get:

The term appearing on the right hand side has a name — we call it the log odds. The term “odds” is the same as you would use colloquially, as in, if you had a $1/3$ chance of winning a game, your odds are 1:2. If you plug in $p(x) = 1/3$ into the formula above, the log odds on the right hand side is $\log(\frac{1}{2})$. Thus, a 1 unit increase in $\beta_0$ implies that there is a 1 unit increase of the log odds of the event happening. Similarly, a 1 unit increase in $x$ implies that there is a $\beta_1$ unit increase on the scale of log odds. If you’re a “show me the math” type person, suppose we have two observations of someone eating 55 fruit snacks, and someone else eating 56 fruit snacks, $x_1 = 55, x_2 = 56$ (a one unit difference in $x$). If we take the difference of these two scenarios,

We say that $\beta_1$ is the slope on the log odds scale, and $\beta_0$ is the intercept on the log odds scale when $x = 0$. In our scenario of $x_1 = 55, x_2 = 56$ with 1 fruit snack difference (where the denominator is zero), $\beta_1$ is directly interpreted as the difference of log odds.

To be frank, interpreting the coefficients on a log odds scale is quite futile — as shown above in the logistic regression function, changes on the log odds scale depends on where you are on the scale, with larger changes in probability when you are close to 0. The most I interpret, is if the log odds are positive, this means the probability is greater than half.

Even exponentiating both sides to get something on the “odds” scale is quite dubious. I suspect the reason odds are used so often in the gambling industry is because they prey on the fact that odds are fraught with misinterpretations. Seeing an odds of 1:2, you’ll likely jump to thinking about a probability of 1/2, yet there’s quite a large difference between winning probabilities of 1/3 and 1/2. Inversely, when talking about “payouts”, for example on a craps table, you’ll commonly see that an Any 7 bet pays “5 for 1”, whose true meaning is a 4:1 payout, you get 4$ for every 1$ you put down. When it comes to taking your money, the casino will make it seem it’s less likely you’ll lose money (odds), but for paying you out, they’ll make it seem that you get more money (proportion)! Mixing the usage of this is intentionally confusing, and that’s not even discussing “difference” of odds, which is what $\exp(\beta_1)$ would mean, and still difficult to interpret. Avoid odds if you can.³

My recommendation is to use the function to translate things onto the probability scale, you’ll likely have a more intuitive understanding of the values (even better yet, is to plot the entire probability function as a function of your variable). However, one of my professors pointed me toward a Car Talk riddle about potatoes that highlights some misconceptions about the probability scale as well. Be careful with interpretation, it’s very easy to get tripped up!

The Likelihood Principle

Likelihood is a frequently encountered principle for estimating parameters given some data. It is different from the concept of “probability”, even though commonly we’d use those words interchangeably. In fact, the fundamental idea of likelihood is that it’s kind of the “opposite” of probability. If we think about our data model as the statistical world, and our observations as the real world, the probability function is how we can move from our statistical model to the observations. In turn, the likelihood function uses what we observe in the real world to inform us of better parameters for our data model in the statistical world. We feed parameters into our probability function, and we feed data into our likelihood function.

Let’s skip thinking about the regression $p(x_i)$ for now, and just think about observations from a data model with a single parameter $p$. The definition of likelihood for a single observation from Bernoulli distribution is:

For multiple independent observations, we multiply the probabilities,

What?! Likelihood is defined by the probability function?! So does this mean likelihood can be interpreted as the probability that that $p$ is a particular value? Unfortunately no — and I get why this is confusing at first. In frequentist statistics, parameters are thought of as fixed values of the population, i.e. there is no randomness associated with parameters. Hence, it does not make sense to discuss the probability that $p$ takes a certain value.⁴ If you pause and think about what the definition says though, it will make sense because there must be a correspondence between the data model and our observations. It’s a two way road, but it’s just a matter of perspective which direction you’re currently moving in.

Substituting in our pmf for a Bernoulli distribution above, we get

Higher values of $\mathcal{L(p)}$ mean that value of $p$ is more likely to have generated the data that we fed in. Relative values of likelihood are also meaningful, in fact, this forms the basis of the likelihood ratio test. In order to find the most likely value of $p$ given our data, we would try maximizing the function. We call the value of $p$ that maximizes this function the maximum likelihood estimate (MLE) of $p$. In order to denote estimates of parameters, you’ll commonly see little hats on the parameter, $\widehat p_{MLE}$. Since the objective is to maximize this function, we normally take the log of the likelihood function, which doesn’t change the maximum since log is a monotonic function. The log likelihood allows us to avoid dealing with the exponential terms, and the expression is,

In this situation, it’s not difficult to calculate a closed form solution for the maximum likelihood estimate of the function. We’ll take the derivative and set it equal to 0, and solve for $p$ to get:

There’s a good reason the likelihood principle is so prominent — the likelihood function is very useful for guiding us toward good estimates of parameters for our model. However, finding an estimate for $\beta_0$ and $\beta_1$ is not as straightforward when we try to associate data with our logistic function, and no closed form solution exists like the simplified case.

Likelihood for Regression

The likelihood function for logistic regression is instead a function of the parameters $\beta_0, \beta_1$. Instead of just the number of 1’s and 0’s we observe in the data $\bm y = y_1, y_2, \ldots, y_n$, we must also take into consideration the values of $\bm x = x_1, x_2, \ldots, x_n$ that we observe. In a sense, the likelihood function transfers along to the next stage of modeling since $p_i$ itself is subsequently a function of $\beta_0, \beta_1$ and $x_i$.The likelihood is given by substituting the logistic function for $p$ in the binary likelihood function. If we simplify the expression, we get

It’s certainly not a friendly expression, and we are now dealing with two parameters $\beta_0, \beta_1$ instead of just a single $p$, so we are optimizing over a surface instead of just a curve, but we can still visualize the likelihood.

Hauck-Donner Effect

Above we saw that the logistic function “disappears” when we separate the data. Can you rationalize why this happens as you play around with the likelihood plot and the logistic regression function? The maximum likelihood no longer exists (or not numerically stable) in these situations because the likelihood surface becomes very flat, and the coefficients escape off to infinity. We can make the magnitude of $\beta_1$ arbitrarily large, because the logistic function will become sharper and there will also be many values that $\beta_0$ can take on. This phenomenon is known as the Hauck-Donner effect, and happens often in logistic regression, especially in higher dimensional spaces. It’s not the only way this can happen, but it’s important to keep in mind as you diagnose possible issues with your model. One sign is that the standard errors of your $\beta$ coefficients will be very large.

There are various ways of dealing with this issue, philosophical by nature:

• If you’re using logistic regression for classification, you can choose any of the values and get the same classification. You can also consider using a simpler model, like a threshold on $x$.
• Change your model by including different or excluding variables. It’s possible that with this separation, you were asking an ill-posed question in the first place, and a different formulation should be considered.
• You can also adjust the likelihood function with a penalty for large values of coefficients, commonly a ridge or lasso penalty. This is the approach of python’s scikit-learn implementation.
• Changing to a bayesian framework, if we add a prior distribution on the coefficients, we can get a well defined posterior distribution. bayesglm in R package “arm” uses this approach.
• Other solutions are listed on this stack exchange post.

Optimization Topics

Given that this is a relatively standard nonlinear (convex) optimization problem, you have a choice of many methods that will give you the same optimal answer. I won’t get too much into the nitty gritty, as chances are this will be mostly done already by a computing program. Probably the most common implementations are:

• Gradient Descent — A general purpose descent algorithm that is normally used for simplicity. It’s not as efficient as other methods.
• Newton Raphson — A faster algorithm that uses the second derivative of the likelihood function to find the optimal solution. However looking at the surface of the likelihood function, a quadratic approximation is not always reliable especially as the likelihood function becomes flatter.
• Iteratively Reweighted Least Squares — Used for models in generalized linear models because each iteration is solving a least squares problem, which also contains useful statistical information by default. This method is the same as “Fisher-Scoring”, which is similar to Newton-Raphson, but instead uses the expected value of the second derivative of the likelihood. This is the method used in R’s glm.fit function.

If you come from machine learning or engineering disciplines, it’s common instead to think about minimizing loss functions instead of maximizing likelihood functions. There is a correspondence. In the case of logistic regression, maximizing the likelihood is equivalent to minimizing the cross-entropy of the fitted probabilities of our model and the data values. There is an additional factor of $\frac{1}{n}$ that does not change the results of the optimization.