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). 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.
bc_obs <- c(230, 334, 423, 990, 1009, 1510)
bc_cen <- c(rep(50, 288), rep(150, 148), rep(250, 124), rep(350, 111), rep(450, 106),
rep(550, 99), rep(650, 110), rep(750, 114), rep(850, 119), rep(950, 127),
rep(1050, 123), rep(1150, 93), rep(1250, 47), rep(1350, 41), rep(1450, 27),
rep(1550, 11), rep(1650, 6), rep(1850, 1), rep(2050, 2))
bc <- c(bc_obs, bc_cen)
cen <- c(rep(1, length(bc_obs)), rep(0, length(bc_cen)))
print(length(bc))
#> [1] 1703
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).
distns <- c("cauchy", "gamma", "llogis", "logis", "norm", "weibull")
par(mfrow = c(2, 3)) # plot in a 2-row, 3-column grid
for (i in 1:length(distns)) {
try(crplot(dataset = bc, alpha = 0.1, distn = distns[i], cen = cen,
sf = c(0, 2), ylas = 1, main = distns[i]))
}
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.
x <- crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, jumpinfo = TRUE, showjump = TRUE,
main = "with repairs")
#> [1] "Confidence region plot complete; made using 276 boundary points."
crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, repair = FALSE,
main = "without repairs")
#> [1] "Confidence region plot complete; made using 134 boundary points."
jumppoints <- t(matrix(c(x$repair$q2jumpxy, x$repair$q2jumpL, x$repair$q2jumpR,
x$repair$q4jumpxy, x$repair$q4jumpL, x$repair$q4jumpR), nrow = 2))
points(jumppoints, pch = 24, bg = rep(c("red", "green", "yellow")))
print(labels(x$repair))
#> [1] "q2jumpuphill" "q4jumpuphill" "q2jumpshift" "q4jumpshift"
#> [5] "q2jumpxy" "q4jumpxy" "q2jumpL" "q4jumpL"
#> [9] "q2jumpR" "q4jumpR" "q2gaptype" "q4gaptype"
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:
q1, q2, q3, q4.
Each label begins with an abbreviation specifying what quadrant (with
respect to the MLE) the jump-center repair is made. For example, a
q1 prefix indicates a jump-center repair located in the 1st
quadrant (the jump-center horizontal and vertical values are both
greater than the MLE’s respective values). In our example,
q2 and q4 prefixes indicate jump-center
repairs in the 2nd and 4th quadrant relative to the MLE (to its
upper-left and lower-right, respectively).
Suffixes associated with repair are:
jumpuphill quantifies the amount uphill (with regards
to level of significance) from the confidence region boundary that the
jump-center will locate,jumpshift quantifies the point along the vertical or
horizontal “gap” (as a fractional value between 0 and 1) that the
azimuth from the MLE to the jump-center will cross,jumpxy stores the coordinates for the jump-center,jumpL stores the left (with perspective taken from the
jump-center location towards its inaccessible region) confidence region
boundary point bordering its inaccessible region, andjumpR stores the left (with perspective taken from the
jump-center location towards its inaccessible region) confidence region
boundary point bordering its inaccessible region.
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:
jumpuphilljumpshiftellipse_nmaxcount
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).
crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, jumpuphill = 0.1, showjump = TRUE)
#> [1] "Confidence region plot complete; made using 277 boundary points."
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).
crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, jumpshift = 0.1, showjump = TRUE)
#> [1] "Confidence region plot complete; made using 292 boundary points."
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.
crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, showjump = TRUE, ellipse_n = 80)
#> [1] "Confidence region plot complete; made using 387 boundary points."
\(\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.
crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, maxcount = 50)
#> [1] "Confidence region plot complete; made using 476 boundary points."
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.
x <- crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, showjump = TRUE,
jumpshift = 0.07, jumpuphill = 0.05, jumpinfo = TRUE)
#> [1] "Confidence region plot complete; made using 240 boundary points."
plot(c(x$theta, x$theta[1]), c(x$kappa, x$kappa[1]), type = "l", axes = FALSE,
xlab = expression(theta), ylab = expression(kappa))
axis(side = 1, at = range(x$theta))
axis(side = 2, at = range(x$kappa))
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.
par(mar = c(4, 4, 1.5, 1))
# top-left confidence region extreme (zoomed-in)
plot(x$theta, x$kappa, xlim = c(min(x$theta), x$thetahat), ylim = c(x$kappahat, max(x$kappa)),
type = 'l', xlab = expression(theta), ylab = expression(kappa),
main = "top-left (zoom-in)")
# augment the top-left region graphic with plot points and jump-center information
points(x$theta, x$kappa)
jumppoints <- t(matrix(c(x$repair$q2jumpxy, x$repair$q2jumpL, x$repair$q2jumpR), nrow = 2))
points(jumppoints, pch = 24, bg = rep(c("red", "green", "yellow")))
# bottom-right confidence region extreme (zoomed-in)
plot(x$theta, x$kappa, ylim = c(min(x$kappa), 0.95), type = 'l',
xlab = expression(theta), ylab = expression(kappa), main = "bottom-right (zoom-in)")
points(x$theta, x$kappa)
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.
# top-left confidence region extreme (zoomed-in)
x <- crplot(bc, alpha = 0.1, distn = "gamma", cen = cen, showjump = TRUE, jumpinfo = TRUE,
jumpshift = c(0.5, 0.5, 0.5, 0.07), jumpuphill = c(0.01, 0.25, 0.01, 0.05),
xlim = c(min(x$theta), x$thetahat), ylim = c(x$kappahat, 4.15),
main = "top-left (zoom-in)")
#> [1] "Confidence region plot complete; made using 299 boundary points."
# bottom-right confidence region extreme (zoomed-in)
plot(x$theta, x$kappa, ylim = c(min(x$kappa), 0.95), type = 'l',
xlab = expression(theta), ylab = expression(kappa), main = "bottom-right (zoom-in)")
points(x$theta, x$kappa)
Satisfied with plot accuracy in all regions, the final \(90\%%\) confidence region for the heavily
right-censored bearing fracture data is given:
# final plot, with all necessary repairs complete
plot(c(x$theta, x$theta[1]), c(x$kappa, x$kappa[1]), type = "l", axes = FALSE,
xlab = expression(theta), ylab = expression(kappa), main = "gamma 90% CR for BF data (final)")
axis(side = 1, at = range(x$theta))
axis(side = 2, at = range(x$kappa))
points(x$thetahat, x$kappahat, pch = 3)
