---
title: "Introduction to stagelife: Two-Stage Life Testing"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to stagelife: Two-Stage Life Testing}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4
)
old_opts <- options(digits = 4)
```

```{r setup}
library(stagelife)
```

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

### 1.1 Baseline distributions

```{r}
slt_families()
```

```{r}
fam <- slt_dist("weibull")
fam$p(1, lambda = 1, alpha = 2)
fam$q(0.5, lambda = 1, alpha = 2)
```

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

### 1.2 Parameter-sharing cases

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:

```{r}
slt_cases()
```

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`:

```{r}
slt_reduce_params(theta0 = c(0.2, 1.4), theta1 = c(0.5, 1.4), case = "common")
```

## 2. Simulating data

The classical simple-step-stress special case (`prop_transfer = 1`):

```{r}
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

d <- rslt(n = 200, family = "weibull", theta0 = theta0, theta1 = theta1, tau1 = tau1, prop_transfer = 1)
str(d)
```

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):

```{r}
d_palt <- rslt(n = 200, family = "weibull", theta0 = theta0, theta1 = theta1,
                       tau1 = tau1, design = "typeP", prop_transfer = 0.5)
str(d_palt)

d_palt_m <- rslt(n = 200, family = "weibull", theta0 = theta0, theta1 = theta1,
                         tau1 = tau1, design = "typeM", R0_1 = 100)
str(d_palt_m)
```

## 3. Maximum likelihood estimation

```{r}
fit <- fit_slt(d, family = "weibull", case = "shape1", tau1 = tau1)
fit
summary(fit)
```

`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:

```{r}
coef(fit)
vcov(fit)
logLik(fit)
AIC(fit)
BIC(fit)
```

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`:

```{r}
confint(fit)                     # type = "lognormal" (default), always positive
confint(fit, type = "normal")    # can occasionally go negative for small n
```

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

```{r, eval = FALSE}
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.

### 5.1 Reproducing the paper's worked example exactly

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:

```{r}
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"])
# Paper reports: theta0_hat = 32.35, theta1_hat = 27.55
```

### 5.2 Planning a stage life test: `slt_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:

```{r}
slt_prob_stage1(12, 0.5, c(1, 1), design = "typeP", pi1 = 0.25)  # paper: n*=8,  prob=0.9864756680
slt_prob_stage1(12, 0.5, c(1, 1), design = "typeM", R0_1 = 6)    # paper: n*=11, prob=0.9999862301
```

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

## 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:

```{r, eval = FALSE}
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`:

```{r}
fit0 <- fit_slt(d, "weibull", "shape11", tau1 = tau1)
slt_lrt(fit, fit0)
```

For exploratory comparison across any set of family/case fits (not
necessarily nested), use [`slt_model_comparison()`]:

```{r}
fitc <- fit_slt(d, "weibull", "common", tau1 = tau1)
slt_model_comparison(list("shape1" = fit, "shape11" = fit0, "common" = fitc))
```

## 7. Diagnostic plots

```{r, fig.show = "hold"}
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):

```{r}
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)
```

```{r}
best_fit <- fits[[slt_model_comparison(fits)$model[1]]]
summary(best_fit)
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.

```{r, include = FALSE}
options(old_opts)
```