Vignette for the cir Package

Overview

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 (ED50) 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 ED50s 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 α = 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 = (1 − diff(x))/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 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')
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 Plain and ED50 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 ED50, 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.1018371
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.6 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.1005961
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 α = 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.087 (>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)
PSropiEstimate
#   target      point lower83conf upper83conf
# 1    0.5 0.09287671   0.0831183  0.09971887
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 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:

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

Notes

  • The 2007 CIs 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 plain. 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 by local slope around target, which we believe is the more correct 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 bootstrap’s 83% CI, but also outside the 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.

References

  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.↩︎