Introduction to stagelife: Two-Stage Life Testing

1. The model

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

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]

slt_reduce_params(theta0 = c(0.2, 1.4), theta1 = c(0.5, 1.4), case = "common")
#> lambda0 lambda1   alpha 
#>     0.2     0.5     1.4

2. Simulating data

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 88

The 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 100

3. Maximum likelihood estimation

fit <- 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] 1122

One 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.

4. Confidence intervals

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.2660

Bootstrap 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 withdrawal

5. Design theory from the paper: reproducing Laumen & Cramer (2019)

This 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.

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.55

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] 1

6. Goodness of fit and model comparison

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:

slt_gof_ks(fit, B = 1000)

"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.0006428

For 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      TRUE

7. Diagnostic plots

plot_cdf(fit)
plot_survival(fit)
plot_pp(fit)
plot_qq(fit)

8. Real-data example

Four 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      TRUE

best_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.