Overview
I demonstrate the cir package with data from the
anesthesiology experiment of Benhamou et al. (2003). This experiment used
the Up-and-Down (UD) dose-finding design, and was re-analyzed a few
years later by Pace and Stylianou (2007). Here I re-analyze the
dataset yet again, using state-of-the-art Centered Isotonic Regression
(CIR),
the core estimator found in our package.
Up-and-down data analysis was the original motivation for developing
CIR and cir. For a general overview of UD designs, see Oron
et al. 2022. However, the package is compatible with any
binary-outcome data for which a monotone dose-response relationship is
expected - regardless of design. The cir package includes
methods for both “forward” estimation, i.e., the
expected response rate conditional on dose, and
“inverse” estimation, a.k.a. dose-finding. Both
estimation directions have quick, more user-friendly all-in-one
functions, and more detailed in-depth functions.
If your interest is mainly in UD data analysis, then you are
encouraged to use the dedicated UD package upndown. It
includes wrappers for even simpler, single-command plotting and
estimation of UD datasets, using cir package functions with
the defaults already dialed onto the UD context. You can download the
latest version of that package using
remotes::install_github(assaforon/upndown). CRAN also has a
copy but it would always be a couple of steps back, as
upndown is still new and changing a lot.
The Study in our Example
Benhamou et al. (2003) randomized pregnant volunteers in labor into
two groups, receiving as analgesic either ropivacaine or
levobupivacaine, to estimate each agent’s median effective dose (ED\(_{50}\)) for this condition, and to test
whether there is significant difference in equal-concentration efficacy
by comparing the two estimates. The original study used a “traditional”
UD dose-averaging estimation method, concluding that even though
levobupivacaine seemed 19% more potent, the difference between ED\(_{50}\)s was not significant. To derive
inference about this difference, they apparently used the standard
errors of dose averages in a \(t\)-test
like manner.
Pace and Stylianou (2007) re-analyzed the data using isotonic
regression (IR) and \(83\%\) bootstrap
confidence intervals (CIs); under certain assumptions, examining whether
these CIs overlap is equivalent to rejecting the Null hypothesis with
\(\alpha = 0.05.\) They found a larger
difference: levobupivacaine was 37% more potent according to the IR
estimates. Furthermore, the difference was deemed statistically
significant because the \(83\%\) CIs
did not overlap.
Here we revisit this experiment yet again. To reduce some confusion,
we drop the percent sign from description of doses used (they were
given as concentrations in percent in the original).
Experimental Trajectories
We do not have the original data table, nor do the original authors
have access to it anymore; but we can read the sequence of administered
doses off of Benhamou et al.’s Figure 1.
# For brevity, we initially use integers to denote the doses.
xropi = c(11:9,10:8,9,10,9,10:7,8:11,10:12,11:7,8,7:10,9,8,9,8:10,9,10,9,10)
xlevo = c(11,10,11,10,11:9,10:7,8,7,8:5,6:8,7,8:6,7,6,7,6,7:5,6,7,6:12)
The study design being “classical” or median-finding UD, the
responses (\(y\)) can be read off
directly from the doses (\(x\)): a
positive increment indicates no effectiveness (\(y=0\)), and vice versa. Symbolically, \[y = \left( 1-diff(x) \right) / 2.\] We use
this and the DRtrace() constructor utility, to create
objects that store doses (converted to their physical units) and
responses; as we call them here, the experimental “trace” or
trajectory.
library(cir)
bhamou03ropi = DRtrace(x=xropi[-40]/100, y=(1-diff(xropi))/2)
bhamou03levo = DRtrace(x=xlevo[-40]/100, y=(1-diff(xlevo))/2)
In the construction above, we gave up the 40th and last observation
in each arm, because we only know its dose but not the response. Since
UD datasets (and more generally, small-sample dose-response datasets)
are rather compact, as shown above entering the data takes no more than
2-3 lines of code. But you can also use .csv file input if
preferred, with two columns headed x and
y.
DRtrace objects have a native plotting method:
par(mfrow=c(1,2), las=1, mar=c(4,4,4,1)) # image format parameters
doserange = c(5,12)/100
plot(bhamou03ropi, ylim=doserange, ylab="Concentration (%)", main='Ropivacaine Arm')
legend('bottomright',legend=c('Effective','Ineffective'),pch=c(19,1),bty='n')
plot(bhamou03levo, ylim=doserange, ylab="Concentration (%)", main='Levobupivacaine Arm')
The “traditional” estimates in the original articles, just used some
type of average of the doses administered. This is based on the premise
that over time, up-and-down designs concentrate these doses roughly
symmetrically around the target percentile.
Generally speaking, the “traditional” estimates are rather outdated
and lack robustness. We recommend CIR as the standard estimator for
up-and-down experiments (see e.g., Oron
et al. 2022 supplement (pdf link)).
Dose-Response Plot and ED\(_{50}\)
Estimates
To derive CIR estimates, we take the DRtrace trajectory
objects and generate doseResponse dose-rate-count
summaries.
bhamou03ropiRates = doseResponse(bhamou03ropi)
bhamou03levoRates = doseResponse(bhamou03levo)
knitr::kable(bhamou03ropiRates, row.names=FALSE,align='ccr',digits=c(2,4,0))
| 0.07 |
0.0000 |
3 |
| 0.08 |
0.3750 |
8 |
| 0.09 |
0.3846 |
13 |
| 0.10 |
0.8000 |
10 |
| 0.11 |
0.7500 |
4 |
| 0.12 |
1.0000 |
1 |
knitr::kable(bhamou03levoRates, row.names=FALSE,align='ccr',digits=c(2,4,0))
| 0.05 |
0.0000 |
2 |
| 0.06 |
0.2500 |
8 |
| 0.07 |
0.5455 |
11 |
| 0.08 |
0.8333 |
6 |
| 0.09 |
0.3333 |
3 |
| 0.10 |
0.4000 |
5 |
| 0.11 |
0.7500 |
4 |
Let us visualize the response frequencies on the
dose-response plane, which is “where CIR estimation action
happens.” (and likewise for any regression-based estimation
of dose-response data)
Analogously to the plots above, the plots below are native
cir methods for doseResponse objects. Symbol
area is proportional to the number of observations at each dose. To
change to fixed-size, use varsize=FALSE.
par(mfrow=c(1,2), las=1, mar=c(4,4,4,1)) # image format parameters
plot(bhamou03ropiRates, xlab="Concentration (%)",
ylab="Proportion Effective", main='Ropivacaine Arm', ylim=0:1)
# Showing the target response rate
abline( h=0.5, col='purple', lty=2)
plot(bhamou03levoRates, xlab="Concentration (%)",
ylab="Proportion Effective", main='Levobupivacaine Arm', ylim=0:1)
abline( h=0.5, col='purple', lty=2)
Target-dose estimation via regression is an “Inverse
Problem”: we draw a line at the desired \(y\) value (in this case, \(50\%\) or 0.5, marked in purple), and look
for the best-fitting \(x\) value. With
the Ropivacaine arm data (left), there’s a relatively clear transition
between low-response and high-response regions, so the target seems to
lie between concentrations \(0.09\) and
\(0.10\). We cannot be very confident
of that, given the limited amount of information - but at least there’s
a clear candidate.
With Levobupivacaine data (right) the picture is far more murky.
Concentrations \(0.07\) and \(0.08\) go above \(50\%\) response, but \(0.09\) and \(0.10\) go back below that level! As we
visualize further below, CIR helps clear that murky picture via the
simplest pooled response-rate averages that makes the overall
dose-response curve obey monotonicity.
But before going behind the scenes: from the
doseResponse object, CIR estimation of \(F(x)\) is only one step away.
quickIsotone(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE)
# x y lower90conf upper90conf
# 0.07 0.07 0.1250000 0.01388341 0.4566917
# 0.08 0.08 0.3888889 0.17029265 0.6144262
# 0.09 0.09 0.3928571 0.20777116 0.6148574
# 0.1 0.10 0.6721501 0.46544685 0.8286039
# 0.11 0.11 0.8553030 0.55419939 0.9356434
# 0.12 0.12 1.0000000 0.57538267 1.0000000
quickIsotone(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE)
# x y lower90conf upper90conf
# 0.05 0.05 0.1666667 0.01687173 0.4598125
# 0.06 0.06 0.2777778 0.10187230 0.5556723
# 0.07 0.07 0.5416667 0.31190997 0.7406421
# 0.08 0.08 0.5542328 0.34408940 0.7480747
# 0.09 0.09 0.5705255 0.37555944 0.7607139
# 0.1 0.10 0.6352627 0.39780747 0.8410408
# 0.11 0.11 0.7000000 0.42005550 0.9213677
Estimating the ED\(_{50}\), a.k.a.
the target dose, is done very similarly:
ropiTargetCIR = quickInverse(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE)
ropiTargetCIR
# target point lower90conf upper90conf
# 1 0.5 0.09383622 0.07837103 0.101341
levoTargetCIR = quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE)
levoTargetCIR
# target point lower90conf upper90conf
# 1 0.5 0.06842105 0.05904391 0.09233012
You can also do the estimation single-step from the x, y
data; the functions will create a doseResponse object on
the fly.
Notes
- The target-dose estimate appears under the word
point:
0.0938 for ropivacaine and 0.0684 for levobupivacaine.
- We specify the target response-rate as a fraction, via the
target argument.
- Default CIs are 90% rather than 95%, due to the large variability
and limited confidence in these typically small-sample experiments. To
change it use the
conf argument.
- The
adaptiveShrink option performs an empirical
correction of the bias induced by the adaptive design, a bias discovered
and described recently by Flournoy and Oron. This type of bias is
induced not just by UD designs, but also by model-based adaptive designs
such as the Continual Reassessment Method. It is minimal near the target
dose, but flares out to substantial magnitudes in both directions.
Because of this bias, we strongly recommend to not
estimate any target except 0.5 (or a value very close to 0.5)
from these data, except for exploratory/illustrative purposes.
The correction serves mostly to provide better CI coverage for the
target dose estimate.
Let us calculate 83% CIs, to evaluate the evidence for different
potencies via the overlapping-intervals method.
ropi83 = quickInverse(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE, conf=0.83)
ropi83
# target point lower83conf upper83conf
# 1 0.5 0.09383622 0.07966513 0.1000783
levo83 = quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE, conf=0.83)
levo83
# target point lower83conf upper83conf
# 1 0.5 0.06842105 0.06038051 0.0867241
Paradoxically, despite using a point-estimation method very close to
Pace and Stylianou’s (CIR is a minor upgrade of IR), we reach a
conclusion similar to the original authors: the 83% confidence
intervals do overlap, suggesting no evidence for a difference in
potency at the \(\alpha=0.05\)
level.
Indeed, our point estimates are very similar to the Pace-Stylianou
reanalysis, and so is our estimated levo/ropi potency ratio (1.37, quite
a bit higher than Benhamou et al.’s 1.19). Where we part ways
with Pace and Stylianou - dramatically - is indeed the confidence
intervals for which our method is quite different. Pace and
Stylianou used an adaptation of the bootstrap, whereas we use an
analytical approach based on Morris’ theoretical work with intervals for
datasets with monotone dose-response data, the Delta method, and
additional modifications to adapt to typical dose-finding data
limitations.
Looking in more detail at the two confidence limits with the biggest
difference between our method and the boostrap: their \(83\%\) LCL for Ropivacaine is \(0.087\) - only 0.6 of a dose-spacing step
below their point estimate - whereas ours is 0.08, nearly 1.5 spacings
below our point estimate. Their Levobupivacaine \(83\%\) UCL is \(0.081\) (1.3 spacings above point-estimate)
vs. 0.087 (>2 spacings above that estimate). Conceptually, the
bootstrap CIs appear unrealistically narrow given the amount of
information available - only 39 dependent binary observations spread
over several doses.
By the way, we can use the same quickInverse() function
to reproduce the Pace-Stylianou point estimates (but not the
CIs, since cir does not implement the bootstrap
approach):
PSropiEstimate = quickInverse(bhamou03ropiRates, target=0.5, estfun=oldPAVA, conf=0.83)
PSropiEstimate
# target point lower83conf upper83conf
# 1 0.5 0.09287671 0.0831183 0.09791556
PSlevoEstimate = quickInverse(bhamou03levoRates, target=0.5, estfun=oldPAVA, conf=0.83)
PSlevoEstimate
# target point lower83conf upper83conf
# 1 0.5 0.06846154 0.06240652 0.08618716
Note the argument specifying the oldPAVA() estimation
function (PAVA is the name of the leading algorithm to produce IR
estimates; the package default is cirPAVA(), i.e.,
CIR).
Let us visualize more clearly the CIR estimates, and the rationale
supporting wider intervals from this dataset.
Visualizing CIR on the Dose-Response Plane
To construct the full estimated \(F(x)\) curves from IR and CIR, we call
cir’s “long-form” forward-estimation
functions:
ropiCurveIR = oldPAVA(bhamou03ropiRates, full=TRUE)
ropiCurveCIR = cirPAVA(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE, full=TRUE)
levoCurveIR = oldPAVA(bhamou03levoRates, full=TRUE)
levoCurveCIR = cirPAVA(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE, full=TRUE)
With full=TRUE, each of these returns a list of three
doseResponse objects. Here’s one example:
levoCurveCIR
# $output
# x y weight
# 0.05 0.05 0.1666667 2
# 0.06 0.06 0.2777778 8
# 0.07 0.07 0.5416667 11
# 0.08 0.08 0.5542328 6
# 0.09 0.09 0.5705255 3
# 0.1 0.10 0.6352627 5
# 0.11 0.11 0.7000000 4
#
# $input
# x y weight
# 0.05 0.05 0.1666667 2
# 0.06 0.06 0.2777778 8
# 0.07 0.07 0.5416667 11
# 0.08 0.08 0.7857143 6
# 0.09 0.09 0.3750000 3
# 0.1 0.10 0.4166667 5
# 0.11 0.11 0.7000000 4
#
# $shrinkage
# x y weight
# 0.05 0.05000000 0.1666667 2
# 0.06 0.06000000 0.2777778 8
# 0.07 0.07000000 0.5416667 11
# 0.08 0.08928571 0.5659014 14
# 0.11 0.11000000 0.7000000 4
$output is the estimates at the original doses
(which is what quickIsotone() returns sans the
CIs),$input is the incoming data, and$shrinkage is the output at the doses where the
algorithm makes the pooled estimates (then interpolated to generate
$output). IOW, $shrinkage is the
actual regression curve found by the algorithm. With
oldPAVA(), this component will always be identical to
$output.
We end with the dose-response plots, the regression curves and the
target estimates, with the two arms aligned and arranged vertically to
help visualize the CI overlap.
As can be seen below, in cir this takes quite a
bit of coding. In the new upndown package dedicated to
up-and-down designs (currently in beta version on GitHub
assaforon/upndown, and soon also on CRAN), most of
what you see in the plot can be achieved in a single command with the
drplot() utility (which does of course rely on
cir functions in the backend).
par(mfrow=c(2,1), las=1, mar=c(4,4,4,1)) # image format parameters
plot(bhamou03ropiRates, xlab="Concentration (%)",
ylab="Proportion Effective", main='Ropivacaine Arm',
ylim=0:1, xlim=c(0.05, 0.12), dosevals = (5:12)/100)
# Adding IR and CIR lines with the same colors/lines as article’s Fig. 4
lines(y~x, data=ropiCurveIR$output, lty=2)
lines(y~x, data=ropiCurveCIR$shrinkage, col='blue',lwd=2)
# Showing the CIR estimate, and confidence interval as a horizontal line
points(target ~ point, data=ropiTargetCIR, col='purple', pch=19, cex=2)
lines(c(ropi83$lower83conf,ropi83$upper83conf), rep(0.5,2), col='purple', lwd=2)
# The estimate appearing in Pace and Stylianou 2007, nudged 0.01 units down:
points(I(target-0.01) ~ point, data=PSropiEstimate, cex=2)
lines(c(0.087, 0.097), rep(0.49,2))
# Adding legend:
legend('topleft', legend=c("Observed Proportions", 'Isotonic Regression',
'Centered Isotonic Regression', paste(c('CIR', 2007), 'Estimate +/- 83% CI')),
bty='n',pch=c(4,NA,NA,16,1), col=c('black','black','blue','purple', 'black'), lty=c(0,2,1,1,1))
### Now, second plot for Levobupivacaine
plot(bhamou03levoRates, xlab="Concentration (%)",
ylab="Proportion Effective", main='Levobupivacaine Arm',
ylim=0:1, xlim=c(0.05, 0.12), dosevals = (5:12)/100)
lines(y~x, data=levoCurveIR$output, lty=2)
lines(y~x, data=levoCurveCIR$shrinkage, col='blue',lwd=2)
points(target ~ point, data=levoTargetCIR, col='purple', pch=19, cex=2)
lines(c(levo83$lower83conf,levo83$upper83conf), rep(0.5,2), col='purple', lwd=2)
points(I(target-0.01) ~ point, data=PSlevoEstimate, cex=2)
lines(c(0.059, 0.081), rep(0.49,2))
Notes
The 2007 estimates were shifted slightly downward in the plot, to
make them visible.
Referring back to the “murky” picture around target in the
Levobupivacaine dose-response plots, we see that both IR (black dashes)
and CIR (solid blue) pool the observations from \(x = 0.08, 0,09, 0.10\) into a single
weighted average \(y\) value, which is
only slightly higher than the estimate at \(x
= 0.07.\) The main difference is that IR creates a flat “stretch”
with that pooled \(y\) value, whereas
CIR also pools along the \(x\) axis
into a single point on the dose-response plane. So both these
nonparametric algorithms “convert the murkiness” into
concluding (with very limited confidence, given the overall small \(n\)) that \(F(x)\) is very shallow for a long stretch,
right above the \(y = 0.5\) target
response rate.
We see how the CIR intervals stretch to accommodate this
shallowness. Indeed, our confidence-interval method for inverse
estimation is driven directly by local slope around target,
which we believe is the more realistic approach.
Take in particular the Levobupivacaine plot at \(x = 0.10\). The CIR estimate of \(F(x)\) there is 0.635. The IR estimate
which is what Pace and Stylianou used to inform their bootstrap, is
0.571 - even closer to the \(50\%\)
study target-rate, and statistically indistinguishable from it at these
sample sizes. Yet, \(0.10\) lies not
only outside the 2007 article bootstrap’s \(83\%\) CI, but also outside their \(95\%\) CI (which is \(0.058, 0.095\)). This suggests that the
bootstrap intervals are substantially too short.
Prior to cir package version 2.3.0, the estimated
slope informing the CI was the same to the right and left of target, and
hence CIs were usually near-symmetric (and in this particular
example, much shorter). We now estimate different “left”
and “right” slopes. To revert to the single-slope version, use
the argument slopeRefinement = FALSE in your
quickInverse() call.
In this demonstration the CIR curves differ from the ordinary IR
curves not only in the CIR algorithm that “shrinks” horizontal intervals
to single points, but also in the bias correction found by Flournoy and
Oron (2020) and mentioned earlier. In principle this correction is
compatible with both methods, but here we applied it only to CIR,
because it did not exist at the time of Pace and Stylianou’s
article.