STATWI0001.02 Cover Page – page 0 2018-03-23

USDA Center for Veterinary Biologics

Statistics Section

Work Instructions

This document approved for the indicated purposes

Internal use Yes

External distribution Yes

Public web site Yes

Document: STATWI0001.02

Title: Non-parametric Estimation of Median Effective Dose

Author: David Siev

Approved by: David Siev on 2018-03-23

Version history

01 2017-01-24 First version

02 2018-03-23 Added appendix
Nonparametric Estimation of Median Effective Dose

David Siev

April 2011

Contents

1 Median effective dose 1

2 The tolerance distribution 2

3 Non-parametric estimators 3

3.1 Spearman-Karber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.2 Reed-Muench . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.3 Dragstedt-Behrens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.4 The skrmdb package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1 Median effective dose

In dose-response studies, the median effective dose (ED50) is the estimated dose that produces a response in half the population. Particular types of dose-response studies are

associated with specific forms of the ED50. In accute toxicity studies, it is the median of a

latent distribution known as the tolerance distribution (Section 2). The response is usually

death, and it is termed the median lethal dose (LD50). In vaccination-challenge studies,

a median protective dose (P D50) is estimated. It can be thought of as the median of an

immunocompetence distribution, analagous to a tolerance distribution.1

Another ‘dose 50’ that should be mentioned for completeness is the tissue culture infective

dose (T CID50). In virus titrations, the T CID50 is the dose that produces evidence of

infection in half the wells to which it is applied. It is estimated in a similar way as the

preceding measures, but unlike them the T CID50 is not the median of an underlying latent

1A distribution of immunocompentences only makes sense if the response to challenge is constant in

every individual. Since that is often implausible, the distribution is in fact a mixture of immunocompetences

and susceptibilities.

1

STATWI0001.02 page 1 of 11
2 THE TOLERANCE DISTRIBUTION 2

distribution. It is simply the volume of the virus suspension that would contain an average

of 0.7 infective virus particles (e.g. PFU).

In modern times, dose-response studies are usually handled by statistical modeling to estimate the dose-response curve. This has obvious advantages over the non-parametric

measures discussed here. One of the most salient benefits is, of course, the estimation of

the entire dose-response curve itself, rather than a single measure of one of its features, its

center. Dose-response curves were not so easily estimated in the days before computers,

and non-parametric estimators of ED50 were important historically. A glance at this document’s references will show that the ones discussed here were initially published between

1908 and 1938.

Of the three, the Spearman-Karber estimator is still somewhat useful today for its valuable

statistical properties. The Reed-Muench and Dragstedt-Behrens estimators are largely

of historical interest. Their use is still surprisingly widespread, however, so it is worth

being familiar with their mechanics. Users should be warned that they can be worse than

inaccurate, they can produce meaningless results.

2 The tolerance distribution

In bioassays with binary response,2the probability of response is often thought to reflect

an underlying latent distribution known as the tolerance distribution. An individual’s

tolerance is the smallest dose that produces a response. The tolerance distribution describes

the distribution of tolerances in the population.

Consider for example an old fashioned ‘kill-em-and-count-em’ acute toxicicity assay, in

which various doses of a toxin are applied to a test species. The response is death, and

the probability of the response is conditional on the dose. The dose expected to produce a

specified response probability is a quantile of the tolerance distribution.

For individual j let yj denote its response, which is observed, and xjits tolerance, which

is not. Let di be the ith dose. Then the probability of a response is the probability that

the tolerance is no greater than the dose: Pr(yj = 1|di) = Pr(xj ≤ di). The conditional

response distribution is yj|di ∼ BERN(πi). The tolerance distribution is xj ∼ f(μ, σ2),

where f(·) is the PDF of a location-scale distribution.

The expectation of the conditional response distribution is related to the standardized

tolerance distribution by

πi = F

 

di − μ

σ

!

A binomial generalized linear model, g(πi) = α+βdi, connects the response distribution to

2

In the past, the term ‘quantal’ was often used for a binary response.

STATWI0001.02 page 2 of 11
3 NON-PARAMETRIC ESTIMATORS 3

the tolerance distribution through the link function: πi = g

−1

(α + βdi) = F ((di − μ)/σ ).

The estimated mean and variance of the tolerance distribution are then given by the regression parameter estimates, and μ = −α/β, σ = 1/β .

Logit and probit link functions correspond to logistic and normal tolerance distributions,

respectively. Since they are symmetrical, their means and medians are the same:

m = μ = −α/β . The complementary-log-log and log-log link functions correspond to

extreme value distributions, which are asymmetrical. For those distributions, the median

would be m = {log (log(2)) − α}/β .

0 1 2 3 4 5

0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

−log (dose) d

Fraction Positive

●

●

●

●

●

●

●

●

●

●

●

●

2.33

0

10

3

10

4

10

6

9

7

9

10

10

Figure 1: The tolerance distribution.

Figure 1 illustrates this relationship. The data shown at the top of the plot are the number

of mice dying out of a group of mice administered a dose of toxin. Each point is a realization

from a binomial response distribution that is conditional on the dose administered. The

curve, estimated in this case with a logit link, represents the tolerance distribution. The

ED50 is shown by the arrow.

3 Non-parametric estimators

Three of the non-parametric estimators most commonly used to estimate ED50 are described. The Spearman-Karber estimator is an explicit estimator of the mean, not the the

median. The mean corresponds to the ED50 only for symmetrical distributions, of course,

STATWI0001.02 page 3 of 11
3 NON-PARAMETRIC ESTIMATORS 4

and there is rarely enough data to evaluate this assumption in the types of experiments

where it is used — it is an article of faith and hope. As an estimator of the discretized

mean, it is uniform minimum variance unbiased and the maximum likelihood estimator.

The Reed-Muench and Dragstedt-Behrens estimators were intended as estimators of the

median. Miller (1973) points out the surprising result that they are, in fact, asymptotically

equivalent to Spearman-Karber and hence are estimators of the mean.

3.1 Spearman-Karber

The Spearman-Karber method (Spearman 1908; Karber 1931) gives a non-parametric estimate of the mean of a tolerance distribution from its empirical probability mass function

(PMF). The observed data are thought to give an empirical estimate of the cumulative

distribution function (CDF) of the tolerance distribution(Figure 2(a)).3 The empirical

probability mass function (PMF) is derived from the CDF by differencing (Figures 2(b)–

(c)).

The estimator is P(x · f(x)), the usual one for the mean of a discrete distribution, except

that here f(x) is the empirical PMF obtained from the observed data. This estimator depends on the complete distribution, which may not be available in a particular experiment.

If the CDF does not cover the entire support of x, a common practice is to extend it by

assuming the next lower dose would produce zero response and the next higher dose would

produce complete response. Although this is not always a good idea, it is the default in

the function SpearKarb(). While the tolerance distribution is not discrete, it is discretized

by virtue of the interval censoring inherent in experiments of this type.

3.2 Reed-Muench

The Reed-Muench method (Reed and Muench 1938) takes a different approach. It begins

with the belief that there is more information in the experiment than is given by the

observed responses. Instead, we can assume that we know how some of the mice would

have responded had they been given a different dose. A mouse that died at a lower dose

would certainly die at a higher dose, and one that survived a higher dose would certainly

survive a lower dose. This approach can only be considered quasi-statistical, since it treats

observed responses as known constants rather than random variables. There is also the

problem that some of the subjects contribute more information than others.

Here’s how it works. Accumulate the sums in both directions that represent the hypothetical number that at each dose would have died or survived. The actual numbers of observed

responses are shown in columns 1 and 2 of the table in Figure 3(a), and the hypothetical

3Note that the empirical CDF of the tolerance distribution estimated in a dilution experiment of this

type is not the same as an empirical distribution function (EDF).

STATWI0001.02 page 4 of 11
3 NON-PARAMETRIC ESTIMATORS 5

Figure 2: The Spearman-Karber method

(a) Observed CDF of tolerance distribution

0 1 2 3 4 5

0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

−log (dose) d

Fraction Positive

0 3 4 6 7 10

10 7 6 3 2 0

Positive

Negative

(b) Differencing the CDF

0 1 2 3 4 5

0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

−log (dose) d

Fraction Positive

(c) Empirical PMF

0 1 2 3 4 5

0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Fraction Positive

(0.5 × 0.3) + (1.5 × 0.1) +(2.5 × 0.27)+(3.5 × 0.11)+(4.5 × 0.22)+ (5.5 × 0) = 2.36

number of known responses are shown in columns 3 and 4. (When the group sizes are

unequal, as they are in this example, an adjustment is made in the cumulative sums that

effectively averages the group sizes, as shown in columns 5 and 6.)

Next, find the doses that bracket the ED50, i.e. the one where fewer than half of the

hypothetical responses are positive and the one where more than half the hypothetical

responses are positive. The ED50 is found by interpolation between the bracketing doses

to find the dose at which the hypothetical responses would be equal. It is given by the

intersection of the line that connects the hypothetical positive responses and the line that

connects the hypothetical negative responses at the bracketing doses (Figure 3(b)).

STATWI0001.02 page 5 of 11
3 NON-PARAMETRIC ESTIMATORS 6

Figure 3: The Reed-Muench method

(a) Observed responses (Pos, Neg), hypothetical responses (CumPos, CumNeg),

adjusted hypothetical response (AdjCumPos, AdjCumNeg)

Pos Neg CumPos CumNeg AdjCumPos AdjCumNeg

0 10 0 28 0.0 27.6

3 7 3 18 2.9 17.9

4 6 7 11 6.8 11.2

6 3 13 5 13.2 5.4

7 2 20 2 20.7 2.1

10 0 30 0 30.4 0.0

(b) Hypothetical responses

●

●

●

●

●

●

●

●

●

●

●

 0 ●

 5

10

15

20

25

30

0 1 2 3 4 5

2.36

−log (dose) d

Number of Subjects

0 3 7 13 20 30

28 18 11 5 2 0

30 37 40 40 36 28

Positive

Negative

???????

3.3 Dragstedt-Behrens

The Dragstedt-Behrens method (Dragstedt and Lang 1928; Behrens 1929) is very similar

to the Reed-Muench method and is based on the same cumulative sums. For some reason,

most microbiology textbooks present the Dragstedt-Behrens method under the name ReedMuench method.4

Instead of working with the cumulative sums directly, the Dragstedt-Behrens method uses

the fraction of the cumulative sums that are positive at each dose (Figure 4). The ED50is

estimated by interpolation on the line that connects the hypothetical fractions of the bracketing doses.

4For years I struggled to figure out why there were two distinct formulations of the Reed-Muench

method, until Don Kolbe, a microbiologist in the CVB Bacteriology lab, pointed out that one of them was

actually the Dragstedt-Behrens method.

STATWI0001.02 page 6 of 11
3 NON-PARAMETRIC ESTIMATORS 7

●

●

●

●

●

●

−log (dose) d

Positive/(Positive+Negative)

0 3 7 13 20 30

28 21 18 18 22 30

30 37 40 40 36 28

Positive

Pos+Neg

???????

2.33

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 1 2 3 4 5

Figure 4: The Dragstedt-Behrens method

3.4 The skrmdb package

The skrmdb package provides functions for the three nonparameteric estimators described

above. (See the package documentation for more details.) It is, admittedly, a bit comical

to use a computer to perform simple calculations that were devised to avoid the need for

difficult calculations. Nevertheless, these estimators persist, and the package is handy.

> require(skrmdb)

> # use data from skrmdb plots.r

> tmp <- data.frame(

+ y=c(0,3,4,6,7,10),

+ n=c(10,10,10,9,9,10),

+ x=c(0:5)

+ )

> #

> # fit the GLM

> fit <- glm(cbind(y,n-y)~x,binomial,tmp)

> ed <- -fit$coef[1]/fit$coef[2]

> #

> # Spearman-Karber estimate

> skx <- SpearKarb(cbind(y,n)~x,tmp)$ed

> #

> # Reed-Muench estimate

STATWI0001.02 page 7 of 11
REFERENCES 8

> rmx <- ReedMuench(cbind(y,n)~x,tmp)$ed

> #

> # Dragstedt-Behrens estimate

> dbx <- DragBehr(cbind(y,n)~x,tmp)$ed

The model fit gives the estimate ˆμ = −α/ˆ βˆ = 2.6066/1.1163 = 2.3350 . The table shows

all the estimates.

ED50

GLM Fit 2.335

Spearman-Karber 2.356

Reed-Muench 2.360

Dragstedt-Behrens 2.333

References

Behrens, B. (1929), “Zur Auswertung der Digitalisblatter im Froschversuch,” Arkiv fur

Experimentelle Pathologie und Pharmakologie, 140, 297–256.

Dragstedt, C. A. and Lang, V. F. (1928), “Respiratory Stimulants in acute poisoning in

rabbits,” Journal of Pharmacology, 32, 215–222.

Karber, G. (1931), “Beitrag zur kollektiven Behandlung Parmakogischer Reihenversuche,”

Archiv fur Experimentelle Pathologie und Pharmakologie, 162, 480–487.

Miller, R. G. (1973), “Nonparametric estimateors of the mean tolerance in bioassay,”

Biometrika, 60, 535–542.

Reed, L. J. and Muench, H. (1938), “A simple method of estimating fifty percent endpoints,”

American Journal of Hygiene, 27, 493–497.

Spearman, C. (1908), “The method of ”right and wrong cases” (”constant stimuli”) without

Gauss’s formulae,” British Journal of Psychology, 2, 227–242.

STATWI0001.02 page 8 of 11
Appendix

ED50 Formulas

skrmdb package

This appendix provides formulas for calculating ED50 estimates by the methods in the

skrmdb package. See the vignette for the principles underlying the methods. Some relevant

discussion may also be found in Finney (1964).

Notation

yj Number of positives at dilution j. Positives are those with increasing

response.1

j Indexes the data from 1 . . . J, where 1 is the dilution with the lowest

response (smallest number of positives), and J is the dilution with the

greatest response.

nj Total number at dilution j.

pj = yj/nj Fraction positive at dilution j.

xj The ‘dose’ at dilution j. Most often it is the log dilution.

dj = xj+1 − xj Difference of log dilution, for j < J.

2

aj =

Pj

k=1 yk Cumulative sum of the positives from the ‘bottom up’

bj =

Pj

k=J nk − yk Cumulative sum of the negatives from the ‘top down’

zj =

aj

aj + bj

Fraction of cumulative sums

1

If the response is decreasing, either use the complementary response (e.g. affected rather than unaffected)

or reverse the order of the data set.

2With a constant dilution factor, the log dilutions are evenly spaced. In that case, the subscript may be

dropped and d is constant.

Appendix page 1

STATWI0001.02 page 9 of 11
1 Dragstedt-Behrens

ED50 is estimated by x ̃DB, defined as the dilution at which z =

1

2 The dilutions bracketing

x ̃DB are xlow and xhigh. Find them by their corresponding z: zlow = max(z

−

≤ 0.5); zhigh =

min(z

−

> 0.5).

And the ED50 is found by interpolation along the z line segment connecting them.

x ̃DB = xlow + dlow

1

2 − zlow

zhigh − zlow

2 Reed-Muench

ED50 estimated by x ̃RM, defined as the dilution at which a = b. Find the bracketing dilutions

as for x ̃DB, but instead of interpolating on the z line segment, find the intersection of the a

and b line segments.

x ̃RM = xlow + dlow

blow − alow

nlow − ylow + yhigh

It is easy to see that x ̃RM and x ̃DB are estimating the same thing, since a = b ⇔ z =

1

2

.

However the estimates often differ slightly, since x ̃DB is calculated by interpolation along a

single line segment, while x ̃RM is calculated by the intersection of two line segments.3

3 Spearman-Karber

The Spearman-Karber method requires that the pj range from no response to complete

response; i.e. p1 = 0% and pJ = 100%.

The ED50 estimate is

xˆSK =

XJ−1

k=1

(pk+1 − pk) {xk − (xk+1 − xk)/2 }

In estimating the mean of the probability mass function, the first term in the summation is

the estimated mass and the second term assigns it to the midpoint of the dilution interval.

Unlike the quasi-statistical Dragstedt-Behrens and Reed-Muench methods, the variance of

the Spearman-Karber estimator can be calculated. It is:

V ar(ˆxSK) = PJ

k=1 n

(xk+1 − xk)

2

pk (1 − pk)

o

/(nk − 1)

3See Figures 3 and 4 in the vignette

Appendix page 2

STATWI0001.02 page 10 of 11
References

Finney DJ (1964). Statistical Method in Biological Assay, Second Edition, New York: Hafner

Publishing. Chapter 20, particularly sections 20.8 (Reed-Muench) and 20.9 (DragstedtBehrens). (And yes, this is the penultimate 2nd edition.)

Appendix page 3

STATWI0001.02 page 11 of 11