7. Diagnostic plots
This package implements the stage life testing (SLT)
model of Laumen & Cramer (2019, Naval Research Logistics
66:632-647), which unifies progressive censoring (with fixed censoring
times) and step-stress accelerated life testing. Units begin on a
baseline stage and, at a fixed stage-change time tau1,
R*_1 of the units still surviving are randomly selected and
moved to an elevated stage; the rest continue failing under the original
stage conditions. [rslt()] is this package’s only
data-generating function, implementing the paper’s own Procedure 2.3.
R*_1 is generated by one of two rules:
design = "typeP" transfers a fixed target
proportion pi1 of survivors:
R*_1 = floor(pi1 * (n - D1)), where D1 is the
number of stage-0 failures observed before tau1.design = "typeM" transfers a fixed target
count R0_1, capped by the number of survivors
actually available.Setting pi1 = 1 (every survivor transferred) collapses
the model to the classical “simple step-stress” design – Remark 2.2(4)
of the paper – which is why no separate function is needed for that
case: it is simply
rslt(..., design = "typeP", prop_transfer = 1). This is
also the structure of a typical published step-stress data set
(see Section 8), so it is the default prop_transfer used
throughout this vignette’s non-Section-5 examples for simplicity, even
though the package’s actual namesake model is the general
random-withdrawal case.
The two stages are linked by the cumulative exposure
model: a unit that has survived to time tau under
the baseline stage is treated, after switching, as if it had already
accumulated an equivalent age
\[ \nu = Q\big(F(\tau;\theta_0);\ \theta_1\big) \]
under the elevated-stage distribution, where
F/Q are the CDF/quantile function of a common
two-parameter baseline family \theta = (\lambda, \alpha)
(lambda a rate/scale-type parameter, alpha a
shape parameter). [slt_equiv_age()] computes
\nu directly, and this exactly reproduces equation (1.1) of
the paper. [dslt()], [pslt()] and
[qslt()] give the density/CDF/quantile function of a single
unit’s marginal lifetime under this model – but note that these are only
the correct pooled reference distribution for the combined
Y/Z sample when
prop_transfer = 1; see the note in
[slt_gof_ks()] for why.
A parameterisation note for anyone comparing directly against
the paper’s tables: the paper parameterises the
Exponential/Weibull distributions by their mean
\theta (F(x) = 1 - exp(-x/\theta)), whereas
this package parameterises every family by a rate-like
lambda (F(x) = 1 - exp(-\lambda x) for the
alpha=1 case), so lambda = 1/\theta. Section 5.1 below
reproduces the paper’s own worked example numerically, using this
correspondence.
slt_families()
#> [1] "weibull" "powerlindley" "explindley" "expexponential"
#> [5] "moexponential" "lomax" "expteissier" "loglogistic"
#> [9] "chen" "lognormal" "gamma" "powermuth"fam <- slt_dist("weibull")
fam$p(1, lambda = 1, alpha = 2)
#> [1] 0.6321
fam$q(0.5, lambda = 1, alpha = 2)
#> [1] 0.8326| Family | alpha = 1 reduces to |
|---|---|
weibull |
Exponential |
expexponential |
Exponential |
moexponential |
Exponential |
gamma |
Exponential |
powerlindley |
Lindley |
explindley |
Lindley (same limit as powerlindley) |
lomax |
(boundary member of the Lomax family) |
expteissier |
baseline Teissier |
powermuth |
baseline Teissier (same limit as expteissier) |
loglogistic |
Log-Logistic, shape 1 |
chen |
Chen, shape 1 |
lognormal |
Log-normal, log-sd 1 |
case = "shape1"/"shape11" (Section 1.2)
fixes the elevated-stage shape parameter to exactly 1 – as the table
shows, this only makes the elevated stage Exponential for
weibull, expexponential,
moexponential, and gamma.
powerlindley, expexponential,
lomax, and expteissier are the paper’s own
baseline-stage extension points; the remaining seven families are this
package’s own additions, covering hazard shapes not available among the
first four (a non-monotonic hazard for loglogistic; a
bathtub-shaped hazard for chen;
gamma/lognormal, two of the most standard
lifetime distributions;
explindley/moexponential, the “exponentiate
the cdf” and “tilt parameter” siblings of
powerlindley/expexponential respectively). See
?slt_dist for the pdf/cdf and literature reference of
each.
With \theta_0 = (\lambda_0, \alpha_0) at the baseline
stage and \theta_1 = (\lambda_1, \alpha_1) at the elevated
stage, four parameter-sharing structures are logically possible. This
package supports three of them:
slt_cases()
#> case n_par parameters
#> 1 common 3 lambda0, lambda1, alpha
#> 2 shape1 3 lambda0, alpha0, lambda1
#> 3 shape11 2 lambda0, lambda1
#> description
#> 1 theta0=(lambda0,alpha), theta1=(lambda1,alpha) [shared shape]
#> 2 theta0=(lambda0,alpha0), theta1=(lambda1,1) [stage-1 shape fixed at 1]
#> 3 theta0=(lambda0,1), theta1=(lambda1,1) [both shapes fixed at 1]The fourth, fully unrestricted structure (\alpha_0 and
\alpha_1 both free, 4 parameters) is not
offered: with a single change-point and no replicated stress levels,
that structure is empirically non-identifiable – the profile
log-likelihood in (\alpha_0, \alpha_1) is extremely flat,
producing unstable estimates and frequent optimiser non-convergence.
"shape1" (baseline shape free, elevated shape fixed at 1)
is the structure most useful as a primary working model;
"shape11" (both shapes fixed at 1) is its natural nested
null, used by [slt_lrt()]; "common" (a single
shared shape, both scales free) is a well-identified alternative.
[slt_reduce_params()] converts a full
(theta0, theta1) pair into a case’s free parameter vector –
mostly useful in simulation studies, to put a known “true” parameter set
into the same format fit_slt() returns, for comparing
against the fitted par:
The classical simple-step-stress special case
(prop_transfer = 1):
set.seed(2026)
theta0 <- c(0.05, 1.5) # (lambda0, alpha0)
theta1 <- c(0.20, 1) # (lambda1, alpha1 = 1) -- consistent with case = "shape1"
tau1 <- slt_tau("weibull", theta0, p = 0.5) # change-point at the stage-0 median
tau1
#> [1] 5.771
d <- rslt(n = 200, family = "weibull", theta0 = theta0, theta1 = theta1, tau1 = tau1, prop_transfer = 1)
str(d)
#> List of 4
#> $ Y : num [1:112] 0.0831 0.192 0.2282 0.242 0.4735 ...
#> $ Z : num [1:88] 5.79 5.83 5.85 5.97 6.01 ...
#> $ n_early : int 112
#> $ n_transferred: int 88The general random-withdrawal case – the paper’s actual SLT model –
uses the same function with a smaller prop_transfer
("typeP", a target proportion) or
design = "typeM" (a target count):
d_palt <- rslt(n = 200, family = "weibull", theta0 = theta0, theta1 = theta1,
tau1 = tau1, design = "typeP", prop_transfer = 0.5)
str(d_palt)
#> List of 4
#> $ Y : num [1:149] 0.233 0.38 0.584 0.592 0.657 ...
#> $ Z : num [1:51] 5.94 5.94 6.02 6.05 6.06 ...
#> $ n_early : int 98
#> $ n_transferred: int 51
d_palt_m <- rslt(n = 200, family = "weibull", theta0 = theta0, theta1 = theta1,
tau1 = tau1, design = "typeM", R0_1 = 100)
str(d_palt_m)
#> List of 4
#> $ Y : num [1:100] 0.233 0.323 0.324 0.337 0.506 ...
#> $ Z : num [1:100] 5.8 5.81 5.83 5.83 5.88 ...
#> $ n_early : int 95
#> $ n_transferred: int 100fit <- fit_slt(d, family = "weibull", case = "shape1", tau1 = tau1)
fit
#> Two-stage lifetime model fit ('slt_fit')
#> Family : Weibull
#> Case : shape1
#> tau1 : 5.771
#> n : 200 (n0 = 112, n1 = 88)
#> logLik : -553 AIC: 1112 BIC: 1122
#>
#> Coefficients (natural scale):
#> Estimate Std. Error
#> lambda0 0.0749 0.0158
#> alpha0 1.3669 0.1182
#> lambda1 0.2200 0.0235
summary(fit)
#> Summary of two-stage lifetime model fit
#> Family : Weibull
#> Case : shape1
#> n = 200 logLik = -553.01 AIC = 1112.0201 BIC = 1121.915
#>
#> Coefficients (natural scale, 95% Wald log-normal CI):
#> Estimate Std. Error Lower Upper
#> lambda0 0.0749 0.0158 0.0495 0.1132
#> alpha0 1.3669 0.1182 1.1537 1.6194
#> lambda1 0.2200 0.0235 0.1786 0.2712
#>
#> theta0 (stage 0) = (0.0749, 1.3669)
#> theta1 (stage 1) = (0.22, 1)fit_slt() runs a Nelder-Mead warm start followed by BFGS
(with a numerical Hessian for inference) from several perturbed starting
values, keeping the converged fit with the highest log-likelihood – this
guards against both local optima and the occasional optimiser failure.
Standard S3 methods work as expected:
coef(fit)
#> lambda0 alpha0 lambda1
#> 0.07485 1.36688 0.22004
vcov(fit)
#> lambda0 alpha0 lambda1
#> lambda0 2.497e-04 -1.670e-03 3.238e-13
#> alpha0 -1.670e-03 1.398e-02 -2.312e-12
#> lambda1 3.238e-13 -2.312e-12 5.502e-04
logLik(fit)
#> 'log Lik.' -553 (df=3)
AIC(fit)
#> [1] 1112
BIC(fit)
#> [1] 1122One family-specific note: lomax’s shape parameter has a
notoriously high-variance MLE at small-to-moderate sample sizes (a well
known feature of the family, not a defect of the optimiser). If you plan
to fit family = "lomax", prefer larger samples and/or the
bootstrap (rather than the asymptotic normal) confidence intervals
below.
Asymptotic (Wald) intervals: type = "lognormal" (the
default) builds the interval on the log scale and exponentiates,
guaranteeing a positive interval; type = "normal" applies
the delta-method standard error directly on the natural scale and can
occasionally produce a negative lower bound for a
positive parameter, especially at small n:
confint(fit) # type = "lognormal" (default), always positive
#> 2.5% 97.5%
#> lambda0 0.04949 0.1132
#> alpha0 1.15375 1.6194
#> lambda1 0.17856 0.2712
confint(fit, type = "normal") # can occasionally go negative for small n
#> 2.5% 97.5%
#> lambda0 0.04388 0.1058
#> alpha0 1.13517 1.5986
#> lambda1 0.17407 0.2660Bootstrap alternatives – non-parametric (case resampling) and
parametric (simulate from the fitted model) – each with
percentile, bootstrap-t, and
log-scale bootstrap-t variants. The plain
bootstrap-t interval has the same occasional-negative-bound
issue as the natural-scale asymptotic interval above, for the same
reason (it is an ordinary linear pivot); the log-scale version applies
the identical pivot construction on log(theta) and
exponentiates, guaranteeing positivity, directly analogous to why
type = "lognormal" is preferred above. The percentile
interval is already always positive (it is built from quantiles of the
always-positive bootstrap estimates themselves).
b <- slt_bootstrap_ci(fit, B = 1000, method = "nonparametric")
b$percentile # always positive
b$boot_t # can occasionally go negative for small n, like type = "normal" above
b$log_boot_t # always positive, like type = "lognormal" above
slt_bootstrap_ci(fit, B = 1000, method = "parametric") # prop_transfer = 1 default
slt_bootstrap_ci(fit, B = 1000, method = "parametric", prop_transfer = 0.5) # genuine partial withdrawalThis section validates the package’s fidelity to the paper directly, and documents where the package’s inferential tools go beyond, or fall short of, the paper’s own results.
Table 2 of Laumen & Cramer (2019) gives a "typeP"
SLT sample (n=16, tau1=15, true
theta0=40, theta1=20, i.e. Exponential
means – see the parameterisation note in Section 1). Entering
that data directly and fitting the pure-Exponential model
(case = "shape11") recovers the paper’s own reported MLEs
exactly:
Y <- c(2.55, 7.77, 9.01, 9.23, 12.34, 18.61, 19.39, 27.54, 39.67, 43.25, 91.46)
Z <- c(24.51, 30.73, 36.88, 55.26, 65.39)
fit_paper <- fit_slt(list(Y = Y, Z = Z), family = "weibull", case = "shape11", tau1 = 15)
c(theta0_hat = 1 / fit_paper$par["lambda0"], theta1_hat = 1 / fit_paper$par["lambda1"])
#> theta0_hat.lambda0 theta1_hat.lambda1
#> 32.35 27.55
# Paper reports: theta0_hat = 32.35, theta1_hat = 27.55slt_prob_stage1()
(Lemma 4.1)The MLE of the elevated-stage parameter only exists when at least one
unit is actually transferred and fails there (R*_1 > 0).
slt_prob_stage1() implements the paper’s Lemma 4.1 for this
probability, and exactly reproduces the paper’s own Table 5:
slt_prob_stage1(12, 0.5, c(1, 1), design = "typeP", pi1 = 0.25) # paper: n*=8, prob=0.9864756680
#> $n
#> [1] 12
#>
#> $tau1
#> [1] 0.5
#>
#> $F0_tau1
#> [1] 0.3935
#>
#> $design
#> [1] "typeP"
#>
#> $pi1
#> [1] 0.25
#>
#> $R0_1
#> NULL
#>
#> $n_star
#> [1] 8
#>
#> $prob
#> [1] 0.9865
slt_prob_stage1(12, 0.5, c(1, 1), design = "typeM", R0_1 = 6) # paper: n*=11, prob=0.9999862301
#> $n
#> [1] 12
#>
#> $tau1
#> [1] 0.5
#>
#> $F0_tau1
#> [1] 0.3935
#>
#> $design
#> [1] "typeM"
#>
#> $pi1
#> NULL
#>
#> $R0_1
#> [1] 6
#>
#> $n_star
#> [1] 11
#>
#> $prob
#> [1] 1This is useful for planning: before running a real experiment, check
that your intended n, tau1, and
pi1/R0_1 give a comfortably high
prob under a plausible guess at theta0.
Because the reference distribution’s parameters were estimated from
the same data, the classical Kolmogorov-Smirnov p-value is
anti-conservative (the well known Lilliefors problem).
[slt_gof_ks()] reports both the naive value and a
parametric-bootstrap-corrected one:
"shape11" is nested in both "shape1" and
"common" (fixing the remaining free shape(s) to 1), giving
a likelihood-ratio test of H0: shape = 1:
fit0 <- fit_slt(d, "weibull", "shape11", tau1 = tau1)
slt_lrt(fit, fit0)
#> Likelihood-ratio test (Weibull): H1 = 'shape1' vs. H0 = 'shape11'
#> LR statistic = 11.6477, df = 1, p-value = 0.0006428For exploratory comparison across any set of family/case fits (not
necessarily nested), use [slt_model_comparison()]:
fitc <- fit_slt(d, "weibull", "common", tau1 = tau1)
slt_model_comparison(list("shape1" = fit, "shape11" = fit0, "common" = fitc))
#> model family case n_par logLik AIC BIC converged
#> 1 shape1 Weibull shape1 3 -553.0 1112 1122 TRUE
#> 2 common Weibull common 3 -553.5 1113 1123 TRUE
#> 3 shape11 Weibull shape11 2 -558.8 1122 1128 TRUEFour real step-stress data sets are bundled with the package (see
?slt_han_kundu_2014 and related help pages). Here we
compare all twelve baseline distributions under
case = "common" on the Han & Kundu (2014) data set
(case = "common" estimates a single shared shape rather
than fixing the elevated-stage shape at exactly 1, which suits an
across-family comparison better – see Section 1.2):
data(slt_han_kundu_2014)
d_real <- list(Y = slt_han_kundu_2014$y0, Z = slt_han_kundu_2014$y1)
tau1_real <- slt_han_kundu_2014$tau1
fits <- lapply(slt_families(), function(fam) {
fit_slt(d_real, family = fam, case = "common", tau1 = tau1_real, n_starts = 5)
})
names(fits) <- slt_families()
slt_model_comparison(fits)
#> model family case n_par logLik AIC BIC
#> 1 expteissier Exponentiated Teissier common 3 -184.3 374.6 378.9
#> 2 powermuth Power Muth (Power Teissier) common 3 -184.7 375.3 379.6
#> 3 weibull Weibull common 3 -185.0 376.1 380.4
#> 4 moexponential Marshall-Olkin Exponential common 3 -185.2 376.3 380.6
#> 5 explindley Exponentiated Lindley common 3 -185.3 376.7 381.0
#> 6 gamma Gamma common 3 -185.4 376.8 381.1
#> 7 expexponential Exponentiated Exponential common 3 -185.5 377.0 381.3
#> 8 powerlindley Power Lindley common 3 -185.6 377.2 381.5
#> 9 chen Chen (2000) common 3 -185.7 377.4 381.7
#> 10 lomax Lomax common 3 -186.3 378.6 382.9
#> 11 loglogistic Log-Logistic (Fisk) common 3 -187.7 381.4 385.7
#> 12 lognormal Log-normal common 3 -188.1 382.2 386.5
#> converged
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 TRUE
#> 5 TRUE
#> 6 TRUE
#> 7 TRUE
#> 8 TRUE
#> 9 TRUE
#> 10 TRUE
#> 11 TRUE
#> 12 TRUEbest_fit <- fits[[slt_model_comparison(fits)$model[1]]]
summary(best_fit)
#> Summary of two-stage lifetime model fit
#> Family : Exponentiated Teissier
#> Case : common
#> n = 31 logLik = -184.3182 AIC = 374.6364 BIC = 378.9383
#>
#> Coefficients (natural scale, 95% Wald log-normal CI):
#> Estimate Std. Error Lower Upper
#> lambda0 0.0015 0.0004 0.0009 0.0024
#> lambda1 0.0162 0.0033 0.0108 0.0243
#> alpha 0.5476 0.1419 0.3295 0.9102
#>
#> theta0 (stage 0) = (0.0015, 0.5476)
#> theta1 (stage 1) = (0.0162, 0.5476)
plot_cdf(best_fit)
Not every family need converge for every case/data set: a fixed
elevated-stage shape of exactly 1
(case = "shape1"/"shape11") can be a poor
match for a family whose shape parameter’s natural scale depends on the
data (e.g. chen, where the shape is an exponent on the
failure time feeding directly into an exponential), in which case
[slt_model_comparison()] reports
converged = FALSE/NA for that family rather
than a misleading result; trying case = "common", as above,
is a reasonable first thing to check.