Why Use Logistic Regression for Case-Control Association?
In a case-control association analysis, the outcome has two values: case or control. A linear regression model can produce fitted values below zero or above one, so it is not an appropriate model for the probability of being a case. Logistic regression instead models the log odds of case status as a linear function of genotype dosage and covariates. Its SNP coefficient is estimated on the log-odds scale.
Here, genotype dosage \(G_i\) is the number, or expected number, of copies of the effect allele for individual \(i\), from 0 to 2. The additive specification treats a one-copy dosage increase as the same change in log odds across the dosage range. It does not require that every person with the same dosage has the same outcome probability.
The model assumes that, conditional on the included predictors, outcomes are independent Bernoulli variables and their probabilities follow the stated logistic form. The covariate effects are linear on the log-odds scale as written; categorical covariates must be represented by appropriate indicator variables, and nonlinear effects require additional terms. For ordinary maximum-likelihood estimation, the design matrix must have full column rank, and the data must not exhibit complete or quasi-complete separation. These conditions allow the coefficients to be identified and a finite maximum-likelihood estimate to exist.
In a case-control study, participants are sampled according to outcome status. If sampling depends on case/control status but not on genotype or covariates within outcome groups, the intercept for the sampled data generally differs from the population intercept. The association slope can still be estimated from the logistic model under this sampling scheme. Thus the fitted intercept should not automatically be read as the disease prevalence in the source population.
From Log Odds to a Fitted Probability
The inverse-logit function converts a linear predictor into a probability between zero and one. Coefficients are additive changes in the log odds when a predictor increases by one unit, holding the other included predictors fixed. This tutorial focuses on how those coefficients are specified and estimated; interpretation of the corresponding per-allele odds ratio is a separate step.
Worked Example: Computing a Fitted Case Probability
Suppose a hypothetical fitted model uses dosage \(G_i\), age in years, and a binary covariate \(S_i\) coded 1 for one category and 0 for the reference category: $$ \operatorname{logit}(\pi_i)=-2.20+0.35G_i+0.040\,\text{Age}_i-0.20S_i. $$ For an individual with dosage 2, age 60 years, and \(S_i=1\), calculate the fitted probability.
First calculate the linear predictor, retaining the units and coding in the model: $$ \eta_i=-2.20+(0.35)(2)+(0.040)(60)-(0.20)(1) =-2.20+0.70+2.40-0.20=0.70. $$ The fitted probability is $$ \hat{\pi}_i=\frac{1}{1+\exp(-0.70)} =\frac{1}{1+0.4966}\approx 0.6682. $$ The model therefore assigns this hypothetical individual a fitted case probability of about 0.668, conditional on the specified predictors and the sampling context. This is a model-based probability, not a claim that the genotype causes disease.
Maximum-Likelihood Estimation
Let \(x_i\) be the row of the design matrix containing an intercept, \(G_i\), and the covariate values. Write the coefficient vector as \(\theta=(\beta_0,\beta_G,\gamma_1,\ldots,\gamma_m)^\mathsf{T}\). Then \(\eta_i=x_i^\mathsf{T}\theta\) and \(\pi_i=1/\{1+\exp(-x_i^\mathsf{T}\theta)\}\). Under the independent Bernoulli model, the likelihood is the product of the individual probabilities of their observed outcomes.
Maximum likelihood chooses the coefficient vector that maximizes \(\ell(\theta)\). In general, the logistic likelihood does not yield a closed-form coefficient formula: the probabilities depend nonlinearly on the coefficients. An iterative numerical algorithm finds the maximum. The score, the vector of first derivatives of the log likelihood, and the Hessian, the matrix of second derivatives, characterize that optimization.
Why the Score Equations Lead to Iteration
For one observation, the log-likelihood contribution is \(Y_i\eta_i-\log\{1+\exp(\eta_i)\}\). Because the derivative of \(\pi_i\) with respect to \(\eta_i\) is \(\pi_i(1-\pi_i)\), differentiating that contribution with respect to \(\eta_i\) gives \(Y_i-\pi_i\). Since the derivative of \(\eta_i\) with respect to \(\theta\) is \(x_i\), the contribution to the coefficient score is \(x_i(Y_i-\pi_i)\). Summing across observations gives \(U(\theta)=X^\mathsf{T}(y-\pi)\). Differentiating once more gives the negative curvature \(-X^\mathsf{T}WX\).
Newton–Raphson uses the local curvature to update the coefficient estimates. When the matrix \(X^\mathsf{T}WX\) is invertible, its update can be written as $$ \theta_{\text{new}}=\theta_{\text{old}}+ (X^\mathsf{T}WX)^{-1}X^\mathsf{T}(y-\pi). $$ In practice, iteratively reweighted least squares is a closely related way of carrying out these updates. The iterations stop when changes in coefficients or in the likelihood are sufficiently small and the score is close to zero. This computational procedure is not a formula that guarantees a sensible fit: model identifiability and convergence still need to be checked.
Worked Example: Estimating the Dosage Coefficient From Grouped Counts
Consider hypothetical case-control counts at dosage values 0, 1, and 2. Suppose the data are:
| Dosage \(G\) | Cases | Controls | Total |
|---|---|---|---|
| 0 | 10 | 30 | 40 |
| 1 | 20 | 30 | 50 |
| 2 | 40 | 30 | 70 |
Fit the model \(\operatorname{logit}(\pi)=\beta_0+\beta_GG\), with no additional covariates. The observed case proportions are \(10/40=0.25\), \(20/50=0.40\), and \(40/70\approx0.5714\). Their logits are $$ \log\left(\frac{0.25}{0.75}\right)=\log(1/3)\approx-1.0986,\quad \log\left(\frac{0.40}{0.60}\right)=\log(2/3)\approx-0.4055,\quad \log\left(\frac{40/70}{30/70}\right)=\log(4/3)\approx0.2877. $$ The logit increases by approximately \(0.6931\) from dosage 0 to 1 and again by approximately \(0.6931\) from dosage 1 to 2. Thus the additive logistic model fits these three group proportions exactly with $$ \hat{\beta}_0=\log(1/3)\approx-1.0986,\qquad \hat{\beta}_G=\log(2)\approx0.6931. $$ This is a maximum-likelihood estimate, not merely a line drawn through the proportions: at these fitted probabilities, the predicted number of cases in each dosage group equals the observed number. Therefore, within each group, the case residual sum is zero. Since the model score sums those residuals and their dosage-weighted values, both score components are zero. The dosage coefficient is expressed in log odds per one-copy increase in effect-allele dosage.
Including Covariates in the Estimation
Adding covariates changes the design matrix and the likelihood optimization, not the basic meaning of the dosage coefficient. The fitted \(\beta_G\) is the model’s dosage association adjusted for the covariates as specified. Adjustment cannot correct for omitted or incorrectly coded variables automatically; it only conditions on the predictors that are actually included.
Worked Example: A Dosage Coefficient Adjusted for a Binary Covariate
Suppose a hypothetical grouped data set has the following case and control counts for each dosage \(G\) and binary covariate \(C\). The fitted model is \(\operatorname{logit}(\pi)=\beta_0+\beta_GG+\gamma C\).
| \(G\) | \(C\) | Cases | Controls |
|---|---|---|---|
| 0 | 0 | 1 | 1 |
| 1 | 0 | 2 | 1 |
| 0 | 1 | 3 | 1 |
| 1 | 1 | 6 | 1 |
The case-to-control odds in the four cells are \(1/1=1\), \(2/1=2\), \(3/1=3\), and \(6/1=6\). Under the additive model, these correspond to $$ \exp(\beta_0)=1,\quad \exp(\beta_0+\beta_G)=2,\quad \exp(\beta_0+\gamma)=3,\quad \exp(\beta_0+\beta_G+\gamma)=6. $$ Solving gives \(\hat{\beta}_0=0\), \(\hat{\beta}_G=\log(2)\approx0.6931\), and \(\hat{\gamma}=\log(3)\approx1.0986\). The probabilities predicted for the four cells are \(1/2\), \(2/3\), \(3/4\), and \(6/7\). For example, in the cell \(G=1,C=1\), the predicted case count is \(7(6/7)=6\), matching the observed six cases. The same match holds in each cell, so the residual score contributions sum to zero for the intercept, dosage, and covariate columns. Each cell contains both cases and controls, and the four predictor combinations identify the intercept and two slopes; the fit has finite estimates. The dosage coefficient is estimated while including \(C\) in the model.
Diagnostics, Failure Modes, and Alternatives
A fitted model should be checked before its coefficients are used. First confirm that the outcome is coded consistently, that dosage uses the intended effect allele, and that covariates have the intended units and reference categories. Then check that the fitting algorithm converged, that the estimated coefficients and standard errors are finite, and that the design matrix is not rank-deficient. Implausibly large coefficients, huge standard errors, or sensitivity to small data changes can signal sparse data or separation. Residual or calibration checks can reveal poor fit, although a good overall fit does not prove that confounding or selection bias is absent.
Complete separation occurs when a linear combination of predictors perfectly distinguishes cases from controls. In that situation, the likelihood can keep increasing as one or more coefficients grow without bound, so a finite maximum-likelihood estimate does not exist. Quasi-complete separation can likewise prevent finite estimates. A numerical optimizer may stop and print a result even when a coefficient is effectively diverging; convergence messages and coefficient magnitudes should therefore be examined rather than treated as proof of a valid fit.
The additive dosage model is a parsimonious choice when the scientific model assumes a constant per-copy change in log odds. A model with genotype indicators instead estimates separate contrasts for dosage 1 and dosage 2 relative to dosage 0, allowing a nonadditive pattern but using more parameters. A linear probability model is computationally simple, but can give fitted values outside the interval from zero to one and does not use the Bernoulli likelihood. These alternatives answer different modeling questions; neither removes the need to specify the outcome, genotype coding, and covariates carefully.
A common mistake is to treat the fitted logistic coefficient as a change in probability. For example, saying that a coefficient of \(0.35\) means the case probability rises by 0.35 for each additional allele is incorrect. The coefficient is a change in log odds, and the probability change depends on the starting probability and the other predictors. Another mistake is to assume that adjusting for a covariate proves the association is free of confounding. Adjustment only addresses the modeled covariates under the model’s assumptions; unmeasured or misspecified factors remain possible sources of bias.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In the model \(\operatorname{logit}(\pi_i)=\beta_0+\beta_GG_i+\gamma C_i\), what does a one-unit increase in \(G_i\) change directly?
- For a fitted linear predictor \(\eta=0\), calculate the fitted probability.
- In the grouped dosage example, why do the displayed coefficients satisfy the likelihood score equations?
- Why can a case-control sample have an intercept that does not describe disease prevalence in the source population?
- Name one warning sign of separation and state its consequence for ordinary maximum-likelihood estimation.
Solutions. (1) It changes the model’s log odds by \(\beta_G\), holding \(C_i\) fixed; it does not add a fixed amount to probability. (2) \(\pi=1/\{1+\exp(0)\}=1/2=0.5\). (3) At the displayed coefficients, fitted case counts equal observed case counts in each dosage group, so the residuals sum to zero and their dosage-weighted sum is zero. Those are the intercept and dosage score equations. (4) Sampling by case/control status changes the sampled case fraction, and therefore the sampled intercept; under outcome-dependent sampling, the source-population intercept is not recovered simply by using that sample intercept. (5) A coefficient that grows extremely large, often with a very large standard error, can indicate separation. Complete separation prevents a finite ordinary maximum-likelihood estimate.