Maximum Likelihood Estimation

In plain language

Maximum likelihood estimation (MLE) finds the model parameter values that make your observed data as probable as possible under an assumed statistical model — think of it as asking "which settings would have produced these exact numbers most often if we repeated the study many times?" For a simple example, if 7 of 10 patients had an event, MLE says the best estimate of the true event rate is 70%, because that value makes those exact results more probable than any other. MLE is the estimation engine behind logistic regression, Cox survival models, Poisson models, and virtually every other regression tool used in real-world evidence research, so understanding it unlocks the logic underlying the standard errors, confidence intervals, and model-comparison statistics you see in every analysis output. Its main limitation is that results can become unreliable when data are sparse or when a predictor perfectly separates the outcome groups, in which case specialized corrections such as Firth penalization are needed.

What maximum likelihood estimation is

Maximum likelihood estimation (MLE) answers one question: which parameter values would have made the observed data most probable? The key conceptual move is treating the probability model backwards. Ordinarily a fixed parameter — say, an event rate p — generates data; MLE instead fixes the observed data and asks which p best explains what was seen. The likelihood function L(θ | data) is the probability (or probability density) of the observed data evaluated as a function of the unknown parameter θ. Critically, it is not a probability distribution over θ — it is a surface over parameter space that quantifies how consistent each candidate parameter value is with the data that actually appeared. The MLE θ-hat is the value that maximizes L(θ | data) across all admissible parameter values.

Because likelihoods are products of many small probabilities — one factor per observation — working with the log-likelihood ℓ(θ) = log L(θ | data) is universally preferred. Products become sums, numerical underflow is eliminated, and derivatives are far easier to compute. The MLE is unchanged: the θ that maximizes ℓ(θ) also maximizes L(θ), because log is a strictly monotone transformation.

Score function and Fisher information

The score function S(θ) = ∂ℓ/∂θ is the first derivative of the log-likelihood with respect to the parameter. Setting S(θ) = 0 and solving yields the MLE for regular parametric families where the likelihood is smooth and unimodal. For the binomial family with n trials and k events, ℓ(p) = k log(p) + (n−k) log(1−p) (excluding the constant binomial coefficient). The score is S(p) = k/p − (n−k)/(1−p), and setting it to zero gives k(1−p) = (n−k)p, which simplifies to p-hat = k/n — the sample proportion.

The observed Fisher information I(θ) = −∂²ℓ/∂θ² is the negative second derivative of the log-likelihood evaluated at θ-hat; it measures how sharply peaked the log-likelihood surface is at its maximum. A steep peak (large I) signals that the data are highly informative about θ; a flat peak (small I) means substantial uncertainty persists. The Cramér-Rao lower bound states that the minimum possible asymptotic variance of any unbiased estimator is 1/I(θ₀). MLE achieves this bound in large samples, making it asymptotically efficient. In practice, SE = sqrt(1/I(θ-hat)) is the Wald standard error reported alongside every regression coefficient.

Asymptotic properties: consistency and normality

Under standard regularity conditions — correctly specified model, identically distributed observations, smooth likelihood — the MLE has two celebrated large-sample properties. (1) Consistency: θ-hat converges in probability to the true θ₀ as n grows without bound. (2) Asymptotic normality: sqrt(n)(θ-hat − θ₀) converges in distribution to a normal distribution with mean 0 and variance 1/I(θ₀). These properties underpin Wald confidence intervals: θ-hat ± z * SE, where SE = sqrt(1/I(θ-hat)) and z is the appropriate normal quantile. The quality of this large-sample approximation improves with n but can deteriorate badly at small n, near boundary values (probability near 0 or 1, variance near 0), or under serious model misspecification.

The trinity of tests: Wald, likelihood-ratio, and score

Three hypothesis tests emerge naturally from the likelihood framework and are collectively known as the "trinity" of classical tests.

Wald test: (θ-hat − θ₀) / SE(θ-hat) ~ N(0,1) under H₀. Reported as the z-statistic (or t-statistic with df adjustment) in nearly every regression output table. Simple to compute once the MLE and SE are available. Wald confidence intervals are symmetric: θ-hat ± z * SE.

Likelihood-ratio (LR) test: 2[ℓ(θ-hat) − ℓ(θ₀)] ~ χ²(df) under H₀, where df equals the number of parameters being tested. Compares the log-likelihood at the unconstrained MLE to its value constrained to the null. The LR test is generally preferred to the Wald test because it uses the full shape of the log-likelihood surface rather than only local curvature at the peak. Near boundary values and with sparse data, the log-likelihood surface is asymmetric, and the normal approximation underlying the Wald test deteriorates; the LR test remains better calibrated. Deviance differences between nested GLMs are LR statistics.

Score (Rao) test: evaluates the slope S(θ₀) at the null parameter value without computing the unrestricted MLE — advantageous when fitting the full model is computationally expensive.

Profile likelihood confidence intervals invert the LR test: they include all θ values for which 2[ℓ(θ-hat) − ℓ(θ)] is below the chi-squared critical value. Profile CIs respect the natural curvature and asymmetry of the log-likelihood, producing intervals that need not be symmetric. The Clopper-Pearson exact interval for a binomial proportion is the textbook profile CI. Near boundary values, with sparse data, or when separation warnings appear in logistic regression, profile CIs should replace Wald CIs, and Firth penalized likelihood extends this logic further (→ firth-penalized-regression-rwe).

Deviance, AIC, and BIC: likelihood-based model comparison

The deviance D = −2ℓ(θ-hat) is the scaled negative log-likelihood. The difference in deviance between two nested models — 2[ℓ(larger) − ℓ(smaller)] — is a likelihood-ratio statistic distributed approximately as χ²(Δdf), providing a formal test of whether the additional parameters improve fit beyond chance. This is the standard test for adding a covariate or interaction term to a GLM or survival model.

The Akaike Information Criterion AIC = −2ℓ + 2k penalizes the number of estimated parameters k to guard against overfitting. The Bayesian Information Criterion BIC = −2ℓ + k log(n) applies a stronger complexity penalty that grows with sample size. Both are directly comparable across models fit to the same dataset by maximum likelihood; lower values indicate a better penalized fit. AIC and BIC cannot be meaningfully compared across datasets with different n, different response distributions, or transformed responses. AIC tends to favor more complex models and is preferred for prediction; BIC penalizes complexity more strongly and is preferred for model identification.

Where MLE quietly powers the catalog

MLE is the estimation engine behind virtually every model used in RWE and HEOR. Logistic regression maximizes the binomial log-likelihood over a logit-linear predictor. Poisson regression maximizes the Poisson log-likelihood; negative-binomial models add an overdispersion parameter estimated jointly by MLE. The Cox proportional hazards model uses a partial likelihood — a MLE variant that conditions out the unspecified baseline hazard — to estimate log-hazard ratios without parametrizing the hazard shape. Accelerated failure time (AFT) models maximize parametric likelihoods (Weibull, log-normal, log-logistic). Linear mixed models and GLMMs integrate over random effects to form a marginal likelihood, which is then maximized by REML or full ML. Understanding MLE therefore unifies interpretation of standard errors, confidence intervals, model comparison statistics, and test procedures across all these methods (→ generalized-linear-models for the GLM subfamily).

Failure modes and protections

Complete separation: in logistic regression, when a predictor perfectly predicts the outcome in one direction, the score equation has no finite solution — the log-likelihood surface is monotone and the MLE diverges to ±∞. Software may silently report implausibly large coefficients accompanied by enormous standard errors, or may fail to converge. Firth penalized likelihood (→ firth-penalized-regression-rwe) adds a Jeffreys prior term that regularizes the log-likelihood surface, producing finite bias-reduced estimates.

Boundary parameters: variance components estimated at zero in a mixed model, or a probability estimated at exactly 0 or 1, signal near-identification failure. Profile CIs and exact methods (→ exact-methods-sparse-data-rwe) are the recommended safeguards.

Model misspecification: if the assumed distribution is wrong — Poisson when overdispersion is present, or binomial when observations are correlated — the MLE remains consistent for the working-model parameter, but standard likelihood-based standard errors are invalid because they assume correct model specification. The sandwich (robust / HC) standard error replaces the inverse information matrix with a "meat-bread" estimator that is consistent under misspecification. Quasi-likelihood extends this idea, allowing the mean-variance relationship to be estimated without specifying a full distributional family. The sandwich SE is the first safeguard; a correctly specified model — negative-binomial for overdispersed counts, a mixed model for clustered data — is always preferable when the data structure is understood (→ regression-diagnostics for AIC/BIC-based selection of distributional family).

Pros, cons, and trade-offs

Pros: MLE is asymptotically efficient (minimum-variance among consistent estimators under correct specification), universally implemented across statistical software, and provides a coherent unified basis for point estimation, the full trinity of hypothesis tests, confidence intervals (Wald or profile), and likelihood-based model comparison (AIC/BIC/deviance). The framework scales from simple closed-form cases (binomial proportion, normal mean) to complex numerical optimization problems (mixed models, survival models, custom likelihoods).

Cons: MLE requires correct model specification to be unbiased in finite samples. Asymptotic guarantees may not hold with small n, sparse cells, or boundary parameter values. Wald standard errors and symmetric CIs are unreliable near boundaries; separation causes divergence in logistic models. MLE ignores prior information (contrast with Bayesian MAP estimation, which incorporates a prior and does not diverge under separation).

Key trade-offs: Profile likelihood CIs are more accurate than Wald CIs near boundaries but require repeated model fitting; sandwich SEs protect against misspecification at some cost of efficiency when the model is correctly specified; the LR test outperforms the Wald test at small n and near boundaries but requires fitting both restricted and unrestricted models.

When to use

MLE is the default estimation method for all parametric regression and survival models in RWE/HEOR: logistic regression for binary outcomes, Poisson or negative-binomial for counts, Cox partial likelihood for time-to-event outcomes, GLMs for any exponential-family outcome, and mixed models for clustered or longitudinal data. Use likelihood-ratio tests for nested model comparison — adding a covariate, testing an interaction, or comparing distributional families. Use profile likelihood CIs when Wald CIs may be unreliable: sparse cells, small n, parameters near 0 or 1, or any setting where a separation or convergence warning appears. Use AIC and BIC for non-nested model selection across link functions, distributional families, or covariate subsets.

When NOT to use

Do not rely on Wald tests and symmetric Wald CIs when parameters are near the boundary of their parameter space — log-odds near ±∞, variance near 0, probability near 0 or 1 — use profile likelihood CIs or exact methods instead. Do not ignore separation warnings in logistic regression: a very large estimated coefficient paired with a very large standard error signals that the MLE has not converged to a finite value; apply Firth penalization rather than reporting the diverged estimate. Do not use standard likelihood-based standard errors when the model is known to be misspecified — clustered data analyzed without cluster-level terms, overdispersed count data fit with a plain Poisson model — report sandwich SEs or re-specify the model. For very small samples (fewer than approximately 20 observations), consider exact methods or a Bayesian analysis with informative priors rather than relying on the large-sample normal approximation underlying Wald inference.

Interpreting the output

In the worked example, seven events are observed in ten independent Bernoulli trials. The MLE for the event probability is p-hat = 7/10 = 0.7. The binomial likelihood at p = 0.7 is approximately 0.267, versus approximately 0.117 at p = 0.5 — the observed data are more than twice as probable under p = 0.7 as under p = 0.5. The Wald variance is 0.7 * 0.3 / 10 = 0.021, giving SE approximately 0.145 and a Wald 95% CI of approximately [0.42, 0.98]. With n = 10, the profile likelihood (Clopper-Pearson exact) CI of approximately [0.35, 0.93] is more reliable because the Wald interval assumes a symmetric normal approximation for p-hat that is questionable at this sample size.

(1) Formal interpretation. p-hat = 0.7 is the value of the event probability that maximizes the binomial log-likelihood given 7 events in 10 trials. The Wald 95% CI [0.42, 0.98] is a repeated-sampling interval: in approximately 95% of repeated experiments of the same design, a CI constructed by this procedure would contain the true parameter value. The specific interval [0.42, 0.98] is one realization of this procedure — it does not mean the true p lies within this interval with 95% probability. The likelihood-ratio statistic comparing p-hat = 0.7 to the null p = 0.5 is approximately 1.64 on 1 degree of freedom, well below the chi-squared critical value of 3.84, so p = 0.5 cannot be excluded at alpha = 0.05 given only 10 trials. The MLE is a conditional statement about the parameter given the observed data; causal language is not appropriate unless the study design supports a causal estimand.

(2) Practical interpretation. If 7 of 10 pilot patients responded to a treatment, the best single estimate of the true response rate is 70%. But the data are sparse: the 95% CI spans from roughly 42% to 98%, meaning the true rate could plausibly be near a coin flip or near certainty. A decision-maker should not treat 70% as a precise characterization of treatment efficacy — the sample is too small. A study of approximately 100 or more patients would be needed to narrow the CI enough for pricing or coverage decisions, and profile likelihood CIs should be used throughout because the normal approximation is unreliable at n = 10.

Worked example

Scenario

A clinical trialist runs a small pilot study of 10 independent patients. Each patient either has the primary event (outcome = 1) or does not (outcome = 0). Seven patients have the event. The trialist wants to estimate the true event probability p using maximum likelihood, confirm that p = 0.7 is the MLE, compare the likelihood at p = 0.7 versus the null hypothesis value p = 0.5, and compute a Wald 95% confidence interval for the estimate.

Dataset

Raw binary outcomes for 10 independent trial participants. Seven have the event (outcome = 1); three do not (outcome = 0). MLE is computed on these 10 observations.

Steps

• Count events and trials: 7 events observed in 10 independent trials. The analytic MLE for the binomial probability parameter is the sample proportion: p_hat = 7/10 = 0.7.

• Confirm by calculus. The log-likelihood (excluding the constant log(C(10,7))) is l(p) = 7log(p) + 3log(1-p). The score is dl/dp = 7/p - 3/(1-p). Setting the score to zero gives 7(1-p) = 3p, so 7 - 7p = 3p, so 7 = 10*p, confirming p_hat = 7/10 = 0.7 as the exact maximum of the log-likelihood.

• Compare likelihoods. The binomial coefficient C(10,7) = 120 appears in L(p) = 120 p^7 (1-p)^3 at every value of p and cancels in any ratio, so it does not affect which p is the MLE. At p=0.7 the likelihood L(0.7) is approximately 0.267; at p=0.5 the likelihood L(0.5) is approximately 0.117. The observed data are more than twice as probable under p=0.7 as under p=0.5.

• Compute the Wald variance and standard error. The variance of the MLE p-hat is p_hat(1-p_hat)/n. Variance = 0.7 0.3 / 10 = 0.021. SE = sqrt(0.021) approximately 0.145.

• Construct the Wald 95% CI: p_hat +/- 1.96 * SE approximately 0.7 +/- 0.284, giving approximately [0.42, 0.98]. Because n=10 is small and p-hat=0.7 is not close to 0.5, a profile likelihood (Clopper-Pearson exact) CI of approximately [0.35, 0.93] is more reliable here — it correctly captures the asymmetry of the log-likelihood at this sample size and does not extend above 1.0 as Wald can when p-hat is large.

Result

p_hat = 7/10 = 0.7. Variance of p_hat = 0.7 * 0.3 / 10 = 0.021. Wald SE approximately 0.145. Wald 95% CI approximately [0.42, 0.98]. L(0.7) approximately 0.267 versus L(0.5) approximately 0.117: p=0.7 makes the observed data more than twice as probable as p=0.5. The LR statistic comparing p=0.7 to p=0.5 is approximately 1.64 on 1 degree of freedom, below the chi-squared critical value of 3.84, so p=0.5 cannot be excluded at alpha=0.05 with n=10. The Clopper-Pearson profile CI of approximately [0.35, 0.93] is preferred.

Runnable example

import math
import numpy as np
from scipy import optimize, stats
from scipy.stats import beta as beta_dist
import statsmodels.formula.api as smf
import pandas as pd

# ── Part 1: Binomial MLE — matches worked example (7 events in 10 trials) ──
n_trials, n_events = 10, 7

def neg_log_lik(p):
    """Negative binomial log-likelihood; C(n,k) constant omitted (doesn't change MLE)."""
    if p <= 0 or p >= 1:
        return 1e10
    return -(n_events * math.log(p) + (n_trials - n_events) * math.log(1 - p))

result = optimize.minimize_scalar(neg_log_lik, bounds=(1e-9, 1 - 1e-9), method="bounded")
p_hat = result.x
print(f"MLE p-hat = {p_hat:.4f}  (analytic: {n_events}/{n_trials} = {n_events/n_trials})")

# Likelihood comparison at p=0.7 vs p=0.5 (C(10,7)=120 appears in both; cancels in ratio)
binom_coef = math.comb(10, 7)                    # = 120
L_07 = binom_coef * (0.7 ** 7) * (0.3 ** 3)
L_05 = binom_coef * (0.5 ** 10)
print(f"L(0.7) = {L_07:.4f},  L(0.5) = {L_05:.4f},  ratio = {L_07/L_05:.2f}")

# Wald SE and 95% CI
wald_se = math.sqrt(p_hat * (1 - p_hat) / n_trials)
wald_lo, wald_hi = p_hat - 1.96 * wald_se, p_hat + 1.96 * wald_se
print(f"Wald SE = {wald_se:.4f},  Wald 95% CI = [{wald_lo:.3f}, {wald_hi:.3f}]")

# Profile likelihood (Clopper-Pearson exact) CI — preferred at small n or near boundaries
cp_lo = beta_dist.ppf(0.025, n_events,     n_trials - n_events + 1)
cp_hi = beta_dist.ppf(0.975, n_events + 1, n_trials - n_events)
print(f"Clopper-Pearson exact CI = [{cp_lo:.3f}, {cp_hi:.3f}]")

# LR statistic vs H0: p = 0.5 (chi-sq 1 df)
ll_07   = n_events * math.log(0.7) + (n_trials - n_events) * math.log(0.3)
ll_05   = n_events * math.log(0.5) + (n_trials - n_events) * math.log(0.5)
lr_stat = 2 * (ll_07 - ll_05)
lr_p    = stats.chi2.sf(lr_stat, df=1)
print(f"LR statistic = {lr_stat:.3f},  p = {lr_p:.4f}  (H0: p=0.5)")

# SE via observed Fisher information from Hessian
opt_h   = optimize.minimize(neg_log_lik, x0=[0.5], method="L-BFGS-B",
                            bounds=[(1e-9, 1-1e-9)], options={"gtol": 1e-12})
# Note: scipy finite-diff Hessian; use numdifftools for analytic Hessian in production
fisher_info = 1 / (p_hat * (1 - p_hat) / n_trials)    # analytic Fisher info for binomial
print(f"SE from analytic Fisher information: {math.sqrt(1/fisher_info):.4f}")

# ── Part 2: Logistic regression via statsmodels (MLE via IRLS / Newton-Raphson) ──
np.random.seed(42)
df = pd.DataFrame({"x": np.random.normal(size=200),
                   "y": np.random.binomial(1, 0.4, size=200)})
fit = smf.logit("y ~ x", data=df).fit(disp=False)
print(f"\nLog-likelihood at MLE:  {fit.llf:.4f}")
print(f"Null log-likelihood:    {fit.llnull:.4f}")
print(f"LR statistic (2*delta): {-2*(fit.llnull - fit.llf):.4f}")
print(f"AIC = {fit.aic:.4f},  BIC = {fit.bic:.4f}")
# Profile likelihood CIs — more accurate than Wald near boundaries
print("\nProfile likelihood 95% CIs:")
print(fit.conf_int(method="profile"))
print("\nWald 95% CIs (symmetric, may be unreliable near |coef| >> 0):")
print(fit.conf_int(method="normal"))

# ── Part 1: Binomial MLE (matches worked example: 7 events in 10 trials) ──
n_trials <- 10L; n_events <- 7L

# Negative log-likelihood (excluding constant log(C(n,k)))
neg_log_lik <- function(p) {
  if (p <= 0 || p >= 1) return(1e10)
  -(n_events * log(p) + (n_trials - n_events) * log(1 - p))
}

res   <- optimize(neg_log_lik, interval = c(1e-9, 1 - 1e-9))
p_hat <- res$minimum
cat(sprintf("MLE p-hat = %.4f  (analytic: %d/%d = %.4f)\n",
            p_hat, n_events, n_trials, n_events / n_trials))

# Likelihood at p=0.7 vs p=0.5
L_07 <- choose(10, 7) * 0.7^7 * 0.3^3
L_05 <- choose(10, 7) * 0.5^10
cat(sprintf("L(0.7) = %.4f,  L(0.5) = %.4f,  ratio = %.2f\n", L_07, L_05, L_07/L_05))

# Wald SE and 95% CI
wald_se <- sqrt(p_hat * (1 - p_hat) / n_trials)
cat(sprintf("Wald SE = %.4f,  95%% CI = [%.3f, %.3f]\n",
            wald_se, p_hat - 1.96*wald_se, p_hat + 1.96*wald_se))

# Clopper-Pearson exact (profile likelihood) CI — preferred at small n
cp_lo <- qbeta(0.025, n_events,     n_trials - n_events + 1)
cp_hi <- qbeta(0.975, n_events + 1, n_trials - n_events)
cat(sprintf("Clopper-Pearson exact CI = [%.3f, %.3f]\n", cp_lo, cp_hi))

# LR statistic vs H0: p = 0.5
ll_07   <- n_events * log(0.7) + (n_trials - n_events) * log(0.3)
ll_05   <- n_events * log(0.5) + (n_trials - n_events) * log(0.5)
lr_stat <- 2 * (ll_07 - ll_05)
lr_p    <- pchisq(lr_stat, df = 1, lower.tail = FALSE)
cat(sprintf("LR statistic = %.3f,  p = %.4f  (H0: p=0.5)\n", lr_stat, lr_p))

# SE from observed Fisher information via Hessian (optim with hessian=TRUE)
opt_h    <- optim(par = 0.5, fn = neg_log_lik, method = "L-BFGS-B",
                  lower = 1e-9, upper = 1 - 1e-9, hessian = TRUE)
se_hess  <- sqrt(1 / opt_h$hessian[1, 1])
cat(sprintf("SE from observed Fisher information (Hessian): %.4f\n", se_hess))

# ── Part 2: Logistic regression via glm() (MLE by IRLS) ──
set.seed(42)
df   <- data.frame(x = rnorm(200), y = rbinom(200, 1, 0.4))
fit  <- glm(y ~ x, data = df, family = binomial(link = "logit"))
fit0 <- glm(y ~ 1, data = df, family = binomial)

cat("\nLogistic regression coefficients:\n")
print(summary(fit)$coefficients)
cat(sprintf("Log-likelihood at MLE: %.4f\n", logLik(fit)))
cat(sprintf("AIC = %.4f,  BIC = %.4f\n", AIC(fit), BIC(fit)))

# LR test via anova() — preferred over Wald z-test for adding / dropping a covariate
lr_glm <- anova(fit0, fit, test = "Chisq")
print(lr_glm)
cat("Note: anova() LR test preferred over Wald z for nested model comparison.\n")

# Profile likelihood CIs via MASS::confint.glm (inverts LR test)
# Uncomment if MASS is available:
# library(MASS); confint(fit)   # profile CIs — asymmetric; more reliable near boundaries