From b5253b8acc493d8d3a9b6baeb78cb020983883d3 Mon Sep 17 00:00:00 2001 From: PharmCat <13901158+PharmCat@users.noreply.github.com> Date: Wed, 25 Mar 2020 03:31:32 +0300 Subject: [PATCH] Dev (#7) * callback * linear algebra * memory optimization * isposdef * memory optimization * minor fix, code cosmitics * svd -> svdvals * LCL' gradient optimizations * lclgf deprecate * optimizations, drop iterim data, reml function unification * deprecation, cosmetics * result structure, drop typeiii computing by default * update docs, remove structure description (not included in public API) * public API fix, deprecated functions removed * v1.0.6 * inv cholesky * update * update * update * update * update * update * update * initvar optim * simulation update * minor optimization * rbe settings init * optim * update * generalized simulation * update * update * update * update * update * cosmitics * MvNormal fix * update * update * update * fix Co-authored-by: PharmCat --- .travis.yml | 3 +- Project.toml | 12 +-- README.md | 1 + chagelog.md | 12 +++ examples/simulation.jl | 48 ++++++++++- src/ReplicateBE.jl | 3 +- src/algebra.jl | 65 ++++++++++----- src/deprecated.jl | 12 +++ src/generalfunc.jl | 180 +++++++++++++++++++-------------------- src/randrbeds.jl | 185 +++++++++++++++++++++++++++++++++++++++-- src/rbe.jl | 49 +++++++---- src/rbesettings.jl | 31 +++++++ src/utils.jl | 5 ++ test/test.jl | 15 ++++ 14 files changed, 476 insertions(+), 145 deletions(-) create mode 100644 src/rbesettings.jl diff --git a/.travis.yml b/.travis.yml index daf05bd..6a59491 100644 --- a/.travis.yml +++ b/.travis.yml @@ -9,6 +9,7 @@ julia: - 1.0 - 1.2 - 1.3 + - 1.4 branches: only: @@ -16,7 +17,7 @@ branches: matrix: allow_failures: - - julia: 1.3 + - julia: 1.4 notifications: email: false diff --git a/Project.toml b/Project.toml index 31c563e..9d25ba2 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,7 @@ name = "ReplicateBE" uuid = "671d9d50-c343-11e9-1a9c-fdd992682823" keywords = ["bioequivalence", "mixedmodel"] desc = "Mixed model solution for replicate designed bioequivalence study." -version = "1.0.6" +version = "1.0.7" [deps] Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -19,6 +19,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Optim = "429524aa-4258-5aef-a3af-852621145aeb" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extras] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" @@ -30,12 +31,13 @@ test = ["CSV", "Test", "StatsBase"] [compat] CategoricalArrays = "0.7" -julia = "1.0, 1.1, 1.2" +julia = "1.0, 1.1, 1.2, 1.3, 1.4" StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32" -DataFrames = "0.19" +DataFrames = "0.19, 0.20" StatsModels = "0.6" -Distributions = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21" +Distributions = "0.20, 0.21, 0.22, 0.23" PDMats = "0.9" ForwardDiff = "0.10" -Optim = "0.19" +Optim = "0.19, 0.20" LineSearches = "7.0" +StaticArrays = "0.11, 0.12" diff --git a/README.md b/README.md index f158819..3d716cd 100644 --- a/README.md +++ b/README.md @@ -11,6 +11,7 @@ Mixed model solution for replicate designed bioequivalence study. This can be us [![codecov](https://codecov.io/gh/PharmCat/ReplicateBE.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/PharmCat/ReplicateBE.jl) [![Latest docs](https://img.shields.io/badge/docs-latest-blue.svg)](https://pharmcat.github.io/ReplicateBE.jl/latest/) [![doi](https://img.shields.io/badge/doi-10.13140%2FRG.2.2.27418.39363-blue)](https://doi.org/10.13140/RG.2.2.27418.39363) +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3543823.svg)](https://doi.org/10.5281/zenodo.3543823) Install: ``` diff --git a/chagelog.md b/chagelog.md index 871f717..9775822 100644 --- a/chagelog.md +++ b/chagelog.md @@ -1,3 +1,15 @@ + +- v1.0.7 + * settings: initial + * randrbeds "light" generation + * optimizations + * Distributions bump + * "generalized" simulation + * Optim 0.20 + * Julia 1.3 + + + - v1.0.6 * many optimizations * linear algebra custom functions diff --git a/examples/simulation.jl b/examples/simulation.jl index adc53b4..b1ab065 100644 --- a/examples/simulation.jl +++ b/examples/simulation.jl @@ -1,6 +1,6 @@ using ReplicateBE -task = ReplicateBE.RandRBEDS(;n=24, +task = ReplicateBE.randrbetask(;n=24, sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], inter=[0.05, 0.04, 0.6], intra=[0.02, 0.04], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], @@ -15,7 +15,7 @@ Number: 100000 Result: 0.05078660225829358 =# -task = ReplicateBE.RandRBEDS(;n=24, +task = ReplicateBE.randrbetask(;n=24, sequence=[1,2], design = ["T" "R" "T"; "R" "T" "R"], inter=[0.05, 0.04, 0.4], intra=[0.02, 0.04], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0], @@ -30,7 +30,7 @@ Result: 0.04997148203368122 #E63 =# -task = ReplicateBE.RandRBEDS(;n=24, +task = ReplicateBE.randrbetask(;n=24, sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], inter=[0.05, 0.04, 0.3], intra=[0.02, 0.04], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], @@ -38,3 +38,45 @@ formcoef = [0.0, log(0.8)], seed = 0, dropobs = 10) result = ReplicateBE.simulation(task; num = 2000, seed = 6564683774737, verbose = true) + +task = ReplicateBE.randrbetask(;n=24, +sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], +inter=[0.01, 0.01, 0.9], intra=[0.09, 0.09], +intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], +formcoef = [log(0.95), 0.0], seed = 0) + +rds = ReplicateBE.randrbeds(task) +be = ReplicateBE.rbe!(rds, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence) + +result = ReplicateBE.simulation(task; num = 1002, verbose = true, rsabe = false) + + + +task = ReplicateBE.randrbetask(;n=24, +sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], +inter=2000.0, intra=[0.09, 0.09], +intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], +formcoef = [log(0.95), 0.0], seed = 0) + +rds = ReplicateBE.randrbeds(task) +result = ReplicateBE.simulation(task; num = 1002, verbose = true, rsabe = true) + + + + +task = ReplicateBE.randrbetask(;n=24, +sequence=[1,1], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], +inter=[0.04, 0.04, 1.0], intra=[0.02, 0.02], +intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], +formcoef = [0.0, log(0.8)], seed = 0, dropobs = 0) + +out = zeros(Float64, 0) +function simfunc!(out, be) + push!(out, ReplicateBE.theta(be)[5]) +end +result = ReplicateBE.simulation!(task, out, simfunc!; num = 10, seed = 10, verbose = false) +m = mean(result) + +using Plots + +histogram(result) diff --git a/src/ReplicateBE.jl b/src/ReplicateBE.jl index c474b94..af12f3f 100644 --- a/src/ReplicateBE.jl +++ b/src/ReplicateBE.jl @@ -5,7 +5,7 @@ module ReplicateBE -using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches, CategoricalArrays, Printf +using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches, CategoricalArrays, Printf, StaticArrays export rbe, rbe!, reml2, nobs, coef, stderror, dof, coefnum, fixed, theta, typeiii, design, show, confint, contrast, estimate, optstat, randrbeds, randrbetask import Base.show @@ -16,6 +16,7 @@ using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlge const LOG2PI = log(2π) const MEMOPT = true +include("rbesettings.jl") include("rbetable.jl") include("memalloc.jl") include("design.jl") diff --git a/src/algebra.jl b/src/algebra.jl index fc07ab7..e3543b2 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -57,7 +57,6 @@ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::Vector, c) end @inbounds β[i] += c[n] * C[n] end - for n = 1:p for m = 1:q @inbounds θ[i, n] += A[m, n] * c[m] @@ -77,7 +76,7 @@ function mulall(A::Vector, B::AbstractMatrix) for n = 1:q for m = 1:q #@inbounds c[n] += B[m, n] * A[m] - θ += B[m, n] * A[m] * A[n] + @inbounds θ += B[m, n] * A[m] * A[n] end #@inbounds θ += A[n] * c[n] end @@ -101,7 +100,7 @@ function mulθ₃(y::Vector, X::AbstractMatrix, β::Vector, V::AbstractMatrix, c end for n = 1:q for m = 1:q - θ += V[m, n] * (y[m] - c[m]) * (y[n] - c[n]) + @inbounds θ += V[m, n] * (y[m] - c[m]) * (y[n] - c[n]) end end return θ @@ -120,12 +119,12 @@ function mulαβαt(A::AbstractMatrix, B::AbstractMatrix) c .= 0 for n = 1:q for m = 1:q - c[n] += A[i, m] * B[n, m] + @inbounds c[n] += A[i, m] * B[n, m] end end for n = 1:p for m = 1:q - mx[i, n] += A[n, m] * c[m] + @inbounds mx[i, n] += A[n, m] * c[m] end end end @@ -143,16 +142,17 @@ function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) c .= 0 for n = 1:q for m = 1:q - c[n] += A[i, m] * B[n, m] + @inbounds c[n] += A[i, m] * B[n, m] end end for n = 1:p for m = 1:q - mx[i, n] += A[n, m] * c[m] + @inbounds mx[i, n] += A[n, m] * c[m] end end end - mx .+ C + mx .+= C + end function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c::AbstractVector) q = size(B, 1) @@ -163,21 +163,18 @@ function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c: c .= 0 for n = 1:q for m = 1:q - try - c[n] += A[i, m] * B[n, m] - catch - println("$q, $p, $i, $n, $m" ) - error() - end + @inbounds c[n] += A[i, m] * B[n, m] end end for n = 1:p for m = 1:q - mx[i, n] += A[n, m] * c[m] + @inbounds mx[i, n] += A[n, m] * c[m] end end end - mx .+ C + mx .+= C + #SMatrix{p,p,eltype(mx)}(mx) + end """ A * B * A' + C -> O @@ -191,12 +188,12 @@ function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C c .= 0 for n = 1:q for m = 1:q - c[n] += A[i, m] * B[n, m] + @inbounds c[n] += A[i, m] * B[n, m] end end for n = 1:p for m = 1:q - O[i, n] += A[n, m] * c[m] + @inbounds O[i, n] += A[n, m] * c[m] end end end @@ -211,14 +208,42 @@ function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C c .= 0 for n = 1:q for m = 1:q - c[n] += A[i, m] * B[n, m] + @inbounds c[n] += A[i, m] * B[n, m] end end for n = 1:p for m = 1:q - O[i, n] += A[n, m] * c[m] + @inbounds O[i, n] += A[n, m] * c[m] end end end O .+= C end + +function invchol(M) + q = size(M, 1) + v = zeros(eltype(M), q, q) + if q == 1 + v[1,1] = 1/M[1,1] + return v + end + il = inv(cholesky(M).U) + for n = 1:q + for m = n:q + @inbounds v[n, n] += il[n, m]^2 + end + end + for n = 1:(q-1) + for m = (n+1):q + for i = m:q + @inbounds v[n, m] += il[n, i] * il[m, i] + end + end + end + for n = 1:q-1 + for m = 2:q + @inbounds v[m, n] = v[n, m] + end + end + return v +end diff --git a/src/deprecated.jl b/src/deprecated.jl index ad29eda..ad9c2e9 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -1,4 +1,16 @@ +""" + Feel M with zero feeled matrices +""" +function matvecz!(M, Zv) + n = length(Zv) + for i = 1:n + M[i] = zeros(size(Zv[i])[1], size(Zv[i])[1]) + end + return +end + + struct MemAlloc{T <: AbstractFloat} mem1::Vector{Matrix{T}} mem2::Vector{Matrix{T}} diff --git a/src/generalfunc.jl b/src/generalfunc.jl index a425ad0..5b969d4 100644 --- a/src/generalfunc.jl +++ b/src/generalfunc.jl @@ -78,27 +78,28 @@ end return end function mvmat(G::AbstractMatrix, σ::Vector, Z::Matrix, cache)::Matrix - h = hash(tuple(σ, Z)) - if h in keys(cache) - return cache[h] + #h = hash(tuple(σ, Z)) + if Z in keys(cache) + return cache[Z] else - #V = mulαβαtc(Z, G, Diagonal(Z*σ), mem) - V = Z * G * Z' + Diagonal(Z*σ) - cache[h] = V + V = mulαβαtc(Z, G, Diagonal(Z*σ)) + #V = Z * G * Z' + Diagonal(Z*σ) + cache[Z] = V return V end end function mvmat(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache)::Matrix - h = hash(tuple(σ, Z)) - if h in keys(cache) - return cache[h] + #h = hash(tuple(σ, Z)) + if Z in keys(cache) + return cache[Z] else V = mulαβαtc(Z, G, Diagonal(Z*σ), mem) #V = Z * G * Z' + Diagonal(Z*σ) - cache[h] = V + cache[Z] = V return V end end + function mvmatall(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache) #h = hash(Z) #if h in keys(cache) @@ -108,41 +109,47 @@ function mvmatall(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache) else V = mulαβαtc(Z, G, Diagonal(Z*σ), mem) #V = Z * G * Z' + Diagonal(Z*σ) - iV = inv(V) - ldV = logdet(V) - cache[Z] = (V, iV, ldV) - return V, iV, ldV - end -end -""" - Feel M with zero feeled matrices -""" -function matvecz!(M, Zv) - n = length(Zv) - for i = 1:n - M[i] = zeros(size(Zv[i])[1], size(Zv[i])[1]) + V⁻¹ = nothing + try + V⁻¹ = invchol(V) + catch e + if typeof(e) <: PosDefException + V⁻¹ = inv(V) + else + throw(e) + end + end + #V⁻¹ = inv(V) + log│V│ = logdet(V) + cache[Z] = (V, V⁻¹, log│V│) + return V, V⁻¹, log│V│ end - return end #println("θ₁: ", θ1, " θ₂: ", θ2, " θ₃: ", θ3) -function minv(M::Matrix, cache::Dict)::Matrix - h = hash(M) - if h in keys(cache) - return cache[h] +function minv(G::AbstractMatrix, σ::Vector, Z::Matrix, cache::Dict)::Matrix + #h = hash(M) + if Z in keys(cache) + #return cache[h] + return cache[Z] else - iM = inv(M) - cache[h] = iM - return iM + V = mulαβαtc(Z, G, Diagonal(Z*σ)) + #V = Z * G * Z' + Diagonal(Z*σ) + V⁻¹ = invchol(V) + #V⁻¹ = inv(V) + #ldV = logdet(V) + cache[Z] = V⁻¹ + return V⁻¹ end end + function mlogdet(M::Matrix, cache::Dict) - h = hash(M) - if h in keys(cache) - return cache[h] + #h = hash(M) + if M in keys(cache) + return cache[M] else iM = logdet(M) - cache[h] = iM + cache[M] = iM return iM end end @@ -152,84 +159,83 @@ end """ -2 REML function for ForwardDiff """ -function reml2(data::RBEDataStructure, θvec::Vector, β::Vector; memopt::Bool = true) - +function reml2(data::RBEDataStructure, θ, β::Vector; memopt::Bool = true) #memory optimizations to reduse allocations (cache rebuild) #empty!(data.mem.dict) - rebuildcache(data, promote_type(eltype(data.yv[1]), eltype(θvec))) + rebuildcache(data, promote_type(eltype(data.yv[1]), eltype(θ))) cache = Dict{Matrix, Tuple{Matrix, Matrix, Number}}() #cache = data.mem.dict #--------------------------------------------------------------------------- - G = gmat(θvec[3:5]) - θ1 = 0 - θ2 = zeros(promote_type(eltype(data.yv[1]), eltype(θvec)), data.p, data.p) - θ3 = 0 - iV = nothing + G = gmat(θ[3:5]) + θ₁ = 0 + θ₂ = zeros(promote_type(eltype(data.yv[1]), eltype(θ)), data.p, data.p) + θ₃ = 0 + V⁻¹ = nothing for i = 1:data.n if MEMOPT && memopt - V, iV, ldV = mvmatall(G, θvec[1:2], data.Zv[i], first(data.mem.svec), cache) - θ1 += ldV + V, V⁻¹, log│V│ = mvmatall(G, θ[1:2], data.Zv[i], first(data.mem.svec), cache) + θ₁ += log│V│ else - @inbounds V = vmat(G, rmat(θvec[1:2], data.Zv[i]), data.Zv[i]) - iV = inv(V) - θ1 += logdet(V) + @inbounds V = vmat(G, rmat(θ[1:2], data.Zv[i]), data.Zv[i]) + V⁻¹ = invchol(V) + θ₁ += logdet(V) end #----------------------------------------------------------------------- #θ2 += Xv[i]'*iV*Xv[i] - @inbounds mulαtβαinc!(θ2, data.Xv[i], iV, first(data.mem.svec)) + @inbounds mulαtβαinc!(θ₂, data.Xv[i], V⁻¹, first(data.mem.svec)) #----------------------------------------------------------------------- #r = yv[i] - Xv[i] * β #θ3 += r' * iV * r - @inbounds θ3 += mulθ₃(data.yv[i], data.Xv[i], β, iV, first(data.mem.svec)) + @inbounds θ₃ += mulθ₃(data.yv[i], data.Xv[i], β, V⁻¹, first(data.mem.svec)) end - return θ1 + logdet(θ2) + θ3 + data.remlc + return θ₁ + logdet(θ₂) + θ₃ + data.remlc + end """ -2 REML estimation with β recalculation for ForwardDiff """ -function reml2bfd(data::RBEDataStructure, θvec::Vector; memopt::Bool = true) - return reml2b(data, θvec; memopt = memopt)[1] +function reml2bfd(data::RBEDataStructure, θ; memopt::Bool = true) + return reml2b(data, θ; memopt = memopt)[1] end -function reml2b(data::RBEDataStructure, θvec::Vector; memopt::Bool = true) +function reml2b(data::RBEDataStructure, θ; memopt::Bool = true) - rebuildcache(data, promote_type(eltype(data.yv[1]), eltype(θvec))) + rebuildcache(data, promote_type(eltype(data.yv[1]), eltype(θ))) cache = Dict() #--------------------------------------------------------------------------- - G = gmat(θvec[3:5]) - iVv = Array{Array{eltype(θvec), 2}, 1}(undef, data.n) + G = gmat(θ[3:5]) + V⁻¹ = Vector{Matrix{eltype(θ)}}(undef, data.n) # Vector of V⁻¹ matrices V = nothing - ldV = nothing - θ1 = 0 - θ2 = zeros(promote_type(eltype(first(data.yv)), eltype(θvec)), data.p, data.p) - θ3 = 0 - βm = zeros(promote_type(eltype(first(data.yv)), eltype(θvec)), data.p) - β = zeros(promote_type(eltype(first(data.yv)), eltype(θvec)), data.p) + log│V│ = nothing # Vector log determinant of V matrix + θ₁ = 0 + θ₂ = zeros(promote_type(eltype(first(data.yv)), eltype(θ)), data.p, data.p) + θ₃ = 0 + βm = zeros(promote_type(eltype(first(data.yv)), eltype(θ)), data.p) + β = zeros(promote_type(eltype(first(data.yv)), eltype(θ)), data.p) for i = 1:data.n if MEMOPT && memopt - - @inbounds V, iVv[i], ldV = mvmatall(G, θvec[1:2], data.Zv[i], first(data.mem.svec), cache) - θ1 += ldV + @inbounds V, V⁻¹[i], log│V│ = mvmatall(G, θ[1:2], data.Zv[i], first(data.mem.svec), cache) + θ₁ += log│V│ else - @inbounds R = rmat(θvec[1:2], data.Zv[i]) + @inbounds R = rmat(θ[1:2], data.Zv[i]) @inbounds V = vmat(G, R, data.Zv[i]) - @inbounds iVv[i] = inv(V) - θ1 += logdet(V) + @inbounds V⁻¹[i] = invchol(V) + θ₁ += logdet(V) end #----------------------------------------------------------------------- #θ2 += Xv[i]'*iVv[i]*Xv[i] #βm += Xv[i]'*iVv[i]*yv[i] - mulθβinc!(θ2, βm, data.Xv[i], iVv[i], data.yv[i], first(data.mem.svec)) + mulθβinc!(θ₂, βm, data.Xv[i], V⁻¹[i], data.yv[i], first(data.mem.svec)) #----------------------------------------------------------------------- end - mul!(β, inv(θ2), βm) + mul!(β, inv(θ₂), βm) for i = 1:data.n # r = yv[i] - Xv[i] * β # θ3 += r' * iVv[i] * r - @inbounds θ3 += mulθ₃(data.yv[i], data.Xv[i], β, iVv[i], first(data.mem.svec)) + @inbounds θ₃ += mulθ₃(data.yv[i], data.Xv[i], β, V⁻¹[i], first(data.mem.svec)) end - return θ1 + logdet(θ2) + θ3 + data.remlc, β, θ2 + return θ₁ + logdet(θ₂) + θ₃ + data.remlc, β, θ₂ end #------------------------------------------------------------------------------- @@ -241,7 +247,7 @@ non inverted C matrix gradient function """ function cmatgf(Xv::Vector, Zv::Vector, θ::Vector; memopt::Bool = true) p = size(Xv[1], 2) - jC = ForwardDiff.jacobian(x -> cmatvec(Xv, Zv, x; memopt = memopt), θ) + jC = ForwardDiff.jacobian(x -> cmatvec(Xv, Zv, x; memopt = memopt), SVector{length(θ), eltype(θ)}(θ)) result = Vector{Matrix}(undef, 0) for i in 1:length(θ) push!(result, reshape(jC[:,i], p, p)) @@ -251,21 +257,21 @@ end """ non inverted C matrix in vector form for gradient """ -function cmatvec(Xv::Vector, Zv::Vector, θ::Vector; memopt::Bool = true) +function cmatvec(Xv::Vector, Zv::Vector, θ; memopt::Bool = true) p = size(Xv[1], 2) G = gmat(θ[3:5]) C = zeros(promote_type(eltype(Zv[1]), eltype(θ)), p, p) cache = Dict() - cachem = Dict() + #cachem = Dict() for i = 1:length(Xv) if memopt - iV = minv(mvmat(G, θ[1:2], Zv[i], cachem), cache) + V⁻¹ = minv(G, θ[1:2], Zv[i], cache) else R = rmat(θ[1:2], Zv[i]) - iV = inv(vmat(G, R, Zv[i])) + V⁻¹ = invchol(vmat(G, R, Zv[i])) end - #C += Xv[i]' * iV * Xv[i] - mulαtβαinc!(C, Xv[i], iV) + #C += Xv[i]' * V⁻¹ * Xv[i] + mulαtβαinc!(C, Xv[i], V⁻¹) end return C[:] end @@ -344,23 +350,17 @@ end #------------------------------------------------------------------------------- #------------------------------------------------------------------------------- #------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- """ Initial variance computation """ -function initvar(df::DataFrame, dv::Symbol, fac::Symbol, sbj::Symbol)::Vector - u = unique(df[:, sbj]) +function initvar(df::DataFrame, dv::Symbol, fac::Symbol)::Vector f = unique(df[:, fac]) - fv = Array{eltype(df[!, dv]), 1}(undef, 0) - sb = Array{eltype(df[!, dv]), 1}(undef, 0) for i in f push!(fv, var(df[df[!, fac] .== i, dv])) end - for i in u - sv = var(df[df[!, sbj] .== i, dv]) - if sv > 0 push!(sb, sv) end - end - push!(fv, mean(sb)) return fv end #------------------------------------------------------------------------------- @@ -396,7 +396,7 @@ function rholinksigmoid2r(ρ, m) return tan(ρ) end -function varlinkmap(θ::Vector, r1::Union{Int, UnitRange}, r2::Union{Int, UnitRange}, f1::Function, f2::Function) +function varlinkmap(θ, r1::Union{Int, UnitRange}, r2::Union{Int, UnitRange}, f1::Function, f2::Function) θl = similar(θ) θl[r1] = f1.(θ[r1]) θl[r2] = f2.(θ[r2]) diff --git a/src/randrbeds.jl b/src/randrbeds.jl index d6801e1..ddfdb17 100644 --- a/src/randrbeds.jl +++ b/src/randrbeds.jl @@ -22,7 +22,7 @@ mutable struct RandRBEDS n::Int sequence::Vector design::Matrix - inter::Vector + inter::Union{Vector, Real} intra::Vector intercept::Real seqcoef::Vector @@ -33,7 +33,7 @@ mutable struct RandRBEDS seed function RandRBEDS(n::Int, sequence::Vector, design::Matrix, - θinter::Vector, θintra::Vector, + θinter, θintra::Vector, intercept::Real, seqcoef::Vector, periodcoef::Vector, formcoef::Vector, dropobs::Int, seed) new(n, sequence, @@ -188,7 +188,8 @@ function randrbeds(n::Int, sequence::Vector, subjmx[:, 2] = design[i,:] subjmx[:, 3] = collect(1:pnum) subjmx[:, 4] .= sqname[i] - subjmx[:, 5] = rand(rng, MvNormal(PDMat(Vv[i]))) + Mv[i] + #subjmx[:, 5] = rand(rng, MvNormal(PDMat(Vv[i]))) + Mv[i] + subjmx[:, 5] = rand(rng, MvNormal(Vv[i])) + Mv[i] subj += 1 for c = 1:pnum push!(subjds, subjmx[c, :]) @@ -241,7 +242,17 @@ Count successful BE outcomes. * ``seed = 0`` - seed for random number generator (0 - random seed) """ -function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l = log(0.8), u = log(1.25), seed = 0) +function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l = log(0.8), u = log(1.25), seed = 0, rsabe = false, rsabeconst = 0.760, reference = "R", alpha = 0.05) + tl = l + tu = u + if rsabe + if reference ∉ task.design + rsabe = false + @warn "Reference value not found in design. RSABE set false." + end + end + + #max range 69,84-143,19 task.seed = 0 #rng = MersenneTwister() if isa(seed, Array) @@ -269,6 +280,12 @@ function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l println(io, "Simulation seed: $(seed)") end println(io, "Task hash: $(hash(task))") + println(io, "Alpha: $(alpha)") + println(io, "RSABE: $(rsabe)") + if rsabe + println(io, "Regulatory const: $(rsabeconst)") + println(io, "Reference formulation: $(reference)") + end end for i = 1:num @@ -277,9 +294,11 @@ function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l try be = rbe(rds, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence) - q = quantile(TDist(be.fixed.df[end]), 0.95) - ll = be.fixed.est[end] - q*be.fixed.se[end] - ul = be.fixed.est[end] + q*be.fixed.se[end] + #q = quantile(TDist(be.fixed.df[end]), 1.0 - alpha) + #ll = be.fixed.est[end] - q*be.fixed.se[end] + #ul = be.fixed.est[end] + q*be.fixed.se[end] + ll, ul = confint(be, 2 * alpha)[end] + #println("CI: $(exp(ll)) - $(exp(ul)) ") #! if verbose if !optstat(be) printstyled(io, "Iteration: ", i, ", seed ", seeds[i], ": unconverged! \n"; color = :yellow) end @@ -287,7 +306,27 @@ function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l printstyled(io, "Iteration: ", i, ", seed ", seeds[i], ": Hessian not positive defined! \n"; color = :yellow) end end - if ll > l && ul < u + # If RSABE is true - calculating CI limits + if rsabe + σ² = intravar(be)[reference] + if geocv(σ²) > 0.30 + bconst = rsabeconst * sqrt(σ²) + tl = -bconst + tu = bconst + if tu > log(1.4319) || tl < log(0.6984) + tu = log(1.4319) + tl = log(0.6984) + end + if verbose + println("Reference scaled: $(exp(tl)*100) - $(exp(tu)*100)") + end + else + tl = l + tu = u + end + + end + if ll > tl && ul < tu cnt += 1 end #! @@ -309,6 +348,64 @@ function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l return RBEDSSimResult(seed, num, seeds, cnt/(num - err), err) end + + +""" +```julia +simulation!(task::RandRBEDS, out, simfunc!::Function; io = stdout, verbose = false, num = 100, seed = 0) +``` + +Generalized simulation method. + +""" +function simulation!(task::RandRBEDS, out, simfunc!::Function; io = stdout, verbose = false, num = 100, seed = 0) + task.seed = 0 + if isa(seed, Array) + seeds = seed + else + if seed != 0 + rng = MersenneTwister(seed) + else + rng = MersenneTwister() + end + seeds = Array{UInt32, 1}(undef, num) + for i = 1:num + seeds[i] = rand(rng, UInt32) + end + end + n = 0 + err = 0 + if verbose + printstyled(io, "Custom simulation start...\n"; color = :green) + if isa(seed, Array) + println(io, "Simulation seed: Array") + else + println(io, "Simulation seed: $(seed)") + end + println(io, "Task hash: $(hash(task))") + end + for i = 1:num + task.seed = seeds[i] + rds = randrbeds(task) + try + be = rbe(rds, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence) + simfunc!(out, be) + if n > 1000 + println(io, "Iteration: $i") + println(io, "Mem: $(Sys.free_memory()/2^20)") + println(io, "Pow: $(cnt/i)") + println(io, "-------------------------------") + n = 0 + end + n += 1 + catch + err += 1 + printstyled(io, "Iteration: $i, seed $(seeds[i]): $(err): ERROR! \n"; color = :red) + end + end + return out +end + function Base.show(io::IO, obj::RBEDSSimResult) if isa(obj.seed, Array) println(io, "Simulation seed: Array") @@ -320,3 +417,75 @@ function Base.show(io::IO, obj::RBEDSSimResult) println(io, "Errors: $(obj.errn)") println(io, "Result: $(obj.result)") end + + +function randrbeds(n::Int, sequence::Vector, + design::Matrix, + inter::Real, intra::Vector, + intercept::Real, seqcoef::Vector, periodcoef::Vector, formcoef::Vector, + dropobs::Int, seed) + if seed != 0 + rng = MersenneTwister(seed) + else + rng = MersenneTwister() + end + + r = n/sum(sequence) + sn = Array{Int, 1}(undef, length(sequence)) + for i = 1:(length(sequence)-1) + sn[i] = round(r*sequence[i]) + end + sn[length(sequence)] = n - sum(sn[1:(length(sequence)-1)]) + + u = unique(design) + sqname = Array{String, 1}(undef,size(design)[1]) + sqnum = size(design)[1] + pnum = size(design)[2] + for i = 1:sqnum + sqname[i] = join(design[i,:]) + end + Zv = Array{Matrix, 1}(undef, sqnum) + Vv = Array{Vector, 1}(undef, sqnum) + for i = 1:size(design)[1] + Z = Array{Int, 2}(undef, pnum, length(u)) + for c = 1:pnum + for uc = 1:length(u) + if design[i, c] == u[uc] Z[c, uc] = 1 else Z[c, uc] = 0 end + end + end + Zv[i] = Z + Vv[i] = Z * intra + end + Mv = Array{Array{Float64, 1}, 1}(undef, sqnum) + for i = 1:sqnum + Mv[i] = zeros(pnum) .+ intercept .+ seqcoef[i] + periodcoef + Zv[i]*formcoef + end + ndist = Normal() + subjds = DataFrame(subject = Int[], formulation = String[], period = Int[], sequence = String[], var = Float64[]) + subj = 1 + subjmx = Array{Any, 2}(undef, pnum, 5) + for i = 1:sqnum + for sis = 1:sn[i] + subjmx[:, 1] .= subj + subjmx[:, 2] = design[i,:] + subjmx[:, 3] = collect(1:pnum) + subjmx[:, 4] .= sqname[i] + subjmx[:, 5] .= 0 + subjmx[:, 5] .+= rand(rng, ndist)*sqrt(inter) + subj += 1 + for c = 1:pnum + subjmx[c, 5] += Mv[i][c] + rand(rng, ndist)*sqrt(Vv[i][c]) + push!(subjds, subjmx[c, :]) + end + end + end + if dropobs > 0 && dropobs < size(subjds, 1) + dellist = sample(rng, 1:size(subjds, 1), dropobs, replace = false) + deleterows!(subjds, sort!(dellist)) + end + categorical!(subjds, :subject); + categorical!(subjds, :formulation); + categorical!(subjds, :period); + categorical!(subjds, :sequence); + return subjds +end diff --git a/src/rbe.jl b/src/rbe.jl index db75052..411b7e4 100644 --- a/src/rbe.jl +++ b/src/rbe.jl @@ -3,21 +3,24 @@ # """ ```julia -struct RBE{T <: AbstractFloat} +mutable struct RBE{T <: AbstractFloat} model::ModelFrame # Model frame rmodel::ModelFrame # Random effect model design::Design - factors::Vector{Symbol} # Factor list θ0::Vector{T} # Initial variance paramethers vlm::T θ::Tuple{Vararg{T}} # Final variance paramethers reml::T # -2REML fixed::EffectTable - typeiii::ContrastTable - G::Matrix{T} # G matrix + G::AbstractMatrix{T} # G matrix C::Matrix{T} # C var(β) p×p variance-covariance matrix A::Matrix{T} # asymptotic variance-covariance matrix ofb θ H::Matrix{T} # Hessian matrix + X::Matrix{T} # Matrix for fixed effects + Z::Matrix{T} # Matrix for random effects + data::RBEDataStructure # Data for model fitting + result::RBEResults + detH::T # Hessian determinant preoptim::Union{Optim.MultivariateOptimizationResults, Nothing} # Pre-optimization result object optim::Optim.MultivariateOptimizationResults # Optimization result object @@ -50,6 +53,7 @@ mutable struct RBE{T <: AbstractFloat} optim::Optim.MultivariateOptimizationResults # Optimization result object end + """ ```julia rbe(df; dvar::Symbol, @@ -118,7 +122,7 @@ function rbe(df; dvar::Symbol, #Check if any(x -> x ∉ names(df), [subject, formulation, period, sequence]) throw(ArgumentError("Names not found in DataFrame!")) end if !(eltype(df[!,dvar]) <: AbstractFloat) - @warn "Responce variable ∉ Float!" + @warn "Responce variable ∉ AbstractFloat!" end if !(typeof(df[!,subject]) <: CategoricalArray) @warn "Subject variable not Categorical, use rbe!()!" @@ -168,9 +172,10 @@ function rbe(df; dvar::Symbol, if length(init) == 5 θvec0 = init else - iv = initvar(df, dvar, formulation, subject) + intra = sum(replace!(var.(yv) .* (length.(yv) .- 1), NaN => 0))/(sum(length.(yv))-1) + iv = initvar(df, dvar, formulation) iv = iv .+ eps() - θvec0 = [iv[3], iv[3], iv[1], iv[2], 0.5] + θvec0 = [intra, intra, iv[1], iv[2], 0.5] end #Variance link function @@ -184,8 +189,7 @@ function rbe(df; dvar::Symbol, varlink = (x, y) -> varlinkmap(x, 1:4, 5, vlink, z -> rholinksigmoid2(z, y)) rvarlink = (x, y) -> varlinkmap(x, 1:4, 5, vlinkr, z -> rholinksigmoid2r(z, y)) else - varlink = (x, y) -> varlinkmap(x, 1:4, 5, vlink, z -> rholinkpsigmoid(z, y)) - rvarlink = (x, y) -> varlinkmap(x, 1:4, 5, vlinkr, z -> rholinkpsigmoidr(z, y)) + throw(ArgumentError("Unknown link function! Check rholink!")) end θvec0 = rvarlink(θvec0, vlm) @@ -265,9 +269,10 @@ function rbe(df; dvar::Symbol, lcl = L*C*Lt #lcl = L*C*L' lclr = rank(lcl) se[i] = sqrt((lcl)[1]) - F[i] = β'*L'*inv(lcl)*L*β/lclr #F[i] = (L*β)'*inv(L*C*L')*(L*β)/lclr - g = lclg(gradc, L) - df[i] = max(1, 2*((lcl)[1])^2/(g'*(A)*g)) + F[i] = β'*L'*inv(lcl)*L*β/lclr + df[i] = sattdf(data, result, L, lcl) #F[i] = (L*β)'*inv(L*C*L')*(L*β)/lclr + #g = lclg(gradc, L) + #df[i] = max(1, 2*((lcl)[1])^2/(g'*(A)*g)) t[i] = ((L*β)/se[i])[1] pval[i] = ccdf(TDist(df[i]), abs(t[i]))*2 end @@ -319,13 +324,23 @@ function rbe!(df; dvar::Symbol, categorical!(df, sequence); end sort!(df, [subject, formulation, period]) - return rbe(df, dvar=dvar, subject=subject, formulation=formulation, period=period, sequence=sequence, - g_tol=g_tol, x_tol=x_tol, f_tol=f_tol, iterations=iterations, - store_trace=store_trace, extended_trace=extended_trace, show_trace=show_trace, - memopt=memopt, init=init, postopt=postopt, vlm = vlm, maxopttry = maxopttry, rhoadjstep = rhoadjstep, + return rbe(df, dvar = dvar, subject = subject, formulation = formulation, period = period, sequence = sequence, + g_tol = g_tol, x_tol = x_tol, f_tol = f_tol, iterations = iterations, + store_trace = store_trace, extended_trace = extended_trace, show_trace = show_trace, + memopt = memopt, init = init, postopt = postopt, vlm = vlm, maxopttry = maxopttry, rhoadjstep = rhoadjstep, rholink = rholink, singlim = singlim) +end +function fit!(rbe::RBE) end + +function rbe(df, settings; dvar::Symbol, + subject::Symbol, + formulation::Symbol, + period::Symbol, + sequence::Symbol) +end + #------------------------------------------------------------------------------- #returm -2REML """ @@ -534,7 +549,7 @@ end function Base.show(io::IO, rbe::RBE) rcoef = coefnames(rbe.rmodel); θ = theta(rbe) - println(io, "Bioequivalence Linear Mixed Effect Model (status: $(Optim.converged(rbe.optim) ? "converged" : printstyled(io, "not converged"; color = :red)))") + print(io, "Bioequivalence Linear Mixed Effect Model (status:"); Optim.converged(rbe.optim) ? print(io,"converged") : printstyled(io, "not converged"; color = :red); println(io, ")") if !isposdef(Symmetric(rbe.H)) printstyled(io, "Hessian not positive defined!"; color = :yellow) println(io, "") diff --git a/src/rbesettings.jl b/src/rbesettings.jl new file mode 100644 index 0000000..6626678 --- /dev/null +++ b/src/rbesettings.jl @@ -0,0 +1,31 @@ +struct RBESettings{T} + g_tol::Real + x_tol::Real + f_tol::Real + iterations::Int + store_trace::Bool + extended_trace::Bool + show_trace::Bool + memopt::Bool + init::Vector{T} + postopt::Bool + vlm::Real + maxopttry::Int + rholink::Symbol + singlim::Real + function RBESettings(;g_tol::Real = 1e-12, x_tol::Real = 0.0, f_tol::Real = 0.0, iterations::Int = 100, + store_trace = false, extended_trace = false, show_trace = false, + memopt = true, + init = zeros(0), + postopt = false, vlm = 1.0, maxopttry = 100, + rholink = :psigmoid, + singlim = 1e-10) + new{typeof(init)}(g_tol, x_tol, f_tol, iterations, + store_trace, extended_trace, show_trace, + memopt, + init, + postopt, vlm, maxopttry, + rholink, + singlim) + end +end diff --git a/src/utils.jl b/src/utils.jl index f66f0cb..9fcc462 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -258,3 +258,8 @@ function sbjnbyf(df, subj, fac, f) return length(sbj) end #------------------------------------------------------------------------------- +function intravar(rbe::RBE) + terms = rbe.rmodel.f.rhs.terms[2].contrasts.termnames + θ = theta(rbe) + return Dict(terms .=> θ[1:2]) +end diff --git a/test/test.jl b/test/test.jl index edabf05..a2c7cae 100644 --- a/test/test.jl +++ b/test/test.jl @@ -500,6 +500,7 @@ end end @testset " # Simulation " begin + #Power simulation by task io = IOBuffer() task = ReplicateBE.randrbetask( ; @@ -517,4 +518,18 @@ end ) pow = ReplicateBE.simulation(task; io = io, num = 100, seed = 1234, verbose = true) @test pow.result == 0.04 + + #Custom simulation + + task = ReplicateBE.randrbetask(;n=24, + sequence=[1,1], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], + inter=[0.04, 0.04, 1.0], intra=[0.02, 0.02], + intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], + formcoef = [0.0, log(0.8)], seed = 0, dropobs = 0) + out = zeros(Float64, 0) + function simfunc!(out, be) + push!(out, ReplicateBE.theta(be)[5]) + end + result = ReplicateBE.simulation!(task, out, simfunc!; num = 10, seed = 10, verbose = false) + @test mean(result) ≈ 0.8904498245891423 end