Methods and validation

Nonlinear least squares (Levenberg-Marquardt)

Used by both tools

Every curve on this site is fitted by Levenberg-Marquardt: an iterative method that solves the damped normal equations (JᵀWJ + λ·diag(JᵀWJ))δ = JᵀWr for a parameter step, accepting steps that reduce the weighted sum of squares and increasing the damping when they do not. Scaling the damping by the diagonal of JᵀWJ rather than adding λI is Marquardt’s own refinement, and it matters for exactly these models: a plateau parameter measured in absorbance units and a dimensionless Hill slope have diagonal entries six orders of magnitude apart, and undamped-in-one-direction is indistinguishable from not converging.

Parameter standard errors are the square roots of the diagonal of σ²(JᵀWJ)⁻¹ evaluated at the optimum, where σ² = SSE/(n − p). This is the asymptotic covariance: it is exact only for a model that is linear in its parameters, which none of these are. Where that approximation is known to be poor — the confidence interval on an IC50 from a curve with no upper plateau — a profile-likelihood interval is computed and reported alongside it, and the figure states which method produced which number.

The solver is checked against the NIST Statistical Reference Datasets for nonlinear regression, whose values are certified: computed by NIST in extended precision and rounded to eleven significant figures, so they are the true answer rather than another program’s answer. Both logistic datasets are used, from NIST’s deliberately distant starting values, and the certified parameter standard deviations are checked as well as the estimates — which tests the covariance path and not merely the optimiser.

The four-parameter logistic (4PL)

Used by Dose-response (IC50)

Y = Bottom + (Top − Bottom) / (1 + 10^((LogIC50 − X) · HillSlope))     X = log₁₀(concentration)

This is the equation GraphPad Prism calls “log(inhibitor) vs. response — Variable slope”, in Prism’s own parameterisation and with Prism’s own parameter names. That choice is deliberate: it means a number produced here can be compared with a number produced there without anybody having to work out a sign convention or a change of base first. The R package drc writes the same curve with a natural logarithm and with the slope’s sign reversed, so its b is −HillSlope·ln(10) and its e is the IC50 — a mapping worth knowing if you are checking against a paper that used drc.

The fit is performed on LogIC50 rather than on IC50. Not for convenience: a concentration cannot be negative, so fitting IC50 directly means either a constrained optimiser or a fit that can wander into meaningless territory, and the log parameter is unbounded. It also makes the likelihood closer to quadratic, which is what makes the asymptotic standard error a reasonable description of the uncertainty. The consequence, which is correct and often surprises people, is that the confidence interval on the IC50 is asymmetric: it is symmetric on the log scale and back-transformed.

The IC50 confidence interval, by profile likelihood

Used by Dose-response (IC50)

the interval is {θ : SSE(θ) ≤ SSE_min · (1 + F(1, n−p, 0.95)/(n−p))}

Two intervals are reported for the IC50 and they usually disagree, so it is worth being explicit about which is which. The asymptotic interval is estimate ± t·SE on the log scale, back-transformed. It assumes the sum-of-squares surface is a paraboloid near the optimum. For a well-sampled sigmoid with plateaus at both ends, it is close to right.

The profile-likelihood interval makes no such assumption. LogIC50 is stepped away from its fitted value; at each step the OTHER parameters are re-fitted with LogIC50 held fixed; and the interval is the range over which the resulting sum of squares stays inside the F-distributed threshold above. Where the surface is genuinely curved — which is what a truncated dose-response curve produces — the two intervals differ substantially, and the profile one is the one to quote. It can also be one-sided, correctly: a curve with no upper plateau has no finite upper confidence limit for its IC50, and reporting “not determined” for that bound is more useful than reporting the value at which the optimiser happened to stop.

The Kaplan-Meier product-limit estimator

Used by Kaplan-Meier survival

Ŝ(t) = ∏ over tᵢ ≤ t of (nᵢ − dᵢ) / nᵢ

The survival curve is the running product, over every time at which an event occurred, of the fraction of the risk set that survived it. nᵢ is the number still at risk immediately before tᵢ and dᵢ the number of events at tᵢ. Censored subjects leave the risk set without producing a step, which is the entire point of the estimator: somebody lost to follow-up at 14 months contributes their fourteen months of survival and nothing about what happened afterwards.

The convention for a tie between an event and a censoring at the same recorded time is that the censoring is treated as occurring immediately after. So a subject censored at t is still counted in nᵢ for an event at t. This matters more often than it sounds — study data is routinely recorded to the day or the week, so exact ties are common rather than exceptional — and the opposite convention shifts the curve visibly.

Greenwood’s variance, and the confidence band

Used by Kaplan-Meier survival

Var(Ŝ(t)) = Ŝ(t)² · Σ over tᵢ ≤ t of dᵢ / (nᵢ(nᵢ − dᵢ))

Greenwood’s formula is the standard variance estimator for the product-limit estimate, and the standard error it gives is what the confidence band is built from. The band is not computed on the survival scale directly, though, because an interval of Ŝ ± 1.96·SE runs outside [0, 1] wherever the curve is near either end — and a survival probability of 1.04 is not a defensible thing to draw. The interval is therefore computed on a transformed scale and mapped back, which keeps it inside the unit interval by construction.

A confidence band around a survival curve is pointwise: it is a 95% interval for the survival probability at each individual time, not a 95% region containing the whole true curve. The probability that the true curve stays inside the band everywhere is considerably less than 95%. This is a widely-made reading error and it is why the figure labels the band rather than leaving it to be inferred.

The log-rank test

Used by Kaplan-Meier survival

At each event time the test compares the events observed in each group with the number expected if survival were identical across groups, given the sizes of the risk sets at that moment. Summing those differences and standardising by their covariance gives a statistic that is chi-square distributed with one fewer degree of freedom than there are groups.

The test asks one question — are these curves the same? — and answers only that. It gives no effect size, and it has most of its power when the hazard ratio is roughly constant over time. Two curves that cross can produce a thoroughly non-significant log-rank p-value while being obviously and importantly different, because the early and late differences cancel in the sum. A p-value from this test should always be read next to the curves, never instead of them.

The Gehan-Breslow-Wilcoxon test

Used by Kaplan-Meier survival

U = Σ over event times of nₜ·(d₁ₜ − dₜ·n₁ₜ/nₜ)

The same accumulation as the log-rank, with each event time weighted by the number of subjects still at risk at that moment rather than weighted equally. Early times, when the risk set is large, therefore count for more — which makes this the more powerful test when the curves separate early and converge later, and the less powerful one when a constant hazard ratio holds throughout. Which of the two to use is a decision about where the difference is expected to be, and it has to be made before either p-value is seen; both are printed on the figure precisely so that neither can be quietly swapped in after the fact.

That weight is not an arbitrary choice of emphasis. With wₜ = nₜ the score collapses to an exact identity: it becomes Gehan’s W statistic, the sum over every pair of subjects — one from each group — of +1 when the first demonstrably outlived the second, −1 for the reverse, and 0 when the censoring pattern leaves the pair unorderable. On uncensored data that is the ordinary Wilcoxon rank-sum test, which is where the name comes from. The identity is what this implementation is validated against: Gehan’s 1965 paper works his test through on the 6-MP data used elsewhere on this page and prints W = 335 − 64 = 271, which the score reproduces exactly.

The variance is a different matter, and the distinction matters enough to state plainly. Gehan’s original refers W to a permutation distribution: the variance of W over all the ways the observed pattern of failures and censorings could have been split between two groups of the given sizes. What this tool computes — and what “Gehan-Breslow-Wilcoxon” denotes in Prism and in SAS PROC LIFETEST — is the sum of per-event-time hypergeometric variances, the same covariance the log-rank uses with the weights squared into it. The two are different null distributions for the same statistic, not two attempts at one number, and they do not agree to the printed digits: on Gehan’s own worked example his variance gives chi-square 13.01 and this one gives 13.46.