Experimental Examples
In selecting examples, I sought relatively recent open-access publications, preferably from different fields. The publications also had to share clearly both the design and the raw data.
Material Strength Testing
We begin with a study testing the single-tooth bending failure (STBF) load of steel helicopter gears, published by an Italian team in 2017.5 The article summarized a long-running series of STBF experiments and simulations, and can perhaps be seen as a meta-analysis. The experiments were 9 separate median-finding, or as we call them “Classical” UDD runs on gears made of 9 different types of steel. In each such run, each tested gear was subjected to up to 10 million high-frequency cycles at a specific load strength, presumably approximating an entire gear-lifespan worth of loads.
- If the gear survived these cycles without breaking a tooth, it was considered a success and the next gear was subject to a load \(1kN\) higher. (“up”).
- If a tooth broke, it was a failure and the next gear received \(1kN\) lower load (“down”).
Classical UDD was first developed during World War II, to estimate the median effective dropping height of explosives.6 Independently, Nobel-Prize winning audiologist von Bekesy developed a near-identical design to estimate the hearing threshold.7
Gorla et al. visualize the raw data of their 9 experiments (Tables 4-12) in the very same manner that Dixon and Mood visualized an explosive-testing dataset in their 1948 article: as a text graph with “x” and “o” symbols. We re-plot the data in the somewhat more modern fashion encountered in other fields. Here are two of the experiments:
library(upndown)
# From Gorla et al. (2017) Table 6
gorla751x = 39 + c(3:0, 1, 2, 1:3, 2, 3, 2, 3)
# With UD data, one can usually discern the y's (except the last one) from the x's:
gorla751y = c( (1 - diff(gorla751x) ) / 2, 1)
udplot(x=gorla751x, y=gorla751y, main = "Gorla et al. 2017, Material 751",
ytitle = "Load (kN)", las = 1)
legend('bottomright', legend = c(expression('Survived 10'^7*' cycles'),
'Failed'), pch=c(1, 19), bty = 'n')
- This run had \(n=13,\) which for live-subject experiments would be on the low side. Given these gears are industrial products which likely undergo tight process control, such a small \(n\) might suffice.
- That said, the random walk nature can be seen in how the first half of the run is a bit different from the second. After 6-7 observations, one might jump to conclude that the target ( = median failure threshold) is surely below 41 kN, probably around 40 kN. But the last 6 tell a rather different story, of the target tightly locked between 41 and 42 kN.
Which half is correct? We don’t know. Each half provides a ridiculously small amount of information, and even combined they leave a lot of uncertainty (if the uncertainty is properly accounted for).
Here’s the dose-response plot, and below it the numerical target estimate:
drplot(x=gorla751x, y=gorla751y, addest = TRUE, target = 0.5, addcurve = TRUE,
percents = TRUE, main = "Gorla et al. 2017, Material 751",
xtitle = "Load (kN)", ytitle = "Percent Failure", las = 1)
legend('bottomright', legend = c('CIR target estimate and 90% CI',
'CIR load-failure curve estimate'), lty = 1,
col = c('purple', 'blue'), bty = 'n')
udest(gorla751x, gorla751y, target = 0.5)
# target point lower90conf upper90conf
# 1 0.5 41.17241 39.57807 41.7665The 'X' marks denote raw response frequencies, with
symbol area proportional to sample size.
Note that indeed the confidence interval is fairly wide.
It is also asymmetric, mostly because cir
dose-finding functions (starting version 2.3.0) separately estimate the
local slope of \(F(x)\) to the right
and left of target. Note also that upndown default (and
cir default too) is a \(90\%\) CI rather than \(95\%\). You hopefully agree that with 13
dependent binary observations, pretending to know \(95\%\) of anything is a bit overly
ambitious.
CIR stands for Centered Isotonic Regression, an adaptation of that longstanding nonparametric method to the dose-response and dose-finding realms. It implicitly assumes that the dose-response (in this particular case, load-failure) curve is smooth and strictly monotone, and therefore both interpolates between observations, and shrinks the characteristic flat stretches produced by isotonic regression towards single points.
Gorla et al. used the second-oldest UD estimator, suggested by
Brownlee et al. in 1953.8 Like most early UD estimators, it is some
average of the sequence of doses, and the first one to deploy the trick
of adding the \(n+1\)st dose, which is
pre-determined by the \(n\)th dose and
response. In addition, the estimator excludes the first sequence of
doses with identical responses - in UD parlance, the doses before
the first reversal. This yields 40.91 kN.
upndown includes a generic implementation of
reversal-driven estimates, so this number can be replicated by choosing
the right arguments:
# Note the manual addition of "dose n+1"
reversmean(c(gorla751x, 41), gorla751y, rstart = 1, all = TRUE)
# [1] 40.90909Side Note: Observation Bias and Empirical Correction
Estimating the target with CIR is as simple as reading off the dose-response plot, and seeing where the CIR curve crosses \(y=target.\)
There is one catch though: the up-and-down process induces dependence between consecutive doses and responses. Therefore, even though each response is independent conditional upon the dose it received, responses are not marginally independent and they do not behave according to Binomial assumptions - contrary to common misconceptions in the dose-finding field (of which yours truly had also been guilty in the past).
In particular: the dependence induces a bias upon observed response rates, a bias that tends to “flare out” away from target.9 Near target the bias is minimal, and this is why up-and-down and other dose-finding designs still work - but using observed frequencies away from target at face value is not recommended.
The blue curve does include an empirical bias fix, based on that “flaring out” property, which is why it doesn’t cross the middle of the ‘X’ marks denoting observed frequencies. The main motivation for putting that fix in as default, is not to get the curve right - that would be a bit of a fool’s errand, because dose-finding designs are really about estimating a point, not an entire curve - but to get the confidence intervals more realistic. The “flaring-out” bias makes \(F(x)\)’s apparent slope too steep. The bias correction mitigates that steepness and hence widens the CI.
What if we do try to estimate away from target? Here is the (relatively) proper way to do it. Assume we’re interested in knowing a “safe” load, under which only \(5\%\) of gears are expected to fail:
udest(x=gorla751x, y=gorla751y, target = 0.05, balancePt=.5)
# Warning in udest(x = gorla751x, y = gorla751y, target = 0.05, balancePt = 0.5): We strongly advise against estimating targets this far from the design's balance point.
# target point lower90conf upper90conf
# 1 0.05 39.13333 37.58147 39.31827We change the target argument to the desired percentile.
But we must also inform udest() that the experiment was
actually designated to revolve around another percentile, which we call
the balance point.10 The bias would flare out from there, not
from where we want to estimate now. As you might note,
udest() kindly reminds us that this is not recommended.
(If balancePt is left empty, udest()
assumes it is equal to target.)
Want to find out that safe load? Design an experiment that targets a lower percentile (we’ll discuss that soon). Or an experiment that estimates an entire stretch of the load-failure curve. That would not be an UDD though.
On to the second experiment, presented below in one fell swoop.
op <- par(no.readonly = TRUE)
par(mfrow=1:2, mar=c(4,4,4,1))
# From Gorla et al. (2017) Table 9
gorla951x = 35 + c(1:0, 1:4, 3:2, 3:0, 1, 2, 1)
gorla951y = c( (1 - diff(gorla951x) ) / 2, 1)
udplot(x=gorla951x, y=gorla951y, main = "Gorla et al. 2017, Material 951",
ytitle = "Load (kN)", las = 1)
drplot(x=gorla951x, y=gorla951y, addest = TRUE, target = 0.5, addcurve = TRUE,
percents = TRUE, main = "Gorla et al. 2017, Material 951",
xtitle = "Load (kN)", ytitle = "Percent Failure", las = 1)
udest(gorla951x, gorla951y, target = 0.5)
# target point lower90conf upper90conf
# 1 0.5 36.26829 35.28684 38.65569
par(op) # Back to business as usual ;)This run was a bit less well-behaved than 751, and produced a monotonicity violation in observed response frequencies. CIR (and isotonic regression in general) was designed to deal with that, but as a result the confidence interval is even wider than the 751 estimate’s CI, despite having 2 (wow!) additional data points.
Curiously, Gorla et al. mixed and matched two estimators in their study. For this run, they used the original Dixon-Mood estimator (sometimes called “Dixon-Massey”), which also a weighted average of doses but in fact a more sophisticated than the Brownlee et al. one. Using that, they got 36.63 kN with a standard error of 0.62 kN, which translates to a \(90\%\) CI width of 2.0 kN (since Dixon and Mood recommend using Normal tables rather than, e.g., \(t\)). Our CI is \(\sim 1.7x\) wider.
We have a version of the Dixon-Mood in our package as well.
dixonmood(x=gorla951x, y=gorla951y)
# [1] 36.78571The value is a bit different from Gorla et al.’s, because there’s a myriad variations on how that estimator is implemented. Our package uses the formula as provided in the original article. However, we consider that estimator rather obsolete, and provide it only for comparisons like this.
\(ED_{90}\) of Vasopressor during Caesarean Section
Whew! We move from the “Wild West” of engineering UDD applications, to a field where the design is used with substantially more recent methodology. Anesthesiology is probably where one finds nowadays the richest and best-informed variety of UDD application articles, with much thanks to an influential outreach/tutorial article published in 2007 by Pace and Stylianou.11
Pace and Stylianou have successfully nudged anesthesiologists towards using isotonic regression for estimation instead of mid-20th-Century improvisations, and also introduced to that field a UDD that can target percentiles other than the median. Anesthesiologists find much use for estimating a high percentile - the 80th, 90th or the 95th - to identify a dose that would work most of the time, but would not be excessive.
The Up-and-Down Biased Coin Design (BCD) (not to be confused with BCDs in other applications) is part of a 1990s body of work by Durham, Flournoy, and others, the first substantial modern revisitation of UDDs after a generation-long gap.12 BCD can target any percentile \(p.\) For targets above the median -
- After a negative response, move up
- After a positive response, “toss a biased coin” and with probability \((1-p)/p\) move down.
- If the toss “fails”, the next patient/specimen/etc. will receive the same dose.
For a target below the median, flip the coin-tossing to be carried out after a negative response, and invert the toss probability to \(p/(1-p).\)
The anesthesiology study we chose was run by an international team, including a leading popularizer of UDDs in Anesthesiology, Malachy Columb.13 They sought to estimate the \(ED_{90}\) of phenylephrine for treatment of hypotension during Cesarian birth.
Unlike machined helicopter gears, people do not undergo industrial process control; we come in a broad diversity shapes, sizes, ages, and sensitivities. As seen in the previous examples, even with an industrially-controlled sample, \(n\) of 13-15 seems very questionable for dose-finding. With a sample of humans or other animals, such a sample size is not advised. Generally, anesthesiologists have been rather responsible in using more realistic sample sizes for UDD studies; typical samples encountered in that field are \(n=30-40\) for median-finding, and larger samples for extreme percentiles. The George et al. study described here ended up with \(n=45\).
The context was Cesarean delivery, during which many mothers experience hypotension. Phenylephrine is one of the recommended treatments, and the team wanted to estimate the \(ED_{90}.\) Inspired by Pace and Stylianou’s article, they chose BCD, albeit to simplify the “toss” probability they rounded it down from \(1/9\) to \(1/10.\) Since the treatment would be administered only to mothers experiencing hypotension, they initially enrolled 65 mothers, 45 of whom were treated for hypotension. The starting dose was the one most commonly used in practice: \(100\ \mu g.\)
george10x = 80 + 20 * c(1, rep(2, 5), 1, 1, 0, 0, rep(1, 7), 0:2, 2, 2, rep(1, 4), 2, 1, 1, 2, 2,
rep(3, 5), 4, 5, 5, rep(4, 6))
george10y = c(ifelse(diff(george10x) > 0, 0, 1), 1)
udplot(x=george10x, y=george10y, ytitle = dlabel )
The study’s UDD balance point was \(p=10/11,\) even as the experimental target - the percentile to be estimated - was still \(0.9\). This is certainly close enough to not worry about that pesky adaptive-response bias, but we should use the actual balance point to get our numbers right.
On to the dose-response curve and the estimates!
The CI veering off further to the right rather than to the left (which could seem plausible looking at the ‘X’ marks) might raise questions. Here are some insights:
- Remember \(y\) cannot exceed 1. Therefore, there’s a limit on how much \(F(x)\) can move to the left while remaining reasonably close to the data and below 1. Equivalently, the dose of \(100\ \mu g\) which remains outside the CI despite having a 13 of 17 (\(76\%\)) success rate is both “buffered” from the point estimates by the doses of 120 and 140 having even higher observed rates, and itself having relatively many observations increasing the confidence that indeed, the rate at that dose is below target.
- (aside: to avoid misleading the reader, our CI process does not simulate various curves; the limited ability of curves to “move left” is represented via the forward (\(F(x)\)) CIs having to remain below 1)
- Conversely, to the right of the point estimate, there’s more “room” for the curve to move (i.e., the lower forward bounds are much further down from the point estimate), and the shallowness of the estimated CIR curve suggests the ED\(_{90}\) estimate can indeed move quite a bit rightward (and/or become even shallower) while remaining close to the data.
Generally with high/low targets, with our CI method one might expect the “external” confidence bound to be further away. We feel this does represent genuine uncertainty regarding how \(F(x)\) levels off towards the \(y=0\) or \(y=1\) boundary. In many cases, only a much larger sample size can probably narrow that down considerably.