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.

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

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

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

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

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

- The estimate appears under the word
`point`

: 0.0938% for ropivacaine and 0.0684% for levobupivacaine. - We specify the target as a fraction, via the
`target`

argument. - Default CIs are 90% rather than 95%, due to the large variabily in
these typical small-sample experiments. To change it use the
`conf`

argument. - The
`adaptiveShrink`

option performs an empirical correction the bias induced by the adaptive design, a bias discovered and described recently by Flournoy and Oron.^{5}It is induced not just by UD designs, but also by model-based designs such as the Continual Reassessment Method.**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 from these data, except for exploratory/illustrative purposes. - In an experiment with non-median targets such as the ED\(_{90}\), we recommend adding the option
`adaptiveCurve=TRUE`

to the command above. This broadens the confidence intervals a bit to ensure proper coverage, which is more challenging the closer the target gets to the edge of the dose-response curve.

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

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

Each of these returns a list of three `doseResponse`

objects. Here’s one of the lists:

```
ropiCurveCIR# $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'),
bty='n',pch=c(4,rep(NA,2),16),col=c('black','black','blue','purple'),lty=c(0,2,1,1))
### 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)
```

- The ‘x’ symbols showing the observed response rates come in varying
sizes: the area they occupy is proportional to the number of
observations. This is the
`plot.doseResponse`

default; to change it use`varsize=FALSE`

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

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