Vignette for the cir Package

Using Up-and-Down data from Benhamou et al. (2003)

Assaf Oron



I demonstrate the cir package with data from the anesthesiology experiment of Benhamou et al. (2003)1. This experiment used the Up-and-Down (UD) dose-finding design, and was re-analyzed a few years later by Pace and Stylianou (2007).2 Here I re-analyze the dataset again, using state-of-the-art Centered Isotonic Regression (CIR),3 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.4 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 I have a dedicated UD package called upndown in “advanced Beta” stage (as of Spring 2023). 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 that package using remotes::install_github(assaforon/upndown).

The Study in our Example

Benhamou et al. (2003) randomized women 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 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.\) 5 They found a larger difference: levobupivacaine was 37% more potent according to the IR estimates. Furthermore, they diffrence 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.

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 adding 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')
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 Plain 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))
x y weight
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))
x y weight
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 near the end, 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)
#   target      point lower90conf upper90conf
# 1    0.5 0.09383622  0.07837103   0.1020148
levoTargetCIR = quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE)
#   target      point lower90conf upper90conf
# 1    0.5 0.06842105   0.0586975   0.0985161

You can also do the estimation single-step from the x, y data; the functions will create a doseResponse object on the fly.


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)
#   target      point lower83conf upper83conf
# 1    0.5 0.09383622  0.07966513   0.1007468
levo83 = quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE, conf=0.83)
#   target      point lower83conf upper83conf
# 1    0.5 0.06842105  0.06006152   0.0911622

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,7 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.091 (>2 spacings above that estimate). Conceptually, the bootstrap CIs appear unrealistically narrow given the amount of information available - only 39 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)
#   target      point lower83conf upper83conf
# 1    0.5 0.09287671   0.0831183  0.09971887
PSlevoEstimate =quickInverse(bhamou03levoRates, target=0.5, estfun=oldPAVA, conf=0.83)
#   target      point lower83conf upper83conf
# 1    0.5 0.06846154  0.06239488  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 argument for 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:

# $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

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 background).

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))



  1. Benhamou D, Ghosh C, Mercier FJ. A Randomized Sequential Allocation Study to Determine the Minimum Effective Analgesic Concentration of Levobupivacaine and Ropivacaine in Patients Receiving Epidural Analgesia for Labor. Anesthesiology. 2003;99(6):1383-1386.↩︎

  2. Pace NL, Stylianou MP. Advances in and Limitations of Up-and-down Methodology: A Précis of Clinical Use, Study Design, and Dose Estimation in Anesthesia Research. Anesthesiology. 2007;107(1):144-152.↩︎

  3. Oron AP, Flournoy N. Centered Isotonic Regression: Point and Interval Estimation for Dose-Response Studies. Stat Biopharm Res. 2017;9(3):258-267.↩︎

  4. Oron AP, Souter MJ, Flournoy N. Understanding Research Methods: Up-and-down Designs for Dose-finding. Anesthesiology 2022; 137:137–50. See also the online supplement.↩︎

  5. Payton ME, Greenstone MH, Schenker N. Overlapping confidence intervals or standard error intervals: What do they mean in terms of statistical significance? J Insect Sci. 2003;3(1).↩︎

  6. Flournoy N, Oron AP. Bias induced by adaptive dose-finding designs. J Appl Stat. 2020;47(13-15):2431-2442.↩︎

  7. Morris MD. Small-Sample Confidence Limits for Parameters under Inequality Constraints with Application to Quantal Bioassay. Biometrics. 1988;44:1083-1092.↩︎