From ce6ec127403b48ecbd0c0b9e8d50006f0b954ea9 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sat, 30 Oct 2021 00:37:56 +0300 Subject: [PATCH] v0.1 --- .gitignore | 20 +++++++ LICENSE | 21 ++++++++ Project.toml | 21 ++++++++ README0.md | 57 ++++++++++++++++++++ change.log | 3 ++ src/SeqBounds.jl | 11 ++++ src/functions.jl | 135 +++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 7 +++ 8 files changed, 275 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 Project.toml create mode 100644 README0.md create mode 100644 change.log create mode 100644 src/SeqBounds.jl create mode 100644 src/functions.jl create mode 100644 test/runtests.jl diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..59b8824 --- /dev/null +++ b/.gitignore @@ -0,0 +1,20 @@ +# Build and Release Folders +bin-debug/ +bin-release/ +[Oo]bj/ +[Bb]in/ + +# Other files and folders +.settings/ +./docs/src/pdf +./docs/src/pdf/* + +# Executables +*.swf +*.air +*.ipa +*.apk + +# Project files, i.e. `.project`, `.actionScriptProperties` and `.flexProperties` +# should NOT be excluded as they contain compiler settings and other important +# information for Eclipse / Flash Builder. diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..9fc4096 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 PharmCat + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..f9e4b71 --- /dev/null +++ b/Project.toml @@ -0,0 +1,21 @@ +name = "SeqBounds" +uuid = "343dbc49-a00e-4db9-80e1-8930f4861f86" +authors = ["PharmCat "] +version = "0.1.0" + +[deps] +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" + +[compat] +PrettyTables = "1" +Distributions = "0.20, 0.21, 0.22, 0.23, 0.24" +Roots = "1" +julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README0.md b/README0.md new file mode 100644 index 0000000..38c3421 --- /dev/null +++ b/README0.md @@ -0,0 +1,57 @@ +# SeqBounds + +Group sequential design bounds. + +## Install + +``` +import Pkg; Pkg.add("SeqBounds") +``` + +## Using + +``` +using SeqBounds + +bounds([0.25, 0.5, 0.75, 1.0], 0.05; h = 0.05) +``` + +Result: + +``` +julia> bounds([0.25, 0.5, 0.75, 1.0], 0.05; h = 0.05) +One-sided group sequential design +Alpha spending function: O'Brien-Fleming, Alpha = 0.05 +┌─────────┬────────────────┬────────────┬─────────┬────────────┐ +│ Portion │ Function value │ Spend │ Z │ Nominal p │ +├─────────┼────────────────┼────────────┼─────────┼────────────┤ +│ 0.25 │ 8.85754e-5 │ 8.85754e-5 │ 3.74955 │ 8.85754e-5 │ +│ 0.5 │ 0.0055746 │ 0.00548602 │ 2.53993 │ 0.00554366 │ +│ 0.75 │ 0.0236251 │ 0.0180505 │ 2.01604 │ 0.0218979 │ +│ 1.0 │ 0.05 │ 0.0263749 │ 1.72014 │ 0.0427037 │ +└─────────┴────────────────┴────────────┴─────────┴────────────┘ +``` + +## API +``` + bounds(v::Vector, alpha::Float64; h::Float64 = 0.05) +``` + +Where: + +* `v` - vector of information portion for each interim analysis; +* `alpha` - total alpha level; +* `h` - grid step multiplier, default 0.05, use 0.025 for better precision. + +Now: + +* Only O'Brien-Fleming alpha spending function implemented. +* Only one-sided bounds implemented. + +## Reference + +* Reboussin, D. M., DeMets, D. L., Kim, K., & Lan, K. K. G. (2000). Computations for Group Sequential Boundaries Using the Lan-DeMets Spending Function Method. Controlled Clinical Trials, 21(3), 190–207. doi:10.1016/s0197-2456(00)00057-x +* Lan KKG, DeMets DL. Discrete sequential boundaries for clinical trials. Biometrika. 1983; 70:659-63. +* DeMets DL, Lan G. “The alpha spending function approach to interim data analyses” in Recent Advances in Clinical Trial Design and Analysis, ed. PF Thall. Kluwer Academic Publishers, Boston; 1995. +* Armitage P, McPherson CK, Rowe BC. Repeated significance tests on accumulating data. Journal of the Royal Statistical Society. 1969; 132:235-44 +* ldbounds: Lan-DeMets Method for Group Sequential Boundaries - https://CRAN.R-project.org/package=ldbounds (comparation) diff --git a/change.log b/change.log new file mode 100644 index 0000000..7ad4c53 --- /dev/null +++ b/change.log @@ -0,0 +1,3 @@ + +0.1.0 + - Initial diff --git a/src/SeqBounds.jl b/src/SeqBounds.jl new file mode 100644 index 0000000..b3be85d --- /dev/null +++ b/src/SeqBounds.jl @@ -0,0 +1,11 @@ +module SeqBounds + + using Distributions, Roots, PrettyTables + + import Base:show + + export bounds + + include("functions.jl") + +end # module diff --git a/src/functions.jl b/src/functions.jl new file mode 100644 index 0000000..da83506 --- /dev/null +++ b/src/functions.jl @@ -0,0 +1,135 @@ + +#SeqBounds.jl +struct SeqBoundsResult + sf::String + v::Vector + alpha::Float64 + zb::Vector + asfv::Vector + asfvdiff::Vector + stdv::Vector + sdproc::Vector +end + +function obf(p, alpha) + 2(1 - cdf(Normal(), quantile(Normal(), 1-alpha/2)/sqrt(p))) +end + +function bounds(v::Vector, alpha::Float64; h::Float64 = 0.05) + if alpha ≥ 1 || alpha ≤ 0 error("alpha should be in range: 0 < alpha < 1!") end + if v[end] > 1 error("Last value of v shoul be ≤ 1!") end + if length(v) > 1 + for i = 2:length(v) + if !(v[i] > v[i-1]) error("v[i] shoul be > v[i-1]") end + end + end + zninf = -8.0 + side = 1 + asfv = obf.(v, alpha) # Alpha spending function values + asfvdiff = similar(asfv) + vdiff = similar(v) + asfvdiff[1] = asfv[1] + vdiff[1] = v[1] + for i = 2:length(v) + vdiff[i] = v[i] - v[i-1] + asfvdiff[i] = asfv[i] - asfv[i-1] + end + stdv = sqrt.(vdiff) + sdproc = sqrt.(v) + za = similar(asfv) + zb = similar(asfv) + ya = similar(asfv) + yb = similar(asfv) + nints = zeros(Int, length(v)) + grids = Vector{StepRangeLen}(undef, length(v)) + # + za[1] = zninf # For one side + zb[1] = quantile(Normal(), 1-asfvdiff[1]/side) + ya[1] = za[1]*stdv[1] + yb[1] = zb[1]*stdv[1] + + nints[1] = ceil((yb[1] - ya[1])/(h*stdv[1])) + grids[1] = range(ya[1], yb[1], length=nints[1] + 1) + last = pdf.(Normal(0, stdv[1]), grids[1]) + + for i = 2:length(v) + if i > 2 + # For next step + nints[i-1] = ceil((yb[i-1]-ya[i-1])/(h*stdv[i-1])) + grids[i-1] = range(ya[i-1], yb[i-1], length=nints[i-1]+1) + last = highordgrid(last, grids[i-1], grids[i-2], stdv[i-1]) + end + + za[i] = zninf + ya[i] = zninf*sdproc[i] + zerof = x -> asfvdiff[i] - integrategrid(x, last, grids[i-1], stdv[i]) + + yb[i] = find_zero(zerof, (0.0, zb[i-1])) + zb[i] = yb[i]/sdproc[i] + end + SeqBoundsResult("O'Brien-Fleming", + v, + alpha, + zb, + asfv, + asfvdiff, + stdv, + sdproc) +end + +function integrategrid(x, last::Vector{T}, grid, stdv) where T + hlast = step(grid) + DIST = Normal(x, stdv) + fst = cdf(DIST, grid[1]) * last[1] + lst = cdf(DIST, grid[end]) * last[end] + sfn = zero(T) + for i = 2:length(grid)-1 + sfn += cdf(DIST, grid[i]) * last[i] + end + hlast*(sfn + fst/2 + lst/2) +end + +function highordgrid(last::Vector{T}, x, grid, stdv) where T + h = step(grid) + rvec = zeros(T, length(x)) + for n = 1:length(x) + fst = last[1] * pdf(Normal(x[n], stdv), grid[1]) + lst = last[end] * pdf(Normal(x[n], stdv), grid[end]) + for m = 2:length(grid)-1 + rvec[n] += last[m] * pdf(Normal(x[n], stdv), grid[m]) + end + rvec[n] += fst/2 + lst/2 + end + rvec * h +end + +function Base.show(io::IO, obj::SeqBoundsResult) + println(io, """One-sided group sequential design + Alpha spending function: $(obj.sf), Alpha = $(obj.alpha)""") + pretty_table(io, hcat(obj.v, obj.asfv, obj.asfvdiff, obj.zb, ccdf.(Normal(), obj.zb)); header = ["Portion", "Function value", "Spend", "Z", "Nominal p"]) +end +################################################################################ +#= +Dounds can be found by hcubature with HCubature.jl but it takes more time and memory +=# +#= +f1i(x) = quadgk(ab -> pdf(Normal(0, stdv[1]), ab), ya[1], yb[1], rtol=1e-8)[1] +f1i(yb[1]) +f1(x) = pdf(Normal(0, stdv[1]), x) + +f21(x) = quadgk(u -> cdf(Normal(x, stdv[2]), u)*pdf(Normal(0, stdv[1]), u), ya[1], yb[1], rtol=1e-8)[1] +f21(yb[2]) + +# OR + +f22(x) = quadgk(u -> cdf(Normal(x, stdv[2]), u)*f1(u), ya[1], yb[1], rtol=1e-8)[1] +zerof = x ->asfvdiff[2] - f23(x) +fz = find_zero(zerof, (0.0, zb[1])) +zbf = fz/sdproc[2] + +f3(x) = hcubature(u -> cdf(Normal(x, stdv[i]), u[1]) * pdf(Normal(u[1], stdv[i-1]), u[2]) * pdf(Normal(0, stdv[1]), u[2]) , [ya[i-1], ya[i-2]], [yb[i-1],yb[i-2]], rtol=1e-8)[1] +f3(yb[3]) +zerof = x ->asfvdiff[3] - f3(x) +fz = find_zero(zerof, (0.0, zb[2])) +zbf = fz/sdproc[3] +=# diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..2ae1c88 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,7 @@ +using SeqBounds +using Test + +@testset "SeqBounds.jl" begin + + +end