Merge pull request #47 from PharmCat/dev

v0.3.0

    #44 , #45 , #46
    update test
    Move alpha from CTask to Hypothesis
    besim method for bioequivalence simulation
    checkhyp changed for simulations methods
    minor changes in text output
    theo test dataset; theo :luld AUClast test
    TaskResult changes: field result -> res, t.result returns t.res[:result]
    simulation hypothesis generalization
    minor changes in power formulas: predict negative power
    simulation bug: fix Equivalence checkhyp
    pval: experimental function
    changes in Show
    documenting
    code cosmetics

.travis.yml

@@ -2,19 +2,19 @@ language: julia
 os:
   - linux
   - osx
+  - windows
 
 julia:
   - 1.0
   - 1.1
   - 1.2
   - 1.3
   - 1.4
+  - 1.5
 
 branches:
   only:
   - master
-  - dev
-
 
 notifications:
   email: false
@@ -29,6 +29,12 @@ after_success:
   - julia -e 'using Pkg; cd(Pkg.dir("ClinicalTrialUtilities")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
 
 jobs:
+  allow_failures:
+    - os:
+        - windows
+    - julia:
+        - 1.4
+        - 1.5
   include:
     - stage: "Documentation"
       julia: 1.0

Project.toml

@@ -1,18 +1,18 @@
 authors = ["Vladimir Arnautov ([email protected])"]
 name = "ClinicalTrialUtilities"
 uuid = "535c2557-d7d0-564d-8ff9-4ae146c18cfe"
-version = "0.2.7"
+version = "0.3.0"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 
 [extras]
@@ -27,10 +27,10 @@ test = ["CSV", "Test", "Plots"]
 julia = "1.0, 1.1, 1.2, 1.3, 1.4"
 Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23"
 StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33"
-DataFrames = "0.19, 0.20"
-QuadGK = "2.0, 2.1, 2.2, 2.3"
+QuadGK = "2.0, 2.1, 2.2, 2.3, 2.4"
 SpecialFunctions = "0.8, 0.9, 0.10"
 Roots = "0.7, 0.8, 1.0"
-CSV = "0.5, 0.6"
 RecipesBase = "0.7, 0.8, 1.0"
 Reexport = "0.1, 0.2"
+DataFrames = "0.19, 0.20"
+CSV = "0.5, 0.6"

README.md

@@ -8,6 +8,10 @@
 [![Coverage Status](https://coveralls.io/repos/github/PharmCat/ClinicalTrialUtilities.jl/badge.svg?branch=master)](https://coveralls.io/github/PharmCat/ClinicalTrialUtilities.jl?branch=master)
 [![Latest docs](https://img.shields.io/badge/docs-latest-blue.svg)](https://pharmcat.github.io/ClinicalTrialUtilities.jl/dev/)
 
+## Difference since v0.2.7
+
+  Changes in Equivalence Hypothesis alpha level! See [here](https://github.com/PharmCat/ClinicalTrialUtilities.jl/issues/46)!
+
 ## Description
 
 The package is designed to perform calculations related to the planning and analysis of the results of clinical trials. The package includes the basic functions described below, as well as a few modules to perform specific calculations.

cange.log

@@ -1,3 +1,20 @@
+v0.3.0
+  - #44 , #45 , #46
+  - update test
+  - Move alpha from CTask to Hypothesis
+  - besim method for bioequivalence simulation
+  - checkhyp changed for simulations methods
+  - minor changes in text output
+  - theo test dataset; theo :luld AUClast test
+  - TaskResult changes: field result -> res, t.result returns t.res[:result]
+  - simulation hypothesis generalization
+  - minor changes in power formulas: predict negative power
+  - simulation bug: fix Equivalence checkhyp
+  - pval: experimental function
+  - changes in Show
+  - documenting
+  - code cosmetics
+
 v0.2.7
   - setkelauto!
   - getkelauto

docs/src/ci.md

@@ -2,6 +2,24 @@
 
 ## Proportions
 
+### One proportion
+
+P = n / x
+
+where n - number of outcomes; x - number of observations.
+
+### Absolute risk difference
+
+Diff(𝛿) = P₁ - P₂ = n₁ / x₁ - n₂ / x₂
+
+### Risk Ratio
+
+RR = P₁ / P₂ = (n₁ / x₁) - (n₂ / x₂)
+
+### Odd Ratio
+
+OR = (n₁ / (x₁ - n₁)) - (n₂ / (x₂ - n₂))
+
 ### propci
 ```@docs
 ClinicalTrialUtilities.propci
@@ -29,7 +47,7 @@ ClinicalTrialUtilities.orpropci
 ClinicalTrialUtilities.meanci
 ```
 
-### .diffmeanci
+### diffmeanci
 ```@docs
 ClinicalTrialUtilities.diffmeanci
 ```

docs/src/pkplot.md

@@ -8,25 +8,23 @@ ClinicalTrialUtilities.pkplot
 #### Examples
 
 
-```@example 1
+```@example pkplots
 using ClinicalTrialUtilities, DataFrames, CSV, Plots
-pkdatapath = joinpath(dirname(pathof(ClinicalTrialUtilities)))*"\\..\\test\\csv\\pkdata2.csv"
+pkdatapath = joinpath(dirname(pathof(ClinicalTrialUtilities)))*"/../test/csv/pkdata2.csv"
 pkdata     = CSV.File(pkdatapath) |> DataFrame
 pkds       = ClinicalTrialUtilities.pkimport(pkdata, [:Subject, :Formulation]; conc = :Concentration, time = :Time)
 plot1      = pkplot(pkds[1], legend = false);
-#savefig("pkplot1.svg"); # hide
 plot2      = pkplot(pkds;  pagesort = [:Formulation], typesort = [:Subject])[1];
-#savefig("pkplot2.svg"); nothing # hide
 ```
 
 Plot for subject:
 
-```@example 1
+```@example pkplots
 plot1
 ```
 
 First plot for DataSet:
 
-```@example 1
+```@example pkplots
 plot2
 ```

docs/src/samplesize.md

@@ -2,6 +2,50 @@
 
 Sample size estimation.
 
+## Introduction
+
+### Hypothesis
+
+#### Equality
+
+Two-sided hypothesis:
+
+H₀: A - B = 0
+Hₐ: A - B ≠ 0
+
+100\*(1 - α)% two-sided confidence interval used to demonstrate it.
+
+#### Equivalence | Bioequivalence
+
+Two-sided hypothesis or two one-sided hypothesis (TOST):
+
+H₀: |A − B| ≥ δ
+Hₐ: |A − B| < δ
+
+or
+
+H₀: A - B ≤ δˡ || A - B ≥ δᵘ
+Hₐ: δˡ < A - B < δᵘ
+
+100\*(1 - 2\*α)% two-sided confidence interval used to demonstrate it.
+
+See Walker&Nowacki, 2011, Understanding Equivalence and Noninferiority Testing.
+
+*Chow at al. Sample Size Calculations in Clinical Research 3-rd ed. for sample size calculation used alpha level for testing corresponding 100\*(1 - α)% two-sided confidence interval*
+
+For bioequivalence clinical trial 90% *(100\*(1 - 2\*α)%)* two-sided confidence interval is used (α = 0.05), but for therapeutic equivalence 95% two-sided confidence interval should be used (in this case use α = 0.025).
+
+#### Non-inferiority | Superiority
+
+One-sided hypothesis:
+
+H₀: A − B ≤ δ
+Hₐ: A − B > δ
+
+100\*(1 - α)% one-sided confidence interval used to demonstrate it.
+
+In clinical research α = 0.025 or less should be used, or 100*(1 - α)% two-sided confidence interval should be used.
+
 ### ctsamplen
 ```@docs
 ClinicalTrialUtilities.ctsamplen

src/ClinicalTrialUtilities.jl

@@ -21,12 +21,23 @@ using Distributions, Random, Roots, QuadGK, RecipesBase, Reexport
 import SpecialFunctions
 import Base: show, findfirst, getproperty, showerror, getindex, length, in, iterate, eltype, deleteat!, findall
 import StatsBase.confint
-import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows!, rename!
+import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows!, rename!, AbstractDataFrame
 
-function lgamma(x)
-    return SpecialFunctions.logabsgamma(x)[1]
+try
+    methods(SpecialFunctions.logabsgamma)
+    global lgamma(x) = SpecialFunctions.logabsgamma(x)[1]
+catch
+    global lgamma(x) = SpecialFunctions.lgamma(x)
 end
 
+try
+    if collect(methods(DataFrames.delete!, (AbstractDataFrame, Any)))[1].file == Symbol("deprecated.jl")
+        DataFrames.delete!(df::AbstractDataFrame, inds) = deleterows!(df, inds)
+    end
+catch
+end
+
+
 const ZDIST  = Normal()
 const LOG2   = log(2)
 const PI2    = π * 2.0
@@ -41,11 +52,11 @@ include("proportion.jl")
 
 include("means.jl")
 
+include("hypothesis.jl")
+
 #Confidence interval calculation
 include("ci.jl")
 
-include("hypothesis.jl")
-
 #Owen function calc: owensQ, owensQo, ifun1, owensTint2, owensT, tfn
 include("owensq.jl")
 #powerTOST calc: powerTOST, powertostint, powerTOSTOwenQ, approxPowerTOST, power1TOST, approx2PowerTOST, cv2se, designProp
@@ -129,6 +140,7 @@ setkelauto!,
 getkelauto,
 # Simulation
 ctsim,
+besim,
 #Plots
 pkplot,
 pkplot!

src/ci.jl

@@ -14,12 +14,16 @@ struct ConfInt
     upper::Real
     estimate::Real
     alpha::Real
+    #method::Symbol
     function ConfInt(lower, upper, estimate)
         new(lower, upper, estimate, NaN)::ConfInt
     end
     function ConfInt(lower, upper, estimate, alpha)
         new(lower, upper, estimate, alpha)::ConfInt
     end
+    #function ConfInt(lower, upper, estimate, alpha, method)
+    #    new(lower, upper, estimate, alpha, method)::ConfInt
+    #end
 end
 
 function getindex(a::ConfInt, b::Int)
@@ -68,6 +72,43 @@ end
 function StatsBase.confint(param::DiffMean{true}; level = 0.95, method = :default)::ConfInt where T <: Mean
     diffmeanci(param.a.m, param.a.sd^2, param.a.n, param.b.m, param.b.sd^2, param.b.n; alpha = 1 - level, method = method)
 end
+
+"""
+
+P value.
+"""
+function pval(param::DiffProportion{Proportion{true}, Proportion{true}}, hyp::Equivalence; method = :default, atol = 1E-6)
+    basepval(param, hyp; method = method, atol = atol)
+end
+function pval(param::DiffProportion{Proportion{true}, Proportion{true}}, hyp::Superiority; method = :default, atol = 1E-6)
+    basepval(param, hyp; method = method, atol = atol)
+end
+function pval(param::DiffProportion{Proportion{true}, Proportion{true}}, hyp::Equality; method = :default, atol = 1E-6)
+    basepval(param, hyp; method = method, atol = atol)
+end
+function basepval(param::T, hyp::Equivalence; method = :default, atol = 1E-6) where T <: AbstractParameter
+    fx = x -> hyp.llim - confint(param; level = 1 - x, method = method).lower
+    p1 = find_zero(fx, (1-eps(), eps()), atol = atol)/2
+    fx = x -> hyp.ulim - confint(param; level = 1 - x, method = method).upper
+    p2 = find_zero(fx, (1-eps(), eps()), atol = atol)/2
+    max(p1, p2)
+end
+function basepval(param::T, hyp::Superiority; method = :default, atol = 1E-6) where T <: AbstractParameter
+    if confint(param; method = method).estimate > hyp.diff
+        fx = x -> hyp.diff - confint(param; level = 1 - x, method = method).lower
+        return find_zero(fx, (1-eps(), eps()), atol = atol)/2
+    else return 1.0 end
+end
+function basepval(param::T, hyp::Equality; method = :default, atol = 1E-6) where T <: AbstractParameter
+    if confint(param; method = method).estimate > hyp.val
+        fx = x -> hyp.val - confint(param; level = 1 - x, method = method).lower
+        p = find_zero(fx, (1-eps(), eps()), atol = atol)
+    else
+        fx = x -> hyp.val - confint(param; level = 1 - x, method = method).upper
+        p = find_zero(fx, (1-eps(), eps()), atol = atol)
+    end
+    p
+end
 """
     propci(x::Int, n::Int; alpha=0.05, method = :default)::ConfInt
 
@@ -534,28 +575,28 @@ end
     #Method of Mee 1984 with Miettinen and Nurminen modification n / (n - 1) Newcombe 1998
     #Score intervals for the difference of two binomial proportions
     #6
-    function propdiffmnci(x1::Int, n1::Int, x2::Int, n2::Int, alpha::Real)::ConfInt
+    function propdiffmnci(x1::Int, n1::Int, x2::Int, n2::Int, alpha::Real; atol::Float64 = 1E-6)::ConfInt
         ci       = propdiffnhsccci(x1, n1, x2, n2, alpha) #approx zero bounds
         p1       = x1/n1
         p2       = x2/n2
         estimate = p1 - p2
         z        = quantile(Chisq(1), 1-alpha)
         fmnd(x)  = mndiffval(p1, n1, p2, n2, estimate, x) - z
-        if fmnd(ci.lower)*fmnd(estimate-1e-8) < 0.0
+        if fmnd(ci.lower)*fmnd(estimate-eps()) < 0.0
             ll = ci.lower
-            lu = estimate-1e-8
+            lu = estimate-eps()
         else
-            ll = -1.0+1e-8
+            ll = -1.0+eps()
             lu = ci.lower
         end
-        if fmnd(ci.upper)*fmnd(estimate+1e-8) < 0.0
-            ul = estimate+1e-8
+        if fmnd(ci.upper)*fmnd(estimate+eps()) < 0.0
+            ul = estimate+eps()
             uu = ci.upper
         else
             ul = ci.upper
-            uu = 1.0-1e-8
+            uu = 1.0-eps()
         end
-        return ConfInt(find_zero(fmnd, (ll, lu), atol=1E-6), find_zero(fmnd, (ul, uu), atol=1E-6), estimate, alpha)
+        return ConfInt(find_zero(fmnd, (ll, lu), atol = atol), find_zero(fmnd, (ul, uu), atol = atol), estimate, alpha)
     end
     @inline function mndiffval(p1::Real, n1::Int, p2::Real, n2::Int, estimate::Real, Δ::Real)::AbstractFloat
         return (estimate-Δ)^2/((n1+n2)/(n1+n2-1)*mlemndiff(p1, n1, p2, n2, Δ))
@@ -700,9 +741,9 @@ end
     function meandiffev(m1::Real, σ²1::Real, n1::Real, m2::Real, σ²2::Real, n2::Real, alpha::Real)::ConfInt
         diff   = float(m1 - m2)
         stddev = sqrt(((n1 - 1) * σ²1 + (n2 - 1) * σ²2) / (n1 + n2 - 2))
-        stderr = stddev * sqrt(1/n1 + 1/n2)
-        d      = stderr*quantile(TDist(n1+n2-2), 1-alpha/2)
-        return ConfInt(diff-d, diff+d, diff, alpha)
+        stderr = stddev * sqrt(1 / n1 + 1 / n2)
+        d      = stderr * quantile(TDist(n1 + n2 - 2), 1 - alpha / 2)
+        return ConfInt(diff - d, diff + d, diff , alpha)
     end
     function meandiffev(a1::AbstractVector{T}, a2::AbstractVector{S}, alpha::Real)::ConfInt where {T<:Real,S<:Real}
         return meandiffev(mean(a1), var(a1), length(a1), mean(a2), var(a2), length(a2), alpha)
@@ -714,8 +755,8 @@ end
         diff   = float(m1 - m2)
         v      = (σ²1 / n1 + σ²2 / n2)^2 /( σ²1^2 / n1^2 / (n1 - 1) + σ²2^2 / n2^2 / (n2 - 1))
         stderr = sqrt(σ²1 / n1 + σ²2 / n2)
-        d      = stderr*quantile(TDist(v), 1-alpha/2)
-        return ConfInt(diff-d, diff+d, diff, alpha)
+        d      = stderr * quantile(TDist(v), 1 - alpha / 2)
+        return ConfInt(diff - d, diff + d, diff, alpha)
     end
     function meandiffuv(a1::AbstractVector{T}, a2::AbstractVector{S}, alpha::Real)::ConfInt where {T<:Real,S<:Real}
         return meandiffuv(mean(a1), var(a1), length(a1), mean(a2), var(a2), length(a2), alpha)