crplot Advanced Options

Christopher Weld ceweld@email.wm.edu, Lawrence Leemis leemis@math.wm.edu

2022-05-20

Introduction

This is the second of two vignettes on the crplot function. It details advanced features which are helpful when troubleshooting plot issues, focusing on its jump-center algorithm parameters. For a basic comprehensive summary of crplot usage, please reference the first vignette titled crplot (users should familiarize with it prior to viewing this vignette). Both vignettes are available via links on the conf package webpage.

The crplot function generates a two-dimensional confidence region plot for the specified two-parameter parametric distribution, fitted to a dataset. Details of the crplot algorithm are available in its corresponding publication1.

Jump-Center Components

A jump-center is necessary to plot confidence region boundary points that are inaccessible via the radial profile technique (plotting algorithm of crplot) because a line joining them to the MLE crosses the confidence region boundary multiple times. The jump-center algorithm “repairs” the inaccessible region by providing an alternate point-of-reference (or alternate “center-point”, hence its name) from which to execute the radial profile plotting technique. See crplot vignette 1 of 2 for additional details.

Many advanced crplot options involve customizing its jump-center parameters. Before exploring those methods, it is first important to familiarize with jump-center components. The plot annotations below will aid in those definitions.

(note: select word coloring below is chosen to correspond with the appropriate reference points above)

Plot errors

Plot issues using crplot are relatively uncommon, but possible. Certain datasets and/or distributions pose greater plot challenges than others. In particular, small datasets and/or heavily censored datasets can result in unusual (strong deviations from elliptical, even to the point of non-convex) shapes, making their boundary difficult to map using the radial profile log likelihood technique2 employed by crplot.

When the crplot algorithm is stretched beyond its default limitations, adjusting a series of optional arguments can further extend its practical reach. Optional arguments in the crplot function enable algorithm customization, which is often enough to overcome plot issues when they arise.

Example

This heavily right-censored example produces a particularly challenging confidence region shape, enabling us to illustrate a suite of optional arguments available to its user to customize the crplot algorithm. Understanding the impact of optional argument adjustments will facilitate their use elsewhere when appropriate.

Heavily Right-Censored Dataset

Bearing cage fracture data is taken from Abernethy et. al. (1983)3. Among 1703 samples there are six observed failures. The remaining 1697 right-censored samples occur over a variety of right-censored values between 50–2050 hours. Both data subsets are given below with the variables bc_obs and bc_cen, respectively.

Default Confidence Region Results

Confidence regions corresponding to a variety of parametric distributions for this heavily right-censored dataset are given below. Each plot uses a level of significance of alpha = 0.1 (a 90% confidence region).

Investigating Plot Issues

The gamma distribution confidence region plot for the bearing cage fracture dataset takes a suspect shape. Nearer the bottom of its plot there are long segments that are seemingly missing confidence region boundary points. Relatively sharp vertex angles flank each of these line segments, a characteristic indicative of a plot issue.

To assess what is happening, jump-center reference points are plot using showjump = TRUE. Additionally, its plot information is returned using the optional argument jumpinfo = TRUE. The jumpinfo argument is used in lieu of its alternative info = TRUE because we want jump-center repair information returned in addition to the final confidence region plot information. It is then available for subsequent use, as shown by plotting those points on the “without repairs” plot below.

Setting jumpinfo = TRUE returns a repair attribute, which is a list relevant to all jump-center repairs. The repair labels are output to the console above, and given in greater detail next.

Prefixes associated with repair are:

Suffixes associated with repair are:

Seeing its jump-center location enables a targeted strategy—using the optional arguments given next—to improve the plot when necessary.

Adjusting Optional Arguments

Optional arguments whose adjustment can refine crplot results include:

Optional argument adjustments can be applied separately or in combination, each serving a unique purpose. Combinations of jumpshift and/or jumpuphill are the recommended first option to address plot regions that remain inaccessible to its jump-center.

The fourth quadrant repairs (with respect to the MLE) of the aforementioned example will illustrate optional argument uses below since they are more easily visible than its second quadrant repairs.

\(\circ\) jumpuphill

Increasing this optional argument from its default value, jumpuphill = min(alpha, 0.01), moves the jump-center nearer the MLE and further from its confidence region boundary. Doing so in our example will enable the fourth quadrant jump-center to achieve a better line-of-sight to the bottom portion of its confidence region (where a perceived lack of plot points appears problematic).

The confidence region’s bottom-most curve benefits substantially from this jumpuphill adjustment, but a problematic area (lacking plot points) remains in the upper curve at higher \(\theta\) values.

\(\circ\) jumpshift

Adjusting this optional argument from its default value, jumpshift = 0.5, shifts the radial azimuth (with respect to the MLE) to its jump-center. Decreasing its value generally “pushes” the jump-center location further along its confidence region boundary, nearer the inaccessible region. Doing so in our example enables the fourth quadrant jump-center to achieve a better line-of-sight to the upper portion of its confidence region curve at higher \(\theta\) values (where a perceived lack of plot points appears problematic).

The confidence region’s upper curve (at higher \(\theta\) values) benefits substantially from this jumpshift adjustment, but does so at the expense of its bottom curve. The bottom curve (at minimal \(\kappa\) values) plot difficulties are compounded by the perspective taken by this new jump-center location.

Caution: although decreasing jumpshift generally improves access to inaccessible regions, inevitably a jumpshift value too near zero will result in multiple solutions for its jump-center at the specified jumpuphill value. Under those circumstance, the point nearer its MLE is typically returned, resulting in a jump-center whose location is not between the inaccessible boundary points, and is therefore unable to plot its inaccessible region. Using jumpshift = 0.01 rather than jumpshift = 0.1 with the syntax above demonstrates this undesirable paradigm, and is left as an exercise to the reader. Another example of a jumpshift value that is prohibitively low for its circumstance is given later in this vignette, in the Optional Argument Combinations section.

\(\circ\) ellipse_n

Adjusting this optional argument from its default value, ellipse_n = 4, combines the elliptically-oriented plot technique with the default smoothing-search algorithm. Details of the elliptically-oriented plot algorithm are given in the first crplot vignette. It increases the final number of plot points, albeit at a computational expense. Used here, it improves but does not entirely correct our example’s plot concerns.

\(\circ\) maxcount

The maxcount parameter specifies the number of smoothing-search iterations crplot performs before terminating its algorithm. The default value (30) is typically sufficient to plot all accessible regions to its maximum degree tolerance constraint (maxdeg), so reaching this limit is indicative of an inaccessible region plot issue. On rare occasion however, increasing maxcount may improvement select plot regions, albeit at a noteworthy computational cost. This is true for our example plot, whose bottom region continues to (slowly) populate with relevant points as maxcount increases.

Increasing maxcount is effective at identifying additional relevant points here due to poor circumstance surrounding its fourth quadrant jump-center position. The smoothing search algorithm iteratively locates additional plot points based on azimuths (or angles) associated with the mid-point of its existing points (wherever maxdeg angle constraints are not yet met). The mid-point between adjacent points along its bottom curve are so offset and spread from the jump-center that, despite aiming towards the mid-point of its long vacant region, additional points are successively added much nearer its very far-right existing point.

\(\circ\) Optional argument combinations

Combinations of the aforementioned optional arguments are typically capable of correcting the most challenging plot circumstance. For the bearing cage fracture data example, jumpuphill and jumpshift parameter adjustments—influencing the jump-center distance and direction from the MLE, respectively—prove sufficient to adequately repair plot issues in its fourth quadrant. One solution, using jumpshift = 0.07 and jumpuphill = 0.05, is given below.

A closer review of its lower-right and upper-left plot regions is appropriate to ensure those extremes are adequately captured. Those plots reveal that the quadrant four repairs were successful, but applying its customizations uniformly for all jump-center locations results in an inaccurate plot region for quadrant two. It encounters an issue discussed as “Caution: …” in the jumpshift section above. Specifically, its decreased jumpshift value causes the MLE-to-jump-center azimuth to pass too near its inaccessible region right boundary. It subsequently fails to locate its jump-center by the “far” confidence region boundary, and instead crosses another solution for the jumpuphill parameter before its inaccessible region.

Different jumpshift and jumpuphill values for quadrant two and four repairs are necessary to adequately plot each of its extreme regions to a high level of detail. Itemizing customizations in this manner is possible, and detailed in the next section.

Varying jumpshift and jumpuphill values by quadrant.

A plot with multiple jump-center repairs may warrant a different jumpshift and jumpuphill value for each respective jump-center. Given scalar jumpshift and/or jumpuphill arguments, the crplot algorithm assumes its value applies regardless of the jump-center repair location (or quadrant). To customize these parameters to jump-center location-specific values, assign jumpshift and/or jumpuphill a vector of length four; its values then correspond uniquely to the four respective quadrants surrounding the MLE.

When providing a quadrant-specific vector of jumpshift or jumpuphill values, any placeholder value is acceptable for quadrants with no jump-center. For our example, quadrants one and three contain no jump-center repairs and we can therefore arbitrarily assign the default values of 0.5 and 0.01 for jumpshift and jumpuphill in those respective locations.

Through separate analysis (analogous to that for quadrant four), we determine that the quadrant two jump-center adequately plots its associated region using jumpuphill = 0.25 (keeping the default value for jumpshift = 0.5). Notice how the maximum \(\kappa\) value increases as a result. This solution is combined with the quadrant four solution (jumpshift = 0.07 and jumpuphill = 0.05) below.

Satisfied with plot accuracy in all regions, the final \(90\%%\) confidence region for the heavily right-censored bearing fracture data is given:

Trouble-shooting R Uniroot Errors

Albeit uncommon, if crplot is unable to return a confidence region plot the following message prints to the console:

#> [1] "--------------------------------------------------------------------------------------"   
#> [1] "R uniroot failure searching for confidence region boundary---challenging parameters or shape."   
#> [1] "Unable to produce a confidence region for the given sample and/or parameterization.   "   
#> [1] "--------------------------------------------------------------------------------------"

The R uniroot function numerically solves for confidence region boundary points within crplot. No plot returns when this numeric method fails. Failures can occur due-to, or despite uniroot parameters upper and lower properly bracketing the confidence region (uniroot solves using a bisection method). This typically occurs under the most challenging plot circumstance such as heavily censored and/or small sample sizes, or scales that approach the limits of R’s finite precision arithmetic.

Good options to troubleshoot this error include: increasing alpha, setting repair = FALSE, and adjusting the jumpshift and/or jumpuphill parameters. Increasing alpha is effective because it reduces the confidence region size, and this is can remedy uniroot upper-bound issues. The latter two are effective options because both influence the jump-center repairs algorithm. The jump-center repairs algorithm is usually the source of uniroot errors since it addresses the most “extreme” corners of the confidence region plot, where uniroot numeric issues tend to occur. Setting repair = FALSE may enable return of a working-copy of the confidence region (complete with exception of its radially inaccessible regions). To attain a complete confidence region, adjusting jumpshift and/or jumpuphill inherently will assign its jump-center a different location, from where its uniroot calculations may not experience similar numeric difficulties.

The example below demonstrates circumstance resulting in a plot error, and attainable solutions using the strategy described above.

# crplot is unable to plot this 98% confidence region 
crplot(dataset = c(1.5, 2), alpha = 0.01, distn = "invgauss")
#> [1] "--------------------------------------------------------------------------------------"   
#> [1] "R uniroot failure searching for confidence region boundary---challenging parameters and/or shape."   
#> [1] "Unable to produce a confidence region for the given sample and/or parameterization.   "   
#> [1] "--------------------------------------------------------------------------------------"

The three troubleshooting options given above are illustrated next. For this example, each successfully produces a confidence region plot.

# a plot without jump-center repairs is attainable, but its 3rd and 4th quadrants, 
# relative to the MLE, are in need of jump-center repairs
crplot(dataset = c(1.5, 2), alpha = 0.01, distn = "invgauss", repair = FALSE, 
       sf = c(3, 3), main = "without repairs")
#> [1] "99% confidence region plot complete; made using 131 boundary points."
# a complete plot is returned by increasing alpha to values >= 0.03 
crplot(dataset = c(1.5, 2), alpha = 0.05, distn = "invgauss", main = "95% CR", 
       sf = c(3, 3))
#> [1] "95% confidence region plot complete; made using 276 boundary points."
# adjusting jumpshift and jumpuphill parameters
x <- crplot(dataset = c(1.5, 2), alpha = 0.02, "invgauss", jumpinfo = TRUE, 
       sf = c(3, 3), jumpshift = 0.1, jumpuphill = 0.001, main = "98% CR")
#> [1] "98% confidence region plot complete; made using 234 boundary points."
plot(c(x$mu, x$mu[1]), c(x$lambda, x$lambda[1]), type = "l", axes = FALSE, xlab = expression(mu), 
     ylab = expression(lambda), main = "98% CR")
axis(side = 1, at = format(range(x$mu), digits = 3))
axis(side = 2, at = format(range(x$lambda), digits = 3))

  1. Weld, C., Loh, A., and Leemis, L. (in press), “Plotting Likelihood-Ratio Based Confidence Regions for Two-Parameter Univariate Probability Models”, The American Statistician↩

  2. Jaeger, A. (2016), “Computation of Two- and Three-Dimensional Confidence Regions With the Likelihood Ratio”, The American Statistician, 70, 395–398↩

  3. Abernethy, R. B., Breneman, J. E., Medlin, C. H., and Reinman, G. L. (1983), “Weibull Analysis Handbook”, Aero Propulsion Laboratory, 43–47↩