Vignette for ‘cir’ Package

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

Assaf Oron



We demonstrate the cir package with data from the anesthesiology experiment of Benhamou et al.1. This experiment used the Up-and-Down (UD) dose-finding design, and was re-analyzed a few years later by Pace and Stylianou.2 Here we re-analyze them again using state-of-the-art Centered Isotonic Regression (CIR),3 the core estimator found in our package.

The study randomized women 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 test whether they are different. The original study used a “traditional” dose-averaging estimation method, concluding that even though levobupivacaine seemed 19% more potent, the difference between ED\(_{50}\)’s, and hence the difference in potency, was not significant. To derive inference about this difference, they apparently used the standard errors of the dose averages in a \(t\)-test like manner.

Pace and Stylianou 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. 4 They found a larger and statistically significant difference between the two ED\(_{50}\)’s. Per their estimates, levobupivacaine was 37% more potent.

Here we revisit this experiment yet again.

Experimental Trajectories

Since UD datasets are rather compact, 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”. 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.

The figure sizes have been customised so that you can easily put two images side-by-side.

# For brevity, we initially use integers to denote the doses. 
# We make use of R’s shorthand for consecutive sequences, 
# e.g., 1:3 is really 1,2,3
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.

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

ED\(_{50}\) Estimates

To derive IR and 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

From here, ED\(_{50}\) estimation is one step away:

ropiTargetCIR=quickInverse(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE)
#   target      point lower90conf upper90conf
# 1    0.5 0.09383622  0.08090006   0.1060014
levoTargetCIR=quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE)
#   target      point lower90conf upper90conf
# 1    0.5 0.06842105  0.05395897  0.08389953


Let us calculate 83% CIs, to evaluate the evidence for different potencies via the overlapping-intervals method.

quickInverse(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE, conf=0.83)
#   target      point lower83conf upper83conf
# 1    0.5 0.09383622   0.0828185   0.1043032
quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE, conf=0.83)
#   target      point lower83conf upper83conf
# 1    0.5 0.06842105  0.05602481  0.08173084

The lower 83% confidence limit for ropivacaine’s ED\(_{50}\) is almost touching the upper 83% limit for levobupivacaine, but there’s no overlap. One may conclude that there is moderate evidence for a difference in potency at the \(\alpha=0.05\) level. Thus, our conclusion is similar to Pace and Stylianou’s reanalysis, but not as strong. Our point estimates are also very similar to theirs, and so is the levo/ropi potency ratio (1.37, quite a bit higher than Benhamou et al.’s 1.19). This should not be too surprising, because CIR is a direct upgrade of IR which Pace and Stylianou used, whereas Benhamou et al. used dose averaging.

Looking more closely for the source of discrepancy, the difference is not in the ropivacaine estimates (0.092 vs. 0.093, and 0.094, by Benhamou et al., Pace-Stylianou, and us, respectively) - but rather in levobupivacaine (0.077, 0.068, 0.068). Looking even deeper, the relatively high dose-avaeraging estimate for levobupivacaine in the original article is mostly due to the experiment’s meandering (rather than straight) path downward to the region where most of the doses were given. There is reason to suspect that the starting-point effect was not mitigated well enough. The original dose-averaging estimate probably excluded only the first 1-2 doses. If one excludes, say, the first 10 doses, the average gets closer to IR/CIR (0.072), and the difference may become significant.

Speaking of Pace and Stylianou’s point estimates, we can use the same quickInverse() function to reproduce them:

quickInverse(bhamou03ropiRates, target=0.5, estfun=oldPAVA, conf=0.83)
#   target      point lower83conf upper83conf
# 1    0.5 0.09287671  0.08456467   0.1013001
quickInverse(bhamou03levoRates, target=0.5, estfun=oldPAVA, conf=0.83)
#   target      point lower83conf upper83conf
# 1    0.5 0.06846154  0.05851263   0.0796707

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

Our confidence intervals are still different from Pace and Stylianou’s, even when using the same point-estimation method. For intervals they used the bootstrap, whereas we use an analytical approach based on Morris’ theoretical work and the Delta method.6 Their levobupivacaine 83% CI is slightly wider than ours (0.059-0.081 vs. 0.059-0.080), but their ropivacaine interval is substantially narrower (0.087-0.097 vs. 0.085-0.101). The latter is the more typical difference between the two approaches. In fact, their ropivacaine 83% bootstrap CI is unrealistically narrow given the amount of information in the experiment, taking up only a single dose spacing (0.01).

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

Each of these returns a list of three doseResponse objects. Here’s one of the lists:

# $output
#         x         y weight
# 0.07 0.07 0.1250000      3
# 0.08 0.08 0.3888889      8
# 0.09 0.09 0.3928571     13
# 0.1  0.10 0.6721501     10
# 0.11 0.11 0.8553030      4
# 0.12 0.12 1.0000000      1
# $input
#         x         y weight
# 0.07 0.07 0.1250000      3
# 0.08 0.08 0.3888889      8
# 0.09 0.09 0.3928571     13
# 0.1  0.10 0.7727273     10
# 0.11 0.11 0.7000000      4
# 0.12 0.12 1.0000000      1
# $shrinkage
#              x         y weight
# 0.07 0.0700000 0.1250000      3
# 0.08 0.0800000 0.3888889      8
# 0.09 0.0900000 0.3928571     13
# 0.1  0.1028571 0.7519480     14
# 0.12 0.1200000 1.0000000      1

$input is the incoming data, $output the estimates at the same doses, and $shrinkage the doses where the algorithm makes the estimates (then interpolated to generate $output). With oldPAVA() this will always be identical to $output.

Note: To obtain a quick, “short-form” forward estimate – the forward analogue of quickInverse() used above – use quickIsotone(). You will get an output similar to the $output object, together with confidence intervals for each dose. That function’s arguments are similar to quickInverse().

We end with the dose-response plots.

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')
# 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(ropiTargetCIR$lower90conf,ropiTargetCIR$upper90conf), rep(0.5,2), col='purple', lwd=2)

# Adding legend:
legend('bottomright', legend=c("Observed Proportions",'Isotonic Regression',
                              'Centered Isotonic Regression','Estimate +/- 90% CI'),

### Now, second plot for Levobupivacaine
plot(bhamou03levoRates, xlab="Concentration (%)", 
ylab="Proportion Effective", main='Levobupivacaine Arm', ylim=0:1)
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(levoTargetCIR$lower90conf,levoTargetCIR$upper90conf), rep(0.5,2), col='purple', lwd=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. 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).↩︎

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

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