Logistic Regression: with MLE and more

The following post outlines the concept of logistic regression with the well-known Generalized Linear Model or GLM. Detailed maths is out of scope and we'll just make use of R to aid our learning. To begin, logistic regression is used to model and predict data with the following properties:
     • A categorical dependent variable, e.g. cure vs sick, old vs young (binary values)
     • A continuous, quantitative independent variable

In the discussion herewith, the response variable has 2 possible values, e.g. male/female, success/failure, dead/alive, etc. Our main goal is to predict the probability of either value occurring. A series of such binary responses (Bernoulli process) follows a binomial distribution. Probability density function (pdf) tells us the probability of a particular value for a fixed or known model or distribution, e.g. binomial distribution. Suppose the probability of having Y is denoted as P(Y) or p, where 0 ≤ p ≤ 1. The ratio of having Y and not having Y, p/(1 – p) is called odds. Thus, instead of simply predicting Y from X as in the linear regression, we predict the odds of Y for a given X. 

Logistic Regression as a Special Case of GLM

A simple logistic regression related to P(Y) can be expressed in the form of Equation (1). The word 'logistic' comes from the logistic function of probability shown on the left-hand side. This function can be transformed to the natural logarithmic function of odds called logit as shown on the right-hand side. The transformed model now takes logit as the dependent variable, which reminds us of the linear regression equation. Although the assumption of linearity (and others e.g. normality, homogeneity of variance) does not hold for logistic regression, the linearity between independent variable and logit has to be maintained! An extended version with multiple predictors can be seen in Equation (2).

Fig-1: (Left) An example of a logistic function σ(t) that takes any value between 0, 1. The function has a unique sigmoidal shape. (Right) An example of logistic regression with different values of betas [figure taken from ips8e, supp-Chapter 14, Macmillan].

When we talk about the simple linear regression, we assume X and Y to be linearly related more or less, the random error or residual follows a normal distribution, and the assumption of independence holds true. What happens when some of the assumptions collapse? "Generalized Linear Model" attempts to provide a more general form of regression when such assumptions are not fulfilled. It's sometimes confused with the general linear model found in the ANOVA and neuroimaging analyses. Now, there are three components of GLM:
     • An exponential family distribution for residuals.
     • A set of linear predictors that represents a systematic component.
     • A function that links or connects the expected value or mean of (1) and (2).

In GLM, the outcome variable takes up an exponential family distribution, e.g. normal or Gaussian, beta, gamma, exponential, chi-squared, Poisson, Bernoulli, and Dirichlet distribution. In the case of the ordinary linear regression, it is the Gaussian distribution. Suppose we predict the expected value of the outcome variable as shown by its mean μ. GLM allows a function of μ rather than the mean itself. This function, called the link function, is denoted as g(μ) because it links the mean of the outcome variable to the set of k predictors or explanatory variables.

      In summary, GLM formula states that:    g(μ)  =  a +  b₁X₁  + b₂X₂  + b₃X₃ + .... +  b_kX_k

In the case of linear regression, the link function is the identity link where g(μ) = μ. For binary data such as success/failure, a logit link is used instead. This is the case of the logistic regression. In R, logistic regression can be done by using glm command, instead of lm( ). Now, of so many exponential families, is there a unified way we estimate the parameters in GLM? Yeah, using a technique described in the next section, MLE.

Maximum Likelihood Estimation (MLE)
How do we estimate b? Whereas linear regression uses the least-squares method to minimize the errors, logistic regression uses the maximum likelihood method to arrive at the solution. Key to MLE is really just the "likelihood". Here, we first assume there is a model and this model is valid/correct. If it's correct, we can find the most likely estimate that explains the observed data.

Probabiliy means given a model of an event, what's the distribution of the outcome? Likelihood means given a set of observed outcomes, what are model parameters? 

Undergraduate stats course teaches us binomial probability. For example, suppose we toss an unfair coin 20 times and trials or tosses are independent of each other. In each trial, the probability of a head is p = 0.6, so the probability of not getting a head is 0.4. We try to describe the probability of getting heads r times out of n tosses. In fact, we may obtain a total of r = 0, 1, 2, ..., up to 20 heads from tossing that coin 20 times. We compute,  ₀C²⁰p⁰(1 – p)²⁰ ,  ₁C²⁰p¹(1 – p)¹⁹  , ... up to ₂₀C²⁰p²⁰(1 – p)⁰. Refer to Fig-2. The plot on the left is the binomial probability density function. The shaded area under the curve equals 1.

Fig-2: (Left) Discrete probability density function of tossing a coin 20 times, with the probability of obtaining head in each toss p = 0.6. The X-axis is the number of possible heads from 0 to 20. Obviously, if observing a head in each toss is 0.6 then I would expect 0.6 x 20 = 12 heads in 20 tosses. (Right) Suppose we were to estimate p here based on the observation that there are 13 heads out of 20 tosses. The best estimate of p is the one that has highest likelihood, in this case, the peak is at p = 0.65.





The step shown above is actually generating a possible data set from known parameters (n, p). What if the process is reversed? Twist your mind! We want to estimate p for a given series of n observations or outcomes. A likelihood function is exactly this, i.e. the opposite of probability density function. So, if we observed 13 heads in tossing the coin 20 times, what is the probability of a head in each coin toss? We try to estimate the value of p such that our observation is most likely to happen. Theoretically, if our model is correct, the most likely value of p = 13/20 = 0.65. The right panel Fig-2 shows this is the case. The highest likelihood is when the X-axis or p equals 0.65.

L(data | model) is a function that describes the probability of a model given our observed data set, p(model | data). If independence among trials is assumed, the likelihood function equals the joint-probability density which is none other than the product of each and every p(y_i | x_i), where y_i takes up a binary value.

As the end result of multiplication may be too small to manage, people take the natural logarithm of the likelihood which is then called log-likelihood (LL). Twist your mind again, taking the log of multiplication yields a summation of logs. Look at Equation (3). Now, since the numerical optimizer is designed to find the minimum, not the maximum, we multiply this by –1. Equation (4) is the final equation we want to minimize, just like minimizing sum of squares residuals in linear regression! Numerical optimization is used as there is no closed-form analytical solution. In R, this can be achieved by using optim() function.

Fitting a model to our data in logistic regression means to estimate the coefficient b₀ and b. The first is the intercept, the value of Y when X = 0. But how do we interpret b? The logit part is our dependent variable in the logistic regression. The logit is also easily converted back into the odds, where odds = exp(b₀ + bX). Suppose we have X₁ and X₂ = X₁ + 1; we compare both odds by finding the ratio of odds(X₂) = exp(b₀ + bX₁ + b) and odds(X₁) = exp(b₀ + bX₁), which is just equal to exp(b). Hence, the coefficient b can be interpreted as the log of the relative increase in the odds for a unit increase in score X. 

Assessing the model fit

There are a few common measures of model fit in the logistic regression. The first is the log-likelihood statistics as mentioned previously, i.e. summing the probabilities of the predicted and actual outcomes. This is analogous to the residual sum of squares (SS_R) in multiple regression in the sense that it is an indicator of how much unexplained information there is after the model has been fitted. The larger the value, the poorer the fitting is. A similar measure called deviance statistics is more popular, where deviance = –2 LL. Again, a higher number means a bad fit. We'll talk more about deviance later.

The second measure is Pearson's goodness of fit test. This is a statistical test to investigate how close values predicted by the model are to the observed values. The null hypothesis here is that our data follow the logistic regression model. We first find the chi-squared (c²) value and then the degree of freedom to find the p-value. The third performance assessment is by R², which is related to both the deviance and z² (Wald statistics). This measure is similar to the one we saw in the linear regression but with a different formula. R².can be unnecessarily inflated when we add more predictors. Like in multiple linear regression, we can use BIC or AIC to compute the parsimony measure of R², penalizing for useless extra predictors in the model. AIC/BIC is useful for choosing a model, with the preferred model is the one with a lower AIC/BIC index.

Software packages such as SPSS and R make use of the deviance statistics in the report. More importantly, there are two types of deviance. The null deviance shows how well the response is predicted by a reduced or baseline model. In the case of the simple logistic regression, the model assumes that the outcome can be predicted by nothing but a constant value, that is, the intercept b₀. In other words, it agrees with the null hypothesis that the estimated b = 0. The residual deviance shows how well the response is predicted by the full model when all predictors are included, implying that this is against the null hypothesis. Each deviance carries its own # dof.

The difference between the null and residual deviance is called likelihood ratio, another goodness of fit that also has a chi-squared distribution. The degree of freedom is defined as the difference between two deviances. We then easily get the p-value that suggests whether the model is better than chance at predicting the outcome.

How to assess the importance of the individual predictor in the logistic model? Similar to linear regression, it is assessed by carrying out statistical tests of the significance of the coefficient. Whereas t-test is used in linear regression, Wald test (z²) is used to evaluate the statistical significance of each predictor. It is calculated by taking the ratio of the square of the regression coefficient to the square of the standard error of the coefficient. The null hypothesis is similar, i.e. whether the estimated b is different from zero. As a side note: Wald test has been shown to be less reliable for small sample sizes!

Fig-3: Two R outputs from a logisitc model operation where you have one predictor (intervention) and two (intervention + duration) to predict the same dependent variable. Click to enlarge!


Let's take a look at Fig-3. Suppose there is a data frame containing three fields: whether the disease is cured (dependent variable), whether any intervention is performed (first categorical predictor), and the sickness duration (another, but quantitative, predictor). The model on the left fits the outcome variable (cured vs. not-cured) with whether or not treatment has been performed. What does R show us? First, the coefficients portion suggests there is a significant intervention predictor with estimated b = 1.23, z = 3.07, p < 0.005. This is Wald statistics. A z-value that is sufficiently far from 0 means that the estimate is both precise enough to be statistically different from 0 and large to have an effect on the response.

Next, the model fit. The null deviance = 154.08 and residual deviance = 144.16. The fact that residual deviance is lower suggests that one more predictor is better than just a constant in reducing the residual error. So adding "intervention" improves the fit by 9.93 unit. The difference in # dof of 1 so we can find the p-value of the chi-squared. Accordingly, c²(1) has p-value = 0.00163. The model on the right introduces another predictor, i.e. total duration. Based on the deviance values, there is no benefit of adding one more predictor. The book says that anova() can be used to compare both model fits. The test yields no significant difference between the left and right models.

Special example: Psychometric curve
A fundamental concept in psychophysics, the psychometric function relates a parameter of a sensory stimulus to a subjective response of a participant. We learnt previously how the psychometric curve has a sigmoidal shape and it's getting much clearer why. Take an example of a two-alternative-forced-choice (2AFC) task. What we're interested in is how the probability of responding with one of the two choices varies with the stimulus. Here, we will consider a situation where the varying sensory stimulus is the movement direction θ of a participant's arm to the left or right of the body midline. This becomes the continuous independent variable X. The desired response is binary, either left or right. This response is our dependent variable Y.

The model suitable for this is the simple logistic regression shown in Equation (1). Let p = P(y = right | X) is the probability of responding "right" given a certain direction X. So, 1 – p is equivalent to P(y = left | X). We typically represent the response in numerical format, 1 = right and 0 = left. The value X in which p = 0.5 corresponds to the perceptual threshold that in ideal case is the body midline, but one may have a perceptual bias. From our experiment, we obtain an array of direction X (in degree) and an array of binary responses R (in number 0, 1). First construct our cost function, name it NLL. This is then fed into a numerical optimizer to estimate b₀ and b using MLE. An initial guess for each estimate has to be given, e.g. (0.1, 0.1).

NLL <- function(B,X,R) {
          y = B[1] + B[2]*X
          p = 1/(1+exp(-y))
          NLL = -sum(log(p[R==1])) - sum(log(1-p[R==0]))
  }

out = optim(par=c(-.1,.1),NLL,X=X,R=R)    #Numerical optimization! 
Bfit = out$par    #Retrieve parameter estimates of our model

#Let's construct our model with Bfits above and plot the data...
Xp = seq(-15,15,1)
myModel = logistic(Bfit[1] + Bfit[2]*Xp)

PS: I have recently found another nice R package called "quickpsy" for computing with psychometric functions.


Main references
(1) Andy Field, et al, "Discovering Statistics using R" (2012).
(2) Gribble's note on MLE, Psychology_9041B course.