/
minimize.R
345 lines (302 loc) · 10.8 KB
/
minimize.R
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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
#' Find optimal two-stage design by constraint minimization
#'
#' \code{minimize} takes an unconditional score and
#' a constraint set (or no constraint) and solves the corresponding
#' minimization problem using
#' \href{https://cran.r-project.org/package=nloptr}{\code{nloptr}}
#' (using COBYLA by default).
#' An initial design has to be defined. It is also possible to define
#' lower- and upper-boundary designs. If this is not done, the boundaries are
#' determined automatically heuristically.
#'
#' @param objective objective function
#' @param subject_to constraint collection
#' @param initial_design initial guess (x0 for nloptr)
#' @param lower_boundary_design design specifying the lower boundary.
#' @param upper_boundary_design design specifying the upper boundary
#' @param opts options list passed to nloptr
#' @param ... further optional arguments passed to \code{\link{nloptr}}
#'
#' @return a list with elements:
#' \item{design}{ The resulting optimal design}
#' \item{nloptr_return}{ Output of the corresponding nloptr call}
#' \item{call_args}{ The arguments given to the optimization call}
#'
#' @examples
#' # Define Type one error rate
#' toer <- Power(Normal(), PointMassPrior(0.0, 1))
#'
#' # Define Power at delta = 0.4
#' pow <- Power(Normal(), PointMassPrior(0.4, 1))
#'
#' # Define expected sample size at delta = 0.4
#' ess <- ExpectedSampleSize(Normal(), PointMassPrior(0.4, 1))
#'
#' # Compute design minimizing ess subject to power and toer constraints
#' \dontrun{
#' minimize(
#'
#' ess,
#'
#' subject_to(
#' toer <= 0.025,
#' pow >= 0.9
#' ),
#'
#' initial_design = TwoStageDesign(50, .0, 2.0, 60.0, 2.0, 5L)
#'
#' )
#' }
#'
#' @export
minimize <- function(
objective,
subject_to,
initial_design,
lower_boundary_design = get_lower_boundary_design(initial_design),
upper_boundary_design = get_upper_boundary_design(initial_design),
opts = list(
algorithm = "NLOPT_LN_COBYLA",
xtol_rel = 1e-5,
maxeval = 10000
),
...
) {
args <- c(as.list(environment()), list(...))
f_obj <- function(params) {
evaluate(
objective,
update(initial_design, params),
optimization = TRUE # evaluate in optimization context!
)
}
g_cnstr <- function(params) {
design <- update(initial_design, params)
user_cnstr <- evaluate(subject_to, design, optimization = TRUE)
return(c(
user_cnstr,
design@c1f - design@c1e + ifelse( # ensure c1e > c1f if not one-stage
is(initial_design, "OneStageDesign"), 0, .1)
))
}
res <- nloptr::nloptr(
x0 = tunable_parameters(initial_design),
lb = tunable_parameters(lower_boundary_design),
ub = tunable_parameters(upper_boundary_design),
eval_f = f_obj,
eval_g_ineq = g_cnstr,
opts = opts,
...
)
if (res$status == 5 | res$status == 6)
warning(res$message)
res <- list(
design = update(initial_design, res$solution),
nloptr_return = res,
call_args = args
)
class(res) <- c("adoptrOptimizationResult", class(res))
return(res)
}
#' @rawNamespace S3method(print, adoptrOptimizationResult)
print.adoptrOptimizationResult <- function(x, ...) {
cat(design2str(x$design, TRUE), "\n")
}
#' Initial design
#'
#' The optimization method \code{\link{minimize}} requires an initial
#' design for optimization.
#' The function \code{get_initial_design} provides an initial guess based on a
#' fixed design that fulfills constraints on type I error rate and power.
#' Note that a situation-specific initial design may be much more efficient.
#'
#' @param theta the alternative effect size in the normal case, the
#' rate difference under the alternative in the binomial case
#' @param alpha maximal type I error rate
#' @param beta maximale type II error rate
#' @param type is a two-stage, group-sequential, or one-stage design requried?
#' @param dist distribution of the test statistic
#' @param order desired integration order
#' @template dotdotdot
#'
#' @details
#' The distribution of the test statistic is specified by \code{dist}.
#' The default assumes a two-armed z-test.
#'
#' @examples
#' init <- get_initial_design(
#' theta = 0.3,
#' alpha = 0.025,
#' beta = 0.2,
#' type = "two-stage",
#' dist = Normal(two_armed = FALSE),
#' order = 7L
#' )
#'
#' @export
get_initial_design <- function(theta, alpha, beta,
type = c("two-stage", "group-sequential", "one-stage"),
dist = Normal(), order = 7L, ...) {
type <- match.arg(type)
if (alpha <= 0 || alpha >= 1 || beta <= 0 || beta >= 1)
stop("alpha and beta must be in (0, 1)!")
if (is(dist, "Binomial")) {
p_0 <- dist@rate_control + theta / 2
theta <- theta / sqrt(p_0 * (1 - p_0))
}
c <- quantile(dist, 1 - alpha, 1, 0)
n <- floor(2 * (c + quantile(dist, 1 - beta, 1, 0))^2 / theta^2)
if (type == "one-stage")
return(OneStageDesign(n, c))
else if (type == "group-sequential")
return(GroupSequentialDesign(n/2, 0, quantile(dist, 1 - alpha/2, 1, 0), n/2, c, order))
else if (type == "two-stage")
return(TwoStageDesign(n/2, 0, quantile(dist, 1 - alpha/2, 1, 0), n/2, c, order))
}
#' Boundary designs
#'
#' The optimization method \code{\link{minimize}} is based on the package
#' \code{nloptr}. This requires upper and lower boundaries for optimization.
#' Such boundaries can be computed via \code{lower_boundary_design}
#' respectively \code{upper_boundary_design}.
#' They are implemented by default in \code{\link{minimize}}.
#' Note that \code{\link{minimize}} allows the user to define its own
#' boundary designs, too.
#'
#' @param initial_design The initial design
#' @param n1 bound for the first-stage sample size n1
#' @param n2_pivots bound for the second-stage sample size n2
#' @param c1_buffer shift of the early-stopping boundaries from the initial ones
#' @param c2_buffer shift of the final decision boundary from the initial one
#' @param ... optional arguments
#'
#' The values \code{c1f} and \code{c1e} from the initial design are shifted
#' to \code{c1f - c1_buffer} and \code{c1e - c1_buffer} in
#' \code{get_lower_boundary_design}, respectively, to \cr
#' \code{c1f + c1_buffer} and \code{c1e + c1_buffer} in
#' \code{get_upper_boundary_design}.
#' This is handled analogously with \code{c2_pivots} and \code{c2_buffer}.
#'
#' @examples
#' initial_design <- TwoStageDesign(
#' n1 = 25,
#' c1f = 0,
#' c1e = 2.5,
#' n2 = 50,
#' c2 = 1.96,
#' order = 7L
#' )
#' get_lower_boundary_design(initial_design)
#'
#' @rdname boundary-designs
#' @export
setGeneric("get_lower_boundary_design",
function(initial_design, ...) standardGeneric("get_lower_boundary_design"))
#' @rdname boundary-designs
#' @export
setGeneric("get_upper_boundary_design",
function(initial_design, ...) standardGeneric("get_upper_boundary_design"))
#' @rdname boundary-designs
#' @export
setMethod("get_lower_boundary_design", signature("OneStageDesign"),
function(initial_design, n1 = 1, c1_buffer = 2, ...) {
lb_design <- OneStageDesign(n1, min(0, initial_design@c1f - c1_buffer))
lb_design@tunable <- initial_design@tunable
return(lb_design)
})
#' @rdname boundary-designs
#' @export
setMethod("get_lower_boundary_design", signature("GroupSequentialDesign"),
function(
initial_design,
n1 = 1,
n2_pivots = 1,
c1_buffer = 2,
c2_buffer = 2,
...
) {
lb_design <- GroupSequentialDesign(
n1,
initial_design@c1f - c1_buffer,
initial_design@c1e - c1_buffer,
n2_pivots,
initial_design@c2_pivots - c2_buffer,
order = length(initial_design@c2_pivots)
)
lb_design@tunable <- initial_design@tunable
return(lb_design)
})
#' @rdname boundary-designs
#' @export
setMethod("get_lower_boundary_design", signature("TwoStageDesign"),
function(
initial_design,
n1 = 1,
n2_pivots = 1,
c1_buffer = 2,
c2_buffer = 2,
...
) {
lb_design <- TwoStageDesign(
n1,
initial_design@c1f - c1_buffer,
initial_design@c1e - c1_buffer,
rep(n2_pivots, length(initial_design@c2_pivots)),
initial_design@c2_pivots - c2_buffer,
order = length(initial_design@c2_pivots)
)
lb_design@tunable <- initial_design@tunable
return(lb_design)
})
#' @rdname boundary-designs
#' @export
setMethod("get_upper_boundary_design", signature("OneStageDesign"),
function(initial_design, n1 = 5 * initial_design@n1, c1_buffer = 2, ...) {
ub_design <- OneStageDesign(n1, initial_design@c1f + c1_buffer)
ub_design@tunable <- initial_design@tunable
return(ub_design)
})
#' @rdname boundary-designs
#' @export
setMethod("get_upper_boundary_design", signature("GroupSequentialDesign"),
function(
initial_design,
n1 = 5 * initial_design@n1,
n2_pivots = 5 * initial_design@n2_pivots,
c1_buffer = 2,
c2_buffer = 2,
...
) {
ub_design <- GroupSequentialDesign(
n1,
initial_design@c1f + c1_buffer,
initial_design@c1e + c1_buffer,
n2_pivots,
initial_design@c2_pivots + c2_buffer,
order = length(initial_design@c2_pivots)
)
ub_design@tunable <- initial_design@tunable
return(ub_design)
})
#' @rdname boundary-designs
#' @export
setMethod("get_upper_boundary_design", signature("TwoStageDesign"),
function(
initial_design,
n1 = 5 * initial_design@n1,
n2_pivots = 5 * initial_design@n2_pivots,
c1_buffer = 2,
c2_buffer = 2,
...
) {
ub_design <- TwoStageDesign(
n1,
initial_design@c1f + c1_buffer,
initial_design@c1e + c1_buffer,
n2_pivots,
initial_design@c2_pivots + c2_buffer,
order = length(initial_design@c2_pivots)
)
ub_design@tunable <- initial_design@tunable
return(ub_design)
})