---
title: "Unified IP-VOI algorithm in R"
format: html
execute:
  echo: true
  warning: false
  message: false
---

This file implements the branching algorithm in the manuscript. Replace the toy `sample_inputs()` and `run_model()` functions with the user's decision model. The branch labels match the manuscript: A = lower-expectation EVPI, B = lower-expectation EVPPI/EVSI, C = fixed-measure VOI envelope, D = p-box expectation bounds, followed by a stability step. The code uses base R only.

## 1. User-supplied model

```{r}
set.seed(123)

strategies <- c("A", "B")
K <- length(strategies)

# Each row of H represents one admissible probability measure P_h.
H <- expand.grid(mean_effect = c(-0.05, 0, 0.05),
                 sd_effect = c(0.08, 0.12),
                 corr = c(-0.3, 0, 0.3))
M <- nrow(H)
N <- 2000

sample_inputs <- function(h, N) {
  # Toy two-input model. Replace with model-specific input simulation under P_h.
  effect <- rnorm(N, mean = h$mean_effect, sd = h$sd_effect)
  cost_penalty <- rlnorm(N, meanlog = log(100), sdlog = 0.25)
  data.frame(effect = effect, cost_penalty = cost_penalty)
}

run_model <- function(theta) {
  # Return an N x K matrix of net benefits. Larger values are preferred.
  nb_A <- rep(0, nrow(theta))
  nb_B <- 5000 * theta$effect - theta$cost_penalty
  out <- cbind(A = nb_A, B = nb_B)
  out
}
```

## 2. Outer simulation over admissible probability measures

```{r}
simulate_outer <- function(H, N) {
  lapply(seq_len(nrow(H)), function(m) {
    theta <- sample_inputs(H[m, , drop = FALSE], N)
    B <- run_model(theta)
    list(h = H[m, , drop = FALSE], theta = theta, B = B)
  })
}

outer_sims <- simulate_outer(H, N)
```


## 3. Branch A: lower-expectation EVPI

This branch uses the lower-expectation decision rule. It estimates the current guaranteed value and the guaranteed value under perfect information.

```{r}
mu <- do.call(rbind, lapply(outer_sims, function(o) colMeans(o$B)))
colnames(mu) <- strategies
eta <- sapply(outer_sims, function(o) mean(apply(o$B, 1, max)))

V0_LE <- max(apply(mu, 2, min))
Vtheta_LE <- min(eta)
EVPI_LE <- Vtheta_LE - V0_LE

data.frame(V0_LE = V0_LE,
           Vtheta_LE = Vtheta_LE,
           EVPI_LE = EVPI_LE,
           current_strategy_LE = names(which.max(apply(mu, 2, min))))
```

## 4. Branch B: lower-expectation EVPPI or EVSI by policy search

Let `Z` denote information observed before the final strategy is chosen. For EVPPI, `Z` can be a parameter group. For EVSI, `Z` can be a simulated future-data summary. Here we use a toy binary information state.

```{r}
make_information_state <- function(theta) {
  ifelse(theta$effect > 0, "positive", "nonpositive")
}

policy_value <- function(B, z, policy) {
  chosen <- policy[as.character(z)]
  idx <- match(chosen, colnames(B))
  mean(B[cbind(seq_len(nrow(B)), idx)])
}

states <- c("nonpositive", "positive")
all_policies <- expand.grid(nonpositive = strategies, positive = strategies,
                            stringsAsFactors = FALSE)

policy_matrix <- matrix(NA_real_, nrow = nrow(all_policies), ncol = M)
for (p in seq_len(nrow(all_policies))) {
  policy <- unlist(all_policies[p, ])
  names(policy) <- states
  for (m in seq_len(M)) {
    z <- make_information_state(outer_sims[[m]]$theta)
    policy_matrix[p, m] <- policy_value(outer_sims[[m]]$B, z, policy)
  }
}

worst_policy_value <- apply(policy_matrix, 1, min)
best_policy_id <- which.max(worst_policy_value)
VZ_LE <- worst_policy_value[best_policy_id]
VOI_Z_LE <- VZ_LE - V0_LE

data.frame(best_policy = best_policy_id,
           VZ_LE = VZ_LE,
           VOI_Z_LE = VOI_Z_LE,
           all_policies[best_policy_id, , drop = FALSE])
```

## 5. Branch C: fixed-measure VOI envelope

This branch computes a conventional VOI value under each admissible precise probability measure and takes the range.

```{r}
evpi_one_measure <- function(B) {
  mean(apply(B, 1, max)) - max(colMeans(B))
}

T_evpi <- sapply(outer_sims, function(o) evpi_one_measure(o$B))
envelope_evpi <- range(T_evpi)
reference_index <- 1
reference_evpi <- T_evpi[reference_index]

data.frame(reference_EVPI = reference_evpi,
           lower_endpoint = envelope_evpi[1],
           upper_endpoint = envelope_evpi[2],
           missed_value = max(0, envelope_evpi[2] - reference_evpi),
           nonguaranteed_value = max(0, reference_evpi - envelope_evpi[1]))
```

## 6. Branch D: p-box expectation bounds

If a scalar output has support `[ell, u]` and lower/upper CDF bounds are available on a grid, lower and upper expectations for the full p-box band can be approximated by numerical integration. These are outer bounds for the original admissible set unless the full p-box band is the intended credal set.

```{r}
pbox_expectation_bounds <- function(x_grid, F_lower, F_upper) {
  stopifnot(length(x_grid) == length(F_lower), length(x_grid) == length(F_upper))
  ord <- order(x_grid)
  x <- x_grid[ord]
  FL <- pmin(pmax(F_lower[ord], 0), 1)
  FU <- pmin(pmax(F_upper[ord], 0), 1)
  dx <- diff(x)
  lower_integrand <- 1 - FU
  upper_integrand <- 1 - FL
  lower_E <- min(x) + sum(dx * (head(lower_integrand, -1) + tail(lower_integrand, -1)) / 2)
  upper_E <- min(x) + sum(dx * (head(upper_integrand, -1) + tail(upper_integrand, -1)) / 2)
  c(lower = lower_E, upper = upper_E)
}

x_grid <- seq(-1000, 1000, length.out = 501)
# CDF height is lower for the larger-mean normal distribution, so F_lower uses mean = 100.
F_lower <- pnorm(x_grid, mean = 100, sd = 300)
F_upper <- pnorm(x_grid, mean = -100, sd = 300)
pbox_expectation_bounds(x_grid, F_lower, F_upper)
```

## 7. Stability step: threshold classification

```{r}
classify_threshold <- function(lower, upper, threshold) {
  if (lower > threshold) return("stable positive")
  if (upper <= threshold) return("stable negative")
  "assumption-dependent"
}

threshold <- 100
classify_threshold(envelope_evpi[1], envelope_evpi[2], threshold)
```

## 8. Practical output table

```{r}
results <- data.frame(
  target = c("Fixed-measure EVPI envelope", "Lower-expectation EVPI"),
  estimate = c(NA, EVPI_LE),
  lower = c(envelope_evpi[1], NA),
  upper = c(envelope_evpi[2], NA),
  reference = c(reference_evpi, NA),
  threshold_classification = c(classify_threshold(envelope_evpi[1], envelope_evpi[2], threshold), NA)
)
results
```