kkmann/adoptr

Adaptive optimal two-stage designs for clinical trials made easy.

adoptr's Issues

two-arm vs single arm

haha, that's a fun one - I think we currently work with a single arm z test which is kinda boring. Extend this via the data distribution interface to two-stages (equal allocation for 0.1, everything else can come later ;))!

finalize package vignette

We need to decide what should go in the final package vignette.

I'd say let's just do the standard case + one design showing off Bayesian prior + conditional power constraint.

We also need to talk about prior capabilities + score/constraint arithmetic/construction.

Could do all in one file (slightly in favor of that) or split prior and score/constraints to different files (currently just dummies)

Regularize n1

with the new AverageN2 regularizer, we do not penalize n1 at all. I am not saying that this is necessary, but maybe one would like to apply the same amount of regularization to n2 and n1 to avoid changing the shape of a design too much. This can be easily done by adding a regularizer just penalizing n1. Non, critical, I am not assigning a milestone but we might want to keep this in the backlog.

plot n() properly

the n() function should be plotted as non-continuous on the stopping boundaries, cf.

https://github.com/kkmann/otsd/blob/7dc292730e9c1ced854c150c8a8ba67428f80002/R/TwoStageDesign.R#L281-L303

revise tests to postprocess = FALSE

since we change the post processing default to TRUE in #98, we might want to manually switch to non-post-processing in most test cases (except those checking post processing ;))

postprocessing of solutions

we need to round the n1/n2 part after optimization and then refit the remaining parameters. This is a heuristic, but probably the best one can do...

design proper namespace

We currently export everything. That's not really necessary. We should go through each item and decide whether or not this actually needs to be exported (e.g., users should not need to call AffineConditionalScore directly - ever). This also relieves some documentation burden. Important but not 0.1.0-important.

proper status information in minimize()

we need to make sure that we return the entire information from the optimization call

https://github.com/kkmann/otsd/blob/7dc292730e9c1ced854c150c8a8ba67428f80002/R/minimize.R#L104

not just the optimal desing

Choose only some parameters optimally

I think we should make it possible to fix $n_1$ or even the complete first stage, as it was also mentioned in the paper's review.
Is this straight-forward or do we even need a new class for this?

Increase test coverage and quality

Well, I'd feel a lot more confident if we had something like 90+ % coverage. Test also need to be revised and less copy / paste!

postprocessing does not work for one-stage desings (also check for GS designs)

check for duplicate scores

We have to fix the ToDo in AffineScoreFunction regarding duplicated scores.

https://github.com/kkmann/otsd/blob/7b9218b85f325a0bd84cec83d6de01c6aa2a48e1/R/AffineScoreFunction.R#L77

implement test vignette for scenario I of our valdation concept

This should be the template for all other scenarios. Work out how exactly we are describing everything, the structure and how we handle test cases (where to throw errors). Also need to check how we get this build and deployed on Travis-CI without running into time-limit problems!

Implement one-stage designs

Now that we have such a wonderfully general interface - why not allow optimizing simple one-stage designs as well? Should be a piece of cake - only two parameter n and c, right?

In practice, there will be a few locations in the code where I implicitly assumed TwoStageDesign - those might be a little bit of a headache but in principle should it should be straight-forward to do.

improve generic integration rule

The generic evaluation currently ignores uncontinuous jumps at continuation area vs rest. should also integrate in three parts (also faster!). @MatheMax That's an easy one...

https://github.com/kkmann/otsd/blob/83de4d296e7a3057c450fb70710e42ccc7b90d79/R/Scores.R#L83-L106

fix function names etc

We might need to rename some stuff or at least discuss it. Feel free to collect any weirdnesses here ;)

implement summary + show for our classes

We should have a generic way of getting summaries for design objects.

summary(design, unconditionalscore1, ..., unconditionalscorek)

should return a data.frame,

show(design)

needs to be implemented...

Come up with a better name

OTSD is as interesting and memorable as this ... other thing I recently heard about and cannot remember now.

I am scratching my head but cannot come up with anything better than BAdaSS for Bayesian Adaptive Sample Size - a controversial choice. I am not quite sure whether using 'BAdaSS designs' would be a pro or con in a grant application...

review documentation

We need to review the entire documentation - see where it can be simplified etc.

Vectorize functions

Check what functions are used by stats::integrate and vectorize them all!

determine final suggested integration order

We should come up with a final suggest default for the order of integration.

Criteria are 'nicety' of the resulting n2/c2 functions (smooth enough?) and run time since I suspect that our crappy optimizer will converge a lot faster with fewer parameters.

Goes hand-in-hand with the default relative tolerance and maximum number of iterations.

This is not a priority, correctness and feature completeness are more important but it should be investigated systematically.

generalize constraint generation syntax

we currently only allow score <= / >= value, we should also add value <= / >= score and score >= /<= score for sake of completeness. boring and non-crucial but should be done sometime.

implement Gaussian quadrature instead of Newton Cotes

cf, #2

revise minimize function

minimize() has become quite cluttered over time - any chance of simplifying that? Maybe we should have a separate method for one stage designs?

Distributions/Scores interface

@MatheMax We finally want to support scores of the form s(design, x, theta) and any form of expectation with respect to a joint distribution over (X, Theta). I.e., allowing to integrate over either X, Theta or both.
Currently a ConditionalScore is only allowed to depend on x and we do not really model joint distributions. This needs to be generalized for full flexibility.

A possible interface could be

# specify joint distribution of outcomes and parameters by specifying
# families for X | theta, n and theta
# i.e. in a normal model with uniform prior over theta (standardized effect size)
cp <- ConditionalPower(
    x ~ Normal(), 
    theta ~ ContinuousPrior(function(x) = 1, c(0, 1))
)

evaluate conditional on both:

eval(design, cp, x = 1.55, theta = .0) # conditional type one error given z = 1.55

or only conditional on observed data

eval(design, cp, x = 1.55) # integrate over theta | x

or just the expected value (probability of rejecting):

eval(design, cp)

or condition before integration:

eval(design, cp, theta = c(.3, 1), x = 1.55)

btw 'power' would then be

eval(design, cp, theta = .4)

Anyway, looks much more consistent than the mixture we have right now. Easyly possible, just a bit 'fiddly'.

Means that people only really need to care about s(design, x, theta) and all 'integral scores' are
truely derived consitently.

show() methods for our custom classes

We need to implement a reasonable show() method for all custom classes. Otherwise, the confusing default list-like method is used. In almost all cases, I'd just print the class name since more detailed output will be too complex/confusing on the terminal, @MatheMax what's your take on this?

get rid of (allmost) all ifelse!

okay, just found another ifelse() that somehow ended up in the code - think that's due to the rounded/not rounded distinction.

There is VERY rarely a reason to use ifelse() when using it you have to be aware of it's quirks!

e.g. here

https://github.com/kkmann/otsd/blob/c40797087fb0d99b64c355b9e070035673224df7/R/TwoStageDesign.R#L154-L159

ifelse leads to breaking the vectorization of n2() since it's output dimension depends on the FIRST argument onyl (-> docs).

@MatheMax please make sure u remove/double check all uses of ifelse. Btw, these problems only show that we invested to little time in testing properly ;)

implement plot()

We should provide sensible methods for plot(). Should show sample size and critical value vs x1 and optionally an arbitrary list of conditional scores

plot(design, conditionalpower, ...)

would be a 1x3 plot with n/n2/conditional power.

Only use base plot, no ggplot2, to avoid additional dependency.

Figure out how to implement SmoothnessPenalty

Hey,

Smoothness_n2 is an ugly mess at the moment. I hacked it together to avoid having to specify a prior etc. but it would be much nicer to fit it (at least approximately) in our system.
We simply need to specify a wide uniform default prior and it should work fine.
Also, do we want this for n2 and c2? seems to work fine for only n2.

simulate() function

Would be nice to have a simulate(design, prior, n) function returning respective z1/z2 values (can be NA), n2, c2, and final decision. This allows for easy cross-checking of power / toer.

allow detailed specification of conditional constraints

we currently simply apply condition scores to the continuation set. That is a reasonable default but it should be possible to pass something any subset of list("futility", "continuation", "efficacy") and apply it accordingly.

How/where do we best put that?

https://github.com/kkmann/otsd/blob/7b9218b85f325a0bd84cec83d6de01c6aa2a48e1/R/constraints.R#L90-L104

Is our .evaluate correct?

Simple question: is the current implementation good enough to handle effect size weighted sample size?
I think yes, but @MatheMax have you thought about this?

https://github.com/kkmann/otsd/blob/83de4d296e7a3057c450fb70710e42ccc7b90d79/R/GQDesign.R#L68-L71

eyecandy for negative scores

The syntax -score should work for both conditional and unconditional scores. We currently need to write -1 * score

examine alternative optimizers

Following #47, it might be worthwhile playing around with other optimizers or at least think about ways to implement that. Stochastic global optimizers might be an option but I am not sure as to how well they handle non-linearly constraint problems.

Support of continuous prior

Unfortunately, it turns out that the implemented integration by computing the posterior does not always work.

I think the problem is the following:
When computing the renormalizing constant, this constant can become very small.
This yields to quite large function values of the posterior density (on a small part of the support) and integrating the posterior does not necessarily return "1" anymore.

For now, I see 3 possible ways:

1. Do not change anything and look for a suitable support for any problem
2. Compute the posterior by using Gaussian quadrature. This should work but will make the minimization very slow.
3. Start creating subclasses of "Continuous prior". This could at least help in the case of a normal prior when we know the posterior distribution.

@kkmann Do you have any suggestions?

plot does not work with one-stage designs

well, probably should raise an error since there is not really anything to plot, right?

quality checks for gaussquad

• implements some tests to see how trustworthy gaussquad actually is (just integrating random stuff)
• implement some check comparing .evaluate with the non-generic variant to see if GQ works fine
• implement tests for different orders of integration (very preliminary now but we should at some point come up with a good pick for the order... more gives smoother approximation but there will be an upper limit that we can do without subdivisions...)

implement more examples to check plausibility

We need (all with toer on boundary of null hypothesis):

• point alternative:
  • large value (0.4)
    • ESS under null
    • ESS under 0.4
  • small value (0.2)
    • ESS under null
    • ESS under 0.2
• continuous normal prior
  • ESS under null
  • ESS under normal prior
• uniform prior on small interval ([.3, .5]?)
  • ESS under null
  • ESS under normal prior
• all of the above with expected power using the alternative
• all of the above (only expected power) with constraint on conditional power

Those should all eventually end up as test cases (ideally these would not be copy paste but
programmatically generated from matrix of hypotheses / conditions).

Constraints and final ESS should be checked usig simulate().

I would star, however, by picking 4-8 most important examples and implement them in a vignette. This vignette 'examples' automatically ends up in the online documentation (pkgdown).
Once we are happy with the results, these can be 'frozen' by using test_known_value again...

Implement group-sequential designs

Well, exactly that. Defining GS as fixed n2, same as before - just that n2 is parameterized by a single parameter.

Optimization framework

Since score arithmetics are preliminarily available, we could think about sketching up an optimization framework.

An interface could look like this:

problem <- minimize(objective,  ...)
optimal_design <- solve(problem, ...) # rather optimize but solve(problem) was just too much of a temptation ...

for this to work, we need a class MinimizationProblem (for which minimize is essentially a convenient constructor). objective can be any kind of unconditional score ... hold arbitrary number of constraints.

problem would then internally store everything needed to pass to nloptr and solve essentially calls nloptr and converts the result back to a design object (and accepts arbitrary optional arguments to pass to nloptr).

We also still need a class for conditional/unconditional constraints for this interface to work.
Thats just a thin wrapper around scores of the form

constraint(score, rhs)

which must also be evaluatable (essentially its just score - constant right?).
I will then implement a NSE version allowing constraint(score <= a) to be a valid expression which is definitely cool, right ;)

So ultimately you could have something like

minimize(expected_sample_size,
    constraint(power >= .8),
    constraint(toer <= .05),
    constraint(cp >= .7)
)

That should do.

Validation strategy

We need to come up with a good validation strategy.

Besides proper unittesting, this should also involve an automatically generated documents with several test cases which can be checked visually for correctness and a detailed explanation of the overall testing strategy. In particular, we would probably include a simple function allowing to sample from a design - we could then validate our numerical solutions by simulation. This is probably the best we can do. Before we start doing that we need to define objectives etc. I.e. to what level should this software be usable in a real clinical trial? I feel that it should be possible to use it, maybe with an intermediate validation step required. We could even go as far as providing a SAS script for validation...

Normal prior

Since for a normal prior the posterior is known analytically, we could implement a subclass of ContinuousPrior and implement the normal distribution case.

Create a hex-sticker

well, yeah.

test evaluate ~= .evaluate for some scores

I think we are currently lacking tests to ensure that the internal version and the version used during optimization return roughly the same values, right?

Assessing convergence

Okay, @MatheMax next major insight (cf. PR #46 ): Turns out that you get an almost perfect solution even for order 9, it just takes a while. I can get rid of the ugly bits by just running the optimizer longer with higher precision (stopped at 10k iterations).

A few random thoughts:

1. The results are much closer to the 'analytical' ones you get from solving the Euler Lagrange equations locally - the function is even almost perfectly concave!
2. That is good, extremely good even (our code is not so crappy after all and we are on the right track;))!
3. Our optimizer is extremely ineffective. That is not surprising and my initial feeling to implement this in julia instead was exactly due to concerns in this direction. The problem seems to be that the COBYLA algorithm (which is the only constraint, gradient-free optimizer in nloptr) is not only super inefficient (well, does not use gradients...) but also converges 'non-uniformly', i.e., it first concentrates on a few parameters before moving to the rest. The advantage in julia would be that one could use source-to-source automatic differentiation with gradient descent to (hopefully) get a more efficient but still generic optimizer.
4. We will stick with R since julia is not that known with the target audience (but keep in mind that this would be an excellent exercise demonstrating the power of julia ;)).
5. We need to accept MUCH longer run time (probably less costly with order 7).
6. We should note in the manual that the default run time might not be safe and that it should be increased when weird things happen ;)
7. We need to be careful when implementing test cases since we do not want the tests to run for more than say 30 minutes, at least not during development.
8. It might be interesting to see the progress during optimization. This can be realized relatively easily in R by using functionals. Couldn't find a good R reference but the idea is the same as here (Python) https://stackoverflow.com/questions/16739065/how-to-display-progress-of-scipy-optimize-function#18471343
9. Ideally, one would also allows early termination but that's optional (shiny app, much later ?)
10. A practical approach to making things usable again might be regularization. If COBYLA reacts well to stronger gradient information, we might get stable results faster without deviating too much from optimality. It is nice to know however, that you get proper convergence (even for relatively large order) if you just wait long enough. Regularization can thus be traded off against speed.

revise test suit

before we let this loose we should ensure that the coding style in our test suite is consistent, we use enough whitespaces and all tests are sufficiently well documented.

Smoothness term

Figure out if absolute value (l1-norm) works better for smoothness term

Resolve non-clean import of gaussquad dependency

@MatheMax as of 83de4d2 gaussquad is imported via the Depends: field in DESCRIPTION. This is unclean since it attaches the gaussquad namespace.

gaussquad functions seem to be unusable without attaching the namespace though - complains about legendre.recurrence not being exported.

This should probably be fixed upstream in the gausquad package - is the maintainer still active? Do we have alternatives?

Implement simpler regularizer

Following #47, we might get faster and more stable convergence by using suitable regularizers.

Regularizing c2() is extremely difficult - at least I have no idea about how to chose a penalty parameter there. The scale is just too abstract. Regularizing n() is much easier.
I would suggest trying to directly penalize the the integral over n2 on the continuation region (computed via Gauss quadrature) divided by the lenght of that interval. @MatheMax that should be easy to implement and could help with convergence without changing the shape of the solution as much as the smoothness term.

check correct vectorization when using expectation(dist,f)

We need to check that we always pass correctly vectorized functions since expectation assumes that!