---
title: "'reslife' Package Demonstration"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{'reslife' Package Demonstration}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

The mathematical derivations, integral forms, and distribution specifications implemented in this package are detailed in the companion paper: [Wang et al. (2025)](https://arxiv.org/pdf/2308.07410.pdf).

```{r, include = FALSE}
knitr::opts_chunk$set(
  error = TRUE,
  collapse = TRUE,
  comment = "#>"
)
```

## Load the Package

```{r setup}
library(reslife)
library(flexsurv)
library(survival)
```

---

## Demonstration of the `reslifefsr` Function

The `reslifefsr()` function calculates residual life values directly from a fitted `flexsurvreg` object. It supports both baseline (intercept-only) models and models with patient-level covariate profiles via `newdata`.

### Baseline Verification via Numerical Integration

We begin by fitting an intercept-only generalized gamma model (`gengamma.orig`) to the breast cancer dataset `bc` from `flexsurv`. We then evaluate the mean residual life at an elapsed time of 1 year (`life = 1`).

```{r}
fs1 <- flexsurvreg(Surv(recyrs, censrec) ~ 1, data = bc, dist = "gengamma.orig")
r <- reslifefsr(obj = fs1, life = 1, type = "mean")
r
```

We can verify the accuracy of this analytical derivation against direct numerical integration of the survival function:
$$\text{MRL}(x) = \frac{1}{S(x)} \int_x^\infty S(t) \, dt$$

```{r}
# Extract parameters on their natural scale from the fitted model
b <- as.numeric(fs1$res["shape", "est"])
a <- as.numeric(fs1$res["scale", "est"])
k <- as.numeric(fs1$res["k", "est"])

# Define survival function
surv_gengamma <- function(t) {
  pgengamma.orig(t, shape = b, scale = a, k = k, lower.tail = FALSE)
}

# Numerical integration
integral_val <- integrate(surv_gengamma, lower = 1, upper = Inf)$value
sx_val <- surv_gengamma(1)
mx_integral <- integral_val / sx_val

all.equal(r, mx_integral)
```

The analytical formula in `reslifefsr()` matches numerical integration.

---

### Covariate Models with Multiple Predictors

When a model includes covariates, `reslifefsr()` evaluates residual life across patient profiles. We illustrate this using the `bc` dataset with an added covariate `age`:

```{r}
set.seed(42)
newbc <- bc
newbc$age <- rnorm(nrow(bc), mean = 65 - scale(newbc$recyrs, scale = FALSE), sd = 5)

fsr2 <- flexsurvreg(Surv(recyrs, censrec) ~ group + age, data = newbc, dist = "weibull")
```

Evaluating residual life across the fitted observations returns a vector with one value per observation:

```{r}
head(reslifefsr(obj = fsr2, life = 4, type = "mean"), 5)
```

We can also inspect the median residual life:

```{r}
head(reslifefsr(obj = fsr2, life = 4, type = "median"), 5)
```

Or extract the full summary table (`mean`, `median`, and a custom `quantile`):

```{r}
head(reslifefsr(obj = fsr2, life = 4, p = 0.8, type = "all"), 5)
```

---

### Predicting on New Covariate Profiles (`newdata`)

We can evaluate specific patient profiles by passing a `data.frame` to `newdata`:

```{r}
newdata <- data.frame(
  group = c("Medium", "Good", "Poor"),
  age   = c(43, 35, 39)
)

reslifefsr(obj = fsr2, life = 4, p = 0.6, type = "all", newdata = newdata)
```

`reslifefsr()` is invariant to column order:

```{r}
# Reordered columns produce identical results
newdata_reorder <- data.frame(age = c(43, 35, 39), group = c("Medium", "Good", "Poor"))
reslifefsr(obj = fsr2, life = 4, p = 0.6, type = "all", newdata = newdata_reorder)
```

Extra unmodeled columns are safely ignored:

```{r}
newdata_extra <- cbind(newdata, extra = c(100, 100, 100))
reslifefsr(obj = fsr2, life = 4, p = 0.6, type = "all", newdata = newdata_extra)
```

---

### Error Handling & Input Guards

If required covariates are missing, `reslifefsr()` informs the user:

```{r}
newdata_missing <- data.frame(group = c("Medium", "Good", "Poor"))
reslifefsr(obj = fsr2, life = 4, type = "all", newdata = newdata_missing)
```

Similarly, if an invalid factor level is supplied:

```{r}
newdata_bad_level <- data.frame(group = c("Medium", "Good", "InvalidLevel"), age = c(43, 35, 39))
reslifefsr(obj = fsr2, life = 4, type = "all", newdata = newdata_bad_level)
```

---

### Plotting Mean Residual Life over Time

We can evaluate how residual life changes over elapsed time $x$ for a specific patient profile:

```{r, fig.width=7, fig.height=4.5}
patient1 <- newbc[1, ]
time_grid <- 1:10

mrl_vals <- sapply(time_grid, function(t0) {
  reslifefsr(obj = fsr2, life = t0, type = "mean", newdata = patient1)
})

plot(time_grid, mrl_vals, type = "b", pch = 19, col = "steelblue",
     xlab = "Elapsed Survival Time (Years)",
     ylab = "Mean Residual Life (Years)",
     main = "Mean Residual Life Function for Patient Profile")
grid()
```

---

## Demonstration of the `residlife` Function

The standalone function `residlife()` computes residual life directly from distribution names and parameter vectors without requiring a `flexsurvreg` object.

Supported distributions include:
* `"exponential"` / `"exp"`
* `"weibull"`
* `"gamma"`
* `"gompertz"`
* `"lognormal"` / `"lnorm"`
* `"llogis"`
* `"gengamma"`
* `"gengamma.orig"`
* `"genf"`
* `"genf.orig"`

### Basic Usage with Vectorized Values

`residlife()` accepts a numeric vector of elapsed times:

```{r}
t_seq <- seq(1, 10, by = 1)
pars_weibull <- c(shape = 1.272, scale = 6.191)

residlife(values = t_seq, distribution = "weibull", parameters = pars_weibull, type = "mean")
```

### Flexible Parameter Naming & Ordering

Parameter order does not matter as long as names match:

```{r}
# Unordered names work identically
residlife(values = 5, distribution = "weibull", parameters = c(scale = 6.191, shape = 1.272))
```

For distributions with alternative parameterizations like Gamma (rate vs. scale), `residlife()` seamlessly converts between them:

```{r}
# Parameterized by scale
residlife(values = 5, distribution = "gamma", parameters = c(shape = 1.272, scale = 6.191))

# Parameterized by rate (rate = 1 / scale)
residlife(values = 5, distribution = "gamma", parameters = c(shape = 1.272, rate = 1 / 6.191))
```

### Summary Table for All Metrics

```{r}
residlife(values = c(1, 3, 5), distribution = "weibull", parameters = pars_weibull, p = 0.75, type = "all")
```

### Plotting MRL Directly from `residlife`

Because `residlife()` natively vectorizes over `values`, generating an MRL curve requires no manual loops:

```{r, fig.width=7, fig.height=4.5}
times <- seq(0.1, 10, by = 0.2)
mrl_curve <- residlife(values = times, distribution = "weibull", parameters = pars_weibull, type = "mean")

plot(times, mrl_curve, type = "l", lwd = 2, col = "darkred",
     xlab = "Elapsed Time (x)",
     ylab = "Mean Residual Life e(x)",
     main = "Weibull Mean Residual Life Curve")
grid()
```

---

## References

To cite the `reslife` package or the theoretical derivations in publications, please use:

* Wang, Z., Crawford, A., Lee, K. L., Lu, L., & Jaganathan, S. (2025). *reslife: Residual Lifetime Analysis Tool in R*. arXiv preprint [arXiv:2308.07410](https://arxiv.org/abs/2308.07410).

```bibtex
@article{wang2025reslife,
  title   = {reslife: Residual Lifetime Analysis Tool in R},
  author  = {Wang, Zekai and Crawford, Andrew and Lee, Ka Lok and Lu, Lin and Jaganathan, Srihari},
  journal = {arXiv preprint arXiv:2308.07410},
  year    = {2025},
  url     = {https://arxiv.org/abs/2308.07410}
}
```

---

## Session Information

```{r session-info}
sessionInfo()
```