---
title: "Tools for Experimenters"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Tools for Experimenters}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

Beyond fitting, `exnexSurv` provides a small set of inference and
model-selection utilities that operate directly on the posterior draws already
stored in a fitted `exnex_surv` object. This vignette walks through each of
them using a small simulated example.

## A quick simulated example

We first generate a basket-trial dataset with `simulate_data()`. By default it
creates `K = 9` baskets with `n = 30` patients each, a single covariate, and
lets us flag some baskets as "resistant" (their location is shifted away from
the healthy population) — exactly the situation the EXNEX model is designed to
handle.

```{r}
library(exnexSurv)
library(survival)

set.seed(42)
d <- simulate_data(
  n = 20,
  beta = 0.5,
  sigma = 1.1,
  outlier_baskets = c(2, 8),
  resist_delta = -1.0,
  censoring_rate = 0.3,
  seed = 1
)

head(d)
```

The dataset stores the true generating parameters as attributes:

```{r}
attr(d, "true_theta")
attr(d, "true_beta")
attr(d, "true_sigma")
```

Now fit the model. We keep the chain short so the vignette builds quickly; for
real work use more iterations and multiple chains.

```{r}
fit <- exnex_surv(
  Surv(time, event) ~ group + x1,
  data = d,
  iter = 2000,
  warmup = 1000,
  chains = 2,
  seed = 7
)
print(fit, show_trace = FALSE)
```

## Survival curves

`survival_curves()` evaluates the posterior survival function
$S(t) = \Pr(T > t)$ on a time grid, returning a data frame with the posterior
median and a credible band:

```{r}
curves <- survival_curves(fit)
head(curves)
```

`plot.survival_exnex()` draws the curves (requires `ggplot2`):

```{r fig.show = "hold"}
plot(curves)
```

By default a single curve is produced with covariates fixed at zero. Passing
`newdata` evaluates one row per subject with its own covariate values and group:

```{r}
nd <- data.frame(group = d$group[1:3], x1 = c(0, 0.2, -0.1))
curves_nd <- survival_curves(fit, newdata = nd)
table(curves_nd$group)
```

## Posterior median survival time

For the log-normal AFT model the median survival time of a linear predictor
$\eta$ is simply $\exp(\eta)$. `median_survival()` reports its posterior
quantiles with credible intervals:

```{r}
median_survival(fit)
```

## Restricted mean survival time (RMST)

`rmst()` computes $\mathrm{RMST}(t_{max}) = \int_0^{t_{max}} S(t)\,dt$ for each
posterior draw and summarises the distribution:

```{r}
rmst(fit, tmax = 10)
```

## Model comparison with WAIC

`compute_waic()` evaluates the pointwise log-likelihood of the observed data
(here the censoring contribution is handled explicitly), giving WAIC, its
standard error, the log pointwise predictive density (`lpd`) and the penalty
`p_waic`:

```{r}
w <- compute_waic(fit)
str(w[c("waic", "se_elpd_waic", "lpd", "p_waic", "elpd_waic")])
```

`compare_waic()` contrasts several fits (e.g. a full model vs. a group-only
model) and returns them ordered by ascending WAIC:

```{r}
fit_group <- exnex_surv(
  Surv(time, event) ~ group,
  data = d,
  iter = 2000, warmup = 1000, chains = 2, seed = 7
)

compare_waic(without_covariate = fit_group, with_covariate = fit)
```

## Probability that one group beats another

`probability_superiority()` estimates $\Pr(\mathrm{summ}_a > \mathrm{summ}_b)$
draw-by-draw for one of three summaries: median survival, survival probability
at a fixed time, or RMST.

```{r}
# Median survival
probability_superiority(fit, a = 1, b = 3, function_of = "median")

# Survival probability at t = 3
probability_superiority(fit, a = 1, b = 2, function_of = "survival", times = 3)

# RMST up to t = 10
probability_superiority(fit, a = 1, b = 2, function_of = "rmst", tmax = 10)
```

In this example baskets 1 and 3 are both healthy (only 2 and 8 are resistant),
so `probability_superiority()` comparing them has no reason to prefer one over
the other and returns a probability close to 0.5. The exact numbers will vary
with the simulated data and the MCMC run.