Worked Likelihood Inference for a Generalized Model

From a Model to an Inference

Likelihood inference is most useful when its steps are connected: specify a probability model, maximize its likelihood, quantify the curvature around the maximum, and compare the fitted likelihoods required by a scientific question. This tutorial carries out those steps for grouped-binomial logistic regression. The calculations also show why the score equations, information matrix, and likelihood-ratio statistic are different summaries of the same fitted model.

Suppose \(Y_i\) is the number of successes among \(m_i\) independent Bernoulli trials in group \(i\), for \(i=1,\ldots,g\). Assume \(Y_i\sim\operatorname{Binomial}(m_i,p_i)\), independently across groups, and use the logistic model \(\operatorname{logit}(p_i)=\eta_i=x_i^{T}\beta\), where \(x_i\) is a known design vector and \(\beta\) is the coefficient vector. The trials are independent, the design is fixed, and all fitted probabilities considered below are strictly between zero and one.

Definition: For this grouped-binomial logistic model, the log-likelihood, up to no omitted parameter-dependent terms, is \( \ell(\beta)=\sum_{i=1}^{g}\{\log {m_i\choose y_i}+y_i x_i^T\beta-m_i\log[1+\exp(x_i^T\beta)]\} \). The score is \(U(\beta)=\partial\ell(\beta)/\partial\beta\). The maximum-likelihood estimate \(\hat\beta\), when it exists, maximizes \(\ell(\beta)\) and solves \(U(\hat\beta)=0\).

Writing \(p_i=\exp(\eta_i)/[1+\exp(\eta_i)]\), differentiation gives \(U(\beta)=\sum_i x_i(y_i-m_i p_i)\). The expected Fisher information is \(I(\beta)=\sum_i m_i p_i(1-p_i)x_i x_i^T\). These expressions will be calculated directly rather than treated as black-box output.

Worked Logistic-Regression Case Study

Consider two independent groups, with \(Y_0=8\) successes among \(m_0=20\) controls and \(Y_1=14\) successes among \(m_1=20\) treated individuals. Let \(T_i=0\) for controls and \(T_i=1\) for treated individuals, and specify \(\operatorname{logit}(p_i)=\alpha+\beta T_i\). Here \(\beta\) is the treatment log-odds ratio.

Worked Example: Solving the Score Equations

The likelihood factors into the two group binomial likelihoods. Its parameter-dependent log-likelihood is \( \ell(\alpha,\beta)=8\alpha-20\log(1+e^\alpha)+14(\alpha+\beta)-20\log(1+e^{\alpha+\beta}) \), plus the constant log-binomial coefficients.

The score for \(\alpha\) is \(8-20p_0+14-20p_1\), and the score for \(\beta\) is \(14-20p_1\), where \(p_0=\operatorname{logit}^{-1}(\alpha)\) and \(p_1=\operatorname{logit}^{-1}(\alpha+\beta)\). Setting the second score to zero gives \(p_1=14/20=0.7\). Substituting into the first gives \(p_0=8/20=0.4\). Both probabilities are interior, so the corresponding finite coefficient estimates exist.

Thus \(\hat\alpha=\log(0.4/0.6)=\log(2/3)\approx-0.4055\), and \(\hat\alpha+\hat\beta=\log(0.7/0.3)=\log(7/3)\approx0.8473\). Therefore \(\hat\beta=\log(7/3)-\log(2/3)=\log(7/2)\approx1.2528\). Exponentiating gives an estimated odds ratio of \(e^{\hat\beta}=3.5\). This verifies the score equations by substitution: the fitted expected success counts are \(20(0.4)=8\) and \(20(0.7)=14\).

Information and a Standard Error

At the fitted values, the information matrix for \((\alpha,\beta)\) has a particularly transparent form. A control contributes information only about \(\alpha\); a treated observation contributes to both coefficients because its linear predictor is \(\alpha+\beta\).

Worked Example: Information-Based Uncertainty

The control weight is \(20(0.4)(0.6)=4.8\), and the treated weight is \(20(0.7)(0.3)=4.2\). Hence

I(\hat\alpha,\hat\beta)= \begin{pmatrix} 4.8+4.2 & 4.2\\ 4.2 & 4.2 \end{pmatrix} = \begin{pmatrix} 9 & 4.2\\ 4.2 & 4.2 \end{pmatrix}.

The determinant is \(9(4.2)-(4.2)^2=37.8-17.64=20.16\), so \(I^{-1}=(1/20.16)\begin{pmatrix}4.2&-4.2\\-4.2&9\end{pmatrix}\). The estimated variance of \(\hat\beta\) is the lower-right entry, \(9/20.16=0.44643\), and its standard error is \(\sqrt{0.44643}\approx0.6682\).

For illustration, a normal-approximation 95% Wald interval for \(\beta\) is \(1.2528\pm1.96(0.6682)=1.2528\pm1.3097\), or approximately \((-0.057,2.562)\). This interval is a local, curvature-based summary; it is not the profile-likelihood interval constructed later.

Result: If the design matrix has full column rank and all \(p_i\) lie strictly between zero and one, then \(I(\beta)\) is positive definite. Indeed, for every nonzero vector \(v\), \(v^T I(\beta)v=\sum_i m_i p_i(1-p_i)(x_i^T v)^2\). Each weight is positive, and full column rank ensures at least one \(x_i^Tv\ne0\), so the sum is positive.

The inverse information is the usual large-sample covariance approximation under regularity conditions: independent observations, an identifiable interior parameter, a correctly specified smooth model, and nonsingular information. In grouped data, the binomial counts arise from independent trials within groups. The matrix calculation itself is exact at the fitted parameter; interpreting it as a sampling variance is an asymptotic step.

Nested-Model Likelihood-Ratio Test

Compare the two-coefficient model with the nested null model \(H_0:\beta=0\), which assigns the same success probability to both groups. Under the null, the maximum-likelihood probability is the pooled proportion \(\hat p_0=(8+14)/(20+20)=22/40=0.55\). Under the alternative, the fitted probabilities are \(0.4\) and \(0.7\).

Worked Example: Comparing the Nested Models

The likelihood-ratio statistic is twice the alternative log-likelihood minus twice the null log-likelihood. The binomial coefficients cancel because both models use the same observations. Substitution gives

\begin{aligned} G^2 &=2\left[ 8\log\left(\frac{0.4}{0.55}\right) +12\log\left(\frac{0.6}{0.45}\right) +14\log\left(\frac{0.7}{0.55}\right) +6\log\left(\frac{0.3}{0.45}\right) \right]\\ &\approx 2[-2.5476+3.4522+3.3763-2.4328]\\ &\approx 3.696. \end{aligned}

The null imposes one restriction, so under the usual regularity conditions the asymptotic reference distribution is \(\chi^2_1\). The approximate upper-tail probability is \(0.0546\). This is evidence against equal group probabilities, but it is not below the conventional 0.05 threshold. The calculation uses the maximized likelihoods themselves; it does not replace the null fit with the alternative estimate.

The familiar deviance form follows directly from the likelihood. For a binomial group, the log-likelihood difference between a fitted probability \(\tilde p_i\) and the saturated estimate \(y_i/m_i\) is \(y_i\log[\tilde p_i/(y_i/m_i)]+(m_i-y_i)\log[(1-\tilde p_i)/(1-y_i/m_i)]\). Twice the negative of the sum of these differences is the binomial deviance. Thus comparing nested models by a likelihood ratio is equivalent to comparing their deviances, with the larger model's deviance subtracted from the smaller model's deviance.

For the chi-square approximation to be justified, the null parameter must be an interior point of a regular, identifiable model; the parameterization must be smooth near it; and the Fisher information must be nonsingular. The number of restrictions determines the degrees of freedom. Small samples or boundary estimates can make the approximation unreliable, even though the likelihood-ratio statistic remains defined.

Profile Likelihood and a Checked Interval

A profile likelihood for a parameter of interest retains that parameter and maximizes over nuisance parameters. In general, if the full parameter is \((\psi,\lambda)\), define \(\ell_p(\psi)=\max_\lambda\ell(\psi,\lambda)\). Under differentiability and an interior, differentiable nuisance maximizer \(\hat\lambda_\psi\), the chain rule gives \(d\ell_p(\psi)/d\psi=\partial\ell(\psi,\hat\lambda_\psi)/\partial\psi\): the term involving \(d\hat\lambda_\psi/d\psi\) vanishes because the nuisance score is zero at its conditional maximum. This explains how the profile curve follows the best-fitting nuisance values as \(\psi\) varies.

The next calculation uses an intercept-only binomial logistic model. It has no nuisance parameter, so its profile likelihood is the ordinary likelihood as a function of the success probability. This simple case makes the likelihood-ratio cutoff and numerical root check explicit.

Worked Example: A Profile-Likelihood Interval for a Binomial Probability

Suppose \(Y=20\) successes are observed in \(n=80\) independent Bernoulli trials, with \(Y\sim\operatorname{Binomial}(80,p)\). The maximum-likelihood estimate is \(\hat p=20/80=0.25\). The log-likelihood-ratio statistic at a candidate \(p\) is

D(p)=2\left[ 20\log\left(\frac{0.25}{p}\right) +60\log\left(\frac{0.75}{1-p}\right) \right].

For a 95% likelihood-ratio interval, retain values satisfying \(D(p)\leq\chi^2_{1,0.95}\approx3.841\). Numerically solving \(D(p)=3.841\) gives endpoints approximately \(0.16396\) and \(0.35181\). At the lower endpoint, the terms inside the brackets are approximately \(20\log(0.25/0.16396)=8.437\) and \(60\log(0.75/0.83604)=-6.516\). Their sum is about \(1.921\), so \(D(0.16396)\approx3.841\), consistent with the cutoff up to rounding. At \(p=0.3970\), the two terms are approximately \(-9.24951\) and \(13.08936\), giving \(D(0.3970)\approx2(3.83985)=7.67971\)—which is not the cutoff, so this proposed upper value must be checked rather than accepted.

At \(p=0.343\), \(20\log(0.25/0.343)\approx-6.3254\) and \(60\log(0.75/0.657)\approx7.9433\), for a sum of about \(1.6179\), hence \(D\approx3.2359\). This falls below the cutoff. At \(p=0.360\), the terms are approximately \(-7.2929\) and \(9.5163\), giving \(D\approx4.4469\). The upper root is therefore between \(0.343\) and \(0.360\), approximately \(0.35181\). The interval is approximately \((0.164,0.352)\).

A likelihood-ratio interval is obtained by inverting the likelihood-ratio statistic, not by taking a fixed number of standard errors on either side of the estimate. Its endpoints need not be symmetric around the estimate, especially when the likelihood is skewed. Profiling is particularly useful when nuisance parameters are present or when a parameter transformation makes a symmetric approximation on the original scale questionable.

Common Mistakes and a Practical Sequence

Using fitted probabilities from the wrong model. For a nested likelihood-ratio test, calculate each likelihood at that model's own maximum, including the restricted fit under the null.
Confusing observed curvature with expected information. They are related in regular likelihood theory, but their formulas are not interchangeable in every finite sample. State which information matrix supplies the standard error.
Reporting a cutoff without checking roots. Evaluate the likelihood-ratio statistic at each proposed endpoint. Small arithmetic errors can move an endpoint or make an alleged interval inconsistent with its defining inequality.
Applying a chi-square reference distribution automatically. The standard approximation requires regularity, including an interior null value and nonsingular information; boundary or nonexistence issues require additional care.

Specify the model.
Write the distribution, linear predictor, design, and parameter restriction being considered.

Fit each required model.
Derive or solve the score equations and verify the fitted values satisfy them.

Calculate the requested inference.
Use inverse information for a curvature-based standard error, maximized likelihoods for a nested comparison, or a profiled likelihood curve for an interval.

Check numerical and theoretical conditions.
Substitute estimates and endpoints into the defining equations, then confirm that the reference approximation's assumptions are plausible.

Key takeaway: A worked likelihood analysis is a chain of verifiable calculations. Score equations identify the fit, information quantifies local uncertainty, and likelihood ratios compare fits or define intervals through their relative likelihood.

Check Your Understanding

Use the definitions and calculations above to answer the following.

In the two-group logistic model, why does the score equation for \(\beta\) determine the treated fitted probability directly?
Verify the determinant of the information matrix and identify which entry of its inverse estimates the variance of \(\hat\beta\).
Under the null model in the nested comparison, why is the fitted probability \(22/40\) rather than either group's observed proportion?
State two regularity conditions needed for the usual chi-square reference distribution of a likelihood-ratio statistic.
Explain why a candidate profile-likelihood endpoint should be substituted into the likelihood-ratio equation before it is reported.