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
|
## Simulate Bayesian single-arm adaptive trial
## Allow early termination due to futility or efficacy
## Binary outcome
## Beta-binomial:
## p ~ beta(a, b)
## x_i ~ binomial(p) i = 1..n
## p|x ~ beta(a + sum(x), b + n - sum(x))
## Efficacy at interim t if Pr(p > p_0 | x_{(t)}) > gamma_e
## Futility at interim t if Pr(p > p_0 | x_{(t_max)}) < gamma_f
## https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4247348/
library('rmutil') ## for betabinom
## Simulate entire trial
## ptru - true probability of outcome (p)
## pref - reference probability of outcome (p_0)
## nint - sample sizes at which to conduct interim analyses
## efft - efficacy threshold
## futt - futility threshold
## apri - prior beta parameter alpha
## bpri - prior beta parameter beta
simtrial <- function(
ptru = 0.15,
pref = 0.15,
nint = c(10, 13, 16, 19),
efft = 0.95,
futt = 0.05,
apri = 1,
bpri = 1) {
## determine minimum number of 'successes' necessary to
## conclude efficacy if study continues to maximum
## sample size
nmax <- max(nint)
post <- sapply(0:nmax, function(nevt)
1-pbeta(pref, apri + nevt, bpri + nmax - nevt))
nsuc <- min(which(post > efft)-1)
## simulate samples
samp <- rbinom(n = nmax, size = 1, prob = ptru)
## simulate interim analyses
intr <- lapply(nint, function(ncur) {
## compute number of current events
ecur <- sum(samp[1:ncur])
## compute posterior beta parameters
abb <- apri + ecur
bbb <- bpri + ncur - ecur
sbb <- abb + bbb
mbb <- abb/(abb+bbb)
## compute efficacy Pr(p > p_0 | x_{(t)})
effp <- 1-pbeta(pref, abb, bbb)
## return for efficacy
if(effp > efft)
return(list(action='stop',
reason='efficacy',
n = ncur))
## number of events necessary in remainder of
## study to conclude efficacy
erem <- nsuc-ecur
## compute success probability Pr(p > p_0 | x_{(t_max)})
if(erem > nmax-ncur) { ## not enough possible events
sucp <- 0
} else { ## not yet met efficacy threshold
sucp <- 1-pbetabinom(q = erem-1,
size = nmax-ncur, m = mbb, s = sbb)
}
if(sucp < futt)
return(list(action='stop',
reason='futility',
n = ncur))
return(list(action='continue',
reason='',
n = ncur))
})
stpi <- match('stop', sapply(intr, [[, 'action'))
return(intr[[stpi]])
}
## Simulate study with max sample size of 200 where true
## probability is identical to reference (i.e., the null
## hypothesis is true). This type of simulation helps us
## determine the overall type-I error rate.
nint <- c(40,80,120,160,200)
nmax <- max(nint)
res <- do.call(rbind, lapply(1:10000,
function(t) as.data.frame(simtrial(ptru = 0.72,
pref = 0.72,
nint = nint,
efft = 0.975,
futt = 0.20))))
## Prob. early termination (PET) due to Futility
mean(res$reason == 'futility' & res$n < nmax)
## PET Efficacy
mean(res$reason == 'efficacy' & res$n < nmax)
## Pr(conclude efficacy) 'type-I error rate'
mean(res$reason == 'efficacy')
## average and sd sample size
mean(res$n); sd(res$n)
barplot(prop.table(table(res$n)),
xlab='Study Size (N)',
main="No Difference")
## Simulate study where true probability is greater than
## reference (i.e., an alternative hypothesis). This type
## of simulation helps us determine the study power.
res <- do.call(rbind, lapply(1:10000,
function(t) as.data.frame(simtrial(ptru = 0.82,
pref = 0.72,
nint = nint,
efft = 0.975,
futt = 0.20))))
## Prob. early termination (PET) due to Futility
mean(res$reason == 'futility' & res$n < nmax)
## PET Efficacy
mean(res$reason == 'efficacy' & res$n < nmax)
## Pr(conclude efficacy) 'power'
mean(res$reason == 'efficacy')
## average and sd sample size
mean(res$n); sd(res$n)
barplot(prop.table(table(res$n)),
xlab='Study Size (N)',
main="35% Reduction")
|
近期评论