v0.1

.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.

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.

Project.toml

@@ -0,0 +1,21 @@
+name = "SeqBounds"
+uuid = "343dbc49-a00e-4db9-80e1-8930f4861f86"
+authors = ["PharmCat <v.s.arnautov@yandex.ru>"]
+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"]

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)

change.log

@@ -0,0 +1,3 @@
+
+0.1.0
+  - Initial

src/SeqBounds.jl

@@ -0,0 +1,11 @@
+module SeqBounds
+
+    using Distributions, Roots, PrettyTables
+
+    import Base:show
+
+    export bounds
+
+    include("functions.jl")
+
+end # module

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]
+=#

test/runtests.jl

@@ -0,0 +1,7 @@
+using SeqBounds
+using Test
+
+@testset "SeqBounds.jl" begin
+
+
+end