Vignette for the cir Package

Overview

I demonstrate the cir package with data from the anesthesiology experiment of Benhamou et al. (2003)¹, which used the Up-and-Down (UD) dose-finding design. It 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 cir package’s core estimator.

Up-and-down data analysis was the original motivation for developing CIR and cir. If you are not familiar with UD designs, an overview for general scientific audiences was recently published by Oron et al. (2022).⁴ An even more recent overview oriented towards statisticians and data scientists, is available in Oron and Flournoy (2025).⁵

Despite its UD-related history, cir is useful with any binary-outcome data for which a monotone dose-response relationship is expected, regardless of how that data were obtained. Given \(n\) observation pairs consisting of doses \(x_i, i=1, \ldots, n\) and binary responses \(y_i, i=1, \ldots, n\), CIR estimates the dose-response function \(F(x)\), or in short, just \(F\):

\[ F\equiv F(x) \equiv \Pr\left(Y = 1 \mid X = x \right),\] making no assumptions except for strict monotonicity. While the package default (and our description here) assumes an increasing monotone \(F\), the decreasing case is also accommodated, usually by adding the argument \(dec = TRUE\) to function calls.

The package includes methods for both “forward” estimation of \(F\), i.e., the expected response rate conditional on dose, and “inverse” estimation of \(F^{-1}\), a.k.a. dose-finding. Both estimation directions have quick, user-friendly all-in-one functions, as well as more detailed in-depth functions.

If your interest is mainly in UD data analysis rather than more general dose-response data, then you are strongly 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, which while more stable it would also tend to be a couple of steps behind, as upndown is still relatively new and changing a lot.

The Study in our Example

Benhamou et al. (2003) randomized consenting pregnant patients in labor into two arm of \(n=40\) each, to receive as analgesic either ropivacaine or levobupivacaine. The goal was to estimate each agent’s median effective dose (ED\(_{50}\)) for this condition, and to evaluate whether they differ in same-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 such CIs overlap is roughly equivalent to rejecting the Null hypothesis with \(\alpha = 0.05.\) ⁶ They found a larger difference than Benhamou et al.: according to the IR estimates, levobupivacaine was 37% more potent . 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 (in the original, concentrations were described as percents).

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)

Given knowledge of the UD design’s rules, the responses (\(y\)) can be read off directly from the doses (\(x\)). In this specific case, 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 had to give up the 40th and last observation in each arm, because we only know its dose but not the response. 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. By the way, DRtrace() and its plotting method can also accommodate data in which observations were collected in groups or cohorts rather than one-by-one as was done here. See function help for details.

par(mfrow=c(1,2), las=1, mar=c(4,4,4,1)) # image format parameters
doserange = (5:12)/100 # This helps set the same y-axis range and labels

plot(bhamou03ropi, dosevals=doserange, ylab="Concentration (%)", main='Ropivacaine Arm')
legend('topright', legend=c('Effective','Ineffective'), pch=c(19,1), bty='n')
plot(bhamou03levo, dosevals=doserange, ylab="Concentration (%)", main='Levobupivacaine Arm')

The “traditional” estimates in the original articles used some type of average of the doses administered. This is based on the premise that over time, UD designs concentrate these doses roughly symmetrically around the target percentile. Generally speaking, these “traditional” estimates are rather outdated and lack robustness. We recommend CIR as the standard estimator for UD 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) )

knitr::kable(bhamou03levoRates, row.names=FALSE, align='ccr', digits=c(2,4,0) )

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. Indeed, statistically the task is described as estimating \(F^{-1}(0.5)\), or more generally \(F^{-1}(\Gamma)\), with \(100\Gamma\) being the target percentile.

With the Ropivacaine arm data (left), there is 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 clarify that picture via the simplest pooled response-rate averages that makes the overall dose-response curve obey monotonicity.

In any case, CIR estimation of the entire \(F(x)\) curve 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, conf = 0.95)
ropiTargetCIR 
#   target      point lower95conf upper95conf
# 1    0.5 0.09383622  0.07709167   0.1030728
levoTargetCIR = quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE, conf = 0.95)
levoTargetCIR
#   target      point lower95conf upper95conf
# 1    0.5 0.06842105  0.05754698  0.09874464

You can also perform the estimation directly from the raw x, y data; the functions will create a doseResponse object on the fly.

The ratio between our point estimates is 1.37, similar to the Pace-Stylianou IR result. Not too surprising, as the methods are closer to each other than to the original authors’ dose-averaging approach; in particular, both are less affected by starting dose, which played a substantial role in the Levobupivacaine arm.

Comparing the Two Arms

Having examined the methodology behind the CI-overlap approach (which to be prudent, should have been between 84% CIs rather than 83% CIs), we do not recommend it anymore. Besides its somewhat contrived nature, this method is not robust to deviations from distribution assumptions, or even to sample size. Instead, we recommend the more mainstream approach of calculating the (95%) CI for the difference, and observing whether it contains zero.

To obtain the CI of the difference we use the assumption that CI widths are roughly proportional to the square-roots of the variances, as well as independence between arms, while allowing for asymmetry of CI widths to the left and right of the point estimate:

\[ \begin{array}{rcl} LHW_{diff.} & = & \sqrt{LHW_1^2 + RHW_2^2} \\ RHW_{diff.} & = & \sqrt{RHW_1^2 + LHW_2^2}, \end{array} \] where LHW, RHW denote the left- and right-hand widths respectively, and group 2 is the one whose estimates are being subtracted from group 1. Let us put numbers into these formulae:

lhw_1 = ropiTargetCIR$point - ropiTargetCIR$lower95conf
rhw_1 = ropiTargetCIR$upper95conf - ropiTargetCIR$point 
lhw_2 = levoTargetCIR$point - levoTargetCIR$lower95conf
rhw_2 = levoTargetCIR$upper95conf - levoTargetCIR$point 

lhw = sqrt(lhw_1^2 + rhw_2^2)
rhw = sqrt(lhw_2^2 + rhw_1^2)

## Lower 95% bound of difference 
ropiTargetCIR$point - levoTargetCIR$point - lhw
# [1] -0.009224404
## Lower 95% bound of difference 
ropiTargetCIR$point - levoTargetCIR$point + rhw
# [1] 0.03968259

Somewhat surprisingly, we reach a conclusion similar to the original authors: the 95% CI for the difference contains zero, suggesting no evidence for a difference in potency at the \(\alpha=0.05\) level.

Why the difference conclusions from Pace and Stylianou? Our CIs are based on an entirely different approach. Per our repeated examinations over the years, the Pace-Stylianou approach substantially lacks in coverage. That approach relied upon the bootstrap, whereas we use an analytical approach based on Morris’ theoretical work with intervals for datasets with monotone dose-response data,⁸ with some modifications.

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.95)
PSropiEstimate
#   target      point lower95conf upper95conf
# 1    0.5 0.09287671  0.07845822   0.1036113
PSlevoEstimate = quickInverse(bhamou03levoRates, target=0.5, estfun=oldPAVA, conf=0.95)
PSlevoEstimate
#   target      point lower95conf upper95conf
# 1    0.5 0.06846154  0.05999665  0.09607577

Note the argument that chooses oldPAVA() as estimation function (PAVA is the name of the leading algorithm to produce IR estimates; the default is cirPAVA(), i.e., CIR).

Incidentally, plugging the 95% (rather than 83%) CIs calculated by Pace and Stylianou (2007) into the CI-of-difference formula above, we obtain (-0.005, 0.037). Somewhat narrower than the CIR-derived interval, but yielding the same inference - in contrast with the conclusion they actually reached using the 83% CI overlap approach. Looking more closely at their CIs, both \(LHW_1\) and \(RHW_2\) (which in this case are the two half-widths that determine whether there is an overlap) shrank in the 83% CI version to less than half of their 95% CI magnitude. Under a Normal distribution assumption, the 83% CI would be \(\sim 0.7x\) as wide; only rather unusual distributions reach a \(<0.5x\) ratio.

To end the vignette, we visualize more clearly the CIR (and IR) 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, 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(ropiTargetCIR$lower95conf,ropiTargetCIR$upper95conf), 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.08, 0.1), rep(0.49,2))

# Adding legend:
legend('topleft', legend=c("Observed Proportions", 'Isotonic Regression',
                              'Centered Isotonic Regression', paste(c('CIR', 2007), 'Estimate +/- 95% 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,  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(levoTargetCIR$lower95conf,levoTargetCIR$upper95conf), rep(0.5,2), col='purple', lwd=2)
points(I(target-0.01) ~ point, data=PSlevoEstimate, cex=2)
lines(c(0.058, 0.095), rep(0.49,2))

References