`cir`

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

.

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

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.
= 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)
xropi = 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) xlevo
```

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)
= DRtrace(x=xropi[-40]/100, y=(1-diff(xropi))/2)
bhamou03ropi = DRtrace(x=xlevo[-40]/100, y=(1-diff(xlevo))/2) bhamou03levo
```

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
= c(5,12)/100
doserange
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)).

To derive CIR estimates, we take the `DRtrace`

trajectory
objects and generate `doseResponse`

dose-rate-count
summaries.

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

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 |

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

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:

```
= quickInverse(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE)
ropiTargetCIR
ropiTargetCIR # target point lower90conf upper90conf
# 1 0.5 0.09383622 0.07837103 0.1020148
= quickInverse(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE)
levoTargetCIR
levoTargetCIR# 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.

- 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**The correction serves mostly to provide better CI coverage for the target dose estimate.*not*estimate any target except 0.5*(or a value very close to 0.5)*from these data, except for exploratory/illustrative purposes.

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

```
= quickInverse(bhamou03ropiRates, target=0.5, estfun=oldPAVA, conf=0.83)
PSropiEstimate
PSropiEstimate# target point lower83conf upper83conf
# 1 0.5 0.09287671 0.0831183 0.09971887
=quickInverse(bhamou03levoRates, target=0.5, estfun=oldPAVA, conf=0.83)
PSlevoEstimate
PSlevoEstimate# 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.

To construct the full estimated \(F(x)\) curves from IR and CIR, we call
`cir`

’s *“long-form”* forward-estimation
functions:

```
= oldPAVA(bhamou03ropiRates, full=TRUE)
ropiCurveIR = cirPAVA(bhamou03ropiRates, target=0.5, adaptiveShrink=TRUE, full=TRUE)
ropiCurveCIR = oldPAVA(bhamou03levoRates, full=TRUE)
levoCurveIR = cirPAVA(bhamou03levoRates, target=0.5, adaptiveShrink=TRUE, full=TRUE) levoCurveCIR
```

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,With`$shrinkage`

is the actual regression curve found by the algorithm.`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 ** (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))
```

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.

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

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

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

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

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

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