reversmean: Reversal-anchored averaging estimators for Up-and-Down

View source: R/averaging.r

reversmeanR Documentation

Reversal-anchored averaging estimators for Up-and-Down

Description

Dose-averaging target estimation for Up-and-Down experiments, historically the most popular approach, but not recommended as primary nowadays. Provided for completeness.

Usage

reversmean(
  x,
  y,
  rstart = 3,
  all = TRUE,
  before = FALSE,
  conf = 0.9,
  maxExclude = NULL,
  full = FALSE,
  weth66revs = TRUE,
  evenrevs = !all,
  ...
)

reversals(y, x = NULL, directional = TRUE, evenrevs = TRUE)

Arguments

x

numeric vector: sequence of administered doses, treatments, stimuli, etc.

y

numeric vector: sequence of observed responses. Must be same length as x or shorter by 1, and must be coded TRUE/FALSE or 0/1. dynamean() only uses y for bootstrap confidence intervals.

rstart

the reversal point from which the averaging begins. Default 3, considered a good compromise between performance and robustness. See Details.

all

logical: from the cutoff point onwards, should all values of x be used (TRUE, default), or only reversal points as in the Wetherill et al. approach? If set to FALSE, then the before flag also defaults to FALSE regardless of user choice.

before

logical: whether to start the averaging one observation earlier than the cutoff point. Default FALSE.

conf

the CI's confidence level, as a fraction in (0,1). To skip CI calculation set conf = NULL.

maxExclude

a fraction in 0,1 indicating the maximum initial fraction of the vector x to exclude from averaging, in case the algorithm-identified transition point occurs late in the experiment. Default 1/2.

full

logical: should more detailed information be returned, or only the estimate? (default FALSE)

weth66revs

(in reversmean()) logical: identical to directional. The argument name used in reversmean() stems from the article that switched to this definition of reversals when introducing reversal-averaging.

evenrevs

logical: should only an even number of reversals be used, meaning that if the total number is odd the last one is discarded? Default TRUE per common practice. However, when setting all=TRUE it makes sense to also set evenrevs=FALSE.

...

Additional parameters, mostly ones passed on to dfboot if conf is not NULL. See that function's help for details.

directional

(in reversals()) logical: should reversals be defined as change in direction (i.e., x; TRUE which is default), or in response (y)?

Details

Up-and-Down designs (UDDs) allocate doses in a random walk centered nearly symmetrically around a balance point. Therefore, a modified average of allocated doses could be a plausible estimate of the balance point's location.

During UDDs' first generation, a variety of dose-averaging estimators was developed, with the one proposed by Wetherill et al. (1966) eventually becoming the most popular. This estimator uses only doses observed at reversal points: points with a negative response following a positive one, or vice versa. More recent research (Kershaw 1985, 1987; Oron et al. 2022, supplement) strongly indicates that in fact it is better to use all doses beginning from some cutoff point, rather than skip most of them and choose only reversals.

The reversals() utility identifies reversal points, whereas reversmean() produces a dose-averaging estimate whose cutoff point (which should perhaps be called the 'cut-on' point) is determined by a reversal. User can choose whether to use all doses from that cut-point onwards, or only the reversals as in the older approaches. A few additional options make the estimate even more flexible.

Starting version 0.2.0, a bootstrap confidence interval (CI) is also provided. See dfboot, dfsim for additional parameters to pass to the bootstrap routine via ..., beyond the confidence level conf. For the "Classical" median-finding UDD, use ⁠design = krow, desArgs = list(k=1)⁠. To skip CI estimation, set conf = NULL. reversmean() is compatible mostly with median-targeting UDDs such as the "Classical" (traditional) design of Dixon and Mood. For general UDD target estimation, particularly off-median targeting designs, we recommend using centered isotonic regression, available via udest, an up-and-down adapted wrapper to cir::quickInverse(). See Oron et al. 2022 (both article and supplement) for further information, as well as the cir package vignette.

Value

For reversals(), the indices of reversal points. For reversmean(), if full=FALSE returns the point estimate and otherwise returns a data frame with the estimate, as well as the index of the cutoff point used to start the averaging.

Author(s)

Assaf P. Oron <assaf.oron.at.gmail.com>

References

  • Kershaw CD: A comparison of the estimators of the ED50 in up-and-down experiments. J Stat Comput Simul 1987; 27:175–84.

  • Oron AP, Souter MJ, Flournoy N. Understanding Research Methods: Up-and-down Designs for Dose-finding. Anesthesiology 2022; 137:137–50. See in particular the open-access Supplement.

  • Wetherill GB, Chen H, Vasudeva RB: Sequential estimation of quantal response curves: A new method of estimation. Biometrika 1966; 53:439–54

See Also

  • udest, the recommended estimation method for up-and-down targets.

  • dynamean, an unpublished but arguably better approach to dose-averaging (this is not the recommended method though; that would be udest referenced above).

Examples

#'  **An up-and-down experiment that has generated some controversy**
#'  
#' Van Elstraete, AC et al. The Median Effective Dose of Preemptive Gabapentin 
#'      on Postoperative Morphine Consumption After Posterior Lumbar Spinal Fusion. 
#'      *Anesthesia & Analgesia* 2008, 106: 305-308.


# It was a classical median-finding up-and-down study.

doses = c(4:7, 6:13, 12:19, 18:21, 20, 19:23, 22, 21:23, 22:19, 20:23, 
          22:24, 23, 22, 23, 22:25, 24:22, rep(23:24,2), 23, 22)
# With U&D, responses (except the last one) can be read off the doses:
responses = c( (1 - sign(diff(doses)))/2, 0 )


### Let us plot the dose-allocation time series.

# Saving current settings as now required by the CRAN powers-that-be :0
op <- par(no.readonly = TRUE)

par(mar=c(4,4,4,1), mgp=c(2.5,0.8,0), cex.axis = 0.7, las = 1)
udplot(doses, responses, main='Van Elstraete et al. 2008 Study', 
       xtitle = "Patient Number", ytitle = 'Gabapentin (mg/kg)') 


#' Overlay the ED50 reported in the article (21.7 mg/kg):
abline(h = 21.7)

#' The authors cite a little-known 1991 article by Dixon as the method source.
#' However, in their author rejoinder they claim to have used the Dixon-Mood (1948) estimate.


# Our package does include the Dixon-Mood point estimate.
#  (w/o the CIs, because we do not endorse this estimation approach)
# Does it reproduce the article estimate?
dixonmood(x = doses, y = responses)

# Not at all! Let us overlay this one in red
abline(h = dixonmood(x = doses, y = responses), col=2)

# We have found that many articles claiming to use Dixon-Mood (or Dixon-Massey) actually
# Do something else. For example, in this article they report that 
#   "it is necessary to reject sequences with three to six identical results".
# Nothing like this appears in the original Dixon-Mood article, where the estimation method
#   involves identifying the less-common response (either 0 or 1), and using only x values
#   associated with these responses; obviating the need to exclude specific sequences.
#
# More generally, these historical estimates have long passed their expiry dates. 
#   Their foundation is not nearly as solid as, e.g., linear regression, 
#      and it's time to stop using them.

# That said, our package does offer two more types of dose-averaging estimates.
# Both are able to take advantage of the "n+1" dose-allocation, which is determined by
#    the last dose and response:
n = length(doses)
dosePlus1 = doses[n] + ifelse(responses[n]==0, 1, -1)
reversmean(c(doses, dosePlus1), responses, conf = NULL)
# Interestingly, in this particular case the answer is very similar to the Dixon-Mood estimate.

# The `reversmean()` default averages all doses from the 3rd reversal point onwards.
# By the way, at what point did the third reversal happen? 
#     It'll be the 3rd number in this vector:
reversals(x = doses, y = responses)

# Far more commonly in literature, particularly in sensory studies, 
#   one encounters the 1960s-era approach (led by Wetherill) of taking *only doses  
#   at reversal points, usually starting from the first one. `reversmean()` can do that too:
wetherill = reversmean(c(doses, dosePlus1), responses, all = FALSE, rstart = 1, conf = NULL)
wetherill
# This one gives an even lower result than the previous ones.
abline(h = wetherill, col = 3)

# There's another approach to dose-averaging, although it is not in use anywhere that we know of.
# It does not require the y values at all. The underlying assumption is that the dose 
#   sequence has done enough meandering around the true balance point, to provide information
#   about where (approximately) the starting-dose effect is neutralized.
# This function now also provides bootstrap CIs, so we need to give it the y values. 
# The default forces the final 2/3 of observations to be included; here in view of the long run-in
#     we are relaxing this
dynamean(c(doses, dosePlus1), responses, maxExclude = 0.5, conf = NULL)
# Again a bit curiously, this relatively recent approach gives a result similar to what
#   the authors reported (but not similar to the original Dixon-Mood).
# This is not too surprising, since here `dynamean()` excludes the first one-third of doses,
#   which is approximately what happened if indeed the authors excluded all those long dose-increase
#   sequences at the start.

# All this shows how dicey dose-averaging, at face value a simple and effective method, can become.
# The sample size here is rather large for up-and-down studies, and yet because of the unlucky
#    choice of starting point (which in many studies, due to safety concerns cannot be evaded)
#    there is really no good option of which observations to exclude.

# This is one reason why we strongly recommend using Centered Isotonic Regression as default. 
# Figure soon to follow.
# But first, have you noted how we keep specifying "conf = NULL"? 
# This is because now at default, these averages calculate
# A bootstrap confidence-interval. 
# These intervals are generally deficient but are the best anyone can do at present.

# If you want to use a confidence interval, you must provide the experiment's target, 
#    or more precisely its balance point, as well as the parameters to use in 
# the bootstrap simulation (which should be the ones generating the original experiment).

# Like this (not run, to avoid violating CRAN's very narrow limits on example runtime):
# dynamean(c(doses, dosePlus1), responses, maxExclude = 0.5, target = 0.5, 
#                  design = krow, desArgs = list(k=1) )


defest = udest(doses, responses, target = 0.5)
abline(h = defest$point, col = 'purple')
# For this dataset, it is the highest of all the estimates.

legend('bottomright', col = c(1:3, 'purple'), 
       legend = c("Article's estimate", 'Dixon-Mood', 'Reversals (Wetherill)', 'Standard (CIR)'), 
       lty = 1, bty='n', cex = 0.8)


par(op) # Back to business as usual ;)

upndown documentation built on April 3, 2025, 10:57 p.m.