deltaInverse: Backend utility to calculate inverse (dose-finding) intervals, using local inversion and the Delta method

Description

Calculate left-bound to right-bound intervals for the dose point estimates, using local slopes at design points (places where observations exist) to invert the forward lower-upper bounds.

Usage

deltaInverse(
  isotPoint,
  target = (1:3)/4,
  intfun = morrisCI,
  conf = 0.9,
  adaptiveCurve = FALSE,
  minslope = 0.01,
  slopeRefinement = TRUE,
  finegrid = 0.05,
  globalCheck = TRUE,
  ...
)

Value

two-column matrix with the left and right bounds, respectively

Arguments

isotPoint: The output of an estimation function such as cirPAVA,doseFind, with the option full=TRUE. Should be at least a list of 3 doseResponse objects named input, output, shrinkage.
target: A vector of target response rate(s), for which the interval is needed. Default (since version 2.3.0) is the 3 quartiles ((1:3) / 4). If changed to NULL, interval will be returned for the $y$ values of isotPoint$output.
intfun: the function to be used for initial (forward) interval estimation. Default morrisCI (see help on that function for additional options).
conf: numeric, the interval's confidence level as a fraction in (0,1). Default 0.9.
adaptiveCurve: logical, should the CIs be expanded by using a parabolic curve between estimation points rather than straight interpolation? Default FALSE. Switch to TRUE recommended when adaptive design was used, and target is outside of (0.4, 0.6).
minslope: minimum local slope (subsequently normalized by the dose-spacing unit) considered positive, passed on to slope. Needed to avoid unrealistically broad intervals. Default 0.01.
slopeRefinement: (new to 2.3.0) logical: whether to allow refinement of the slope estimate, including different slopes to the left and right of target. Default TRUE. See Details.
finegrid: a numerical value used to guide how fine the grid of x values will be during slope estimation. Should be in (0,1) (preferably much less than 1). Default 0.05.
globalCheck: (new to 2.4.0) logical: whether to allow narrowing of the bound, in case the "global" bound (obtained via inverting the forward interval, and generally more conservative) is narrower. Default TRUE.
...: additional arguments passed on to quickIsotone

Details

This function is the "backend engine" for calculating confidence intervals for inverse (dose-finding) estimation. Methodologically this might be the most challenging task in the package. It is expected that most users will not interact with this function directly, but rather indirectly via the convenience wrapper quickInverse.

The Delta method in this application boils down to dividing the distance to the forward (vertical) bounds, by the slope, to get the left/right interval width. Both forward intervals and slopes are calculated across a standard set of $x$ values, then interpolated at horizontal cross-sections determined by target. Slope estimates are performed by slope.

Starting version 2.3.0, by default the slope estimate is different to the right and left of target. The intervals should now better accommodate the sharp slope changes that often happen with discrete dose-response datasets. Operationally, the intervals are first estimated via the single-slope approach described above. Then using a finer grid of $x$ values, weighted-average slopes to the right and left of the point estimate separately are calculated over the first-stage's half-intervals. The weights are hard-coded as quadratic (Epanechnikov).

An alternative and much simpler interval method (dubbed "global") is hard-coded into quickInverse, and can be chosen from there as an option. It is not recommended being far too conservative, and sometimes not existing. It is now also (since version 2.4.0) used in this function as a fallback upper bound on interval width.

Examples

Run this code

# Interesting run (#664) from a simulated up-and-down ensemble:
# (x will be auto-generated as dose levels 1:5)
dat=doseResponse(y=c(1/7,1/8,1/2,1/4,4/17), wt=c(7,24,20,12,17))
# The experiment's goal is to find the 30th percentile
quick1=quickIsotone(dat, adaptiveShrink = TRUE, target = 0.3)

# For inverse confidence intervals "the long way", 
#    we need a full CIR output object:
fwd1=cirPAVA(dat, full=TRUE, adaptiveShrink = TRUE, target = 0.3)
# Inverse intervals. 
# Note: 'target' specifies the y values at which the interval is calculated.
#       They are selected here based on the y range in which there are estimates.
yvals = c(seq(0.15, 0.3, 0.05), 0.33)
invDelta=deltaInverse(fwd1, target = yvals, adaptiveCurve = TRUE)
# stop()
# We added the adaptiveCurve option because the experiment's target is off-center,
#     and inverse-interval coverage tends to be lacking w/o that option.

### Showing the data and the estimates
par(mar=c(3,3,4,1), mgp=c(2,.5,0), tcl=-0.25)
# Following command uses plot.doseResponse()
plot(dat, ylim=c(0.05,0.55), 
     las=1, xlim=c(0,6.5), main="Inverse-Estimation CIs") 

# The true response function; true target is where it crosses the y=0.3 line
lines(seq(0,7,0.1), pweibull(seq(0,7,0.1), shape=1.1615, scale=8.4839), col=4, lwd=1.5)
abline(h=0.3, col=2, lwd=2, lty=3) ### The experiment's official target

# Forward CIs; the "global" inverse interval just draws horizontal lines between them
# To get these "global" intervals calculated for you at specific targets, choose 'delta=FALSE' 
#      when calling quickInverse()
lines(quick1$y, lwd=1.5, col='purple') 
lines(quick1$lower90conf,lty=2,col=3) 
lines(quick1$upper90conf,lty=2,col=3) 
# Note that only the upper forward bounds intersect the horizontal line at y=0.3.
#   Therefore, via the "global" approach there won't be a finite CI for the target estimate.

# Now, the default "local" inverse interval, which is finite for the range of estimated y values.
# In particular, it is finite for y=0.3.
# Note in the plot, how we make it equal to the "global" bound when the latter is narrower.
lines(invDelta[,1], yvals, lty=2, lwd=2)
lines(invDelta[,2], yvals, lty=2, lwd=2)

legend('topleft', pch=c(NA,NA,'X',NA,NA), lty=c(1,1,NA,2,2), 
       col=c(4,'purple', 1,1,3), lwd=c(1.5,1.5,0,2,1),
       legend = c('True Curve', 'CIR Curve', 'Observations', 
                  'Local Interval (default)',
                  'Forward/Global Interval'), bty='n')