multigrain/README.Rmd at main · GSK-Biostatistics/multigrain · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.path = "man/figures/README-",
    out.width = "100%"
)
```


# multigrain

<!-- badges: start -->
[![R-CMD-check](https://github.com/GSK-Biostatistics/multigrain/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/GSK-Biostatistics/multigrain/actions/workflows/R-CMD-check.yaml)
[![covr2gh-coverage](/../covr2gh-storage/badges/main/coverage_badge.svg)](https://github.com/GSK-Biostatistics/multigrain/actions/workflows/covr2gh.yaml)
<!-- badges: end -->

## Overview

**multigrain** (**Mul**tiple **T**esting using **Gr**aphical **A**pproaches; **I**mprove power with **N**umerical Optimisation) is an R package that finds the optimal graphical multiple testing procedure for your clinical trial.

Confirmatory trials routinely test multiple hypotheses — across endpoints, doses, or subpopulations — and must strongly control the family-wise error rate (FWER). [Graphical approaches](https://doi.org/10.1002/sim.3495) are a popular framework for this, but choosing the hypothesis weights and transition matrix that make up the graph is typically done by hand.

**multigrain** automates this: given simulated p-values and a user-defined **trial success measure** (a function that scores each possible pattern of hypothesis rejections), it searches over hypothesis weights and transition weights to find the graph that maximises expected trial success. The result is a valid graphical test that strongly controls FWER.

### Key features

| Step | Function | Purpose |
|------|----------|---------|
| Simulate p-values | `simulate_pvalues()` | Draw p-values from a multivariate normal test-statistic model |
| Define success | `trial_success()` | Specify what "trial success" means (compiled to C++ for speed) |
| Constrain the graph | `graph_constraint()` | Fix weights, edges, or testing hierarchies |
| Optimise | `graph_optimise()` | Find the graph that maximises expected trial success |
| Evaluate | `calc_power_pvals()` | Compute local power, disjunctive/conjunctive power, and custom metrics |

## Installation

**multigrain** is not on CRAN yet. Install the development version from GitHub with:

``` r
# install.packages("pak")
pak::pak("GSK-Biostatistics/multigrain")
```

## Quick start

Below is an end-to-end example with four hypotheses: two primary (H1, H2) and two secondary (H3, H4). Secondaries are gated by their corresponding primaries (H3 can only be tested after H1 is rejected; H4 after H2).

### 1. Simulate p-values

Assume one-sided test statistics are jointly normal. Provide the per-hypothesis nominal power (unadjusted for multiplicity) and a correlation matrix:

``` r
library(multigrain)

power_nominal <- c(0.95, 0.90, 0.85, 0.80)
corr_mat <- matrix(
    c(
        1.0, 0.5, 0.4, 0.2,
        0.5, 1.0, 0.2, 0.4,
        0.4, 0.2, 1.0, 0.5,
        0.2, 0.4, 0.5, 1.0
    ),
    nrow = 4, byrow = TRUE
)

set.seed(1)
pvals <- simulate_pvalues(
    power_nominal = power_nominal,
    corr_matrix = corr_mat,
    nsim = 1e5
)
```

### 2. Define trial success measures

Use `trial_success()` to express what "success" means for the trial. The variables `r1, r2, r3, r4` are binary indicators — `ri = 1` if hypothesis Hi is rejected, 0 otherwise. You can combine them with arithmetic (`+`, `*`) and logical (`&&`, `||`) operators.

``` r
# Average power: mean proportion of hypotheses rejected
avg_power <- trial_success(0.25 * (r1 + r2 + r3 + r4))

# Disjunctive power: at least one rejection
disjunctive <- trial_success(r1 || r2 || r3 || r4)

# Custom: require a primary rejection, with extra credit for the
# gated secondary (H3 adds value only if H1 is also rejected, etc.)
custom_success <- trial_success(
    0.25 * (2 * (r1 || r2) + r1 * r3 + r2 * r4)
)
```

### 3. Constrain the graph (optional)

Use `graph_constraint()` to encode structural rules — for example, that only primaries receive initial alpha, and that secondaries can only receive alpha recycled from their gating primary:

``` r
gc <- graph_constraint(
    hyp_constraint  = c(NA, NA, 0, 0),
    trans_constraint = matrix(
        c(
            0, NA, NA, 0,
            NA, 0, 0, NA,
            0, 1, 0, 0,
            1, 0, 0, 0
        ),
        nrow = 4, byrow = TRUE
    )
)
```

Here `NA` means "free to optimise" and a fixed number means "held constant".

### 4. Optimise

``` r
g_opt <- graph_optimise(
    pvals = pvals,
    trial_success = custom_success,
    graph_constraint = gc
)

summary(g_opt)
plot(g_opt)
```

### 5. Evaluate

Compute multiplicity-adjusted power and any custom trial success measures for the optimised graph:

``` r
calc_power_pvals(
    pvals,
    hyp_weight   = g_opt$hyp_weight,
    trans_matrix = g_opt$trans_matrix,
    custom_power = list(
        AvgPower      = avg_power,
        Disjunctive   = disjunctive,
        CustomSuccess = custom_success
  )
)
```

## Documentation

- [Getting Started](https://gsk-biostatistics.github.io/multigrain/articles/get-started.html) — full walkthrough of the asthma example above
- [Graph Constraints](https://gsk-biostatistics.github.io/multigrain/articles/graph_constraint.html) — specifying and visualising constraints
- [Function Reference](https://gsk-biostatistics.github.io/multigrain/reference/index.html)