From 916abee8a0bde1a762bc0c8b871696e7a86bd33b Mon Sep 17 00:00:00 2001 From: PharmCat Date: Sat, 22 May 2021 21:39:03 +0300 Subject: [PATCH] v0.6.0 --- Project.toml | 2 +- change.log | 14 ++- docs/src/nca.md | 66 +++++++++++ docs/src/ncatable.jl | 1 + src/dataset.jl | 5 + src/descriptives.jl | 29 +++-- src/export.jl | 2 +- src/pk.jl | 226 +++++++++++++++++++++++++------------ src/plots.jl | 103 +++++++++++++---- test/pktest.jl | 175 ++++++++++++++++++++++++++++ test/validation_report.jmd | 142 +++++++++++++++++++++++ test/weave.jl | 5 + 12 files changed, 661 insertions(+), 109 deletions(-) create mode 100644 docs/src/ncatable.jl create mode 100644 test/validation_report.jmd create mode 100644 test/weave.jl diff --git a/Project.toml b/Project.toml index dbf15ec..f1dd7a8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClinicalTrialUtilities" uuid = "535c2557-d7d0-564d-8ff9-4ae146c18cfe" authors = ["Vladimir Arnautov (mail@pharmcat.net)"] -version = "0.5.1" +version = "0.6.0" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" diff --git a/change.log b/change.log index d394ba3..8fd6a4b 100644 --- a/change.log +++ b/change.log @@ -1,9 +1,21 @@ +v0.6.0 + - changes in PK parameters names + - nca! administration settings + - adjust kel range for nca! with adm = :iv + - MRTtauinf used for Vsstau calculation + - changes in plots + - descriptive, CI add + - export fix + - DataSet sort! + - documentation + - nca tests + v0.5.1 - DataFrames 1 - SpecialFunctions 1 v0.5.0 - - StableRNGs test + - StableRNGs test+ - rng keyword for randomtable, randomseq - auc_sparse for PKSubject diff --git a/docs/src/nca.md b/docs/src/nca.md index eb30851..2f30eef 100644 --- a/docs/src/nca.md +++ b/docs/src/nca.md @@ -8,6 +8,72 @@ NCA analysis based on following steps: 4. Exporting to DataFrame; 5. Descriptive statistics / HTML export. +### Description + +```julia +include("ncatable.jl") +``` +#### AUC + +```math +AUC = \sum_{n=1}^N AUC_{n} +``` + +Where AUCn - partial AUC. + +Linear trapezoidal rule + +```math +AUC\mid_{t_1}^{t_2} = \delta t \times \frac{C_1 + C_2}{2} + +AUMC\mid_{t_1}^{t_2} = \delta t \times \frac{t_1 \times C_1 + t_2 \times C_2}{2} +``` + +Logarithmic trapezoidal rule + +```math +AUC\mid_{t_1}^{t_2} = \delta t \times \frac{ C_2 - C_1}{ln(C_2/C_1)} + +AUMC\mid_{t_1}^{t_2} = \delta t \times \frac{t_2 \times C_2 - t_1 \times C_1}{ln(C_2/C_1)} - \delta t^2 \times \frac{ C_2 - C_1}{ln(C_2/C_1)^2} +``` + +Linear interpolation rule + +```math +C_x = C_1 + \frac{(t_x-t_1)\times(C_2 - C_1)}{t_2 - t_1} +``` + +Logarithmic interpolation rule + +```math +C_x = exp\left(ln(C_1) + \frac{(t_x-t_1)\times(ln(C_2) - ln(C_1))}{t_2 - t_1}\right) +``` + +#### \lambda_z - elimination constant + +#### HL + +```math +HL = ln(2) / \lambda_z +``` + +#### AUCinf + +```math +AUCinf = AUClast + \frac{Clast}{\lambda_z} +``` + +#### AUMCinf + +```math +AUMCinf = AUMClast + \frac{tlast\times Clast}{\lambda_z} + \frac{Clast}{\lambda_z^2} +``` + +#### Accumulation index + +```math +Accind = \frac{1}{1 - exp(-\lambda_z \tau)} +``` ### nca! ```@docs diff --git a/docs/src/ncatable.jl b/docs/src/ncatable.jl new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/docs/src/ncatable.jl @@ -0,0 +1 @@ + diff --git a/src/dataset.jl b/src/dataset.jl index 006ca5d..fd60dd3 100644 --- a/src/dataset.jl +++ b/src/dataset.jl @@ -80,6 +80,11 @@ function Base.deleteat!(a::DataSet{T}, inds::Dict) where T return a end +function Base.sort!(a::DataSet; by) + sort!(a.data; by = by) + a +end + ################################################################################ #= function Tables.istable(table::DataSet) diff --git a/src/descriptives.jl b/src/descriptives.jl index bc39f52..e85350b 100644 --- a/src/descriptives.jl +++ b/src/descriptives.jl @@ -128,10 +128,12 @@ Descriptive statistics. function descriptive(data; sort::Union{Symbol, Array{T,1}} = Array{Symbol,1}(undef,0), vars = [], - stats::Union{Symbol, Array{T,1}, Tuple{Vararg{Symbol}}} = :default)::DataSet{Descriptive} where T <: Union{Symbol, String} + stats::Union{Symbol, Array{T,1}, Tuple{Vararg{Symbol}}} = :default, level = 0.95)::DataSet{Descriptive} where T <: Union{Symbol, String} stats = checkstats(stats) - + if isa(vars, UnitRange{Int64}) + vars = names(data)[vars] + end if isa(vars, Array) && length(vars) == 0 vars = filter(x -> x ∉ sort, names(data)) del = [] @@ -158,7 +160,7 @@ function descriptive(data; #pushvardescriptive!(d, vars, Matrix(data[:,vars]), Dict(), stats) mx = Matrix(data[:,vars]) for v = 1:length(vars) #For each variable in list - push!(d, Descriptive(vars[v], nothing, Dict(:Variable => v), descriptive_(mx[:, v], stats))) + push!(d, Descriptive(vars[v], nothing, Dict(:Variable => v), descriptive_(mx[:, v], stats, level))) end return DataSet(d) end @@ -173,15 +175,15 @@ function descriptive(data; for v = 1:length(vars) #For each variable in list dict = Dict{Symbol, Any}(sort .=> sortval) dict[:Variable] = vars[v] - push!(d, Descriptive(vars[v], nothing, dict, descriptive_(mx[:, v], stats))) + push!(d, Descriptive(vars[v], nothing, dict, descriptive_(mx[:, v], stats, level))) end #pushvardescriptive!(d, vars, mx, sortval, stats) #push variable descriprives for mx end return DataSet(d) end -function descriptive(data::Array{T, 1}; stats::Union{Symbol, Vector, Tuple} = :default, var = nothing, varname = nothing, sort = Dict())::Descriptive where T <: Real +function descriptive(data::Array{T, 1}; stats::Union{Symbol, Vector, Tuple} = :default, var = nothing, varname = nothing, sort = Dict(), level = 0.95)::Descriptive where T <: Real stats = checkstats(stats) - return Descriptive(var, varname, sort, descriptive_(data, stats)) + return Descriptive(var, varname, sort, descriptive_(data, stats, level)) end @@ -189,14 +191,16 @@ end Check if all statistics in allstat list. return stats tuple """ @inline function checkstats(stats::Union{Symbol, Array{T,1}, Tuple{Vararg{Symbol}}})::Tuple{Vararg{Symbol}} where T <: Union{Symbol, String} - allstat = (:n, :min, :max, :range, :mean, :var, :sd, :sem, :cv, :harmmean, :geomean, :geovar, :geosd, :geocv, :skew, :ses, :kurt, :sek, :uq, :median, :lq, :iqr, :mode) + allstat = (:n, :min, :max, :range, :mean, :var, :sd, :sem, :cv, :harmmean, :geomean, :geovar, :geosd, :geocv, :skew, :ses, :kurt, :sek, :uq, :median, :lq, :iqr, :mode, :meanci) if isa(stats, Symbol) if stats == :default stats = (:n, :mean, :sd, :sem, :uq, :median, :lq) elseif stats == :all stats = allstat else stats = Tuple(stats) end end stats = Tuple(Symbol.(stats)) - if any(x -> x ∉ allstat, stats) throw(ArgumentError("stats element not in allstat list")) end + if any(x -> x ∉ allstat, stats) + throw(ArgumentError("stats element ($(findall(x -> x ∉ allstat, stats))) not in allstat list")) + end return stats end """ @@ -237,7 +241,7 @@ function notnan(x) return !(x === NaN || x === nothing || x === missing) end -@inline function descriptive_(data::Vector{T}, stats::Union{Tuple{Vararg{Symbol}}, Array{Symbol,1}}) where T <: Real +@inline function descriptive_(data::Vector{T}, stats::Union{Tuple{Vararg{Symbol}}, Array{Symbol,1}}, level) where T <: Real #= dlist = findall(x -> x === NaN || x === nothing || x === missing, data) @@ -395,6 +399,13 @@ end dict[s] = abs(duq-dlq) elseif s == :mode dict[s] = mode(data) + elseif s == :meanci + if dmean === nothing dmean = mean(data) end + if dvar === nothing dvar = var(data) end + if dn === nothing dn = length(data) end + e = sqrt(dvar/dn)*quantile(TDist(dn-1), 1-(1-level)/2) + dict[Symbol(string(s)*"L"*string(level))] = dmean-e + dict[Symbol(string(s)*"U"*string(level))] = dmean+e end end return dict diff --git a/src/export.jl b/src/export.jl index bd95ff9..d07d894 100644 --- a/src/export.jl +++ b/src/export.jl @@ -196,7 +196,7 @@ HTLM export. while s s = false for r = 2:rown - if tablematrix[r,c] !=0 && data[r,c] === data[r-1,c] + if tablematrix[r,c] !=0 && !ismissing(data[r,c]) && !ismissing(data[r-1,c]) && data[r,c] == data[r-1,c] tablematrix[r,c] -= 1; tablematrix[r-1,c] += 1; s = true; diff --git a/src/pk.jl b/src/pk.jl index 33e561d..92ad3d5 100644 --- a/src/pk.jl +++ b/src/pk.jl @@ -220,8 +220,8 @@ function Base.show(io::IO, obj::DataSet{PKPDProfile}) end end function Base.show(io::IO, obj::PKSubject) - println(io, "Pharmacokinetic subject") - println(io, "Observations: $(length(obj))") + println(io, " Pharmacokinetic subject") + println(io, "Observations: $(length(obj)); ", obj.dosetime) println(io, obj.kelrange) println(io, "Time Concentration") for i = 1:length(obj) @@ -249,8 +249,8 @@ end function Base.show(io::IO, obj::ElimRange) print(io, "Elimination range: $(obj.kelstart) - $(obj.kelend) ") if length(obj.kelexcl) > 0 - print(io, "Exclusions: $(obj.kelexcl[i])") - if length(obj.kelexcl) > 1 for i = 1:length(obj.kelexcl) print(io, ", $(obj.kelexcl[i])") end end + print(io, "Exclusions: $(obj.kelexcl[1])") + if length(obj.kelexcl) > 1 for i = 2:length(obj.kelexcl) print(io, ", $(obj.kelexcl[i])") end end print(io, ".") else print(io, "No exclusion.") @@ -316,7 +316,7 @@ end function dropbeforedosetime(time::Vector, obs::Vector, dosetime::DoseTime) s = 0 for i = 1:length(time) - if time[i] < dosetime.time s = i else break end + @inbounds if time[i] < dosetime.time s = i else break end end return time[s+1:end], obs[s+1:end] end @@ -382,7 +382,7 @@ end #Intrapolation #linear prediction bx from ax, a1 < ax < a2 function linpredict(a₁, a₂, ax, b₁, b₂) - return abs((ax - a₁) / (a₂ - a₁))*(b₂ - b₁) + b₁ + return (ax - a₁) / (a₂ - a₁)*(b₂ - b₁) + b₁ end function logtpredict(c₁, c₂, cx, t₁, t₂) @@ -390,7 +390,7 @@ end end function logcpredict(t₁, t₂, tx, c₁, c₂) - return exp(log(c₁) + abs((tx-t₁)/(t₂-t₁))*(log(c₂) - log(c₁))) + return exp(log(c₁) + (tx-t₁)/(t₂-t₁)*(log(c₂) - log(c₁))) end function cpredict(t₁, t₂, tx, c₁, c₂, calcm) if calcm == :lint || c₂ >= c₁ @@ -445,30 +445,38 @@ end #------------------------------------------------------------------------------- """ - nca!(data::PKSubject; calcm = :lint, verbose = false, io::IO = stdout) + nca!(data::PKSubject; adm = :ev, calcm = :lint, intp = :lint, + verbose = false, warn = true, io::IO = stdout) Pharmacokinetics non-compartment analysis for one PK subject. -calcm - calculation method; +`adm` - administration: + +* `:ev` - extravascular; +* `:iv` - intra vascular bolus; + +`calcm` - calculation method; - :lint - Linear trapezoidal everywhere; - :logt - Log-trapezoidat rule after Tmax if c₁ > 0 and c₂ > 0, else Linear trapezoidal used; - :luld - Linear Up - Log Down everywhere if c₁ > c₂ > 0, else Linear trapezoidal used; - :luldt - Linear Up - Log Down after Tmax if c₁ > c₂ > 0, else Linear trapezoidal used; -intp - interpolation rule; +`intp` - interpolation rule; - :lint - linear interpolation; - :logt - log interpolation; -verbose - print to out stream if "true"; +`verbose` - print to out stream if "true"; - true; - false. -io - out stream (stdout by default) +`warn` - warnings enabled; + +`io` - out stream (stdout by default). """ -function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io::IO = stdout) +function nca!(data::PKSubject; adm = :ev, calcm = :lint, intp = :lint, verbose = false, warn = true, io::IO = stdout) result = Dict() #= result = Dict(:Obsnum => 0, :Tmax => 0, :Tmaxn => 0, :Tlast => 0, @@ -479,7 +487,11 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: :Rsq => NaN, :ARsq => NaN, :Rsqn => NaN, :AUCinf => NaN, :AUCpct => NaN) =# - + if warn + if !allunique(data.time) + @warn "Not all time values is unique!" + end + end #if calcm == :logt / :luld / :luldt calculation method used - log interpolation using if calcm == :logt || calcm == :luld || calcm == :luldt intp = :logt end @@ -487,7 +499,7 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: #obs = data.obs time, obs = dropbeforedosetime(data.time, data.obs, data.dosetime) result[:Obsnum] = length(obs) - + rmn = length(data.obs) - length(obs) if result[:Obsnum] < 2 return PKPDProfile(data, result; method = calcm) end @@ -534,16 +546,61 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: result[:Ctaumin] = minimum(view(obs, ncas:ncae)) end + #Elimination + keldata = KelData() + if data.kelauto + if result[:Obsnum] - result[:Tmaxn] >= 3 + logconc = log.(obs) + cmask = trues(result[:Obsnum]) + kelexclrmn = data.kelrange.kelexcl .- rmn + view(cmask, kelexclrmn) .= false + kellast = findlast(x -> x, cmask) + if kellast - result[:Tmaxn] >= 3 + if adm == :iv tmnc = 0 else tmnc = 1 end + for i = result[:Tmaxn] + tmnc:kellast - 2 + cmask .= false + #cmask[i:result[:Obsnum]] .= true + view(cmask, i:kellast) .= true + if length(data.kelrange.kelexcl) > 0 + view(cmask, kelexclrmn) .= false + end + if sum(cmask) >= 3 + sl = slope(time[cmask], logconc[cmask]) + if sl[1] < 0 + push!(keldata, i, kellast + rmn, sl[1], sl[2], sl[3], sl[4]) + end + end + end + end + end + else + logconc = log.(obs) + #cmask = Array{Bool, 1}(undef, result[:Obsnum]) + #cmask .= false + cmask = falses(result[:Obsnum]) + cmask[data.kelrange.kelstart-rmn:data.kelrange.kelend-rmn] .= true + cmask[data.kelrange.kelexcl .- rmn] .= false + sl = slope(view(time, cmask), view(logconc, cmask)) + push!(keldata, data.kelrange.kelstart, data.kelrange.kelend, sl[1], sl[2], sl[3], sl[4]) + end + #Calcalation partial areas (doseaucpart, doseaumcpart) #Dosetime < first time if data.dosetime.time < time[1] - # IV not included!!! - if data.dosetime.tau > 0 - result[:Cdose] = result[:Ctaumin] + if adm == :iv + if obs[1] > obs[2] > 0 + result[:Cdose] = logcpredict(time[1], time[2], data.dosetime.time, obs[1], obs[2]) + else + result[:Cdose] = obs[findfirst(x->x>0, obs)] + end else - result[:Cdose] = 0 + if data.dosetime.tau > 0 + result[:Cdose] = result[:Ctaumin] + else + result[:Cdose] = 0 + end end - doseaucpart, doseaumcpart = aucpart(data.dosetime.time, time[1], result[:Cdose], obs[1], :lint, false) + doseaucpart, doseaumcpart = aucpart(data.dosetime.time, time[1], result[:Cdose], obs[1], :lint, false) #always :lint? elseif data.dosetime.time == time[1] result[:Cdose] = obs[1] else @@ -569,33 +626,11 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: result[:AUMClast] = sum(view(aumcpartl, pmask)) + doseaumcpart result[:MRTlast] = result[:AUMClast] / result[:AUClast] + + result[:Cllast] = data.dosetime.dose / result[:AUClast] #----------------------------------------------------------------------- #----------------------------------------------------------------------- - #Elimination - keldata = KelData() - if data.kelauto - if result[:Obsnum] - result[:Tmaxn] >= 3 - logconc = log.(obs) - cmask = Vector{Bool}(undef, result[:Obsnum]) - for i = result[:Tmaxn] + 1:result[:Obsnum] - 2 - cmask .= false - cmask[i:result[:Obsnum]] .= true - sl = slope(time[cmask], logconc[cmask]) - if sl[1] < 0 * sl[1] - push!(keldata, i, result[:Obsnum], sl[1], sl[2], sl[3], sl[4]) - end - end - end - else - logconc = log.(obs) - #cmask = Array{Bool, 1}(undef, result[:Obsnum]) - #cmask .= false - cmask = falses(result[:Obsnum]) - cmask[data.kelrange.kelstart:data.kelrange.kelend] .= true - cmask[data.kelrange.kelexcl] .= false - sl = slope(view(time, cmask), view(logconc, cmask)) - push!(keldata, data.kelrange.kelstart, data.kelrange.kelend, sl[1], sl[2], sl[3], sl[4]) - end + if length(keldata) > 0 result[:Rsq], result[:Rsqn] = findmax(keldata.ar) data.kelrange.kelstart = keldata.s[result[:Rsqn]] @@ -614,8 +649,11 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: result[:AUCpct] = (result[:AUCinf] - result[:AUClast]) / result[:AUCinf] * 100.0 result[:AUMCinf] = result[:AUMClast] + result[:Tlast] * result[:Clast] / result[:Kel] + result[:Clast] / result[:Kel] ^ 2 result[:MRTinf] = result[:AUMCinf] / result[:AUCinf] - result[:Vzf] = data.dosetime.dose / result[:AUCinf] / result[:Kel] - result[:Clf] = data.dosetime.dose / result[:AUCinf] + result[:Vzlast] = data.dosetime.dose / result[:AUClast] / result[:Kel] + result[:Vzinf] = data.dosetime.dose / result[:AUCinf] / result[:Kel] + result[:Clinf] = data.dosetime.dose / result[:AUCinf] + result[:Vssinf] = result[:Clinf] * result[:MRTinf] + else result[:Kel] = NaN result[:HL] = NaN @@ -661,6 +699,10 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: else result[:Accind] = NaN end + result[:MRTtauinf] = (result[:AUMCtau] + data.dosetime.tau * (result[:AUCinf] - result[:AUCtau])) / result[:AUCtau] + result[:Vztau] = data.dosetime.dose / result[:AUCtau] / result[:Kel] + result[:Cltau] = data.dosetime.dose / result[:AUCtau] + result[:Vsstau] = result[:Cltau] * result[:MRTtauinf] end if verbose aucpartlsum = similar(aucpartl) @@ -696,33 +738,15 @@ function nca!(data::PKSubject; calcm = :lint, intp = :lint, verbose = false, io: return PKPDProfile(data, result; method = calcm) end """ - nca!(data::DataSet{PKSubject}; calcm = :lint, intp = :lint, - verbose = false, io::IO = stdout) - - -Pharmacokinetics non-compartment analysis for PK subjects set. - -calcm - calculation method; - -- :lint - Linear trapezoidal everywhere; -- :logt - Log-trapezoidat rule after Tmax if c₁ > c₂ > 0, else Linear trapezoidal used (same as :luldt - will be deprecated); -- :luld - Linear Up - Log Down everywhere if c₁ > c₂ > 0, else Linear trapezoidal used; - -intp - interpolation rule; -- :lint - linear interpolation; -- :logt - log interpolation; - -verbose - print to out stream if "true"; + nca!(data::DataSet{PKSubject}; adm = :ev, calcm = :lint, intp = :lint, + verbose = false, warn = true, io::IO = stdout) -- true; -- false. - -io - out stream (stdout by default) +Pharmacokinetics non-compartment analysis for PK subjects DataSet. """ -function nca!(data::DataSet{PKSubject}; calcm = :lint, intp = :lint, verbose = false, io::IO = stdout) +function nca!(data::DataSet{PKSubject}; adm = :ev, calcm = :lint, intp = :lint, verbose = false, warn = true, io::IO = stdout) results = Array{PKPDProfile, 1}(undef, 0) for i = 1:length(data.data) - push!(results, nca!(data.data[i]; calcm = calcm, intp = intp, verbose = verbose, io = io)) + push!(results, nca!(data.data[i]; adm = adm, calcm = calcm, intp = intp, verbose = verbose, warn = warn, io = io)) end return DataSet(results) end @@ -993,7 +1017,7 @@ Pharmacodynamics data import. - th - treashold. """ function pdimport(data, sort::Array; resp::Symbol, time::Symbol, bl::Real = 0, th::Real = NaN)::DataSet - №sortlist = unique(data[:, sort]) + sortlist = unique(data[:, sort]) results = Array{PDSubject, 1}(undef, 0) getdatai(data, sort, [resp, time], (x, y) -> push!(results, PDSubject(x[!, time], x[!, resp], bl, th, sort = Dict(sort .=> collect(y)))); sortby = time) return DataSet(results) @@ -1067,29 +1091,71 @@ Set range for elimination parameters calculation for subject. data - PK subject; range - ElimRange object. + +Set kelauto `false` if kelend > kelstart > 0. """ function setelimrange!(data::PKSubject, range::ElimRange) if range.kelend > length(data) throw(ArgumentError("Kel endpoint out of range")) end - setkelauto!(data, false) + if range.kelend > range.kelstart > 0 setkelauto!(data, false) end data.kelrange = range end +""" + setelimrange!(data::DataSet{PKSubject}, range::ElimRange) + +Set range for elimination parameters calculation for DataSet. +data - DataSet of PK subject; +range - ElimRange object. +""" function setelimrange!(data::DataSet{PKSubject}, range::ElimRange) for i = 1:length(data) setelimrange!(data[i], range) end data end +""" + setelimrange!(data::DataSet{PKSubject}, range::ElimRange, subj::Array{Int,1}) + +Set range for elimination parameters calculation for DataSet. + +data - DataSet of PK subject; +range - ElimRange object; +subj - subjects. +""" function setelimrange!(data::DataSet{PKSubject}, range::ElimRange, subj::Array{Int,1}) for i in subj setelimrange!(data[i], range) end data end +""" + setelimrange!(data::DataSet{PKSubject}, range::ElimRange, subj::Int) + +Set range for elimination parameters calculation for DataSet. + +data - DataSet of PK subject; +range - ElimRange object; +subj - subject. +""" function setelimrange!(data::DataSet{PKSubject}, range::ElimRange, subj::Int) setelimrange!(data[subj], range) data end +""" + setdosetime!(data::DataSet{PKSubject}, dosetime::DoseTime, sort::Dict) + +Set range for elimination parameters calculation for DataSet. + +data - DataSet of PK subject; +range - ElimRange object; +sort - subject sort. +""" +function setelimrange!(data::DataSet{PKSubject}, dosetime::DoseTime, sort::Dict) + for i = 1:length(data) + if sort ∈ data[i].sort setelimrange!(data[i], dosetime) end + end + data +end #------------------------------------------------------------------------------- """ setkelauto!(data::PKSubject, kelauto::Bool) @@ -1176,12 +1242,28 @@ function setdosetime!(data::PKSubject, dosetime::DoseTime) data.dosetime = dosetime data end +""" + setdosetime!(data::DataSet{PKSubject}, dosetime::DoseTime) + +Set dosing time and tau for all subjects. + +data - dataset of PK subjects; +dosetime - DoseTime object. +""" function setdosetime!(data::DataSet{PKSubject}, dosetime::DoseTime) for i = 1:length(data) setdosetime!(data[i], dosetime) end data end +""" + setdosetime!(data::DataSet{PKSubject}, dosetime::DoseTime, sort::Dict) + +Set dosing time and tau for subjects in sort dict. + +data - dataset of PK subjects; +dosetime - DoseTime object. +""" function setdosetime!(data::DataSet{PKSubject}, dosetime::DoseTime, sort::Dict) for i = 1:length(data) if sort ∈ data[i].sort setdosetime!(data[i], dosetime) end diff --git a/src/plots.jl b/src/plots.jl index a5b9495..1ce9f68 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -32,15 +32,16 @@ PKPlotStyle(:dot, :purple, :octagon, :purple), PKPlotStyle(:dot, :indigo, :heptagon, :indigo), ] -function randomstyle() - linestyle = ([:solid, :dash, :dot, :dashdot, :dashdotdot])[sample(1:5)] - linecolor = ([:blue, :red, :green, :yellow, :gray, :cyan, :gold, :magenta, :purple, :indigo])[sample(1:10)] - markershape = ([:circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline])[sample(1:20)] - markercolor = ([:blue, :red, :green, :yellow, :gray, :cyan, :gold, :magenta, :purple, :indigo])[sample(1:10)] +function randomstyle(rng) + linestyle = ([:solid, :dash, :dot, :dashdot, :dashdotdot])[sample(rng, 1:5)] + linecolor = ([:blue, :red, :green, :yellow, :gray, :cyan, :gold, :magenta, :purple, :indigo])[sample(rng, 1:10)] + markershape = ([:circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline])[sample(rng, 1:20)] + markercolor = ([:blue, :red, :green, :yellow, :gray, :cyan, :gold, :magenta, :purple, :indigo])[sample(rng, 1:10)] return PKPlotStyle(linestyle, linecolor, markershape, markercolor) end @userplot PKPlot +@userplot PKElimpPlot @recipe function f(subj::PKPlot) x, y = subj.args @@ -61,6 +62,15 @@ end (x, y) end +@recipe function f(subj::PKElimpPlot) + x, y = subj.args + seriestype --> :line + legend --> false + markersize --> 0 + markerstrokealpha --> 0 + (x, y) +end + function plotlabel(d) title = "" if length(d) > 0 @@ -77,11 +87,13 @@ end Plot for subject -subj - subject; -plotstyle - styles for plots. +* subj - subject; +* plotstyle - styles for plots. """ -function pkplot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject +function pkplot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], ls = false, elim = false, kwargs...) where T <: AbstractSubject + time = subj.time + obs = subj.obs kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -97,19 +109,48 @@ function pkplot(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) whe if !(:ylabel in k) kwargs[:ylabel] = "Concentration" end + if :yscale in k + if kwargs[:yscale] in [:ln, :log, :log2, :log10] ls = true end + end + if ls == true + inds = findall(x->x>0, subj.obs) + time = subj.time[inds] + obs = log.(subj.obs[inds]) + end + - p = pkplot(subj.time, subj.obs; + p = pkplot(time, obs; linestyle = plotstyle.linestyle, linecolor = plotstyle.linecolor, markershape = plotstyle.markershape, markercolor = plotstyle.markercolor, kwargs... ) + if elim + if length(subj.keldata) > 0 + rsq, rsqn = findmax(subj.keldata.ar) + lz = subj.keldata.a[rsqn] + lzint = subj.keldata.b[rsqn] + ts = time[subj.kelrange.kelstart] + te = time[subj.kelrange.kelend] + + if ls true + x = [ts,te] + y = lzint .+ lz .* x + else + x = collect(ts:(te-ts)/100:te) + y = exp.(lzint .+ lz .* x) + end + pkelimpplot!(p, x, y) + end + end return p end -function pkplot!(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) where T <: AbstractSubject +function pkplot!(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], ls = false, elim = false, kwargs...) where T <: AbstractSubject + time = subj.time + obs = subj.obs kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -119,8 +160,17 @@ function pkplot!(subj::T; plotstyle::PKPlotStyle = PKPLOTSTYLE[1], kwargs...) wh if !(:legend in k) kwargs[:legend] = true end + if :yscale in k + if kwargs[:yscale] in [:ln, :log, :log2, :log10] ls = true end + end + if ls == true + inds = findall(x->x>0, subj.obs) + time = subj.time[inds] + obs = log.(subj.obs[inds]) + end - p = pkplot!(subj.time, subj.obs; + + p = pkplot!(time, obs; linestyle = plotstyle.linestyle, linecolor = plotstyle.linecolor, markershape = plotstyle.markershape, @@ -139,7 +189,7 @@ function usort(data::DataSet{T}, list) where T <: AbstractData dl end -function pageplot(pagedatatmp, styledict, typesort; kwargs...) +function pageplot(pagedatatmp, styledict, typesort, utypes; kwargs...) kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -156,7 +206,8 @@ function pageplot(pagedatatmp, styledict, typesort; kwargs...) kwargs[:ylabel] = "Concentration" end - utypes = keys(styledict) + #utypes = keys(styledict) + labels = Vector{String}(undef, 0) fst = true p = nothing @@ -187,7 +238,7 @@ function pageplot(pagedatatmp, styledict, typesort; kwargs...) end p end -function pageplot(pagedatatmp; kwargs...) +function pageplot(pagedatatmp; elim = false, kwargs...) kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -206,11 +257,12 @@ function pageplot(pagedatatmp; kwargs...) fst = true p = nothing + if length(pagedatatmp) > 1 elim = false end for d in pagedatatmp kwargs[:label] = "AUTO" if fst - p = pkplot(d; plotstyle = PKPLOTSTYLE[1], kwargs...) + p = pkplot(d; plotstyle = PKPLOTSTYLE[1], elim = elim, kwargs...) fst = false else pkplot!( d; plotstyle = PKPLOTSTYLE[1], kwargs...) @@ -224,11 +276,11 @@ end Plot for subjects in dataset. -data - subjects dataset; -pagesort - subject page groupping; -typesort - subject sorting within page; +* 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 +function pkplot(data::DataSet{T}; pagesort::Union{Nothing, Vector{Symbol}} = nothing, typesort::Union{Nothing, Vector{Symbol}} = nothing, elim = false, kwargs...) where T <: AbstractSubject # Style types kwargs = Dict{Symbol, Any}(kwargs) k = keys(kwargs) @@ -247,6 +299,7 @@ function pkplot(data::DataSet{T}; pagesort::Union{Nothing, Vector{Symbol}} = no end styledict = nothing + local utypes #style = PKPLOTSTYLE[1] if pagesort !== nothing if !isa(pagesort, Array) pagesort = [pagesort] end @@ -259,7 +312,7 @@ function pkplot(data::DataSet{T}; pagesort::Union{Nothing, Vector{Symbol}} = no if i <= 20 styledict[utypes[i]] = PKPLOTSTYLE[i] else - styledict[utypes[i]] = randomstyle() + styledict[utypes[i]] = randomstyle(MersenneTwister(34534)) end end end @@ -282,20 +335,20 @@ function pkplot(data::DataSet{T}; pagesort::Union{Nothing, Vector{Symbol}} = no end end if typesort !== nothing - p = pageplot(pagedatatmp, styledict, typesort; kwargs...) + p = pageplot(pagedatatmp, styledict, typesort, utypes; kwargs...) else - p = pageplot(pagedatatmp; kwargs...) + p = pageplot(pagedatatmp; elim = elim, kwargs...) end push!(plots, p) end else if kwargs[:title] == :AUTO kwargs[:title] = "" end if typesort !== nothing - p = pageplot(data, styledict, typesort; kwargs...) + p = pageplot(data, styledict, typesort, utypes; kwargs...) else - p = pageplot(data; kwargs...) + p = pageplot(data; elim = elim, kwargs...) end - push!(plots, p) + return p end plots end diff --git a/test/pktest.jl b/test/pktest.jl index 9f4bcfe..9db56cd 100644 --- a/test/pktest.jl +++ b/test/pktest.jl @@ -53,12 +53,15 @@ println(" ---------------------------------- ") #AUMClast #MRTlast #Clast_pred + #Clinf / Cl_F_Obs + #Vzinf / Vz_F_Pred # -------------------------------------------------------------------------- # Linear Trapezoidal Linear Interpolation # calcm = :lint, intp = :lint pkds = ClinicalTrialUtilities.pkimport(pkdata2, [:Subject, :Formulation]; conc = :Concentration, time = :Time) + ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(dose = 100, time = 0)) pk = ClinicalTrialUtilities.nca!(pkds) df = DataFrame(pk; unst = true) sort!(df, :Subject) @@ -196,6 +199,93 @@ println(" ---------------------------------- ") 93.761762 38.810857], sigdigits = 6) + @test round.(df[!, :Clinf], sigdigits = 6) == round.([0.0023296437 +0.0061900608 +0.0038422846 +0.0045446122 +0.0038729255 +0.0062493127 +0.0085550864 +0.0064740799 +0.0040216775 +0.012383404], sigdigits = 6) + + @test round.(df[!, :Vzinf], sigdigits = 6) == round.([0.68827768 +0.43881487 +1.1673601 +0.59056646 +0.56843374 +0.8124135 +0.68666158 +0.72497447 +0.71232266 +0.72039519], sigdigits = 6) + +ClinicalTrialUtilities.setdosetime!(pkds, ClinicalTrialUtilities.DoseTime(dose = 120, time = 0, tau = 12)) +pk = ClinicalTrialUtilities.nca!(pkds; adm = :iv) +df = DataFrame(pk; unst = true) +sort!(df, :Subject) + +#Cltau / CLss +@test round.(df[!, :Cltau], sigdigits = 5) == round.([0.07185191 +0.050414459 +0.12240579 +0.070132959 +0.06902661 +0.085106504 +0.083532913 +0.10859036 +0.073251565 +0.092756742], sigdigits = 5) + +#AUMCtau +@test round.(df[!, :AUMCtau], sigdigits = 6) == round.([9984.8168 +14630.069 +6024.4953 +10299.721 +11466.123 +8467.3568 +9003.0193 +6457.0058 +10095.818 +8367.3005], sigdigits = 6) + +#AUCinf +@test round.(df[!, :AUCinf], sigdigits = 6) == round.([42925.019 +16154.93 +26026.183 +22004.078 +25714.393 +16001.76 +11688.953 +15446.21 +24865.246 +8171.1624], sigdigits = 6) + +#MRTtauinf / MRTinf +@test round.(df[!, :MRTtauinf], sigdigits = 6) == round.([302.40303 +75.590599 +312.72083 +148.34069 +172.0933 +130.19061 +91.908297 +161.57402 +176.30461 +70.260736], sigdigits = 6) + + +@test round.(df[!, :Vsstau], sigdigits = 4) == round.([21.728235 +3.8108592 +38.278842 +10.403572 +11.879017 +11.080068 +7.6773678 +17.545381 +12.914588 +6.517157], sigdigits = 4) + # -------------------------------------------------------------------------- # Linear Log Trapezoidal @@ -625,6 +715,91 @@ println(" ---------------------------------- ") 1638.1903 1293.7065], sigdigits = 6) + #Ctau + @test round.(df[!, :Ctau], sigdigits = 6) == round.([144.964 +196.035 +76.027 +132.257 +154.066 +113.751 +123.37 +134.133 +135.58 +106.476], sigdigits = 6) + + #Cavg + @test round.(df[!, :Cavg], sigdigits = 6) == round.([139.17515 +198.35579 +81.695479 +142.58631 +144.87167 +117.49983 +119.71329 +92.089208 +136.51585 +107.80888], sigdigits = 6) + + #Cavg + @test round.(df[!, :Fluctau], sigdigits = 2) .* 100 == round.([32.983619 +32.840987 +35.886931 +53.500928 +10.538983 +34.806007 +24.962976 +4.5542796 +23.269825 +17.629346], sigdigits = 2) + + #Swingtau + @test round.(df[!, :Swingtau], sigdigits = 6) == round.([0.31666483 +0.3322978 +0.38562616 +0.57679367 +0.099100386 +0.35953091 +0.24223069 +0.031267473 +0.23430447 +0.17850032 +], sigdigits = 6) + + #Accind + @test round.(df[!, :Accind], sigdigits = 6) == round.([25.123662 +6.4216193 +25.821565 +11.336753 +12.737742 +11.341063 +7.2010832 +9.8406852 +15.26571 +5.3650315], sigdigits = 6) + + #Cltau / CLss_F + @test round.(df[!, :Cltau], sigdigits = 6) == round.([0.07185191 +0.050414459 +0.12240579 +0.070132959 +0.06902661 +0.085106504 +0.083532913 +0.10859036 +0.073251565 +0.092756742], sigdigits = 6) + +#Vztau / Vz_F +@test round.(df[!, :Vztau], sigdigits = 6) == round.([21.228167 +3.5738929 +37.18924 +9.113687 +10.131115 +11.063884 +6.7046479 +12.160066 +12.974374 +5.3960536], sigdigits = 6) + # -------------------------------------------------------------------------- # Linear-trapezoidal Linear Interpolation # TAU 2 - 10 diff --git a/test/validation_report.jmd b/test/validation_report.jmd new file mode 100644 index 0000000..a84c75c --- /dev/null +++ b/test/validation_report.jmd @@ -0,0 +1,142 @@ +--- +title: ClinicalTrialUtilities:NCA validation report +author: +date: `j import Dates; Dates.Date(Dates.now())` +mainfont: romanuni.ttf +sansfont: NotoSans-Regular.ttf +monofont: NotoSansMono-Regular.ttf +mathfont: texgyredejavu-math.otf +--- + +```julia; echo = false +using Dates +``` + +\pagebreak + +# Introduction and package description + + +Stable documentation: []() + +## Validation purpose and scope + + +## Requirements + + * Julia 1.5.* (or higher) installed + * Julia packages from dependence list installed (see [Project.toml]()) + +## Developer software life cycle + + * Development stage + * Testing procedures development + * Performing testing procedures on local machine + * Push to master branch + * Performing testing procedures with GitHub Actions + * Make pull request to the official registry of general Julia packages (if nessesary) + * Make release (if previous completed) + +### Versions + + * X.Y.Z - patch release (no breaking changes) + * X.Y.# - minor release (may include breaking changes for versions 0.Y.#) + * X.#.# - major release (breaking changes, changes in public API) + * 0.#.# - no stable public API + * ≥1.#.# - stable public API + + +## Build support + +### Tier 1 + + * julia-version: 1.5, 1.6 + * julia-arch: x64 + * os: ubuntu-18.04, macos-10.15, windows-2019 + +\pagebreak + +# Installation + +## System information + + * Julia version: `j Sys.VERSION` + * Current machine: `j Sys.MACHINE` + +## Installation method + + +```julia; eval = false +import Pkg; Pkg.add("") +``` + +## Version check + +The installation process is checking within each testing job via GitHub Actions. +Also GitHub Action [chek](https://github.com/JuliaRegistries/General/blob/master/.github/workflows/automerge.yml) +performed before merging into JuliaRegistries/General repository +(see [Automatic merging of pull requests](https://github.com/JuliaRegistries/General#automatic-merging-of-pull-requests)). + +```julia; echo = false; results = "hidden" +using Pkg, +pkgversion(m::Module) = Pkg.TOML.parsefile(joinpath(dirname(string(first(methods(m.eval)).file)), "..", "Project.toml"))["version"] +ver = pkgversion(Metida) +``` + +Current package version: +```julia; echo = false; results = "tex" +ver +``` + +\pagebreak + +# Operation qualification + +This part of validation based on testing procedures entails running software products under known conditions with defined inputs and +documented outcomes that can be compared to their predefined expectations. All documented public API included in testing procedures and part of +critical internal methods. + +## Coverage + +Code coverage report available on [Codecov.io](). +Test procedures include all public API methods check. + +* Coverage goal: ≥ 90.0% + +## Data + +For operation checks generated data used. For any purpose, this data available in the repository and included in the package. + +Datasets: + + * + +## Testing results + +```julia +Pkg.test("") +``` + +\pagebreak + +# Performance qualification + +Purpose of this testing procedures to demonstrate performance for critical tasks. + + +\pagebreak + +# Glossary + + * Installation qualification (IQ) - Establishing confidence that process equipment and ancillary systems are compliant with appropriate codes and approved design intentions, and that manufacturer's recommendations are suitably considered. + * Operational qualification (OQ) Establishing confidence that process equipment and sub-systems are capable of consistently operating within established limits and tolerances. + * Product performance qualification (PQ) - Establishing confidence through appropriate testing that the finished product produced by a specified process meets all release requirements for functionality and safety. + * Repository - GitHub repository: https://github.com/PharmCat/Metida.jl + * Master branch - main branch on GitHub ([link](https://github.com/PharmCat/Metida.jl/tree/master)). + * Current machine - pc that used for validation report generating. + +# Reference + +* [General Principles of Software Validation; Final Guidance for Industry and FDA Staff](https://www.fda.gov/media/73141/download) +* [Guidance for Industry Process Validation: General Principles and Practices](https://www.fda.gov/files/drugs/published/Process-Validation--General-Principles-and-Practices.pdf) +* [Glossary of Computer System Software Development Terminology](https://www.fda.gov/inspections-compliance-enforcement-and-criminal-investigations/inspection-guides/glossary-computer-system-software-development-terminology-895) diff --git a/test/weave.jl b/test/weave.jl new file mode 100644 index 0000000..b2d3f1c --- /dev/null +++ b/test/weave.jl @@ -0,0 +1,5 @@ +using Weave, ClinicalTrialUtilities +weave(joinpath(dirname(pathof(ClinicalTrialUtilities)), "..", "test", "validation_report.jmd"); +doctype = "pandoc2pdf", +out_path = :pwd, +pandoc_options=["--toc", "-V colorlinks=true" , "-V linkcolor=blue", "-V urlcolor=red", "-V toccolor=gray"])