Merge pull request #22 from PharmCat/dev

Dev

Project.toml

@@ -1,7 +1,7 @@
 authors = ["Vladimir Arnautov (mail@pharmcat.net)"]
 name = "ClinicalTrialUtilities"
 uuid = "535c2557-d7d0-564d-8ff9-4ae146c18cfe"
-version = "0.1.13"
+version = "0.1.14"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -13,6 +13,8 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 
+[compat]
+julia = "1"
 
 [extras]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

README.md

@@ -2,9 +2,6 @@
 
  Clinical trial related calculation: descriptive statistics, power and sample size calculation, power simulations, confidence interval, pharmacokinetics/pharmacodynamics parameters calculation.
 
-Version:0.1.13
-
-
 [![Build Status](https://travis-ci.com/PharmCat/ClinicalTrialUtilities.jl.svg?branch=master)](https://travis-ci.com/PharmCat/ClinicalTrialUtilities.jl)
 [![Build status](https://ci.appveyor.com/api/projects/status/35f8b5vq259sbssg?svg=true)](https://ci.appveyor.com/project/PharmCat/clinicaltrialutilities-jl)
 [![codecov](https://codecov.io/gh/PharmCat/ClinicalTrialUtilities.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/PharmCat/ClinicalTrialUtilities.jl)

cange.log

@@ -1,3 +1,9 @@
+v0.1.14
+  - Add compat to Project.toml
+  - Julian function names (part)
+  - descriptive_!() -> descriptive_()
+  - randomtable() included and add test
+
 v0.1.13
   - MN DIFF CI optimization
   - Index for ConfInt

src/CI.jl

@@ -20,9 +20,9 @@ module CI
             return ConfInt(0,0,x/n)
         end
         if method==:wilson
-            return propWilsonCI(x, n, alpha)
+            return propwilsonci(x, n, alpha)
         elseif method==:wilsoncc
-            return propWilsonCCCI(x, n, alpha)
+            return propwilsonccci(x, n, alpha)
         elseif method==:cp
             return propCPCI(x, n, alpha)
         elseif method==:soc
@@ -105,7 +105,7 @@ module CI
 
     #Wilson’s confidence interval for a single proportion, wilson score
     #Wilson, E.B. (1927) Probable inference, the law of succession, and statistical inferenceJ. Amer.Stat. Assoc22, 209–212
-    function propWilsonCI(x::Int, n::Int, alpha::Float64)::ConfInt
+    function propwilsonci(x::Int, n::Int, alpha::Float64)::ConfInt
         z = abs(quantile(ZDIST, 1-alpha/2))
         p = x/n
         b = z*sqrt((p*(1-p)+(z^2)/(4*n))/n)/(1+(z^2)/n)
@@ -114,7 +114,7 @@ module CI
     end
     #Wilson CC
     #Newcombe, R. G. (1998). "Two-sided confidence intervals for the single proportion: comparison of seven methods". Statistics in Medicine. 17 (8): 857–872. doi:10.1002/(SICI)1097-0258(19980430)17:8<857::AID-SIM777>3.0.CO;2-E. PMID 9595616
-    function propWilsonCCCI(x::Int, n::Int, alpha::Float64)::ConfInt
+    function propwilsonccci(x::Int, n::Int, alpha::Float64)::ConfInt
         z = abs(quantile(ZDIST, 1-alpha/2))
         p = x/n
         l = (2*n*p+z*z-1-z*sqrt(z*z-2-1/n+4*p*(n*(1-p)+1)))/2/(n+z*z)
@@ -301,8 +301,8 @@ module CI
         p2       = (x2/(n2-x2))
         estimate = p1/p2
         Z        = quantile(ZDIST, 1-alpha/2)
-        wilci1   = propWilsonCI(x1, n1, alpha)
-        wilci2   = propWilsonCI(x2, n2, alpha)
+        wilci1   = propwilsonci(x1, n1, alpha)
+        wilci2   = propwilsonci(x2, n2, alpha)
         wilci1   = ConfInt(wilci1.lower/(1-wilci1.lower), wilci1.upper/(1-wilci1.upper), estimate)
         wilci2   = ConfInt(wilci2.lower/(1-wilci2.lower), wilci2.upper/(1-wilci2.upper), estimate)
         lower    = (p1*p2-sqrt((p1*p2)^2 - wilci1.lower*wilci2.upper*(2*p1-wilci1.lower)*(2*p2-wilci2.upper)))/(wilci2.upper*(2*p2 - wilci2.upper))
@@ -347,8 +347,8 @@ module CI
         p2       = x2/n2
         estimate = p1-p2
         z        = quantile(ZDIST, 1 - alpha/2)
-        ci1      = propWilsonCI(x1, n1, alpha)
-        ci2      = propWilsonCI(x2, n2, alpha)
+        ci1      = propwilsonci(x1, n1, alpha)
+        ci2      = propwilsonci(x2, n2, alpha)
         return ConfInt(estimate-z*sqrt(ci1.lower*(1-ci1.lower)/n1+ci2.upper*(1-ci2.upper)/n2),
                        estimate+z*sqrt(ci1.upper*(1-ci1.upper)/n1+ci2.lower*(1-ci2.lower)/n2), estimate)
     end
@@ -358,8 +358,8 @@ module CI
         p2       = x2/n2
         estimate = p1-p2
         z        = quantile(ZDIST, 1 - alpha/2)
-        ci1      = propWilsonCCCI(x1, n1, alpha)
-        ci2      = propWilsonCCCI(x2, n2, alpha)
+        ci1      = propwilsonccci(x1, n1, alpha)
+        ci2      = propwilsonccci(x2, n2, alpha)
         return ConfInt(estimate-sqrt((p1-ci1.lower)^2 + (ci2.upper-p2)^2),
                        estimate+sqrt((ci1.upper - p1)^2 + (p2 - ci2.lower)^2), estimate)
     end
@@ -555,8 +555,8 @@ module CI
         p2       = (x2/n2)
         estimate = (x1/n1)/(x2/n2)
         Z        = quantile(ZDIST, 1-alpha/2)
-        wilci1   = propWilsonCI(x1, n1, alpha)
-        wilci2   = propWilsonCI(x2, n2, alpha)
+        wilci1   = propwilsonci(x1, n1, alpha)
+        wilci2   = propwilsonci(x2, n2, alpha)
         lower    = (p1*p2-sqrt((p1*p2)^2 - wilci1.lower*wilci2.upper*(2*p1-wilci1.lower)*(2*p2-wilci2.upper)))/(wilci2.upper*(2*p2 - wilci2.upper))
         upper    = (p1*p2+sqrt((p1*p2)^2 - wilci1.upper*wilci2.lower*(2*p1-wilci1.upper)*(2*p2-wilci2.lower)))/(wilci2.lower*(2*p2 - wilci2.lower))
         return ConfInt(lower, upper, estimate)

src/ClinicalTrialUtilities.jl

@@ -1,5 +1,4 @@
 # Clinical Trial Utilities
-# Version: 0.1.13
 # Author: Vladimir Arnautov aka PharmCat
 # Copyright © 2019 Vladimir Arnautov aka PharmCat (mail@pharmcat.net)
 # OwensQ/PowerTOST functions rewrited from https://github.com/Detlew/PowerTOST by Detlew Labes, Helmut Schuetz, Benjamin Lang
@@ -14,7 +13,7 @@
 # If you want to check and get R code you can find some here: http://powerandsamplesize.com/Calculators/
 __precompile__(true)
 module ClinicalTrialUtilities
-using Distributions, StatsBase, Statistics
+using Distributions, StatsBase, Statistics, Random
 using QuadGK
 using DataFrames
 import SpecialFunctions.lgamma
@@ -26,9 +25,11 @@ struct CTUException <: Exception
 end
 
 Base.showerror(io::IO, e::CTUException) = print("CTU Exception code: ", e.n, " Message: ", e.var);
-const ZDIST = Normal()
-const LOG2  = log(2)
-const VERSION = "0.1.13"
+const ZDIST  = Normal()
+const LOG2   = log(2)
+const PI2    = π * 2.0
+const PI2INV = 1.0 / (π * 2.0)
+const VERSION = "0.1.14"
 #Exceptions
 
 struct ConfInt
@@ -62,27 +63,29 @@ end
 export CTUException, ConfInt, NCA
 
 #Owen function calc: owensQ, owensQo, ifun1, owensTint2, owensT, tfn
-include("OwensQ.jl")
-#powerTOST calc: powerTOST, powerTOSTint, powerTOSTOwenQ, approxPowerTOST, power1TOST, approx2PowerTOST, cv2se, designProp
-include("PowerTOST.jl")
+include("owensq.jl")
+#powerTOST calc: powerTOST, powertostint, powerTOSTOwenQ, approxPowerTOST, power1TOST, approx2PowerTOST, cv2se, designProp
+include("powertost.jl")
 #Sample sise and Power atomic functions
-include("PowerSampleSize.jl")
+include("powersamplesize.jl")
 #Main sample size and power functions: sampleSize, ctPower, beSampleN
-include("SampleSize.jl")
+include("samplesize.jl")
 #Confidence interval calculation
 include("CI.jl")
 #Simulations
 include("SIM.jl")
 #PK
 include("PK.jl")
 #info function
-include("Info.jl")
+include("info.jl")
 #Descriptive statistics
-include("Descriptives.jl")
+include("descriptives.jl")
 #Frequency
-include("Freque.jl")
+include("freque.jl")
 #Export
 include("Export.jl")
+#Randomization
+include("randomization.jl")
 
 #Sample size
 export ctSampleN, beSampleN

src/Descriptives.jl → src/descriptives.jl

@@ -33,7 +33,7 @@ function descriptives(data::DataFrame;
     else
         vars = Array{Symbol,1}(undef, 0)
         for i in dfnames
-            if eltype(data[i]) <: Real push!(vars, i) end
+            if eltype(data[:, i]) <: Real push!(vars, i) end
         end
         if sort !== nothing filter!(x-> !(x in sort), vars) end
         if length(vars) == 0 error("Real[] columns not found!") end
@@ -69,10 +69,10 @@ function descriptives(data::DataFrame;
     if !isa(sort, Array)               #if no sort
         for v in vars
             sarray = Array{Any,1}(undef, 0)
-            deleteat!(data[:, v], findall(x->x === NaN || x === nothing, data[v]))
+            deleteat!(data[:, v], findall(x->x === NaN || x === nothing, data[:, v]))
             if length(data[:, v]) > 1
                 push!(sarray, v)
-                descriptives_!(sarray, data[:, v], stats)
+                append!(sarray, descriptives_(data[:, v], stats))
                 push!(dfs, sarray)
             end
         end
@@ -90,14 +90,15 @@ function descriptives(data::DataFrame;
                     push!(sdata, data[c, vars])
                 end
             end
-        sarray = Array{Any,1}(undef, 0)
-        push!(sarray, v)
+            sarray = Array{Any,1}(undef, 0)
+            push!(sarray, v)
             for c in sort
                 push!(sarray, sortlist[i, c])
             end
             deleteat!(sdata[:, v], findall(x->x === NaN || x === nothing, sdata[:, v]))
             if length(sdata[:, v]) > 1
-                descriptives_!(sarray, sdata[:, v], stats)
+                append!(sarray, descriptives_(sdata[:, v], stats))
+                #print(sarray)
                 push!(dfs, sarray)
             end
         end
@@ -106,8 +107,8 @@ function descriptives(data::DataFrame;
 end
 
 #Statistics calculation
-function descriptives_!(sarray::Array{Any,1}, data::Array{T,1}, stats::Array{Symbol,1}) where T <: Real
-
+function descriptives_(data::Array{T,1}, stats::Union{Tuple{Vararg{Symbol}}, Array{Symbol,1}}) where T <: Real
+    sarray = Array{Real,1}(undef, 0)
     dn         = nothing
     dmin       = nothing
     dmax       = nothing
@@ -213,4 +214,5 @@ function descriptives_!(sarray::Array{Any,1}, data::Array{T,1}, stats::Array{Sym
             push!(sarray, mode(data))
         end
     end
+    return sarray
 end

src/OwensQ.jl → src/owensq.jl

@@ -10,16 +10,16 @@
 
 #pnorm = cdf(ZDIST,  )
 #dnorm = pdf(ZDIST, )
-function owensQ(nu, t::Float64, delta::Float64, a::Float64, b::Float64)::Float64
+function owensq(nu::Real, t::Float64, delta::Float64, a::Float64, b::Float64)::Float64
     if a < 0 return  throw(CTUException(1011,"owensQ: a can not be < 0")) end
     if a==b return(0) end
     if a > b return throw(CTUException(1012,"owensQ: a can not be > b")) end
-    if a > 0 return owensQ(nu, t, delta, 0, b) - owensQ(nu, t, delta, 0, a)  end #not effective - double integration
+    if a > 0 return owensq(nu, t, delta, 0, b) - owensQ(nu, t, delta, 0, a)  end #not effective - double integration
     if nu < 29 && abs(delta) > 37.62
         if isinf(b)
             return quadgk(x -> ifun1(x, nu, t, delta), 0, 1, rtol=1.0E-8)[1]
         else
-            return owensQo(nu,t,delta,b)
+            return owensqo(nu,t,delta,b)
         end
     else
         if isinf(b)
@@ -43,7 +43,7 @@ end #owensQ
 # Port from PowerTOST - dlabes Mar 2012
 
 #OwensQOwen
-function owensQo(nu,t::Float64,delta::Float64,b::Float64;a::Float64=0.0)::Float64
+function owensqo(nu::Real, t::Float64, delta::Float64, b::Float64; a::Float64=0.0)::Float64
     if nu < 1  throw(CTUException(1001,"owensQo: nu can not be < 1")) end
     if a != 0.0 throw(CTUException(1002,"owensQo: a can not be not 0")) end
     if isinf(b) return cdf(NoncentralT(nu,delta),t) end
@@ -88,9 +88,9 @@ function owensQo(nu,t::Float64,delta::Float64,b::Float64;a::Float64=0.0)::Float6
                 sumt = sumt + M[i+1] + H[i+1]
             end
         end
-        qv = cdf(ZDIST, b) - 2*owensT(b,A-delta/b) -
-                        2*owensT(delta*sB, (delta*A*B-b)/B/delta) +
-                        2*owensT(delta*sB, A) - (delta>=0) + 2*sumt
+        qv = cdf(ZDIST, b) - 2*owenst(b,A-delta/b) -
+                        2*owenst(delta*sB, (delta*A*B-b)/B/delta) +
+                        2*owenst(delta*sB, A) - (delta>=0) + 2*sumt
     else
         if upr>=0
             for i = 0:2:upr
@@ -105,16 +105,16 @@ end #OwensQo
 #-------------------------------------------------------------------------------
 # functions bellow used in owensQ
 #PowerTost owensQ #37 and impl b for #52-54
-@inline function ifun1(x::Float64, nu, t::Float64, delta::Float64; b=0.0::Float64)::Float64
-    return owensTint2(b + x / (1 - x), nu, t, delta) * 1 / (1 - x) ^ 2
+@inline function ifun1(x::Float64, nu::Real, t::Float64, delta::Float64; b=0.0::Float64)::Float64
+    return owenstint2(b + x / (1 - x), nu, t, delta) * 1 / (1 - x) ^ 2
 end
 #.Q.integrand
-@inline function owensTint2(x::Float64, nu, t::Float64, delta::Float64)::Float64
+@inline function owenstint2(x::Float64, nu::Real, t::Float64, delta::Float64)::Float64
     if x == 0 return 0 end
     #Qconst::Float64 = -((nu/2.0)-1.0)*log(2.0)-lgamma(nu/2.0)
     #return sign(x)^(nu-1)*cdf(ZDIST,t*x/sqrt(nu)-delta)*exp((nu-1)*log(abs(x))-0.5*x^2+Qconst)
-    return sign(x) ^ (nu-1) * cdf(ZDIST, t * x / sqrt(nu) - delta) *
-    exp((nu - 1) * log(abs(x)) - 0.5 * x ^ 2 -
+    return sign(x) ^ (nu - 1.0) * cdf(ZDIST, t * x / sqrt(nu) - delta) *
+    exp((nu - 1.0) * log(abs(x)) - 0.5 * x ^ 2 -
     ((nu / 2.0) - 1.0) * LOG2 - lgamma(nu / 2.0))
 end
 #-------------------------------------------------------------------------------
@@ -125,7 +125,7 @@ end
 # https://people.sc.fsu.edu/~jburkardt/m_src/asa076/asa076.html
 # by J. Burkhardt, license GNU LGPL
 # rewrite from PowerTOST
-function owensT(h::Float64, a::Float64)::Float64
+function owenst(h::Float64, a::Float64)::Float64
     #not implemented
     epsilon = eps()
     # special cases
@@ -165,12 +165,12 @@ end #OwensT(h,a)
 # by J. Burkhardt license GNU LGPL
 # is called as tfn(h, a) if a<=1
 # otherwise as tfn(a*h, 1/a)
-@inline function tfn(x::T, fx::T)::Float64 where T <: Real
+@inline function tfn(x::Real, fx::Real)::Float64
     ng  = 5
     r   = Float64[0.1477621, 0.1346334, 0.1095432, 0.0747257, 0.0333357]
     u   = Float64[0.0744372, 0.2166977, 0.3397048, 0.4325317, 0.4869533]
     #tp::Float64  = 1 / 2 / pi
-    tp::Float64  = 0.15915494309189535
+    #tp::Float64  = 0.15915494309189535
     tv1 = eps();
     tv2::Float64 = tv3::Float64 = 15.0
     tv4::Float64 = 1.0E-05
@@ -181,15 +181,6 @@ end #OwensT(h,a)
     if tv3 <= (log(1.0 + fxs) - xs * fxs)
         x1::Float64  = 0.5 * fx
         fxs = 0.25 * fxs
-        #=
-        while true
-            rt  = fxs + 1.0
-            x2  = x1 + (xs * fxs + tv3 - log(rt)) / (2.0 * x1 * (1.0 / rt - xs))
-            fxs = x2 * x2
-            if abs(x2 - x1) < tv4 break end
-            x1 = x2
-        end
-        =#
         while true
             rt  = fxs + 1.0
             x2  = x1 + (xs * fxs + tv3 - log(rt)) / (2.0 * x1 * (1.0 / rt - xs))
@@ -213,23 +204,10 @@ end #OwensT(h,a)
         end
         =#
 
-            r1 = 1.0 .+ fxs .* (0.5 .+ u) .^ 2
-            r2 = 1.0 .+ fxs .* (0.5 .- u) .^ 2
-            rt = r .* (exp.(xs .* r1) ./ r1 .+ exp.(xs .* r2) ./ r2)
+        #Vectorized:
+    r1 = 1.0 .+ fxs .* (0.5 .+ u) .^ 2
+    r2 = 1.0 .+ fxs .* (0.5 .- u) .^ 2
+    rt = r .* (exp.(xs .* r1) ./ r1 .+ exp.(xs .* r2) ./ r2)
 
-    return sum(rt)*x2*tp
+    return sum(rt) * x2 * PI2INV
 end #tfn(x, fx)
-
-"""
-#-------------------------------------------------------------------------------
-#                        DEPRECATED
-#OT_integrand
-function owensTint(x, h)
-    return exp(-0.5*h^2*(1+x^2))/(1+x^2)/2/π
-end
-#OwensT_old
-function owensTo(h, a)
-    return quadgk(x -> owensTint(x, h),0, a)[1]
-end
-#-------------------------------------------------------------------------------
-"""