From 84a2dd3893ef24d90c3d3237256e1af29852a4e0 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sun, 29 Mar 2020 22:24:16 +0300 Subject: [PATCH 1/8] pk plots, docs, type fix, tryfloatparse, test, update --- Project.toml | 2 +- README.md | 29 ++++++-- cange.log | 9 +++ src/ClinicalTrialUtilities.jl | 3 +- src/pk.jl | 98 ++++++++++++++++++++------- src/plots.jl | 121 ++++++++++++++++++++++++++-------- test/pktest.jl | 6 ++ 7 files changed, 209 insertions(+), 59 deletions(-) diff --git a/Project.toml b/Project.toml index 844fa40..53f13a1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ authors = ["Vladimir Arnautov (mail@pharmcat.net)"] name = "ClinicalTrialUtilities" uuid = "535c2557-d7d0-564d-8ff9-4ae146c18cfe" -version = "0.2.5" +version = "0.2.6" [deps] Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/README.md b/README.md index 4ebb31c..696a86d 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,18 @@ The package is designed to perform calculations related to the planning and anal using Pkg; Pkg.add("ClinicalTrialUtilities"); ``` -### Usage +### Main features + +- SampleSize calculation +- Power calculation +- Confidence Interval calculation +- NCA Pharmacokinetics parameters calculation +- Randomization + + +### Examples + +#### SampleSize **NB! Hypothesis types:** @@ -25,9 +36,6 @@ using Pkg; Pkg.add("ClinicalTrialUtilities"); - :ei - Equivalencens: two one-sided hypothesis; - :ns - Non-Inferiority / Superiority: one-sided hypothesis, for some cases you should use two-sided hypothesis for Non-Inferiority/Superiority, you can use alpha/2 for this; -### Examples - -#### SampleSize ``` #Sample size for one proportion equality ctsamplen(param=:prop, type=:ea, group=:one, a=0.3, b=0.5) @@ -93,6 +101,19 @@ pooledcv([0.12, 0.2, 0.25], [14, 22, 32], [:d2x2, :d2x2, :d2x2]) ``` +#### NCA +``` +using CSV, DataFrames, ClinicalTrialUtilities +pkdata2 = CSV.File("pkdata.csv") |> DataFrame +pkds = pkimport(pkdata2, [:Subject, :Formulation]; time = :Time, conc = :Concentration) +pk = nca!(pkds) +ncadf = DataFrame(pk; unst = true) +ds = ClinicalTrialUtilities.descriptive(ncadf, stats = [:n, :mean, :sd], sort = [:Formulation]) +dsdf = ClinicalTrialUtilities.DataFrame(ds; unst = true) + +``` + + ### Copyrights Clinical Trial Utilities diff --git a/cange.log b/cange.log index d09ff5e..7b78d59 100644 --- a/cange.log +++ b/cange.log @@ -1,3 +1,12 @@ +v0.2.6 + + - findfirst + - Documents update + - export applyncarule! + - code cleaning + - improve coverage + - plots kwargs + v0.2.5 - SpecialFunctions 0.10 diff --git a/src/ClinicalTrialUtilities.jl b/src/ClinicalTrialUtilities.jl index 3a43ad3..d50fcf0 100644 --- a/src/ClinicalTrialUtilities.jl +++ b/src/ClinicalTrialUtilities.jl @@ -19,7 +19,7 @@ using Distributions, Random, Roots, QuadGK, RecipesBase, Reexport @reexport using StatsBase import SpecialFunctions -import Base.show +import Base: show, findfirst import Base.showerror import Base.getindex import Base.length @@ -124,6 +124,7 @@ DataFrame, ElimRange, DoseTime, LimitRule, +applyncarule!, keldata, # Simulation ctsim, diff --git a/src/pk.jl b/src/pk.jl index 80249ca..75f110c 100644 --- a/src/pk.jl +++ b/src/pk.jl @@ -16,7 +16,7 @@ function in(a::Pair, b::Dict) return true end #!!! -struct KelData <: AbstractData +struct KelData s::Array{Int, 1} e::Array{Int, 1} a::Array{Number, 1} @@ -143,14 +143,25 @@ mutable struct PKPDProfile{T <: AbstractSubject} <: AbstractData new{T}(subject, method, result, subject.sort) end end +""" + LimitRule(;lloq::Real = NaN, btmax::Real = NaN, atmax::Real = NaN, nan::Real = NaN, rm::Bool = false) + + + LimitRule(lloq::Real, btmax::Real, atmax::Real, nan::Real, rm::Bool) +Rule for PK subject. + +STEP 1 (NaN step): replace all NaN values with nan reyword value (if nan !== NaN); +STEP 2 (LLOQ step): replace values below lloq with btmax value if this value befor Tmax or with atmax if this value after Tmax (if lloq !== NaN); +STEP 3 (remove NaN): rm == true, then remove all NaN values; +""" struct LimitRule lloq::Real btmax::Real atmax::Real nan::Real rm::Bool - function LimitRule(lloq, btmax, atmax, nan, rm) + function LimitRule(lloq::Real, btmax::Real, atmax::Real, nan::Real, rm::Bool) new(lloq, btmax, atmax, nan, rm)::LimitRule end function LimitRule(lloq, btmax, atmax, nan) @@ -159,7 +170,7 @@ struct LimitRule function LimitRule(lloq, btmax, atmax) new(lloq, btmax, atmax, NaN, true)::LimitRule end - function LimitRule(;lloq = NaN, btmax = NaN, atmax = NaN, nan = NaN, rm = false) + function LimitRule(;lloq::Real = NaN, btmax::Real = NaN, atmax::Real = NaN, nan::Real = NaN, rm::Bool = false) new(lloq, btmax, atmax, nan, rm)::LimitRule end end @@ -187,12 +198,14 @@ end =# #------------------------------------------------------------------------------- +#= function obsnum(data::T) where T <: AbstractSubject return length(data.time) end function obsnum(keldata::KelData) return length(keldata.a) end +=# function length(data::T) where T <: AbstractSubject return length(data.time) end @@ -301,7 +314,7 @@ end function ctmax(data::PKSubject) return ctmax(data.time, data.obs) end - +#= function ctmax(data::PKSubject, dosetime) s = 0 for i = 1:length(data) - 1 @@ -314,7 +327,7 @@ end cmax, tmax, tmaxn = ctmax(data.time[mask], data.obs[mask]) if cmax > cpredict return cmax, tmax, tmaxn + s else return cpredict, data.dosetime.time, s end end - +=# function dropbeforedosetime(time::Vector, obs::Vector, dosetime::DoseTime) s = 0 for i = 1:length(time) @@ -545,7 +558,11 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: result[:Cdose] = 0 end doseaucpart, doseaumcpart = aucpart(data.dosetime.time, time[1], result[:Cdose], obs[1], :lint, false) - else #Dosetime >= first time + elseif data.dosetime.time == time[1] + result[:Cdose] = obs[1] + else + error("Some concentration before dose time after filtering!!!") + #= for i = 1:result[:Obsnum] - 1 if time[i] <= data.dosetime.time < time[i+1] if time[i] == data.dosetime.time @@ -554,8 +571,10 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: else if data.dosetime.tau > 0 result[:Cdose] = result[:Ctaumin] + else result[:Cdose] = 0 + end doseaucpart, doseaumcpart = aucpart(data.dosetime.time, time[i + 1], result[:Cdose], obs[i + 1], :lint, false) pmask[i] = false @@ -565,6 +584,7 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: pmask[i] = false end end + =# end #sum all full AUC parts and part before dose @@ -756,7 +776,7 @@ function nca!(data::PDSubject; verbose = false, io::IO = stdout)::PKPDProfile{PD result[:BL] = data.bl result[:RMAX] = maximum(data.obs) result[:RMIN] = maximum(data.obs) - for i = 2:obsnum(data) #BASELINE + for i = 2:length(data) #BASELINE if data.obs[i - 1] <= result[:BL] && data.obs[i] <= result[:BL] result[:TBBL] += data.time[i,] - data.time[i - 1] result[:AUCBBL] += linauc(data.time[i - 1], data.time[i], result[:BL] - data.obs[i - 1], result[:BL] - data.obs[i]) @@ -895,6 +915,29 @@ end #------------------------------------------------------------------------------- +function tryfloatparse!(x) + for i = 1:length(x) + if typeof(x[i]) <: AbstractString + try + x[i] = parse(Float64, x[i]) + catch + x[i] = NaN + end + else + try + if x[i] === missing + x[i] = NaN + else + x[i] = float(x[i]) + end + catch + x[i] = NaN + end + end + end + return convert(Vector{Float64}, x) +end + """ pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, time::Symbol) @@ -909,6 +952,7 @@ Pharmacokinetics data import from DataFrame. """ function pkimport(data::DataFrame, sort::Array, rule::LimitRule; conc::Symbol, time::Symbol)::DataSet results = Array{PKSubject, 1}(undef, 0) + if !(eltype(data[!, conc]) <: Real) throw(ArgumentError("Type of concentration data not <: Real!"))end getdatai(data, sort, [conc, time], (x, y) -> push!(results, PKSubject(copy(x[!, time]), copy(x[!, conc]), sort = Dict(sort .=> collect(y)))); sortby = time) results = DataSet(results) applyncarule!(results, rule) @@ -998,6 +1042,15 @@ function upkimport(data::DataFrame, sort::Array; stime::Symbol, etime::Symbol, c end #------------------------------------------------------------------------------- +""" + applyncarule!(data::PKSubject, rule::LimitRule) + +Apply rule to PK subject . + +STEP 1 (NaN step): replace all NaN values with nan reyword value (if nan !== NaN); +STEP 2 (LLOQ step): replace values below lloq with btmax value if this value befor Tmax or with atmax if this value after Tmax (if lloq !== NaN); +STEP 3 (remove NaN): rm == true, then remove all NaN values; +""" function applyncarule!(data::PKSubject, rule::LimitRule) cmax, tmax, tmaxn = ctmax(data) #NaN Rule @@ -1028,31 +1081,14 @@ function applyncarule!(data::PKSubject, rule::LimitRule) data.obs = data.obs[filterv] end end - #= - function applyncarule(data::DataSet{PKSubject}, rule::LimitRule) - end - - function applyncarule(data::PKPDProfile{PKSubject}, rule::LimitRule) - end - function applyncarule(data::DataSet{PKPDProfile{PKSubject}}, rule::LimitRule) - end - - function applyncarule(data::PKSubject, rule::LimitRule) - end - =# function applyncarule!(data::DataSet{PKSubject}, rule::LimitRule) for i in data applyncarule!(i, rule) end data end - #= - function applyncarule!(data::PKPDProfile{PKSubject}, rule::LimitRule) - end - function applyncarule!(data::DataSet{PKPDProfile{PKSubject}}, rule::LimitRule) - end - =# + #--------------------------------------------------------------------------- """ setelimrange!(data::PKSubject, range::ElimRange) @@ -1175,6 +1211,18 @@ function setdosetime!(data::DataSet{PKSubject}, dosetime::DoseTime, sort::Dict) end data end +""" + findfirst(sort::Dict, data::DataSet{PKSubject}) + +The first item in the dataset that satisfies the condition - sort dictionary. +""" +function findfirst(sort::Dict, data::DataSet{PKSubject}) + for i = 1:length(data) + if sort ∈ data[i].sort + return data[i] + end + end +end """ dosetime(data::PKSubject) diff --git a/src/plots.jl b/src/plots.jl index 8ad9f5c..e3cf887 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -72,36 +72,50 @@ function plotlabel(d) return title end -function plot(subj::T; title = nothing, legend = true, label = "AUTO", ylabel = "Concentration", xlims = nothing, plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject +function plot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject + kwargs = Dict{Symbol, Any}(kwargs) + k = keys(kwargs) - if title === nothing - title = plotlabel(subj.sort) + if !(:title in k) + kwargs[:title] = plotlabel(subj.sort) end - if xlims === nothing xlims = (minimum(subj.time), maximum(subj.time)) end + if !(:xlims in k) + kwargs[:xlims] = (minimum(subj.time), maximum(subj.time)) + end + if !(:legend in k) + kwargs[:legend] = true + end + if !(:ylabel in k) + kwargs[:ylabel] = "Concentration" + end + p = pkplot(subj.time, subj.obs; - title = title, linestyle = plotstyle.linestyle, linecolor = plotstyle.linecolor, markershape = plotstyle.markershape, markercolor = plotstyle.markercolor, - legend = legend, - label = label, - ylabel = ylabel, - xlims = xlims, + kwargs... ) return p end -function plot!(subj::T; legend = true, label = "AUTO", xlims = nothing, plotstyle::PKPlotStyle = PKPLOTSTYLE[1]) where T <: AbstractSubject - if xlims === nothing xlims = (minimum(subj.time), maximum(subj.time)) end +function plot!(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject + kwargs = Dict{Symbol, Any}(kwargs) + k = keys(kwargs) + + if !(:xlims in k) + kwargs[:xlims] = (minimum(subj.time), maximum(subj.time)) + end + if !(:legend in k) + kwargs[:legend] = true + end + p = pkplot!(subj.time, subj.obs; linestyle = plotstyle.linestyle, linecolor = plotstyle.linecolor, markershape = plotstyle.markershape, markercolor = plotstyle.markercolor, - legend = legend, - label = label, - xlims = xlims, + kwargs... ) return p end @@ -115,7 +129,23 @@ function usort(data::DataSet{T}, list) where T <: AbstractData dl end -function pageplot(pagedatatmp, styledict, typesort; title = "", ylabel = "Concentration", legend = true, xlims = nothing) +function pageplot(pagedatatmp, styledict, typesort; kwargs...) + kwargs = Dict{Symbol, Any}(kwargs) + k = keys(kwargs) + + if !(:title in k) + kwargs[:title] = "" + end + if !(:xlims in k) + kwargs[:xlims] = (minimum(subj.time), maximum(subj.time)) + end + if !(:legend in k) + kwargs[:legend] = true + end + if !(:ylabel in k) + kwargs[:ylabel] = "Concentration" + end + utypes = keys(styledict) labels = Vector{String}(undef, 0) fst = true @@ -134,34 +164,68 @@ function pageplot(pagedatatmp, styledict, typesort; title = "", ylabel = "Concen label = "" end end + + kwargs[:label] = label if fst - p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = style, label = label, ylabel = ylabel) + p = plot(d; plotstyle = style, kwargs...) fst = false else - plot!( d; legend = legend, xlims = xlims, plotstyle = style, label = label) + plot!( d; plotstyle = style, kwargs...) end end end end p end -function pageplot(pagedatatmp; title = "", ylabel = "Concentration", legend = true, xlims = nothing) +function pageplot(pagedatatmp; kwargs...) + kwargs = Dict{Symbol, Any}(kwargs) + k = keys(kwargs) + + if !(:title in k) + kwargs[:title] = "" + end + if !(:xlims in k) + kwargs[:xlims] = (minimum(subj.time), maximum(subj.time)) + end + if !(:legend in k) + kwargs[:legend] = true + end + if !(:ylabel in k) + kwargs[:ylabel] = "Concentration" + end + fst = true p = nothing for d in pagedatatmp - label = "AUTO" + + kwargs[:label] = "AUTO" if fst - p = plot(d; title = title, legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label, ylabel = ylabel) + p = plot(d; plotstyle = PKPLOTSTYLE[1], kwargs...) fst = false else - plot!( d; legend = legend, xlims = xlims, plotstyle = PKPLOTSTYLE[1], label = label) + plot!( d; plotstyle = PKPLOTSTYLE[1], kwargs...) end end p end -function plot(data::DataSet{T}; title = nothing, ylabel = "Concentration", legend = true, xlims = nothing, pagesort = nothing, typesort = nothing) where T <: AbstractSubject +function plot(data::DataSet{T}; pagesort = nothing, typesort = nothing, kwargs...) where T <: AbstractSubject # Style types + kwargs = Dict{Symbol, Any}(kwargs) + k = keys(kwargs) + + if !(:title in k) + kwargs[:title] = :AUTO + end + if !(:xlims in k) + kwargs[:xlims] = nothing + end + if !(:legend in k) + kwargs[:legend] = true + end + if !(:ylabel in k) + kwargs[:ylabel] = "Concentration" + end styledict = nothing #style = PKPLOTSTYLE[1] @@ -181,7 +245,7 @@ function plot(data::DataSet{T}; title = nothing, ylabel = "Concentration", legen end end autotitle = false - if title === nothing autotitle = true end + if kwargs[:title] == :AUTO autotitle = true end #--------------------------------------------------------------------------- # PAGES plots = Vector{Any}(undef, 0) @@ -190,7 +254,7 @@ function plot(data::DataSet{T}; title = nothing, ylabel = "Concentration", legen pagedict = usort(data, pagesort) for i in pagedict if autotitle - title = plotlabel(i) + kwargs[:title] = plotlabel(i) end pagedatatmp = Vector{eltype(data)}(undef, 0) for d in data @@ -199,17 +263,18 @@ function plot(data::DataSet{T}; title = nothing, ylabel = "Concentration", legen end end if typesort !== nothing - p = pageplot(pagedatatmp, styledict, typesort; title = title, legend = legend, xlims = xlims, ylabel = ylabel) + p = pageplot(pagedatatmp, styledict, typesort; kwargs...) else - p = pageplot(pagedatatmp; title = title, legend = legend, xlims = xlims, ylabel = ylabel) + p = pageplot(pagedatatmp; kwargs...) end push!(plots, p) end else + if kwargs[:title] == :AUTO kwargs[:title] = "" end if typesort !== nothing - p = pageplot(data, styledict, typesort; title = title, legend = legend, xlims = xlims, ylabel = ylabel) + p = pageplot(data, styledict, typesort; kwargs...) else - p = pageplot(data; title = title, legend = legend, xlims = xlims, ylabel = ylabel) + p = pageplot(data; kwargs...) end push!(plots, p) end diff --git a/test/pktest.jl b/test/pktest.jl index 09b4654..ff8f07c 100644 --- a/test/pktest.jl +++ b/test/pktest.jl @@ -930,6 +930,7 @@ println(" ---------------------------------- ") # NaN PK LimitRule test + nanpkdata.Concentration = ClinicalTrialUtilities.tryfloatparse!(nanpkdata.Concentration) pkds = ClinicalTrialUtilities.pkimport(nanpkdata, [:Subject, :Formulation]; conc = :Concentration, time = :Time) rule = ClinicalTrialUtilities.LimitRule(0, 0, NaN, 0, true) ClinicalTrialUtilities.applyncarule!(pkds, rule) @@ -940,6 +941,11 @@ println(" ---------------------------------- ") rule = ClinicalTrialUtilities.LimitRule(0, 0, NaN, -1, false) ClinicalTrialUtilities.applyncarule!(pkds, rule) @test pkds[3].obs[6] === NaN + + ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(100,0), Dict(:Subject => 1, :Formulation => "R")) + @test pkds[3].dosetime.dose == 100 + + pks = ClinicalTrialUtilities.findfirst(Dict(:Subject => 1, :Formulation => "R"), pkds) end println(" ---------------------------------- ") From d35ecdb8669a622b641589368c52232e4114bdff Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 3 Apr 2020 18:49:54 +0300 Subject: [PATCH 2/8] cmh update --- src/ClinicalTrialUtilities.jl | 12 +- src/ci.jl | 15 ++- src/freque.jl | 243 +++++++++++++++++++++++++++++++++- 3 files changed, 254 insertions(+), 16 deletions(-) diff --git a/src/ClinicalTrialUtilities.jl b/src/ClinicalTrialUtilities.jl index d50fcf0..a29ea37 100644 --- a/src/ClinicalTrialUtilities.jl +++ b/src/ClinicalTrialUtilities.jl @@ -19,13 +19,7 @@ using Distributions, Random, Roots, QuadGK, RecipesBase, Reexport @reexport using StatsBase import SpecialFunctions -import Base: show, findfirst -import Base.showerror -import Base.getindex -import Base.length -import Base.in -import Base.iterate -import Base.eltype +import Base: show, findfirst, getproperty, showerror, getindex, length, in, iterate, eltype import StatsBase.confint import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows! @@ -101,6 +95,7 @@ pooledcv, descriptive, freque, contab, +cmh, htmlexport, #CI confint, @@ -129,7 +124,8 @@ keldata, # Simulation ctsim, #Plots -plot +plot, +plot! diff --git a/src/ci.jl b/src/ci.jl index 6d25ea9..4f4ce0c 100644 --- a/src/ci.jl +++ b/src/ci.jl @@ -880,7 +880,7 @@ function diffmcnmwaldci(a::Int, b::Int, c::Int, d::Int; alpha = 0.05) p1 = (a + b)/n p2 = (a + c)/n estimate = p1 - p2 - se = ((a+d)*(b+c)+4*b*c)/n^3 + se = sqrt(((a+d)*(b+c)+4*b*c)/n^3) z = quantile(ZDIST, 1 - alpha / 2) ConfInt(estimate - z * se, estimate + z * se, estimate , alpha) end @@ -889,8 +889,19 @@ function diffmcnmwaldccci(a::Int, b::Int, c::Int, d::Int; alpha = 0.05) p1 = (a + b)/n p2 = (a + c)/n estimate = p1 - p2 - se = ((a+d)*(b+c)+4*b*c)/n^3 + se = sqrt(((a+d)*(b+c)+4*b*c)/n^3) z = quantile(ZDIST, 1 - alpha / 2) ConfInt(estimate - z * se - 1 / n, estimate + z * se + 1 / n, estimate , alpha) +end +#= +function diffmcnmci(a::Int, b::Int, c::Int, d::Int; alpha = 0.05) + n = a + b + c + d + p1 = (a + b)/n + p2 = (a + c)/n + estimate = p1 - p2 + se = sqrt((p1 * (1 - p1) + p2 * (1 - p2) - 2 * (a * d - b * c) / n^2) / n) + z = quantile(ZDIST, 1 - alpha / 2) + ConfInt(estimate - z * se, estimate + z * se, estimate , alpha) end +=# diff --git a/src/freque.jl b/src/freque.jl index c96b5ac..5752cc8 100644 --- a/src/freque.jl +++ b/src/freque.jl @@ -1,20 +1,60 @@ +abstract type AbstractTest end +abstract type AbstractConTab end - -struct ConTab{Int32, Int32} - tab::Matrix{Int} +struct ConTab{Int32, Int32} <: AbstractConTab + tab::Matrix{Real} row::Vector col::Vector - function ConTab(m::Matrix{Int}) + function ConTab(m::Matrix{T}) where T <: Real row = Vector{Symbol}(undef, size(m, 1)) .= Symbol("") col = Vector{Symbol}(undef, size(m, 2)) .= Symbol("") new{size(m, 1), size(m, 2)}(m, row, col) end - function ConTab(m::Matrix{Int}, row, col) + function ConTab(m::Matrix{T}, row, col) where T <: Real new{size(m, 1), size(m, 2)}(m, row, col) end end +struct McnmConTab <: AbstractConTab + tab::Matrix{Int} + function McnmConTab(m::Matrix{Int}) + new(m) + end +end + +function Base.getproperty(obj::McnmConTab, attr::Symbol) + if attr == :tab + return getfield(obj, :tab) + elseif attr == :a + return getfield(obj, :tab)[1,1] + elseif attr == :b + return getfield(obj, :tab)[1,2] + elseif attr == :c + return getfield(obj, :tab)[2,1] + elseif attr == :d + return getfield(obj, :tab)[2,2] + else + error("No field!") + end +end + +function Base.getproperty(obj::ConTab{2, 2}, attr::Symbol) + if attr == :tab + return getfield(obj, :tab) + elseif attr == :a + return getfield(obj, :tab)[1,1] + elseif attr == :b + return getfield(obj, :tab)[1,2] + elseif attr == :c + return getfield(obj, :tab)[2,1] + elseif attr == :d + return getfield(obj, :tab)[2,2] + else + error("No field!") + end +end + struct Freque val n @@ -87,9 +127,200 @@ end function mcnmtest(a::Matrix{Int}; cc = false) if cc cc = 1 else cc = 0 end (abs(a[1,2] - a[2,1]) - cc) ^ 2 / (a[1,2] + a[2,1]) - end + #= function StatsBase.confint(t::ConTab{2, 2}) end =# + +struct McnmTest + chisq + function McnmTest(ct::McnmConTab) + new(mcnmtest(ct.tab)) + end +end + +struct CMHTest <: AbstractTest + model::Symbol + type::Symbol + esti::Vector + vari::Vector + wts::Vector + est::Real + var::Real + k::Int + q::Real + tausq::Real + chisq::Real + function CMHTest(model::Symbol, t::Symbol, esti::Vector, vari::Vector, wts::Vector, est::Real, var::Real, k::Int, q::Real, tausq, chisq) + new(model, t, esti, vari, wts, est, var, k, q, tausq, chisq) + end +end + +function tabdiff(tab::McnmConTab, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + n = a + b + c + d + return (a + b)/n - (a + c)/n +end +function tabdiffvar(tab::McnmConTab, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + n = a + b + c + d + return ((a + d) * (b + c) + 4 * b * c) / n^3 +end + +function tabdiff(tab::ConTab{2,2}, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + pt = a / (a + b) + pc = c / (c + d) + return pt - pc +end +function tabdiffvar(tab::ConTab{2,2}, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + pt = a / (a + b) + pc = c / (c + d) + return pt * (1 - pt) / (a + b) + pc * (1 - pc) / (c + d) +end + +function tabor(tab::McnmConTab, adj) + b = tab.b + adj + c = tab.c + adj + return b / c +end +function taborvar(tab::McnmConTab, adj) + b = tab.b + adj + c = tab.c + adj + return 1 / b + 1 / c +end + +function tabor(tab::ConTab{2,2}, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + pt = a / (a + b) + pc = c / (c + d) + return (pt / (1 - pt)) / (pc / (1 - pc)) +end +function taborvar(tab::ConTab{2,2}, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + return 1 / a + 1 / b + 1 / c + 1 / d +end + +function tabrr(tab::McnmConTab, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + return (a + b)/(a + c) +end +function tabrrvar(tab::McnmConTab, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + return (b + c) / (a + c) / (a + b) +end + +function tabrr(tab::ConTab{2,2}, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + pt = a / (a + b) + pc = c / (c + d) + return pt / pc +end +function tabrrvar(tab::ConTab{2,2}, adj) + a = tab.a + adj + b = tab.b + adj + c = tab.c + adj + d = tab.d + adj + return 1 / a - 1 / (a + b) + 1 / c - 1 /(c + d) +end +#= +function cmhweight(tab::ConTab{2,2}) + return (tab.a + tab.b) * (tab.c + tab.d) / (tab.a + tab.b + tab.c + tab.d) + +end +=# + + +function cmh(tab::Vector{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real = 0, tau::Symbol =:dl)::CMHTest where T <: AbstractConTab + k = length(tab) + + if type == :diff + θᵢ = tabdiff.(tab, zeroadj) + vᵢ = tabdiffvar.(tab, zeroadj) + elseif type == :or + θᵢ = log.(tabor.(tab, zeroadj)) + vᵢ = taborvar.(tab, zeroadj) + elseif type == :rr + θᵢ = log.(tabrr.(tab, zeroadj)) + vᵢ = tabrrvar.(tab, zeroadj) + else + error("Type keyword unknown!") + end + + wᵢ = 1 ./ vᵢ + v = 1 / sum(wᵢ) + + + θ = sum(wᵢ .* θᵢ) / sum(wᵢ) + chisq = sum(wᵢ .* (θᵢ .^ 2)) + q = sum(wᵢ .* ((θᵢ .- θ) .^ 2)) + + + + #Tau square for random effect + + if tau == :dl + + τ² = max(0, (q - (k - 1))/ (sum(wᵢ) - sum(wᵢ .^ 2) / sum(wᵢ))) + + elseif tau == :ho + + τ² = max(0, 1 / (k - 1) * sum((θᵢ .- mean(θᵢ)) .^ 2) - 1 / k * sum(vᵢ)) + + elseif tau == :hm + + τ² = q ^ 2 / (2 * (k - 1) + q) / (sum(wᵢ) - sum(wᵢ .^ 2) / sum(wᵢ)) + + elseif tau == :ml + #in dev + lnL(m, t, k, vi, yi) = - 1 / 2 * log(2pi) - 1 / 2 * sum(log.(vi .+ t)) - 1 / 2 * sum(((yi .- m) .^ 2)/(vi .+ t)) + + else + error("Tau keyword unknown!") + #= + w = 1 / k * sum(wᵢ) + s² = 1 / (k - 1) * sum( wᵢ .^ 2 .- (k * w ^ 2)) + u = (k - 1) * (w - s² / k / w) + τ² = max(0, (q - k + 1) / u) + =# + end + + + if model == :fixed + return CMHTest(model, type, θᵢ, vᵢ, wᵢ, θ, v, k, q, τ², chisq) + end + + rwᵢ = 1 ./ (vᵢ .+ τ²) + θ = sum(rwᵢ .* θᵢ) / sum(rwᵢ) + v = 1 / sum(rwᵢ) + + return CMHTest(model, type, θᵢ, vᵢ, rwᵢ, θ, v, k, q, τ², chisq) +end From f93c5c22c751269494500d10667736939460682b Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sat, 4 Apr 2020 18:02:19 +0300 Subject: [PATCH 3/8] meta-analysis, test restructure, update readme, docs, code cosmetics, update --- README.md | 2 +- cange.log | 2 + docs/src/ds.md | 15 ++++ src/ClinicalTrialUtilities.jl | 4 +- src/dataset.jl | 31 +++++++ src/descriptives.jl | 24 ------ src/freque.jl | 125 ++++++++++++++++++++++----- test/citest.jl | 2 +- test/csv/meta.csv | 158 ++++++++++++++++++++++++++++++++++ test/ctsamplen.jl | 1 + test/dstest.jl | 2 +- test/freque.jl | 21 +++++ test/internal.jl | 63 ++++++++++++++ test/pktest.jl | 10 +-- test/test.jl | 85 ++---------------- test/testdata.jl | 2 + test/utils.jl | 2 +- 17 files changed, 414 insertions(+), 135 deletions(-) create mode 100644 test/csv/meta.csv create mode 100644 test/freque.jl create mode 100644 test/internal.jl diff --git a/README.md b/README.md index 696a86d..38fb589 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # ClinicalTrialUtilities - Clinical trial related calculation: descriptive statistics, power and sample size calculation, power simulations, confidence interval, pharmacokinetics/pharmacodynamics parameters calculation. + Clinical trial related calculation: descriptive statistics, power and sample size calculation, power simulations, confidence interval, pharmacokinetics/pharmacodynamics parameters calculation. This program comes with absolutely no warranty. No liability is accepted for any loss and risk to public health resulting from use of this software. [![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) diff --git a/cange.log b/cange.log index 7b78d59..5168ad7 100644 --- a/cange.log +++ b/cange.log @@ -6,6 +6,8 @@ v0.2.6 - code cleaning - improve coverage - plots kwargs + - Proportion meta-analysis + - Tests restruct v0.2.5 diff --git a/docs/src/ds.md b/docs/src/ds.md index 5867909..d3b574a 100644 --- a/docs/src/ds.md +++ b/docs/src/ds.md @@ -6,3 +6,18 @@ Calculation descriptive statistics by groups. ```@docs ClinicalTrialUtilities.descriptive ``` + +### freque +```@docs +ClinicalTrialUtilities.freque +``` + +### contab +```@docs +ClinicalTrialUtilities.contab +``` + +### metaprop +```@docs +ClinicalTrialUtilities.metaprop +``` diff --git a/src/ClinicalTrialUtilities.jl b/src/ClinicalTrialUtilities.jl index a29ea37..6fd25bb 100644 --- a/src/ClinicalTrialUtilities.jl +++ b/src/ClinicalTrialUtilities.jl @@ -21,7 +21,7 @@ using Distributions, Random, Roots, QuadGK, RecipesBase, Reexport import SpecialFunctions import Base: show, findfirst, getproperty, showerror, getindex, length, in, iterate, eltype import StatsBase.confint -import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows! +import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows!, rename! function lgamma(x) return SpecialFunctions.logabsgamma(x)[1] @@ -95,7 +95,7 @@ pooledcv, descriptive, freque, contab, -cmh, +metaprop, htmlexport, #CI confint, diff --git a/src/dataset.jl b/src/dataset.jl index 5a7ad7d..0832157 100644 --- a/src/dataset.jl +++ b/src/dataset.jl @@ -1,12 +1,43 @@ struct DataSet{T <: AbstractData} data::Vector{T} end +#length +function Base.length(a::DataSet)::Int + return length(a.data) +end +#Iterate function Base.iterate(iter::DataSet) return Base.iterate(iter.data) end function Base.iterate(iter::DataSet, i::Int) return Base.iterate(iter.data, i) end +#Eltype function Base.eltype(obj::DataSet) return eltype(obj.data) end +#Getindex +function Base.getindex(a::DataSet{T}, i::Int)::T where T + return a.data[i] +end +function Base.getindex(a::DataSet{T}, d::Pair)::T where T + for i = 1:length(a) + if d ∈ a[i].sort return a[i] end + end + throw(ArgumentError("Element not found!")) +end +function Base.getindex(a::DataSet{T}, d::Dict)::T where T + for i = 1:length(a) + if d ∈ a[i].sort return a[i] end + end + throw(ArgumentError("Element not found!")) +end +function Base.getindex(a::DataSet{T}, d::Tuple{Vararg{Pair}})::T where T + return a[Dict(d)] +end + + +function Base.getindex(a::DataSet{T}, i::Int, s::Symbol)::Real where T + return a.data[i].result[s] +end +# diff --git a/src/descriptives.jl b/src/descriptives.jl index 6953752..f103a56 100644 --- a/src/descriptives.jl +++ b/src/descriptives.jl @@ -39,30 +39,6 @@ function Base.show(io::IO, obj::DataSet{Descriptive}) println(io, "____________________") end end -function Base.getindex(a::DataSet{T}, i::Int)::T where T - return a.data[i] -end -function Base.getindex(a::DataSet{T}, i::Int, s::Symbol)::Real where T - return a.data[i].result[s] -end -function Base.getindex(a::DataSet{T}, d::Pair)::T where T - for i = 1:length(a) - if d ∈ a[i].sort return a[i] end - end - throw(ArgumentError("Element not found!")) -end -function Base.getindex(a::DataSet{T}, d::Dict)::T where T - for i = 1:length(a) - if d ∈ a[i].sort return a[i] end - end - throw(ArgumentError("Element not found!")) -end -function Base.getindex(a::DataSet{T}, d::Tuple{Vararg{Pair}})::T where T - return a[Dict(d)] -end -function Base.length(a::DataSet{T})::Int where T <: AbstractData - return length(a.data) -end """ descriptive(data::DataFrame; sort::Union{Symbol, Array{T,1}} = Array{Symbol,1}(undef,0), diff --git a/src/freque.jl b/src/freque.jl index 5752cc8..ce1fdeb 100644 --- a/src/freque.jl +++ b/src/freque.jl @@ -1,25 +1,33 @@ abstract type AbstractTest end -abstract type AbstractConTab end +abstract type AbstractConTab <: AbstractData end struct ConTab{Int32, Int32} <: AbstractConTab tab::Matrix{Real} row::Vector col::Vector + sort::Dict function ConTab(m::Matrix{T}) where T <: Real row = Vector{Symbol}(undef, size(m, 1)) .= Symbol("") col = Vector{Symbol}(undef, size(m, 2)) .= Symbol("") - new{size(m, 1), size(m, 2)}(m, row, col) + new{size(m, 1), size(m, 2)}(m, row, col, Dict()) end function ConTab(m::Matrix{T}, row, col) where T <: Real - new{size(m, 1), size(m, 2)}(m, row, col) + new{size(m, 1), size(m, 2)}(m, row, col, Dict()) + end + function ConTab(m::Matrix{T}, row, col, sort) where T <: Real + new{size(m, 1), size(m, 2)}(m, row, col, sort) end end struct McnmConTab <: AbstractConTab tab::Matrix{Int} + sort::Dict function McnmConTab(m::Matrix{Int}) - new(m) + new(m, Dict()) + end + function McnmConTab(m::Matrix{Int}, sort) + new(m, sort) end end @@ -60,7 +68,11 @@ struct Freque n end +""" + freque(data::DataFrame; vars::Symbol, alpha = 0.05)::DataFrame +Frequencies. +""" function freque(data::DataFrame; vars::Symbol, alpha = 0.05)::DataFrame result = DataFrame(value = Any[], n = Int[], p = Float64[], cil = Float64[], ciu = Float64[]) list = unique(data[:, vars]) @@ -73,8 +85,12 @@ function freque(data::DataFrame; vars::Symbol, alpha = 0.05)::DataFrame end return result end +""" + contab(data::DataFrame; row::Symbol, col::Symbol, sort = Dict())::ConTab -function contab(data::DataFrame; row::Symbol, col::Symbol)::ConTab +Make contingency table. +""" +function contab(data::DataFrame; row::Symbol, col::Symbol, sort = Dict())::ConTab clist = unique(data[:, col]) rlist = unique(data[:, row]) cn = length(clist) @@ -87,7 +103,31 @@ function contab(data::DataFrame; row::Symbol, col::Symbol)::ConTab dfs[ri, ci] = cnt end end - return ConTab(dfs, rlist, clist) + return ConTab(dfs, rlist, clist, sort) +end +""" + contab(data::DataFrame, sort; row::Symbol, col::Symbol) + +Make contingency tables set. +""" +function contab(data::DataFrame, sort; row::Symbol, col::Symbol) + slist = unique(data[:, sort]) + clist = unique(data[:, col]) + rlist = unique(data[:, row]) + rcv = [row, col] + result = Vector{ConTab}(undef, size(slist, 1)) + tempdf = DataFrame(row = Vector{eltype(data[!, row])}(undef, 0), col = Vector{eltype(data[!, col])}(undef, 0)) + rename!(tempdf, rcv) + for si = 1:size(slist, 1) + if size(tempdf, 1) > 0 deleterows!(tempdf, 1:size(tempdf, 1)) end + for i = 1:size(data, 1) + if data[i, sort] == slist[si, :] + push!(tempdf, data[i, rcv]) + end + end + result[si] = contab(tempdf, row = row, col = col, sort = Dict(sort .=> collect(slist[si, :]))) + end + return DataSet(result) end @@ -141,7 +181,7 @@ struct McnmTest end end -struct CMHTest <: AbstractTest +struct MetaProp <: AbstractTest model::Symbol type::Symbol esti::Vector @@ -153,7 +193,7 @@ struct CMHTest <: AbstractTest q::Real tausq::Real chisq::Real - function CMHTest(model::Symbol, t::Symbol, esti::Vector, vari::Vector, wts::Vector, est::Real, var::Real, k::Int, q::Real, tausq, chisq) + function MetaProp(model::Symbol, t::Symbol, esti::Vector, vari::Vector, wts::Vector, est::Real, var::Real, k::Int, q::Real, tausq, chisq) new(model, t, esti, vari, wts, est, var, k, q, tausq, chisq) end end @@ -258,8 +298,31 @@ function cmhweight(tab::ConTab{2,2}) end =# +""" + metaprop(tab::Vector{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real = 0, tau::Symbol =:dl)::MetaProp where T <: AbstractConTab + +Meta-analysis for 2x2 tables. + +tab: vectro of ConTab{2,2} or McnmConTab; + +type - type of measure: +- :rr +- :or/ +- :diff -function cmh(tab::Vector{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real = 0, tau::Symbol =:dl)::CMHTest where T <: AbstractConTab +model: +- :fixed +- :random + +zeroadj - zero adjustment value for all cells; + +tau - τ² calculation method: +- :dl +- :ho +- :hm + +""" +function metaprop(tab::Vector{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real = 0, tau::Symbol =:dl)::MetaProp where T <: AbstractConTab k = length(tab) if type == :diff @@ -278,31 +341,21 @@ function cmh(tab::Vector{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real wᵢ = 1 ./ vᵢ v = 1 / sum(wᵢ) - θ = sum(wᵢ .* θᵢ) / sum(wᵢ) chisq = sum(wᵢ .* (θᵢ .^ 2)) q = sum(wᵢ .* ((θᵢ .- θ) .^ 2)) - - #Tau square for random effect if tau == :dl - τ² = max(0, (q - (k - 1))/ (sum(wᵢ) - sum(wᵢ .^ 2) / sum(wᵢ))) - elseif tau == :ho - τ² = max(0, 1 / (k - 1) * sum((θᵢ .- mean(θᵢ)) .^ 2) - 1 / k * sum(vᵢ)) - elseif tau == :hm - τ² = q ^ 2 / (2 * (k - 1) + q) / (sum(wᵢ) - sum(wᵢ .^ 2) / sum(wᵢ)) - elseif tau == :ml #in dev lnL(m, t, k, vi, yi) = - 1 / 2 * log(2pi) - 1 / 2 * sum(log.(vi .+ t)) - 1 / 2 * sum(((yi .- m) .^ 2)/(vi .+ t)) - else error("Tau keyword unknown!") #= @@ -313,14 +366,40 @@ function cmh(tab::Vector{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real =# end - if model == :fixed - return CMHTest(model, type, θᵢ, vᵢ, wᵢ, θ, v, k, q, τ², chisq) + return MetaProp(model, type, θᵢ, vᵢ, wᵢ, θ, v, k, q, τ², chisq) end - rwᵢ = 1 ./ (vᵢ .+ τ²) θ = sum(rwᵢ .* θᵢ) / sum(rwᵢ) v = 1 / sum(rwᵢ) + return MetaProp(model, type, θᵢ, vᵢ, rwᵢ, θ, v, k, q, τ², chisq) +end +""" + metaprop(data::DataSet{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real = 0, tau::Symbol =:dl)::MetaProp where T <: AbstractConTab + +Meta-analysis for 2x2 tables. + +data - DataSet of ConTab{2,2} or McnmConTab; +""" +function metaprop(data::DataSet{T}; type::Symbol, model::Symbol = :fixed, zeroadj::Real = 0, tau::Symbol =:dl)::MetaProp where T <: AbstractConTab + metaprop(data.data; type = type, model = model, zeroadj = zeroadj, tau = tau) +end + + +function Base.show(io::IO, obj::MetaProp) + println(io, " Meta-analysis") + println(io, " -----------------------") + println(io, " Trial number (k): $(obj.k)") + println(io, " Model: $(obj.model)") + println(io, " Type: $(obj.type)") + if obj.type == :fixed + println(io, " Estimate: $(obj.est)") + println(io, " Variance: $(obj.var)") + else + println(io, " Estimate: $(exp(obj.est))") + println(io, " Variance: $(exp(obj.var))") + end - return CMHTest(model, type, θᵢ, vᵢ, rwᵢ, θ, v, k, q, τ², chisq) + println(io, " Q: $(obj.q)") + print(io, " τ²: $(obj.tausq)") end diff --git a/test/citest.jl b/test/citest.jl index 004c35a..557e2d1 100644 --- a/test/citest.jl +++ b/test/citest.jl @@ -1,5 +1,5 @@ println(" ---------------------------------- ") -@testset " CI Test " begin +@testset "#6 CI Test " begin # ONE PROPORTION ci = ClinicalTrialUtilities.propci(38, 100, alpha=0.05, method=:wald) @test ci.lower ≈ 0.284866005121432 atol=1E-6 diff --git a/test/csv/meta.csv b/test/csv/meta.csv new file mode 100644 index 0000000..8cbf327 --- /dev/null +++ b/test/csv/meta.csv @@ -0,0 +1,158 @@ +trial;group;result +1;1;1 +1;1;1 +1;1;1 +1;1;1 +1;1;1 +1;1;1 +1;1;1 +1;1;2 +1;1;2 +1;1;2 +1;1;2 +1;1;2 +1;1;2 +1;1;2 +1;1;2 +1;1;2 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;1 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +1;2;2 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;1 +2;1;2 +2;1;2 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;1 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +2;2;2 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;1 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;1;2 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;1 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 +3;2;2 diff --git a/test/ctsamplen.jl b/test/ctsamplen.jl index 13968e0..4a97766 100644 --- a/test/ctsamplen.jl +++ b/test/ctsamplen.jl @@ -1,3 +1,4 @@ +println(" ---------------------------------- ") @testset "#1 ctsamplen Test " begin #Sample Size Calculations in Clinical Research, Second Edition, Shein-Chung Chow, Ph.D., 2008, International Standard Book Number‑13: 978‑1‑58488‑982‑3 #1 diff --git a/test/dstest.jl b/test/dstest.jl index 027dc43..1e65545 100644 --- a/test/dstest.jl +++ b/test/dstest.jl @@ -1,5 +1,5 @@ println(" ---------------------------------- ") -@testset " descriptives " begin +@testset "#7 Descriptives " begin df = descriptivedat ds = ClinicalTrialUtilities.descriptive(df, stats = :all, sort = [:C1, :C2], vars=[:P1, :P2]) @test ds[1,:n] ≈ 20 atol=1E-5 diff --git a/test/freque.jl b/test/freque.jl new file mode 100644 index 0000000..bac3f34 --- /dev/null +++ b/test/freque.jl @@ -0,0 +1,21 @@ +println(" ---------------------------------- ") +@testset "#8 Frequency " begin + df = freqdat + ctab = ClinicalTrialUtilities.contab(df, row = :row, col = :col) + @test ctab.tab == [9 8; 5 21] + + frtab = ClinicalTrialUtilities.freque(df; vars=:row, alpha = 0.05) + @test frtab[1,2] == 17 + + + ctds = ClinicalTrialUtilities.contab(metadf, [:trial]; row = :group, col = :result) + cmht = ClinicalTrialUtilities.metaprop(ctds, type = :or, model = :fixed, zeroadj = 0.0, tau = :ho) + @test cmht.est ≈ 0.4164682333774169 atol=1E-5 + @test cmht.var ≈ 0.13107613010773247 atol=1E-5 + @test cmht.tausq ≈ 0.5554302118884364 atol=1E-5 + + cmht = ClinicalTrialUtilities.metaprop(ctds, type = :or, model = :random, zeroadj = 0.0, tau = :dl) + @test cmht.est ≈ 0.46774668790076773 atol=1E-5 + @test cmht.var ≈ 0.2500964678287197 atol=1E-5 + @test cmht.tausq ≈ 0.3284301973957149 atol=1E-5 +end diff --git a/test/internal.jl b/test/internal.jl new file mode 100644 index 0000000..17598df --- /dev/null +++ b/test/internal.jl @@ -0,0 +1,63 @@ + +@testset "Internal functions test " begin +@testset " tfn function " begin + @test ClinicalTrialUtilities.tfn(1.0,2.0) ≈ 0.07846821 atol=1E-8 + @test ClinicalTrialUtilities.tfn(0.1,10.0) > 0 #Not validated with PowerTOST + @test ClinicalTrialUtilities.tfn(0.1,10E20) > 0 #Not validated with PowerTOST +end +@testset " owensTint2 function " begin + @test round(ClinicalTrialUtilities.owenstint2(1.0, 3.0, 20.0, 3.0), digits=7) ≈ 0.4839414 +end +@testset " owensqo function " begin + @test ClinicalTrialUtilities.owensqo(1 ,2.0,1.0,1.0;a=0.0) ≈ 0.321429 atol=1E-6 + @test ClinicalTrialUtilities.owensqo(2 ,1.0,0.5,0.2;a=0.0) ≈ 0.006781741 atol=1E-9 + @test ClinicalTrialUtilities.owensqo(4 ,2.0,1.0,1.0;a=0.0) ≈ 0.03739024 atol=1E-8 + @test ClinicalTrialUtilities.owensqo(7 ,2.0,1.0,1.0;a=0.0) ≈ 0.001888241 atol=1E-9 + @test ClinicalTrialUtilities.owensqo(3 ,2.0,1.0,Inf;a=0.0) ≈ 0.7436299 atol=1E-7 +end +@testset " owensq function " begin + @test ClinicalTrialUtilities.owensq(4 ,100.0,40.0,0.0,Inf) ≈ 0.9584071 atol=1E-7 + @test ClinicalTrialUtilities.owensq(1 ,1.0,1.0,0.0,Inf) ≈ 0.42202 atol=1E-5 + @test ClinicalTrialUtilities.owensq(4 ,100.0,30.0,0.0,0.8) ≈ 0.02702275 atol=1E-8 + @test ClinicalTrialUtilities.owensq(1,100.0,40.0,0.0,1.0) ≈ 0.3718607 atol=1E-7 + @test ClinicalTrialUtilities.owensq(4 ,100.0,40.0,0.0,Inf) ≈ 0.9584071 atol=1E-7 + @test ClinicalTrialUtilities.owensq(1 ,1.0,1.0,0.0,Inf) ≈ 0.42202 atol=1E-5 + @test ClinicalTrialUtilities.owensq(4 ,100.0,30.0,0.0,0.8) ≈ 0.02702275 atol=1E-8 + @test ClinicalTrialUtilities.owensq(1,100.0,40.0,0.0,1.0) ≈ 0.3718607 atol=1E-7 +end +@testset " powerTOSTOwenQ " begin + @test ClinicalTrialUtilities.powertost_owenq(0.05,0.1,0.4,0.05,0.11,23.0) ≈ 0.00147511 atol=1E-8 +end +@testset " approxPowerTOST " begin + @test ClinicalTrialUtilities.powertost_nct(0.05,0.4,0.9,0.05,0.11,23.0) ≈ 1.076964e-06 atol=1E-12 + @test ClinicalTrialUtilities.powertost_nct(0.05,1.0,1.0,0.5,0.2,100.0) == 0 +end +@testset " approx2PowerTOST " begin + @test ClinicalTrialUtilities.powertost_shifted(0.05,0.1,1.0,0.5,0.2,1000.0) ≈ 0.4413917 atol=1E-7 +end +@testset " owenst " begin + @test ClinicalTrialUtilities.owenst(1.0,Inf) ≈ 0.07932763 atol=1E-8 + @test ClinicalTrialUtilities.owenst(-1.0,Inf) ≈ 0.07932763 atol=1E-8 + @test ClinicalTrialUtilities.owenst(1.0,-Inf) ≈ -0.07932763 atol=1E-8 + @test ClinicalTrialUtilities.owenst(-1.0,-Inf) ≈ -0.07932763 atol=1E-8 + @test ClinicalTrialUtilities.owenst(1.0, 0.5) ≈ 0.0430647 atol=1E-8 + @test ClinicalTrialUtilities.owenst(1.0,2.0) ≈ 0.07846819 atol=1E-8 + @test ClinicalTrialUtilities.owenst(Inf, 1.0) == 0 +end + +@testset " designProp " begin + d = ClinicalTrialUtilities.Design(:parallel) + @test d.df(30) ≈ 28 && d.bkni ≈ 1.0 && d.sq ≈ 2 + d = ClinicalTrialUtilities.Design(:d2x2) + @test d.df(31) ≈ 29 && d.bkni ≈ 0.5 && d.sq ≈ 2 + d = ClinicalTrialUtilities.Design(:d2x2x4) + @test d.df(31) ≈ 89 && d.bkni ≈ 0.25 && d.sq ≈ 2 + d = ClinicalTrialUtilities.Design(:d2x4x4) + @test d.df(31) ≈ 89 && d.bkni ≈ 0.0625 && d.sq ≈ 4 + d = ClinicalTrialUtilities.Design(:d2x2x3) + @test d.df(31) ≈ 59 && d.bkni ≈ 0.375 && d.sq ≈ 2 + d = ClinicalTrialUtilities.Design(:d2x3x3) + @test d.df(31) ≈ 59 && d.bkni ≈ 1/6 && d.sq ≈ 3 + +end +end diff --git a/test/pktest.jl b/test/pktest.jl index ff8f07c..80f2aa6 100644 --- a/test/pktest.jl +++ b/test/pktest.jl @@ -1,10 +1,7 @@ println(" ---------------------------------- ") -@testset " PK " begin - +@testset "#10 Pharmacokinetics " begin # Basic tests - - pkds = ClinicalTrialUtilities.pkimport(pkdata, [:Subject, :Formulation]; conc = :Concentration, time = :Time) pk = ClinicalTrialUtilities.nca!(pkds) @test pk[1, :AUCinf] ≈ 1.63205 atol=1E-5 @@ -949,7 +946,7 @@ println(" ---------------------------------- ") end println(" ---------------------------------- ") -@testset " UPK " begin +@testset "#11 Urine PK " begin upk = ClinicalTrialUtilities.upkimport(upkdata, [:subj]; stime = :st, etime = :et, conc = :conc, vol = :vol) unca = ClinicalTrialUtilities.nca!(upk[1]) unca = ClinicalTrialUtilities.nca!(upk) @@ -961,7 +958,8 @@ println(" ---------------------------------- ") end println(" ---------------------------------- ") -@testset " PD " begin +@testset "#12 Pharmacodynamics " begin + #PD pdds = ClinicalTrialUtilities.pdimport(pddata; time=:time, resp=:effect, bl = 3.0) diff --git a/test/test.jl b/test/test.jl index 68bc360..224eb36 100644 --- a/test/test.jl +++ b/test/test.jl @@ -17,6 +17,8 @@ println(" ---------------------------------- ") println(" --------- START TEST --------- ") println(" ---------------------------------- ") println(" ---------------------------------- ") +println("") +include("internal.jl") include("ctsamplen.jl") @@ -186,97 +188,28 @@ end include("utils.jl") -println(" ---------------------------------- ") -@testset "Internal functions test " begin -@testset " tfn function " begin - @test ClinicalTrialUtilities.tfn(1.0,2.0) ≈ 0.07846821 atol=1E-8 - @test ClinicalTrialUtilities.tfn(0.1,10.0) > 0 #Not validated with PowerTOST - @test ClinicalTrialUtilities.tfn(0.1,10E20) > 0 #Not validated with PowerTOST -end -@testset " owensTint2 function " begin - @test round(ClinicalTrialUtilities.owenstint2(1.0, 3.0, 20.0, 3.0), digits=7) ≈ 0.4839414 -end -@testset " owensqo function " begin - @test ClinicalTrialUtilities.owensqo(1 ,2.0,1.0,1.0;a=0.0) ≈ 0.321429 atol=1E-6 - @test ClinicalTrialUtilities.owensqo(2 ,1.0,0.5,0.2;a=0.0) ≈ 0.006781741 atol=1E-9 - @test ClinicalTrialUtilities.owensqo(4 ,2.0,1.0,1.0;a=0.0) ≈ 0.03739024 atol=1E-8 - @test ClinicalTrialUtilities.owensqo(7 ,2.0,1.0,1.0;a=0.0) ≈ 0.001888241 atol=1E-9 - @test ClinicalTrialUtilities.owensqo(3 ,2.0,1.0,Inf;a=0.0) ≈ 0.7436299 atol=1E-7 -end -@testset " owensq function " begin - @test ClinicalTrialUtilities.owensq(4 ,100.0,40.0,0.0,Inf) ≈ 0.9584071 atol=1E-7 - @test ClinicalTrialUtilities.owensq(1 ,1.0,1.0,0.0,Inf) ≈ 0.42202 atol=1E-5 - @test ClinicalTrialUtilities.owensq(4 ,100.0,30.0,0.0,0.8) ≈ 0.02702275 atol=1E-8 - @test ClinicalTrialUtilities.owensq(1,100.0,40.0,0.0,1.0) ≈ 0.3718607 atol=1E-7 - @test ClinicalTrialUtilities.owensq(4 ,100.0,40.0,0.0,Inf) ≈ 0.9584071 atol=1E-7 - @test ClinicalTrialUtilities.owensq(1 ,1.0,1.0,0.0,Inf) ≈ 0.42202 atol=1E-5 - @test ClinicalTrialUtilities.owensq(4 ,100.0,30.0,0.0,0.8) ≈ 0.02702275 atol=1E-8 - @test ClinicalTrialUtilities.owensq(1,100.0,40.0,0.0,1.0) ≈ 0.3718607 atol=1E-7 -end -@testset " powerTOSTOwenQ " begin - @test ClinicalTrialUtilities.powertost_owenq(0.05,0.1,0.4,0.05,0.11,23.0) ≈ 0.00147511 atol=1E-8 -end -@testset " approxPowerTOST " begin - @test ClinicalTrialUtilities.powertost_nct(0.05,0.4,0.9,0.05,0.11,23.0) ≈ 1.076964e-06 atol=1E-12 - @test ClinicalTrialUtilities.powertost_nct(0.05,1.0,1.0,0.5,0.2,100.0) == 0 -end -@testset " approx2PowerTOST " begin - @test ClinicalTrialUtilities.powertost_shifted(0.05,0.1,1.0,0.5,0.2,1000.0) ≈ 0.4413917 atol=1E-7 -end -@testset " owenst " begin - @test ClinicalTrialUtilities.owenst(1.0,Inf) ≈ 0.07932763 atol=1E-8 - @test ClinicalTrialUtilities.owenst(-1.0,Inf) ≈ 0.07932763 atol=1E-8 - @test ClinicalTrialUtilities.owenst(1.0,-Inf) ≈ -0.07932763 atol=1E-8 - @test ClinicalTrialUtilities.owenst(-1.0,-Inf) ≈ -0.07932763 atol=1E-8 - @test ClinicalTrialUtilities.owenst(1.0, 0.5) ≈ 0.0430647 atol=1E-8 - @test ClinicalTrialUtilities.owenst(1.0,2.0) ≈ 0.07846819 atol=1E-8 - @test ClinicalTrialUtilities.owenst(Inf, 1.0) == 0 -end - -@testset " designProp " begin - d = ClinicalTrialUtilities.Design(:parallel) - @test d.df(30) ≈ 28 && d.bkni ≈ 1.0 && d.sq ≈ 2 - d = ClinicalTrialUtilities.Design(:d2x2) - @test d.df(31) ≈ 29 && d.bkni ≈ 0.5 && d.sq ≈ 2 - d = ClinicalTrialUtilities.Design(:d2x2x4) - @test d.df(31) ≈ 89 && d.bkni ≈ 0.25 && d.sq ≈ 2 - d = ClinicalTrialUtilities.Design(:d2x4x4) - @test d.df(31) ≈ 89 && d.bkni ≈ 0.0625 && d.sq ≈ 4 - d = ClinicalTrialUtilities.Design(:d2x2x3) - @test d.df(31) ≈ 59 && d.bkni ≈ 0.375 && d.sq ≈ 2 - d = ClinicalTrialUtilities.Design(:d2x3x3) - @test d.df(31) ≈ 59 && d.bkni ≈ 1/6 && d.sq ≈ 3 - -end -end - - include("citest.jl") -include("sim.jl") - include("dstest.jl") -println(" ---------------------------------- ") -@testset " Frequency " begin - df = freqdat - ctab = ClinicalTrialUtilities.contab(df, row = :row, col = :col) - @test ctab.tab == [9 8; 5 21] +include("freque.jl") - frtab = ClinicalTrialUtilities.freque(df; vars=:row, alpha = 0.05) - @test frtab[1,2] == 17 -end println(" ---------------------------------- ") -@testset " Random " begin +@testset "#9 Random " begin @test ClinicalTrialUtilities.randomseq(seed = 1234) == [1,2,2,1,2,1,1,2,2,1] rdf = ClinicalTrialUtilities.randomtable(seed = 1234) @test rdf[!, :Group] == [1,2,2,1,2,1,1,2,2,1] end + #PK include("pktest.jl") + +include("sim.jl") + + println(" ---------------------------------- ") @testset " Scenario " begin #Pharmacokinetics statistics diff --git a/test/testdata.jl b/test/testdata.jl index 2a1ffaa..dc64136 100644 --- a/test/testdata.jl +++ b/test/testdata.jl @@ -18,3 +18,5 @@ freqdat = CSV.File(path*"/csv/freqdat.csv") |> DataFrame negdat = CSV.File(path*"/csv/negdat.csv") |> DataFrame #Dataset foe descriptive statistics check descriptivedat = CSV.File(path*"/csv/descriptivedat.csv") |> DataFrame +#Meta-analysis +metadf = CSV.File(path*"/csv/meta.csv") |> DataFrame diff --git a/test/utils.jl b/test/utils.jl index 9b2499a..989de0a 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,5 +1,5 @@ println(" ---------------------------------- ") -@testset "#5 utils Test " begin +@testset "#5 Utilities " begin @testset " ci2cv Test " begin #cvms = ClinicalTrialUtilities.cvfromci(;alpha = 0.05, theta1 = 0.9, theta2 = 1.25, n=30, design=:d2x2x4, cvms=true) From ced3a3f8d33ad58589a00701efcacebacd1c0f07 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Tue, 19 May 2020 17:42:50 +0300 Subject: [PATCH 4/8] v0.2.6 --- README.md | 3 ++- cange.log | 2 ++ docs/make.jl | 3 ++- docs/src/pkplot.md | 6 ++++++ src/ClinicalTrialUtilities.jl | 15 ++++++++------ src/dataset.jl | 38 +++++++++++++++++++++++++++++++++++ src/hypothesis.jl | 12 +++++++++++ src/means.jl | 29 ++++++++++++++++++++++++++ src/pk.jl | 17 +--------------- src/plots.jl | 33 +++++++++++++++++++++++------- src/randomization.jl | 2 +- src/samplesize.jl | 32 ++++------------------------- src/sim.jl | 11 +++++----- test/pktest.jl | 12 ++++++++--- 14 files changed, 147 insertions(+), 68 deletions(-) create mode 100644 docs/src/pkplot.md create mode 100644 src/means.jl diff --git a/README.md b/README.md index 38fb589..9054471 100644 --- a/README.md +++ b/README.md @@ -19,10 +19,11 @@ using Pkg; Pkg.add("ClinicalTrialUtilities"); ### Main features -- SampleSize calculation +- Clinical trial sample size calculation - Power calculation - Confidence Interval calculation - NCA Pharmacokinetics parameters calculation +- Descriptive statistics and frequencies - Randomization diff --git a/cange.log b/cange.log index 5168ad7..ea5d9e5 100644 --- a/cange.log +++ b/cange.log @@ -8,6 +8,8 @@ v0.2.6 - plots kwargs - Proportion meta-analysis - Tests restruct + - Hypothesis test for simulation + - Docs v0.2.5 diff --git a/docs/make.jl b/docs/make.jl index b02bfe9..8835996 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -14,7 +14,8 @@ makedocs( "Descriptive statistics" => "ds.md", "NCA" => ["NCA" => "nca.md", "Pharmakokinetics" => "pk.md", - "Pharmacodynamics" => "pd.md"], + "Pharmacodynamics" => "pd.md", + "PK Plots" => "pkplot.md"], "Randomization" => "random.md", "Simulations" => "sim.md", "Utilities" => "utils.md", diff --git a/docs/src/pkplot.md b/docs/src/pkplot.md new file mode 100644 index 0000000..adc17dd --- /dev/null +++ b/docs/src/pkplot.md @@ -0,0 +1,6 @@ +### PK plots + +#### pkplot +```@docs +ClinicalTrialUtilities.pkplot +``` diff --git a/src/ClinicalTrialUtilities.jl b/src/ClinicalTrialUtilities.jl index 6fd25bb..f0784b6 100644 --- a/src/ClinicalTrialUtilities.jl +++ b/src/ClinicalTrialUtilities.jl @@ -19,7 +19,7 @@ using Distributions, Random, Roots, QuadGK, RecipesBase, Reexport @reexport using StatsBase import SpecialFunctions -import Base: show, findfirst, getproperty, showerror, getindex, length, in, iterate, eltype +import Base: show, findfirst, getproperty, showerror, getindex, length, in, iterate, eltype, deleteat!, findall import StatsBase.confint import DataFrames: DataFrame, DataFrames, names!, unstack, deleterows!, rename! @@ -36,10 +36,15 @@ const PI2INV = 0.5 / π include("abstracttypes.jl") include("design.jl") -include("hypothesis.jl") include("proportion.jl") +include("means.jl") + +#Confidence interval calculation +include("ci.jl") + +include("hypothesis.jl") #Owen function calc: owensQ, owensQo, ifun1, owensTint2, owensT, tfn include("owensq.jl") @@ -49,8 +54,6 @@ include("powertost.jl") include("powersamplesize.jl") #Main sample size and power functions: sampleSize, ctpower, besamplen include("samplesize.jl") -#Confidence interval calculation -include("ci.jl") #Simulations include("sim.jl") #info function @@ -124,8 +127,8 @@ keldata, # Simulation ctsim, #Plots -plot, -plot! +pkplot, +pkplot! diff --git a/src/dataset.jl b/src/dataset.jl index 0832157..1ed538b 100644 --- a/src/dataset.jl +++ b/src/dataset.jl @@ -1,6 +1,22 @@ struct DataSet{T <: AbstractData} data::Vector{T} end + +#!!! +function in(a::Dict, b::Dict) + k = collect(keys(a)) + if any(x -> x ∉ collect(keys(b)), k) return false end + for i = 1:length(a) + if a[k[i]] != b[k[i]] return false end + end + return true +end +function in(a::Pair, b::Dict) + if a[1] ∉ collect(keys(b)) return false end + if a[2] != b[a[1]] return false end + return true +end + #length function Base.length(a::DataSet)::Int return length(a.data) @@ -41,3 +57,25 @@ function Base.getindex(a::DataSet{T}, i::Int, s::Symbol)::Real where T return a.data[i].result[s] end # +function Base.findall(a::DataSet{T}, sort::Dict) where T + inds = Vector{Int}(undef, 0) + for i = 1:length(a) + if sort ∈ a.data[i].sort push!(inds, i) end + end + inds +end +# +function Base.deleteat!(a::DataSet{T}, i::Int) where T + deleteat!(a.data, i) + return a +end + +function Base.deleteat!(a::DataSet{T}, inds::AbstractArray) where T + deleteat!(a.data, inds) + return a +end + +function Base.deleteat!(a::DataSet{T}, inds::Dict) where T + deleteat!(a.data, findall(a, inds)) + return a +end diff --git a/src/hypothesis.jl b/src/hypothesis.jl index adcb755..bcc5c29 100644 --- a/src/hypothesis.jl +++ b/src/hypothesis.jl @@ -45,4 +45,16 @@ function refval(h::Superiority) h.llim - h.diff end + +function checkhyp(h::Superiority, ci::ConfInt) + ci.lower > h.diff +end +function checkhyp(h::Equivalence, ci::ConfInt) + ci.lower > h.llim && ci.upper < h.ulim +end +function checkhyp(h::Equality, ci::ConfInt) + ci.lower > h.val || ci.upper < h.val +end + + struct McNemars <: AbstractHypothesis end diff --git a/src/means.jl b/src/means.jl new file mode 100644 index 0000000..8be31d5 --- /dev/null +++ b/src/means.jl @@ -0,0 +1,29 @@ +struct Mean{Bool} <: AbstractMean + m::Real + sd::Real + n::Int + function Mean(m::Real, sd::Real, n::Int) + if sd ≤ 0.0 throw(ArgumentError("SD can't be ≤ 0.0!")) end + new{true}(m, sd, n)::Mean + end + function Mean(m::Real, sd::Real) + if sd ≤ 0.0 throw(ArgumentError("SD can't be ≤ 0.0!")) end + new{false}(m, sd, 0) + end + function Mean(m::Real) + new{false}(m, NaN, 0) + end +end +""" + Mean comparison hypothesis testing + a - Test group value + b - Reference group value {Mean} type, or reference value for one group test {Real} +""" +struct DiffMean{Bool} <: AbstractCompositeMean + a::Mean + b::Mean + function DiffMean(a::Mean, b::Mean) + if typeof(a) <: Mean{false} || typeof(b) <: Mean{false} t = false else t = true end + new{t}(a, b) + end +end diff --git a/src/pk.jl b/src/pk.jl index 75f110c..3b7eaa2 100644 --- a/src/pk.jl +++ b/src/pk.jl @@ -1,22 +1,7 @@ #Pharmacokinetics #Makoid C, Vuchetich J, Banakar V. 1996-1999. Basic Pharmacokinetics. - -#!!! -function in(a::Dict, b::Dict) - k = collect(keys(a)) - if any(x -> x ∉ collect(keys(b)), k) return false end - for i = 1:length(a) - if a[k[i]] != b[k[i]] return false end - end - return true -end -function in(a::Pair, b::Dict) - if a[1] ∉ collect(keys(b)) return false end - if a[2] != b[a[1]] return false end - return true -end #!!! -struct KelData +struct KelData s::Array{Int, 1} e::Array{Int, 1} a::Array{Number, 1} diff --git a/src/plots.jl b/src/plots.jl index e3cf887..3b153c7 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -72,7 +72,16 @@ function plotlabel(d) return title end -function plot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject +""" + pkplot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject + +Plot for subject + +subj - subject; +plotstyle - styles for plots. + +""" +function pkplot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -99,7 +108,8 @@ function plot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where return p end -function plot!(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject + +function pkplot!(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -167,10 +177,10 @@ function pageplot(pagedatatmp, styledict, typesort; kwargs...) kwargs[:label] = label if fst - p = plot(d; plotstyle = style, kwargs...) + p = pkplot(d; plotstyle = style, kwargs...) fst = false else - plot!( d; plotstyle = style, kwargs...) + pkplot!( d; plotstyle = style, kwargs...) end end end @@ -200,16 +210,25 @@ function pageplot(pagedatatmp; kwargs...) kwargs[:label] = "AUTO" if fst - p = plot(d; plotstyle = PKPLOTSTYLE[1], kwargs...) + p = pkplot(d; plotstyle = PKPLOTSTYLE[1], kwargs...) fst = false else - plot!( d; plotstyle = PKPLOTSTYLE[1], kwargs...) + pkplot!( d; plotstyle = PKPLOTSTYLE[1], kwargs...) end end p end -function plot(data::DataSet{T}; pagesort = nothing, typesort = nothing, kwargs...) where T <: AbstractSubject +""" + pkplot(data::DataSet{T}; pagesort::Union{Nothing, Vector{Symbol}} = nothing, typesort::Union{Nothing, Vector{Symbol}} = nothing, kwargs...) where T <: AbstractSubject + +Plot for subjects in dataset. + +data - subjects dataset; +pagesort - subject page groupping; +typesort - subject sorting within page; +""" +function pkplot(data::DataSet{T}; pagesort::Union{Nothing, Vector{Symbol}} = nothing, typesort::Union{Nothing, Vector{Symbol}} = nothing, kwargs...) where T <: AbstractSubject # Style types kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) diff --git a/src/randomization.jl b/src/randomization.jl index 012b091..2a7d607 100644 --- a/src/randomization.jl +++ b/src/randomization.jl @@ -74,7 +74,7 @@ function randomtable(;blocksize::Int = 4, subject::Int = 10, group::Int = 2, rat seqrand = permutedims(hcat(grseq[r]...)) seqrand = hcat(collect(1:length(r)), hcat(r, seqrand)) df = DataFrame(seqrand) - names!(df, [:Subject, :Group, :Sequence]) + rename!(df, [:Subject, :Group, :Sequence]) for i=1:length(grseq[1]) df[!, Symbol("Period_"*string(i))] = getindex.(df[!, :Sequence], i) end diff --git a/src/samplesize.jl b/src/samplesize.jl index 05ae138..6f3a761 100644 --- a/src/samplesize.jl +++ b/src/samplesize.jl @@ -16,35 +16,11 @@ struct SampleSize{T} <: AbstractSampleSize where T <: AbstractFloat end end -struct Mean{Bool} <: AbstractMean - m::Real - sd::Real - n::Int - function Mean(m::Real, sd::Real, n::Int) - if sd ≤ 0.0 throw(ArgumentError("SD can't be ≤ 0.0!")) end - new{true}(m, sd, n)::Mean - end - function Mean(m::Real, sd::Real) - if sd ≤ 0.0 throw(ArgumentError("SD can't be ≤ 0.0!")) end - new{false}(m, sd, 0) - end - function Mean(m::Real) - new{false}(m, NaN, 0) - end -end -""" - Mean comparison hypothesis testing - a - Test group value - b - Reference group value {Mean} type, or reference value for one group test {Real} -""" -struct DiffMean{Bool} <: AbstractCompositeMean - a::Mean - b::Mean - function DiffMean(a::Mean, b::Mean) - if typeof(a) <: Mean{false} || typeof(b) <: Mean{false} t = false else t = true end - new{t}(a, b) - end +function getval(o::T) where T <: AbstractObjective + o.val end + + struct TaskResult{T <: AbstractTask} task::T method::Symbol diff --git a/src/sim.jl b/src/sim.jl index 287dd88..5623d5f 100644 --- a/src/sim.jl +++ b/src/sim.jl @@ -43,15 +43,16 @@ Proportion difference power simulation. """ function ctsim(t::CTask{DiffProportion{P, P}, D, H, Power}; nsim = 100, seed=0, method = :default, dropout = 0.0) where P <: Proportion where D where H <: Superiority if seed != 0 Random.seed!(seed) end - + n₁ = getval(t.objective) + n₂ = Int(ceil(getval(t.objective) / t.k)) pow = 0 - D₁ = Binomial(t.param.a.n, getval(t.param.a)) - D₂ = Binomial(t.param.b.n, getval(t.param.b)) + D₁ = Binomial(n₁, getval(t.param.a)) + D₂ = Binomial(n₂, getval(t.param.b)) for i = 1:nsim n1 = rand(D₁) n2 = rand(D₂) - ci = confint(DiffProportion(Proportion(n1, t.param.a.n), Proportion(n2, t.param.b.n)); level = 1.0 - t.alpha * 2, method = method) - if ci.lower > t.hyp.diff pow += 1 end + ci = confint(DiffProportion(Proportion(n1, getval(t.objective)), Proportion(n2, getval(t.objective))); level = 1.0 - t.alpha * 2, method = method) + if checkhyp(t.hyp, ci) pow += 1 end end return pow/nsim end diff --git a/test/pktest.jl b/test/pktest.jl index 80f2aa6..6c1b757 100644 --- a/test/pktest.jl +++ b/test/pktest.jl @@ -921,10 +921,16 @@ println(" ---------------------------------- ") pkds = ClinicalTrialUtilities.pkimport(glucose2, [:Subject, :Date]; conc = :glucose, time = :Time) pk = ClinicalTrialUtilities.nca!(pkds) df = DataFrame(pk; unst = true) - p = ClinicalTrialUtilities.plot(pkds; pagesort = [:Date], typesort = [:Subject]) + p = ClinicalTrialUtilities.pkplot(pkds; pagesort = [:Date], typesort = [:Subject]) @test length(p) == 2 print(io, pk[1].subject.keldata) - + p = ClinicalTrialUtilities.pkplot(pkds; pagesort = nothing, typesort = [:Subject]) + @test length(p) == 1 + p = ClinicalTrialUtilities.pkplot(pkds; pagesort = nothing, typesort = nothing, legend = false, xlabel = "L1", ylabel = "L2", xlims = (0,12)) + @test length(p) == 1 + p = ClinicalTrialUtilities.pkplot(pkds; pagesort = [:Date], typesort = nothing) + @test length(p) == 2 + p = ClinicalTrialUtilities.pkplot(pkds[1]) # NaN PK LimitRule test nanpkdata.Concentration = ClinicalTrialUtilities.tryfloatparse!(nanpkdata.Concentration) @@ -954,7 +960,7 @@ println(" ---------------------------------- ") @test unca[1][:mTmax] ≈ 1.5 @test unca[1][:ar] ≈ 16.0 @test unca[1][:volume] ≈ 11.0 - p = ClinicalTrialUtilities.plot(upk, ylabel = "Excretion") + p = ClinicalTrialUtilities.pkplot(upk, ylabel = "Excretion") end println(" ---------------------------------- ") From c22a829f11bf83276edc761b2f76afc0ea70d67d Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 29 May 2020 00:17:38 +0300 Subject: [PATCH 5/8] update julia, CSV, RecipeBase --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 53f13a1..62eb2e1 100644 --- a/Project.toml +++ b/Project.toml @@ -24,13 +24,13 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" test = ["CSV", "Test", "Plots"] [compat] -julia = "1.0, 1.1, 1.2, 1.3" +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" SpecialFunctions = "0.8, 0.9, 0.10" Roots = "0.7, 0.8, 1.0" -CSV = "0.5" -RecipesBase = "0.7, 0.8" +CSV = "0.5, 0.6" +RecipesBase = "0.7, 0.8, 1.0" Reexport = "0.1, 0.2" From eec921c8ede0d3306ed172919b7619407849996d Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 29 May 2020 00:17:56 +0300 Subject: [PATCH 6/8] plots fix --- src/plots.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/plots.jl b/src/plots.jl index 3b153c7..a5b9495 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -45,19 +45,19 @@ end @recipe function f(subj::PKPlot) x, y = subj.args - seriestype --> :line - xlabel --> "Time" - link --> :both - legend --> true - grid --> false + seriestype --> :line + xguide --> "Time" + link --> :both + legend --> true + grid --> false #ticks := [nothing :auto nothing] - xlims --> (minimum(x), maximum(x)), - ylims --> (0, maximum(y)*1.1) - seriescolor --> :blue - markershape --> :circle - markersize --> 3 - markercolor --> :match - msalpha --> 0 + xlims --> (minimum(x), maximum(x)), + ylims --> (0, maximum(y)*1.1) + seriescolor --> :blue + markershape --> :circle + markersize --> 3 + markercolor --> :match + markerstrokealpha --> 0 (x, y) end From 647f9f00beed6715ecc3340a009866826846550c Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 29 May 2020 00:18:05 +0300 Subject: [PATCH 7/8] update CI --- README.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/README.md b/README.md index 9054471..ec25260 100644 --- a/README.md +++ b/README.md @@ -102,6 +102,18 @@ pooledcv([0.12, 0.2, 0.25], [14, 22, 32], [:d2x2, :d2x2, :d2x2]) ``` +#### Confidence Intervals +``` +using ClinicalTrialUtilities +ci = propci(38, 100, alpha=0.05, method=:cp) + +ci = orpropci(30, 100, 40, 90; alpha=0.05, method=:mn) + +ci = diffpropci(30, 100, 40, 90; alpha=0.05, method=:wald) + +ci = meanci(30, 10, 30, alpha = 0.05, method=:norm) +``` + #### NCA ``` using CSV, DataFrames, ClinicalTrialUtilities @@ -114,6 +126,12 @@ dsdf = ClinicalTrialUtilities.DataFrame(ds; unst = true) ``` +#### Randomization +``` +using DataFrames, ClinicalTrialUtilities +rt = ClinicalTrialUtilities.randomtable(;blocksize = 4, subject = 32, group = 2, ratio = [1,1], grseq = ["TR", "RT"], seed = 36434654652452) +``` + ### Copyrights From b192ebef04bca64b87bfb5cce46d63855fa6d58a Mon Sep 17 00:00:00 2001 From: PharmCat Date: Fri, 29 May 2020 00:24:21 +0300 Subject: [PATCH 8/8] update --- cange.log | 1 + 1 file changed, 1 insertion(+) diff --git a/cange.log b/cange.log index ea5d9e5..5d740ad 100644 --- a/cange.log +++ b/cange.log @@ -9,6 +9,7 @@ v0.2.6 - Proportion meta-analysis - Tests restruct - Hypothesis test for simulation + - Deps update - Docs v0.2.5