Contents

Logistic Regression

The cover image was generated by ChatGPT.

Introduction

Logistic regression is a type of binary classification model derived from linear regression. Unlike linear regression, which fits a regression line that passes through the data points by minimizing the sum of squared distances, known as the least squares method. The goal of logistic regression is to find a decision boundary that clearly separates the data into two distinct classes.

https://raw.githubusercontent.com/Josh-test-lab/website-assets-repository/refs/heads/main/posts/Logistic%20Regression/linear_logistic.webp
Linear regression and logistic regression.

Odds and Odds Ratio

When discussing logistic regression, it is essential to first understand the concepts of odds and odds ratio.

Odds

If the probability of an event occurring under certain conditions is $p$, then the odds of the event are defined as the ratio of the probability that the event occurs to the probability that it does not:

$$ \text{odds} = \frac{p}{1 - p} $$

In other words, if the odds are $w$, it means the event is $w$ times more likely to occur than not.

Odds Ratio

The odds ratio (OR) is the ratio of two odds, typically used to compare the likelihood of an event occurring between two groups — one with a particular characteristic and one without:

$$ \text{OR} = \frac{\text{odds}_1}{\text{odds}_2} = \frac{p_1 / (1 - p_1)}{p_2 / (1 - p_2)} $$

From the odds ratio, we can infer:

  • If OR $> 1$: the event is more likely to occur under the specific condition.
  • If OR $= 1$: the event is equally likely to occur with or without the condition, no association.
  • If OR $< 1$: the event is less likely to occur under the specific condition.
Example

If the probability of success after receiving a treatment is $0.8$, and for those without the treatment it is $0.5$, then:

  • The odds for the treatment group are $\frac{0.8}{1 - 0.8} = 4$
  • The odds for the control group are $\frac{0.5}{1 - 0.5} = 1$

Thus, the odds ratio OR = $4 / 1 = 4$, meaning the likelihood of success with treatment is 4 times that without treatment.

Logit and the Logit Function

Logit

The logit (or logit function) is closely related to odds. It originates from the logit transformation, which maps a probability value between 0 and 1 to the entire real number line $\mathbb{R}$, and is defined as:

$$ \text{logit}(p) = \ln \left( \frac{p}{1 - p} \right) $$

In other words, the logit is the logarithm of the odds. This transformation enables us to convert a nonlinear probability problem into a linear predictive model.

Furthermore, the difference between the logits of two probabilities $p_1$ and $p_2$ can be shown to equal the logarithm of their odds ratio:

Proof
$$ \begin{align*} \text{logit}(p_1) - \text{logit}(p_2) & = \ln \left( \frac{p_1}{1 - p_1} \right) - \ln \left( \frac{p_2}{1 - p_2} \right) \\ & = \ln \left( \frac{p_1 / (1 - p_1)}{p_2 / (1 - p_2)} \right) \\ & = \ln (\text{OR}) \end{align*} $$

Logistic Function

The logistic function is the core of logistic regression and is also known as the sigmoid function. It is defined as:

$$ \sigma(x) = \frac{1}{1 + e^{-x}} $$

The input value $x$ can be any real number, while the output of the logistic function is constrained to the interval $[0, 1]$.

Proof

The logistic function can be derived from the logit function.

Let $$ \text{logit}(p) = \ln \left(\frac{p}{1 - p}\right) = x $$

We can then solve for $p$:

$$ \begin{align*} & \text{logit}(p) = \ln \left(\frac{p}{1 - p}\right) = x \\ \iff & \frac{p}{1 - p} = e^x \\ \iff & p = e^x (1 - p) \\ \iff & p = e^x - e^x p \\ \iff & p + e^x p = e^x \\ \iff & p(1 + e^x) = e^x \\ \iff & p = \frac{e^x}{1 + e^x} \\ \iff & p = \frac{1}{1 + e^{-x}} = \sigma(x) \end{align*} $$

Thus, we can see that the logistic function is the inverse of the logit transformation:

$$ \sigma(x) = \text{logit}^{-1}(x) $$

Logistic Regression

Now, suppose the dependent variable $Y$ is a binary random variable, where $Y \in [0, 1]$. Here, $1$ represents a successful event, and $0$ represents a failure. Let $Y = f(\boldsymbol{X})$, and model the conditional probability using the logit function:

$$ \begin{align*} f(\boldsymbol{X}) & = \ln \left( \frac{\boldsymbol{p}}{1 - \boldsymbol{p}} \right) \\ & = \begin{bmatrix} \beta_0 + \beta_1 x_{11} + \beta_2 x_{12} + \cdots + \beta_n x_{1n} \\ \beta_0 + \beta_1 x_{21} + \beta_2 x_{22} + \cdots + \beta_n x_{2n} \\ \vdots \\ \beta_0 + \beta_1 x_{m1} + \beta_2 x_{m2} + \cdots + \beta_n x_{mn} \\ \end{bmatrix}_{m \times 1} \\ & = \boldsymbol{X} \boldsymbol{\beta} \end{align*} $$

Here, $\boldsymbol{X}$ is the input variable matrix (also known as the design matrix) that includes the intercept term. Each row represents one observation, with the first column consisting of the constant term $1$, and the remaining columns containing $n$ independent variables (features). $\boldsymbol{\beta} = (\beta_0, \beta_1, \beta_2, \cdots, \beta_n)^\top$ is the coefficient vector.

$$ \boldsymbol{X} = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1n} \\ 1 & x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{m1} & x_{m2} & \cdots & x_{mn} \end{bmatrix}_{m \times (n+1)} $$

And $\boldsymbol{p} = \Pr(Y = 1 \mid \boldsymbol{X}) = (p_1, \cdots, p_m)$ is the probability vector representing the likelihood of success. Based on the previous proof and by applying the logistic function, we obtain the estimated probability vector $\hat{\boldsymbol{p}}$ as:

$$ \hat{\boldsymbol{p}} = \sigma(f(\boldsymbol{X})) = \sigma(\boldsymbol{X} \boldsymbol{\beta}) = \frac{\exp(\boldsymbol{X} \boldsymbol{\beta})}{1 + \exp(\boldsymbol{X} \boldsymbol{\beta})} $$

Similarly, the probability of failure, $1 - \hat{\boldsymbol{p}} = \Pr(Y = 0 \mid \boldsymbol{X})$, can be written as:

$$ 1 - \hat{\boldsymbol{p}} = \frac{1}{1 + \exp(\boldsymbol{X} \boldsymbol{\beta})} $$

Loss Function

For logistic regression, the commonly used loss function is cross-entropy. Since the response variable $Y$ is binary, the discrete version of cross-entropy is used, with the following formula:

$$ \text{entropy} = -y \ln (p) $$

where, $y$ represents the distribution (i.e., the actual label), and $p$ is the predicted probability. Therefore, the loss function for a single variable can be defined as:

$$ \begin{align*} L(y_i, \hat{p_i}) & = \left\{\begin{matrix} -\ln (\hat{p_i}), & \text{if } y_i = 1; \\ -\ln (1 - \hat{p_i}), & \text{if } y_i = 0. \end{matrix}\right. \\ & = -y_i \ln (\hat{p_i}) - (1 - y_i) \ln (1 - \hat{p_i}) \end{align*} $$

where $i = 1, \cdots, m$, $\hat{p}_i = \sigma(X_i^\top \boldsymbol{\beta})$ is the predicted probability for the $i$-th observation, with $X_i = (1, x_{i1}, x_{i2}, \cdots, x_{in})^\top$ and $\boldsymbol{\beta} = (\beta_0, \beta_1, \beta_2, \cdots, \beta_n)^\top$.

For the entire dataset represented by the matrix $\boldsymbol{X}$, the loss function is defined as:

$$ \begin{align*} L(Y, \hat{\boldsymbol{p}}) & = \sum_{i = 1}^m \left[ -y_i \ln (\hat{p_i}) - (1 - y_i) \ln (1 - \hat{p_i}) \right] \\ & = -\sum_{i = 1}^m \left[ y_i \ln(\sigma(X_i^\top \boldsymbol{\beta})) + (1 - y_i) \ln(1 - \sigma(X_i^\top \boldsymbol{\beta})) \right] \\ & = -\sum_{i = 1}^m y_i \ln(\sigma(X_i^\top \boldsymbol{\beta})) -\sum_{i = 1}^m (1 - y_i) \ln(1 - \sigma(X_i^\top \boldsymbol{\beta})) \\ & = -\sum_{y_i = 1} \ln(\sigma(X_i^\top \boldsymbol{\beta})) -\sum_{y_i = 0}\ln(1 - \sigma(X_i^\top \boldsymbol{\beta})) \\ & = -\sum_{y_i = 1} \ln(\hat{p_i}) -\sum_{y_i = 0}\ln(1 - \hat{p_i}) \end{align*} $$

where, $X_i = (1, x_{i1}, x_{i2}, \cdots, x_{in})^\top$, $Y = (y_1, y_2, \cdots, y_m)^\top$, with $i = 1, \cdots, m$.

Likelihood Function

Since $\boldsymbol{X}$ is a fixed and known variable, minimizing the error, that is, making the loss function approach zero, means we are effectively searching for an optimal coefficient vector $\boldsymbol{\beta}$ such that the model’s predicted probabilities $\hat{\boldsymbol{p}}$ are as close as possible to the observed values $Y$.

Assuming each observed sample $y_i$ in $Y$ is conditionally independent and follows a Bernoulli distribution, the probability density function for each sample is:

$$ \Pr(y_i \mid X_i, \boldsymbol{\beta}) = \hat{p}_i^{y_i} (1 - \hat{p_i})^{1 - y_i} $$

Therefore, the likelihood function can be written as:

$$ \mathcal{L}(\boldsymbol{\beta}) = \prod_{i = 1}^m \hat{p}_i^{y_i} (1 - \hat{p_i})^{1 - y_i} $$

Log-Likelihood Function

To find the optimal coefficients $\boldsymbol{\beta}$, we take the first derivative of the likelihood function $\mathcal{L}(\boldsymbol{\beta})$ and identify the coefficients that maximize it. This method is called maximum likelihood estimation (MLE).

However, since the likelihood function involves a product over many terms $\prod$, which is computationally inconvenient, we usually take the logarithm of the likelihood function to convert the product into a sum for easier calculation. This is known as the log-likelihood function, computed as follows:

$$ \begin{align*} \ell(\boldsymbol{\beta}) & = \ln (\mathcal{L}(\boldsymbol{\beta})) \\ & = \ln \left( \prod_{i = 1}^m \hat{p}_i^{y_i} (1 - \hat{p_i})^{1 - y_i} \right) \\ & = \sum_{i=1}^m \left( \ln \left( \hat{p}_i^{y_i} (1 - \hat{p_i})^{1 - y_i} \right) \right) \\ & = \sum_{i=1}^m \left[ \ln \left( \hat{p}_i^{y_i} \right) + \ln \left( (1 - \hat{p_i})^{1 - y_i} \right) \right] \\ & = \sum_{i=1}^m \left[ y_i \ln(\hat{p}_i) + (1 - y_i) \ln(1 - \hat{p}_i) \right] \\ & = \sum_{y_i = 1} \ln(\hat{p_i}) + \sum_{y_i = 0}\ln(1 - \hat{p_i}) \end{align*} $$

Maximum Likelihood Estimation

From the previous calculations, we find that

$$ L(Y, \hat{\boldsymbol{p}}) = -\ell(\boldsymbol{\beta}) $$

Therefore, finding an estimate $\hat{\boldsymbol{\beta}}$ that minimizes the loss function is equivalent to finding the maximum value of the log-likelihood function.

$$ \begin{align*} \hat{\boldsymbol{\beta}} & = \argmin_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} L(Y, \hat{\boldsymbol{p}}) \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \left( - L(Y, \hat{\boldsymbol{p}}) \right) \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \ell(\boldsymbol{\beta}) \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \left( \sum_{y_i = 1} \ln(\hat{p_i}) + \sum_{y_i = 0}\ln(1 - \hat{p_i}) \right) \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \left( \sum_{y_i = 1} \ln \left( \frac{\exp(X_i^\top \boldsymbol{\beta})}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right) + \sum_{y_i = 0}\ln \left( \frac{1}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right) \right) \end{align*} $$

Note

If the logarithm is not used, and the likelihood function is directly maximized to find $\hat{\boldsymbol{\beta}}$, the calculation is as follows:

$$ \begin{align*} \hat{\boldsymbol{\beta}} & = \argmin_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} L(Y, \hat{\boldsymbol{p}}) \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \mathcal{L}(\boldsymbol{\beta}) \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \prod_{i = 1}^m \hat{p}_i^{y_i} (1 - \hat{p_i})^{1 - y_i} \\ & = \argmax_{\boldsymbol{\beta} \in \mathbb{R}^{n + 1}} \prod_{i = 1}^m \left( \frac{\exp(X_i^\top \boldsymbol{\beta})}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right)^{y_i} \left( \frac{1}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right)^{1 - y_i} \end{align*} $$

The computational complexity of this approach is significantly higher compared to directly using the log-likelihood function!

To find the optimal estimate $\hat{\boldsymbol{\beta}}$, we use maximum likelihood estimation (MLE). This involves taking the first derivative of the negative log-likelihood function $\ell(\boldsymbol{\beta})$ with respect to $\boldsymbol{\beta}$, setting it to zero, and solving for the critical points.

First, let’s rewrite the log-likelihood function $\ell(\boldsymbol{\beta})$:

$$ \begin{align*} \ell(\boldsymbol{\beta}) & = \sum_{y_i = 1} \ln \left( \frac{\exp(X_i^\top \boldsymbol{\beta})}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right) + \sum_{y_i = 0}\ln \left( \frac{1}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right) \\ & = \sum_{i=1}^m \left[ y_i \ln \left( \frac{\exp(X_i^\top \boldsymbol{\beta})}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right) + (1 - y_i) \ln \left( \frac{1}{1 + \exp(X_i^\top \boldsymbol{\beta})} \right) \right] \\ & = \sum_{i=1}^m \left[ y_i \ln \left( \exp(X_i^\top \boldsymbol{\beta}) \right) - y_i \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) + (1 - y_i) \ln \left(1 \right) - (1 - y_i) \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right] \\ & = \sum_{i=1}^m \left[ y_i \ln \left( \exp(X_i^\top \boldsymbol{\beta}) \right) - y_i \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) - \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) + y_i \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right] \\ & = \sum_{i=1}^m \left[ y_i \ln \left( \exp(X_i^\top \boldsymbol{\beta}) \right) - \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right] \\ & = \sum_{i=1}^m \left[ y_i X_i^\top \boldsymbol{\beta} - \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right] \\ \end{align*} $$

Its first-order partial derivative is calculated as follows:

$$ \begin{align*} \frac{\partial}{\partial \beta_0} \ell(\boldsymbol{\beta}) & = \frac{\partial}{\partial \beta_0} \ell(\boldsymbol{\beta}) \\ & = \frac{\partial}{\partial \beta_0} \left( \sum_{i=1}^m \left[ y_i X_i^\top \boldsymbol{\beta} - \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right] \right) \\ & = \frac{\partial}{\partial \beta_0} \left( \sum_{i=1}^m y_i X_i^\top \boldsymbol{\beta} - \sum_{i=1}^m \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right) \\ & = \sum_{i=1}^m y_i \frac{\partial}{\partial \beta_0} X_i^\top \boldsymbol{\beta} - \sum_{i=1}^m \frac{\partial}{\partial \beta_0} \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \\ & = \sum_{i=1}^m y_i - \sum_{i=1}^m \hat{p_i} \\ & = \sum_{i=1}^m (y_i - \hat{p_i}) \end{align*} $$

where $\frac{\partial}{\partial \beta_0} X_i^\top \boldsymbol{\beta} = 1$, $\frac{\partial}{\partial \beta_0} \ln(1 + \exp(X_i^\top \boldsymbol{\beta})) = \frac{\exp(X_i^\top \boldsymbol{\beta})}{1 + \exp(X_i^\top \boldsymbol{\beta})} = \sigma(X_i^\top \boldsymbol{\beta}) = \hat{p_i}$.

$$ \begin{align*} \frac{\partial}{\partial \beta_{j, j \neq 0}} \ell(\boldsymbol{\beta}) & = \frac{\partial}{\partial \beta_{j, j \neq 0}} \ell(\boldsymbol{\beta}) \\ & = \frac{\partial}{\partial \beta_{j, j \neq 0}} \left( \sum_{i=1}^m \left[ y_i X_i^\top \boldsymbol{\beta} - \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \right] \right) \\ & = \sum_{i=1}^m y_i \frac{\partial}{\partial \beta_{j, j \neq 0}} X_i^\top \boldsymbol{\beta} - \sum_{i=1}^m \frac{\partial}{\partial \beta_{j, j \neq 0}} \ln \left( 1 + \exp(X_i^\top \boldsymbol{\beta}) \right) \\ & = \sum_{i=1}^m y_i x_{ij} - \sum_{i=1}^m x_{ij} \hat{p_i} \\ & = \sum_{i=1}^m x_{ij} (y_i - \hat{p_i}) \end{align*} $$

where $\frac{\partial}{\partial \beta_{j, j \neq 0}} X_i^\top \boldsymbol{\beta} = x_{ij}$, $\frac{\partial}{\partial \beta_{j, j \neq 0}} \ln(1 + \exp(X_i^\top \boldsymbol{\beta})) = \frac{x_{ij} \exp(X_i^\top \boldsymbol{\beta})}{1 + \exp(X_i^\top \boldsymbol{\beta})} = x_{ij} \sigma(X_i^\top \boldsymbol{\beta}) = x_{ij} \hat{p_i}$, $i = 1, 2, \cdots, m$ , $j = 0, 1, \cdots, n$.

Combining the above two expressions, we obtain:

$$ \frac{\partial}{\partial \beta_{j}} \ell(\boldsymbol{\beta}) = \sum_{i=1}^m x_{ij} (y_i - \hat{p_i}) $$

where $x_{i0} = 1$, $j = 0, 1, \cdots, n$.

Assuming the first-order partial derivative equals zero, we have:

$$ \begin{align*} & \frac{\partial}{\partial \beta_j} \ell(\boldsymbol{\beta}) = 0 \\ \iff & \sum_{i=1}^m x_{ij} (y_i - \hat{p_i}) = 0 \\ \end{align*} $$

At this point, we realize that logistic regression does not have a closed-form solution like linear regression. Therefore, we need to change our approach and use numerical methods to find the solution.

Numerical Methods for Solution

In general, we use gradient descent to minimize the loss function $L(Y, \hat{\boldsymbol{p}}) = L(\boldsymbol{\beta})$, which helps approximate the true coefficient vector $\boldsymbol{\beta}$ by iteratively updating $\hat{\boldsymbol{\beta}}$.

Let the gradient be defined as:

$$ \begin{align*} \nabla L(\boldsymbol{\beta}) & = \frac{\partial}{\partial \boldsymbol{\beta}} L(\boldsymbol{\beta}) \\ & = \left( \frac{\partial L}{\partial \beta_0}, \frac{\partial L}{\partial \beta_1}, \cdots, \frac{\partial L}{\partial \beta_n} \right) \end{align*} $$

From the previous calculations, we know that

$$ \frac{\partial L}{\partial \beta_j} = \frac{\partial}{\partial \beta_j} \left( -\ell(\boldsymbol{\beta}) \right) = \sum_{i=1}^m x_{ij}(\hat{p_i} - y_i) $$

Therefore, the gradient can be expressed as:

$$ \nabla L(\boldsymbol{\beta}) = \boldsymbol{X}^\top (\hat{\boldsymbol{p}} - Y) \in \mathbb{R}^{n+1} $$

The update function is as follows:

$$ \boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - \eta \cdot \nabla L(\boldsymbol{\beta}^{(t)}) $$

Here, $t$ denotes the $t$-th iteration; $\eta$ is the learning rate, which controls the step size of each update. After multiple iterations, $\boldsymbol{\beta}^{(t)}$ will gradually approach the solution that minimizes the loss function. When $\boldsymbol{\beta}^{(t)}$ no longer changes significantly or the maximum number of iterations is reached, we denote it as $\hat{\boldsymbol{\beta}}$, which serves as the final estimated coefficient vector for the model.

Prediction Method

Once the estimate $\hat{\boldsymbol{\beta}}$ is obtained, we can calculate the predicted probabilities as $\hat{\boldsymbol{p}} = \sigma(\boldsymbol{X}\hat{\boldsymbol{\beta}})$, and then convert these probabilities into predicted classes using a threshold.

Since the probability values lie within the interval $[0, 1]$, the threshold is commonly set to 0.5, such that:

$$ \begin{align*} \hat{y_i} = \left\{\begin{matrix} 1, & \text{if \ } \hat{p_i} \geq 0.5 \\ 0, & \text{if \ } \hat{p_i} < 0.5 \end{matrix}\right., && \text{where } i = 1, 2, \cdots, m \end{align*} $$

Hypothesis Testing

After estimating the coefficients $\hat{\boldsymbol{\beta}}$, we want to determine whether these variables have a statistically significant effect on the outcome. This can be tested through hypothesis testing.

Wald Test

For testing individual variables, we can use the Wald statistic to perform the Wald test. The null hypothesis is as follows:

$$ \left\{\begin{align*} & H_0: \theta = \theta_0 \\ & H_a: \theta \neq \theta_0 \end{align*}\right. $$

where, $H_0$ represents the null hypothesis, and $H_a$ represents the alternative hypothesis; $\theta$ is the parameter being tested, and $\theta_0$ is a constant. The Wald test uses the Wald statistic $W$ to perform the test. When the following condition holds, we reject the null hypothesis $H_0: \theta = \theta_0$.

$$ \begin{align*} &W > z_{1 - \alpha / 2} \\ &W < z_{\alpha / 2} \\ &W^2 > \chi^2_{1, 1-\alpha} \end{align*} $$

where, $\alpha$ is the significance level, where $\alpha \in (0, 1)$; $z_\alpha$ is the critical value from the standard normal distribution corresponding to the two-tailed $\alpha$ tail; and $\chi^2_{k, \alpha}$ is the critical value from the chi-square distribution with $k$ degrees of freedom at significance level $\alpha$. The Wald statistic $W$ is calculated as follows:

$$ W = \frac{\hat{\theta}_{MLE} - \theta}{s(\hat{\theta})} \sim N(0, 1) $$

where, $\hat{\theta}_{MLE}$ is the MLE of $\theta$; $N(0, 1)$ denotes the standard normal distribution; and $s(\hat{\theta})$ is the standard error of $\hat{\theta}_{MLE}$, which can be estimated from its variance as follows:

$$ s^2(\hat{\theta}) = \operatorname{Var}(\hat{\theta}) \approx - \frac{1}{\ell’’(\hat{\theta})} $$

where, $\ell’’(\hat{\theta})$ is the second derivative of the log-likelihood function $\ell(\theta)$, which corresponds to the negative inverse of the Fisher information (FI, denoted as $\mathcal{I}$) evaluated at the estimate:

$$ \mathcal{I}(\hat{\theta}) = - \ell’’(\hat{\theta}) $$

Therefore, $s^2(\hat{\theta})$ can also be expressed as:

$$ s^2(\hat{\theta}) = \frac{1}{\mathcal{I}(\hat{\theta})} $$

The standard error can be estimated as follows:

$$ s(\hat{\theta}) = \sqrt{\frac{1}{\mathcal{I}(\hat{\theta})}} = \sqrt{- \frac{1}{\ell’’(\hat{\theta})}} $$

Univariate Hypothesis Testing

We perform the following hypothesis test for the $j$-th variable:

$$ \left\{\begin{align*} & H_0: \beta_j = 0 \\ & H_a: \beta_j \neq 0 \end{align*}\right. $$

where, $j = 0, 1, \cdots, n$. This hypothesis tests whether the $j$-th variable has a significant effect on the prediction target, which can be evaluated using the Wald statistic.

Let $\hat{\boldsymbol{\beta}}$ be the MLE of $\boldsymbol{\beta}$. According to the asymptotic normality, we have:

$$ \hat{\boldsymbol{\beta}} - \boldsymbol{\beta} \sim N(0, s^2(\hat{\boldsymbol{\beta}})) $$

Therefore, the Wald statistic can be defined as follows:

$$ W = \frac{\hat{\beta}_j}{s(\hat{\beta}_j)} $$

where, $s(\hat{\beta}_j)$ is the standard error of the estimate $\hat{\beta}_j$. Based on asymptotic properties, the Wald statistic $W$ approximately follows a standard normal distribution $N(0, 1)$.

When setting a significance level $\alpha$, we reject the null hypothesis $H_0$ in a two-tailed test if:

$$ |W| > z_{1 - \alpha / 2} $$

Alternatively, by squaring the statistic, we can perform the test using the chi-square distribution:

$$ W^2 > \chi^2_{1, 1-\alpha} $$

If we reject $H_0$, it indicates that the $j$-th variable has a statistically significant effect on the target variable; otherwise, there is no significant evidence to suggest that this variable is important.

Multivariate Hypothesis Testing

We commonly use the Likelihood Ratio Test (LRT) to assess whether multiple variables have a significant impact on the model. The principle is to compare the log-likelihood values of two models. The hypotheses for the LRT are as follows:

$$ \left\{\begin{align*} & H_0: \beta_1 = \beta_2 = \cdots = \beta_n = 0 \\ & H_a: \text{At least one of } \beta_j \neq 0 \end{align*}\right. $$

The null hypothesis corresponds to the reduced model, which assumes that none of the independent variables affect the target variable; the alternative hypothesis corresponds to the full model, which assumes that at least one independent variable has a significant effect.

The test statistic is defined as:

$$ \Lambda = -2 \left[ \ell(\boldsymbol{\beta}_{\text{reduced}}) - \ell(\boldsymbol{\beta}_{\text{full}}) \right] $$

where, $\ell(\hat{\boldsymbol{\beta}}_{\text{reduced}})$ is the log-likelihood function of the reduced model; $\ell(\hat{\boldsymbol{\beta}}_{\text{full}})$ is that of the full model; and $\Lambda$ is the likelihood ratio statistic.

Under the null hypothesis $H_0$, $\Lambda$ follows a chi-square distribution with $q$ degrees of freedom ($\chi^2_q$).

We reject the null hypothesis $H_0$ if:

$$ \Lambda > \chi^2_{q, 1 - \alpha} $$

where $\alpha$ is the significance level.

Confidence Interval

For each coefficient $\beta_j$, we can construct a $(1 - \alpha)$% confidence interval in the form:

$$ [\hat{\beta}_j \pm z_{\alpha/2} \cdot s(\hat{\beta}_j)] $$

where $\hat{\beta}_j$ is the estimated coefficient; $z_{\alpha/2}$ is the critical value from the standard normal distribution corresponding to the two-tailed $\frac{\alpha}{2}$ regions; $\alpha$ is the significance level; and $s(\hat{\beta}_j)$ is the standard error of $\hat{\beta}_j$.

This confidence interval can be interpreted as: if the sampling process is repeated many times, there is a $(1 - \alpha)$% probability that the interval contains the true parameter $\beta_j$. If the interval includes 0, it means that at the given significance level, we cannot reject the null hypothesis $\beta_j = 0$, implying the variable may not have a significant effect on the outcome.

If we want to interpret the effect in terms of the odds ratio, the corresponding estimate is:

$$ \widehat{\text{OR}}_j = \exp(\hat{\beta}_j) $$

which represents the multiplicative change in the odds of the event occurring for each one-unit increase in the independent variable vector $(x_{1j}, x_{2j}, \cdots, x_{mj})$. The $(1 - \alpha)$% confidence interval for the odds ratio is:

$$ \left[ \exp\left(\hat{\beta}_j - z_{\alpha/2} \cdot s(\hat{\beta}_j)\right), \exp\left(\hat{\beta}_j + z_{\alpha/2} \cdot s(\hat{\beta}_j)\right) \right] $$

This interval can be interpreted as: if the sampling process is repeated many times, there is a $(1 - \alpha)$% chance that the interval contains the true odds ratio $\text{OR}_j$. If the interval includes 1, it suggests that the variable’s effect may not be statistically significant, since $\exp(0) = 1$ represents no effect.

Python Example

Below is an example using the Iris dataset, where petal width is used as the feature to classify between setosa and versicolor.

First, import the necessary modules.

1
2
3
4
5
6
7
8
# import modules
import numpy as np
import matplotlib.pyplot as plt
import os
from sklearn.datasets import load_iris
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.metrics import confusion_matrix, accuracy_score
import scipy.stats

Load the Iris datasets.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# Iris data
iris = load_iris()  # load data
X = iris.data[:, 3:]  # data
y = iris.target  # target
class_names = iris.target_names  # Iris class names

mask = np.isin(iris.target, [0, 1])
X = X[mask]
y = y[mask]
class_names = class_names[0:2]

print(f'\nclass names: \n{class_names}')
print(f'\nfeature names: \n{iris.feature_names[3:]}')
Execution result reference
1
2
3
4
5
class names: 
['setosa' 'versicolor']

feature names: 
['petal width (cm)']

Set the image path.

1
2
3
# configs
output_dir = 'plots'  # output directory
os.makedirs(output_dir, exist_ok=True)

Perform logistic regression and plot the results.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# logistic regression
log_reg = LogisticRegression()
log_reg.fit(X, y)
y_log_prob = log_reg.predict_proba(x_test)[:, 1]

plt.figure(figsize=(8, 5))
plt.scatter(X[y == 0], y[y == 0], color='red', label=class_names[0])
plt.scatter(X[y == 1], y[y == 1], color='blue', label=class_names[1])
plt.plot(x_test, y_log_prob, label='Logistic Regression')
plt.xlabel('Petal Width (cm)')
plt.ylabel('Probability')
plt.title('Logistic Regression')
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1.0), borderaxespad=0)
plt.tight_layout()
plt.savefig(os.path.join(output_dir, 'logistic.webp'))
plt.show()
plt.close()

https://raw.githubusercontent.com/Josh-test-lab/website-assets-repository/refs/heads/main/posts/Logistic%20Regression/logistic.webp
Logistic regression for classifying setosa and versicolor using petal width.

The regression coefficients $\boldsymbol{\beta}$ are as follows. Since only one independent variable is used, there are only the intercept $\beta_0$ and the coefficient $\beta_1$.

1
2
3
4
# logistic regression coefficients
coef = log_reg.coef_[0][0]
intercept = log_reg.intercept_[0]
print(f'\nCoefficients:\nintercept = {intercept:.4f}\ncoef = {coef:.4f}')
Execution result reference
1
2
3
Coefficients:
intercept = -3.3799
coef = 4.4240

Below are the standard errors, p-values, and 95% confidence intervals for each coefficient.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# logistic regression summary
X_design = np.hstack([np.ones((X.shape[0], 1)), X])
p = log_reg.predict_proba(X)[:, 1]
W = np.diag(p * (1 - p))
Fisher_inv = np.linalg.inv(X_design.T @ W @ X_design)
se = np.sqrt(np.diag(Fisher_inv))

z = np.array([intercept, coef]) / se
p_values = 2 * (1 - scipy.stats.norm.cdf(np.abs(z)))
ci_95 = np.array([
    [intercept, coef] - 1.96 * se,
    [intercept, coef] + 1.96 * se
]).T

print("\nSummary of Logistic Regression Coefficients:")
for i, name in enumerate(['intercept', 'petal width']):
    print(f"{name:12s}: coef = {intercept if i==0 else coef:.4f} | SE = {se[i]:.4f} | p = {p_values[i]} | 95% CI = [{ci_95[i,0]:.4f}, {ci_95[i,1]:.4f}]")
Execution result reference
1
2
3
Summary of Logistic Regression Coefficients:
intercept   : coef = -3.3799 | SE = 0.6242 | p = 6.140579889013509e-08 | 95% CI = [-4.6034, -2.1565]
petal width : coef = 4.4240 | SE = 0.7091 | p = 4.40006697743911e-10 | 95% CI = [3.0342, 5.8138]

Below are the confusion matrix and accuracy score.

1
2
3
4
5
6
# confusion matrix and accuracy
y_log_pred = log_reg.predict(X)
cm = confusion_matrix(y, y_log_pred)
acc = accuracy_score(y, y_log_pred)
print(f'\nConfusion Matrix:\n{cm}')
print(f'\nAccuracy: {acc}')
Execution result reference
1
2
3
4
5
Confusion Matrix:
[[50  0]
 [ 0 50]]

Accuracy: 1.0

Conclusion

Logistic regression is a simple and practical binary classification model that can quickly classify two categories with minimal computation. However, it also has some limitations, such as assuming a linear relationship between the independent variables and the log-odds, and assuming observations are independent. Therefore, when using logistic regression, it is important to perform model diagnostics and proper data preprocessing to ensure the reliability of the analysis results.

As a classic foundational classification model, logistic regression holds an important position in statistical inference and is often used for preliminary modeling and feature selection in machine learning. It remains one of the favorite tools among data analysts.

Further Learning

References