Saddlepoint Calibration for Sparse Binary-Phenotype Tests

Why Sparse Burden Scores Need Tail Calibration

A burden score for a binary phenotype combines genotype information with case-control residuals. When cases are numerous and the score aggregates many small contributions, a normal approximation can be useful. With rare variants, however, only a few people may carry the alleles contributing to the score. The resulting null distribution can be strongly asymmetric, concentrated on a few possible values, or have a tail very different from the normal tail with the same variance. This matters most when a test is compared with a small significance threshold: a modest tail-probability error can change the conclusion.

Saddlepoint calibration uses the full cumulant-generating function of the score under the null model, rather than only its mean and variance. It can improve tail estimates when the score is skewed, provided the Bernoulli model is appropriate and the relevant tail is not beyond the information in the data. The method does not make a small sample behave like a large one; sparse support can still limit its accuracy.

Null Score and Its Cumulant-Generating Function

Condition on the observed genotypes, variant-set definition, weights, and covariates. Under the null logistic model, let $Y_i$ be the binary phenotype for individual $i$, with $Y_i\sim\operatorname{Bernoulli}(\pi_i)$, where $\pi_i$ is the null-model probability given that individual’s covariates. Assume outcomes are independent conditional on the design and fitted null probabilities. For a burden score with fixed coefficient $a_i$ for individual $i$, write $U=\sum_{i=1}^n a_i(Y_i-\pi_i)$. The coefficients can reflect the genotype burden and any covariate adjustment used to construct the score. The null mean is zero, and the variance is $\sum_i a_i^2\pi_i(1-\pi_i)$.

Definition: The cumulant-generating function (CGF) of the centered score $U$ is $K(t)=\log E\{\exp(tU)\}$. Under independent Bernoulli outcomes with fixed $a_i$ and $\pi_i$, $$ K(t)=\sum_{i=1}^n \log\{1-\pi_i+\pi_i\exp(t a_i)\}-t\sum_{i=1}^n a_i\pi_i. $$ Its first derivative $K'(t)$ is the mean under exponential tilting, and its second derivative $K''(t)$ is the variance under that tilted distribution. In particular, $K'(0)=0$ and $K''(0)=\operatorname{Var}(U)$.

For an observed upper-tail score $u>0$, the saddlepoint $\hat t$ solves $K'(\hat t)=u$. If a solution exists, it is positive for an upper-tail observation. The equation locates an exponential tilt under which $u$ is the mean; the CGF then measures how costly that tilt is under the null. A numerical root finder is usually used, with $K'(0)=0$ and the monotonicity of $K'$ providing a useful check.

Tail Approximations and the Lattice Distinction

For a nonlattice score, the Lugannani–Rice approximation to the upper-tail probability uses $w=\operatorname{sign}(\hat t)\sqrt{2\{\hat t u-K(\hat t)\}}$ and $v=\hat t\sqrt{K''(\hat t)}$. For $\hat t>0$, the approximation is

P(U\ge u)\ \approx\ 1-\Phi(w)+\phi(w)\left(\frac{1}{v}-\frac{1}{w}\right),

where $\Phi$ and $\phi$ are the standard normal distribution and density. This is a tail approximation, not an exact probability. A score is lattice-valued only if its support lies in a shifted lattice $b+h\mathbb{Z}$ for some $h>0$; its lattice span is the largest such $h$. Integer coefficients, or coefficients that are integer multiples of a common unit, often give lattice support. Discrete coefficients alone do not guarantee it: coefficients 1 and $\sqrt{2}$ can yield support values whose differences have no common positive span.

For a lattice-valued score and an upper-tail threshold $u$ on its support, a commonly used lattice saddlepoint approximation is

P(U\ge u)\ \approx\ \frac{\exp\{K(\hat t)-\hat t u\}}{\sqrt{2\pi K''(\hat t)}} \frac{h}{1-\exp(-h\hat t)}.

The factor involving $h$ accounts for the spacing of the possible score values in the tail. Use the nonlattice formula when the support is not on a shifted lattice; do not assign a lattice span merely because the score has discrete outcomes. Both approximations rely on a meaningful saddlepoint and can be inaccurate when the support is extremely sparse or the requested tail is driven by just one or two outcomes.

Worked Example: A Sparse Binomial Burden Score

Suppose 100 independent individuals have the same null probability $\pi_i=0.01$, and the score is the centered carrier count $U=\sum_i(Y_i-0.01)=X-1$. The hypothetical observed count is $X=4$, so $u=3$. This is a lattice score with span $h=1$.

The CGF is $K(t)=100\log(0.99+0.01\exp(t))-t$. Its derivative is $K'(t)=100q(t)-1$, where $q(t)=0.01\exp(t)/(0.99+0.01\exp(t))$. Setting $K'(t)=3$ gives $q(t)=0.04$, and therefore $\exp(\hat t)=0.04(0.99)/\{0.01(0.96)\}=4.125$, so $\hat t=\log(4.125)\approx1.417$.

Substitution gives $K(\hat t)=100\log(1.03125)-1.417\approx1.661$, and $K''(\hat t)=100(0.04)(0.96)=3.840$. The lattice approximation is $\exp(1.661-1.417\times3)/\sqrt{2\pi(3.840)}\times1/\{1-\exp(-1.417)\}\approx0.0202$.

For comparison, the exact binomial tail is $P(X\ge4)=1-\sum_{j=0}^3 {100\choose j}(0.01)^j(0.99)^{100-j}\approx0.0184$. A normal approximation with mean 1, variance 0.99, and continuity correction at 3.5 gives an upper tail of about 0.0060. Here saddlepoint calibration is much closer to the exact tail than that normal approximation. The example also shows why the lattice factor matters: the score’s possible values are spaced one unit apart.

Justification and Practical Calculation

The CGF follows directly from independence. For one Bernoulli term, $E\{\exp[t a_i(Y_i-\pi_i)]\}=\exp(-t a_i\pi_i)\{1-\pi_i+\pi_i\exp(t a_i)\}$. Multiplying these expectations over individuals and taking the logarithm gives the stated $K(t)$. Differentiating yields

K'(t)=\sum_i a_i\{q_i(t)-\pi_i\},\qquad K''(t)=\sum_i a_i^2 q_i(t)\{1-q_i(t)\},\qquad q_i(t)=\frac{\pi_i\exp(t a_i)}{1-\pi_i+\pi_i\exp(t a_i)}.

Thus the saddlepoint equation is the score equation under the tilted probabilities $q_i(t)$. Since $K''(t)\ge0$, $K'(t)$ is nondecreasing; if at least one nonzero coefficient has a nondegenerate Bernoulli outcome, it is strictly increasing over finite $t$. This supports a unique finite root whenever the target $u$ lies in the attainable interior range of $K'$. If no such root exists, the requested threshold may be outside the interior of the score’s support, and the usual saddlepoint formula should not be forced.

Specify the null distribution.
Hold the design fixed and state the null probabilities $\pi_i$, coefficients $a_i$, independence assumption, and observed score $u$.

Form the CGF and solve.
Evaluate $K(t)$ and $K'(t)$, then find $\hat t$ satisfying $K'(\hat t)=u$. Check that the root has the expected sign and that $K''(\hat t)>0$.

Choose and assess the tail formula.
Determine whether the support has a common lattice span. Apply the corresponding approximation, then compare with exact enumeration or simulation when feasible.

Worked Example: Computing the Tilt for a Weighted Score

Consider four independent hypothetical Bernoulli outcomes with $\pi_i=0.2$ and fixed coefficients $1,1,2,2$. The centered score has CGF $K(t)=2\log(0.8+0.2\exp(t))+2\log(0.8+0.2\exp(2t))-1.2t$. For the threshold $u=0.8$, solving $K'(t)=0.8$ gives $\hat t\approx0.412$.

At this value, $K(0.412)\approx0.1555$ and $K''(0.412)\approx2.247$. The deviance distance and the signed curvature term used in the nonlattice formula are $w=\sqrt{2\{(0.412)(0.8)-0.1555\}}\approx0.590$ and $v=(0.412)\sqrt{2.247}\approx0.618$.

These calculations are useful checks on the CGF and saddlepoint root, but this particular score is lattice-valued: its uncentered weighted sum lies on the integers, so its centered support has span 1. The appropriate lattice approximation is $\exp(0.1555-(0.412)(0.8))/\sqrt{2\pi(2.247)}\times1/\{1-\exp(-0.412)\}\approx0.661$. The exact upper-tail probability is $P(\text{weighted sum}\ge2)=1-0.8^2+0.8^2(0.2^2)=0.3856$: the weighted sum is at least 2 if either coefficient-2 outcome is 1 (probability $0.36$), or if neither is 1 and both coefficient-1 outcomes are 1 (probability $0.8^2\times0.2^2=0.0256$). The poor approximation is a warning, not a contradiction. With only four Bernoulli outcomes the score has very few support points, so a tail approximation based on local curvature is not reliable. Exact enumeration is preferable here.

Diagnostics, Alternatives, and Limitations

A normal approximation uses only the score’s null mean and variance. For a standardized upper-tail statistic, it maps the observed score to a normal tail using those two quantities. Saddlepoint calibration additionally uses skewness and higher-order features encoded in $K(t)$, which is why it can help for a skewed tail. It is not a replacement for checking the null model, and it is not automatically better in every small sample.

Useful diagnostics include whether $K'(t)=u$ has a finite solution, whether $K''(\hat t)$ is positive and numerically stable, and whether the resulting probability is between 0 and 1. A nonpositive or extreme probability, unstable roots, or a threshold near the endpoint of the support signals that the approximation may be failing. When the support is small, compute the tail by exact enumeration if feasible; for larger problems, a simulation under the fitted null can provide a comparison, with attention to Monte Carlo uncertainty in very small tails.

The Bernoulli-independence calculation conditions on fixed covariates and uses the specified null probabilities. In practice, those probabilities may be estimated from a logistic model, and score coefficients may incorporate nuisance projection. Plugging estimates into $K(t)$ does not make nuisance-parameter uncertainty disappear. Dependence between outcomes, misspecified probabilities, or an incorrect coefficient construction can also invalidate the stated CGF. Calibration should therefore use the same null model and score definition as the test itself.

Common mistake: “The score is discrete, so it must have a lattice span; use the lattice correction.” Discreteness alone does not establish a shifted lattice. A support generated by coefficients such as 1 and $\sqrt{2}$ need not have any positive common span. Check whether all support differences are integer multiples of one $h>0$; use the lattice correction only when that condition holds.

Saddlepoint calibration is particularly useful when the normal tail seems implausible and the CGF can be evaluated reliably. Its limitations are equally important: it is an approximation, can be poor with very few possible score values, and depends on a credible null distribution. Exact or simulation-based calibration can be preferable when computationally practical and when extreme-tail precision is essential.

Check Your Understanding

Answer each question using the null-score framework above. Brief solutions follow the questions.

For independent Bernoulli outcomes, what does $K'(t)$ represent, and what equation defines the saddlepoint for an upper-tail threshold $u$?
A score has possible values formed from independent outcomes with coefficients 1 and $\sqrt{2}$. Is a lattice correction justified solely because the score is discrete?
For the binomial example with $n=100$, $\pi_i=0.01$, and observed $X=4$, state the centered threshold and saddlepoint equation.
Why can the saddlepoint approximation be inaccurate for the four-outcome weighted example even though its CGF and root were calculated correctly?
Name two checks to perform before reporting a saddlepoint tail probability.

Solutions. (1) $K'(t)$ is the mean of the score under exponential tilting; solve $K'(\hat t)=u$. (2) No. The support must lie on a shifted lattice with a common positive span, which arbitrary real coefficients need not produce. (3) $U=X-1$, so $u=3$; solve $100q(t)-1=3$, or $q(t)=0.04$. (4) There are only four Bernoulli outcomes and very few support points, so the local tail approximation can be poor. (5) Check that a finite root exists with the expected sign, that $K''(\hat t)>0$, that the correct lattice or nonlattice formula was selected, and that the resulting probability is valid; compare with exact or simulated calibration when feasible.

Saddlepoint Calibration for Sparse Binary-Phenotype Tests

What You'll Learn