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

## ----setup--------------------------------------------------------------------
library(stagelife)

## -----------------------------------------------------------------------------
slt_families()

## -----------------------------------------------------------------------------
fam <- slt_dist("weibull")
fam$p(1, lambda = 1, alpha = 2)
fam$q(0.5, lambda = 1, alpha = 2)

## -----------------------------------------------------------------------------
slt_cases()

## -----------------------------------------------------------------------------
slt_reduce_params(theta0 = c(0.2, 1.4), theta1 = c(0.5, 1.4), case = "common")

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

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

## -----------------------------------------------------------------------------
fit <- fit_slt(d, family = "weibull", case = "shape1", tau1 = tau1)
fit
summary(fit)

## -----------------------------------------------------------------------------
coef(fit)
vcov(fit)
logLik(fit)
AIC(fit)
BIC(fit)

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

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

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

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

## ----eval = FALSE-------------------------------------------------------------
# slt_gof_ks(fit, B = 1000)

## -----------------------------------------------------------------------------
fit0 <- fit_slt(d, "weibull", "shape11", tau1 = tau1)
slt_lrt(fit, fit0)

## -----------------------------------------------------------------------------
fitc <- fit_slt(d, "weibull", "common", tau1 = tau1)
slt_model_comparison(list("shape1" = fit, "shape11" = fit0, "common" = fitc))

## ----fig.show = "hold"--------------------------------------------------------
plot_cdf(fit)
plot_survival(fit)
plot_pp(fit)
plot_qq(fit)

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

## -----------------------------------------------------------------------------
best_fit <- fits[[slt_model_comparison(fits)$model[1]]]
summary(best_fit)
plot_cdf(best_fit)

## ----include = FALSE----------------------------------------------------------
options(old_opts)