Skip to content

Commit

Permalink
1.0.1 (#2)
Browse files Browse the repository at this point in the history 
- v1.0.1 (hotfix)
    * increase stability
    * simulation stability
    * change default options: g_tol = 1e-12
    * vlm = 1.0 in rbe() & rbe!()
    * initvar() changes
    * theta() - return final estimates
  • Loading branch information
@PharmCat
PharmCat authored Nov 15, 2019
1 parent 7d4e6b5 commit 1199e94
Show file tree
Hide file tree
Showing 8 changed files with 116 additions and 57 deletions.
4 changes: 3 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,10 @@ name = "ReplicateBE"
uuid = "671d9d50-c343-11e9-1a9c-fdd992682823"
keywords = ["bioequivalence", "mixedmodel"]
desc = "Mixed model solution for replicate designed bioequivalence study."
version = "1.0.0"
version = "1.0.1"

[deps]
CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Expand All @@ -27,6 +28,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
test = ["CSV", "Test", "StatsBase"]

[compat]
CategoricalArrays = "0.7"
julia = "1.0, 1.1, 1.2"
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"
Expand Down
8 changes: 8 additions & 0 deletions chagelog.md
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,11 @@
- v1.0.1 (hotfix)
* increase stability
* simulation stability
* change default options: g_tol = 1e-12
* vlm = 1.0 in rbe() & rbe!()
* initvar() changes
* theta() - return final estimates

- v1.0.0 (unstable)

* validation
Expand Down
2 changes: 1 addition & 1 deletion src/ReplicateBE.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

module ReplicateBE

using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches
using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches, CategoricalArrays

export RBE, RandRBEDS, rbe, rbe!, reml2, nobs, coef, stderror, dof, coefnum, fixed, theta, typeiii, design, show, confint, contrast, estimate, optstat, randrbeds
import Base.show
Expand Down
25 changes: 10 additions & 15 deletions src/generalfunc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,7 @@ end
"""
REML estimation with β recalculation
"""
function remlb(yv, Zv, p, Xv, θvec, β; memopt::Bool = true)
function remlb(yv, Zv, p, Xv, θvec, β; memopt::Bool = true, backval = [])
maxobs = maximum(length.(yv))
#some memory optimizations to reduse allocations
mXviV = Array{Array{eltype(θvec), 2}, 1}(undef, maxobs)
Expand Down Expand Up @@ -252,6 +252,10 @@ function remlb(yv, Zv, p, Xv, θvec, β; memopt::Bool = true)
@inbounds r = yv[i] - Xv[i]*βt
@inbounds θ3 += r'*iVv[i]*r
end

if length(backval) == length(θvec)
copy!(backval, θvec)
end
return -(θ1 + logdet(θ2) + θ3 + c)/2
end
"""
Expand Down Expand Up @@ -323,23 +327,14 @@ function initvar(df, dv, fac, sbj)
u = unique(df[:, sbj])
f = unique(df[:, fac])
fv = Array{Float64, 1}(undef, 0)
sb = Array{Float64, 1}(undef, 0)
for i in f
fvd = Array{Float64, 1}(undef, 0)
v = findall(x->x==i, df[:, fac])
for r in v
push!(fvd, df[r, dv])
end
if length(fvd) > 1 push!(fv, var(fvd)) end
push!(fv,var(df[df[!, fac] .== i, dv]))
end
sv = Array{Float64, 1}(undef, 0)
for i in u
fvd = Array{Float64, 1}(undef, 0)
v = findall(x->x==i, df[:, sbj])
for r in v
push!(fvd, df[r, dv])
end
if length(fvd) > 1 push!(sv, var(fvd)) end
sv = var(df[df[!, sbj] .== i, dv])
if sv > 0 push!(sb, sv) end
end
push!(fv, mean(filter!(x -> x !== NaN, sv)))
push!(fv, mean(sb))
return fv
end
69 changes: 50 additions & 19 deletions src/randrbeds.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,9 +126,11 @@ function randrbeds(n::Int, sequence::Vector,
θinter::Vector, θintra::Vector,
intercept::Real, seqcoef::Vector, periodcoef::Vector, formcoef::Vector,
dropsubj::Float64, dropobs::Int, seed)

rng = MersenneTwister()
if seed == 0 Random.seed!(rng) else Random.seed!(seed) end
if seed != 0
rng = MersenneTwister(seed)
else
rng = MersenneTwister()
end

r = n/sum(sequence)
sn = Array{Int, 1}(undef, length(sequence))
Expand Down Expand Up @@ -162,27 +164,30 @@ function randrbeds(n::Int, sequence::Vector,
Mv[i] = zeros(pnum) .+ intercept .+ seqcoef[i] + periodcoef + Zv[i]*formcoef
end

subjds = DataFrame(subject = String[], formulation = String[], period = String[], sequence = String[], var = Float64[])
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
mvnorm = MvNormal(PDMat(Vv[i]))
for sis = 1:sn[i]
subjmx[:, 1] .= string(subj)
subjmx[:, 2] = string.(design[i,:])
subjmx[:, 3] = string.(collect(1:pnum))
subjmx[:, 1] .= subj
subjmx[:, 2] = design[i,:]
subjmx[:, 3] = collect(1:pnum)
subjmx[:, 4] .= sqname[i]
subjmx[:, 5] = rand(MvNormal(PDMat(Vv[i]))) + Mv[i]
subjmx[:, 5] = rand(rng, MvNormal(PDMat(Vv[i]))) + Mv[i]
subj += 1
for c = 1:pnum
push!(subjds, subjmx[c, :])
end
end
end
if dropobs > 0 && dropobs < size(subjds, 1)
dellist = sample(1:size(subjds, 1), dropobs, replace = false)
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
"""
Expand Down Expand Up @@ -222,32 +227,51 @@ Count successful BE outcomes.
"""
function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l = log(0.8), u = log(1.25), seed = 0)
task.seed = 0
rng = MersenneTwister()
if seed == 0 Random.seed!(rng) else Random.seed!(seed) end
seeds = rand(UInt32, num)
#rng = MersenneTwister()
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
cnt = 0
if verbose
printstyled(io, "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
task.seed = seeds[i]
rds = ReplicateBE.randrbeds(task)
be = ReplicateBE.rbe(rds, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence)

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]
#!
if verbose
if !optstat(be) printstyled(io, "Iteration: $i, seed $(seeds[i]): unconverged! \n"; color = :yellow) end
if be.detH <= 0 printstyled(io, "Iteration: $i, seed $(seeds[i]): Hessian not positive defined! \n"; color = :yellow) end
if !optstat(be) printstyled(io, "Iteration: ", i, ", seed ", seeds[i], ": unconverged! \n"; color = :yellow) end
if be.detH <= 0
printstyled(io, "Iteration: ", i, ", seed ", seeds[i], ": Hessian not positive defined! \n"; color = :yellow)
end
end
if ll > l && ul < u
#println(io, "Seed $(task.seed) LL $(ll) UL $(ul)")
cnt += 1
end
#!
Expand All @@ -259,15 +283,22 @@ function simulation(task::RandRBEDS; io = stdout, verbose = false, num = 100, l
n = 0
end
n += 1

catch
err += 1
printstyled(io, "Iteration: $i, seed $(seeds[i]): $(err): ERROR! \n"; color = :red)
end

end
return RBEDSSimResult(seed, num, seeds, cnt/(num - err))
end

function Base.show(io::IO, obj::RBEDSSimResult)
if isa(obj.seed, Array)
println(io, "Simulation seed: Array")
else
println(io, "Seed: $(obj.seed)")
end
println(io, "Seed: $(obj.seed)")
println(io, "Number: $(obj.num)")
println(io, "Result: $(obj.result)")
Expand Down
53 changes: 38 additions & 15 deletions src/rbe.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,16 +119,29 @@ function rbe(df; dvar::Symbol,
formulation::Symbol,
period::Symbol,
sequence::Symbol,
g_tol::Float64 = 1e-8, x_tol::Float64 = 0.0, f_tol::Float64 = 0.0, iterations::Int = 100,
g_tol::Float64 = 1e-12, x_tol::Float64 = 0.0, f_tol::Float64 = 0.0, iterations::Int = 100,
store_trace = false, extended_trace = false, show_trace = false,
memopt = true,
init = [],
postopt = false, vlm = 0.8, maxopttry = 50, rhoadjstep = 0.15)
postopt = false, vlm = 1.0, maxopttry = 50, rhoadjstep = 0.15)
#Check
if any(x -> x ∉ names(df), [subject, formulation, period, sequence]) throw(ArgumentError("Names not found in DataFrame!")) end
if !(eltype(df[!,dvar]) <: Real)
@warn "Responce variable ∉ Real!"
end
if !(typeof(df[!,subject]) <: CategoricalArray)
@warn "Subject variable not Categorical, use rbe!()!"
end
if !(typeof(df[!,formulation]) <: CategoricalArray)
@warn "Formulation variable not Categorical, use rbe!()!"
end
if !(typeof(df[!,period]) <: CategoricalArray)
@warn "Period variable not Categorical, use rbe!()!"
end
if !(typeof(df[!,sequence]) <: CategoricalArray)
@warn "Sequence variable not Categorical, use rbe!()!"
end

#Model
Xf = @eval(@formula($dvar ~ $sequence + $period + $formulation))
Zf = @eval(@formula($dvar ~ 0 + $formulation))
Expand Down Expand Up @@ -165,9 +178,10 @@ function rbe(df; dvar::Symbol,
θvec0 = init
else
iv = initvar(df, dvar, formulation, subject)
if iv[1] < iv[3] || iv[2] < iv[3] iv[1] = iv[2] = 2*iv[3] end
θvec0 = rvarlink([iv[3], iv[3], iv[1]-iv[3], iv[2]-iv[3], 0.05], vlm)
iv = iv .+ eps()
θvec0 = [iv[3], iv[3], iv[1], iv[2], 0.05]
end
θvec0 = rvarlink(θvec0, vlm)
#Prelocatiom for G, R, V, V⁻¹ matrices
G = zeros(2, 2)
Rv = Array{Array{Float64,2}, 1}(undef, n)
Expand All @@ -181,26 +195,29 @@ function rbe(df; dvar::Symbol,
td = TwiceDifferentiable(x -> -2*remlb(yv, Zv, p, Xv, varlink(x, vlm), β; memopt = memopt), θvec0; autodiff = :forward)
opttry = true
optnum = 0
rng = MersenneTwister(hash(θvec0))
while opttry
try
O = optimize(td, θvec0, method=Newton(), g_tol=g_tol, x_tol=x_tol, f_tol=f_tol, allow_f_increases = true, store_trace = store_trace, extended_trace = extended_trace, show_trace = show_trace)
opttry = false
catch
θvec0 = rvarlink(varlink(θvec0, vlm) .+ (rand(rng)-0.5)/10 .* varlink(θvec0, vlm) .+ eps(), vlm)
θvec0[5] = θvec0[5] - rhoadjstep
end
optnum += 1
if optnum > maxopttry
opttry = false
throw(ErrorException("Initial values faild! Iteration $(optnum), θvec0[5] = $(θvec0[5])."))
throw(ErrorException("Optimization faild! Iteration $(optnum), θvec = $(θvec0)"))
end

end
θ = Optim.minimizer(O)
#Get reml
remlv = -reml2b!(yv, Zv, p, n, N, Xv, G, Rv, Vv, iVv, varlink(θ, vlm), β, memalloc)
#Post optimization
if postopt
pO = O
od = OnceDifferentiable(x -> -2*reml(yv, Zv, p, Xv, varlink(x, vlm), β; memopt = memopt), θvec0; autodiff = :forward)
od = OnceDifferentiable(x -> -2*reml(yv, Zv, p, Xv, varlink(x, vlm), β; memopt = memopt), θ; autodiff = :forward)
method = BFGS(linesearch = LineSearches.HagerZhang(), alphaguess = LineSearches.InitialStatic())
O = optimize(od, [-Inf, -Inf, -Inf, -Inf, -Inf], [Inf, Inf, Inf, Inf, Inf], θ, Fminbox(method), Optim.Options(g_tol=g_tol, x_tol=x_tol, f_tol=f_tol))
θ = copy(Optim.minimizer(O))
Expand Down Expand Up @@ -265,7 +282,7 @@ function rbe(df; dvar::Symbol,
termmodelleveln(MF, formulation),
sbf,
p, zxr)
return RBE(MF, RMF, design, fac, θvec0, vlm, Tuple(θ), remlv, fixed, typeiii, Rv, Vv, G, C, A, H, X, Z, Xv, Zv, yv, dH, pO, O)
return RBE(MF, RMF, design, fac, varlink(θvec0, vlm), vlm, Tuple(varlink(θ, vlm)), remlv, fixed, typeiii, Rv, Vv, G, C, A, H, X, Z, Xv, Zv, yv, dH, pO, O)
end #END OF rbe()
"""
This function apply following code for each factor before executing:
Expand All @@ -285,7 +302,7 @@ function rbe!(df; dvar::Symbol,
formulation::Symbol,
period::Symbol,
sequence::Symbol,
g_tol::Float64 = 1e-8, x_tol::Float64 = 0.0, f_tol::Float64 = 0.0, iterations::Int = 100,
g_tol::Float64 = 1e-12, x_tol::Float64 = 0.0, f_tol::Float64 = 0.0, iterations::Int = 100,
store_trace = false, extended_trace = false, show_trace = false,
memopt = true,
init = [],
Expand All @@ -296,13 +313,19 @@ function rbe!(df; dvar::Symbol,
@warn "Responce variable ∉ Real!"
df[!,dvar] = float.(df[!,dvar])
end

categorical!(df, subject);
categorical!(df, formulation);
categorical!(df, period);
categorical!(df, sequence);
if !(typeof(df[!,subject]) <: CategoricalArray)
categorical!(df, subject);
end
if !(typeof(df[!,formulation]) <: CategoricalArray)
categorical!(df, formulation);
end
if !(typeof(df[!,period]) <: CategoricalArray)
categorical!(df, period);
end
if !(typeof(df[!,sequence]) <: CategoricalArray)
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,
Expand Down Expand Up @@ -520,7 +543,7 @@ end
#-------------------------------------------------------------------------------
function Base.show(io::IO, rbe::RBE)
rcoef = coefnames(rbe.rmodel);
θ = varlink(collect(rbe.θ), rbe.vlm)
θ = theta(rbe)
println(io, "Bioequivalence Linear Mixed Effect Model (status: $(Optim.converged(rbe.optim) ? "converged" : printstyled(io, "not converged"; color = :red)))")
if rbe.detH <= 0.0
printstyled(io, "Hessian not positive!"; color = :yellow)
Expand Down
8 changes: 4 additions & 4 deletions src/utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ DF for one-dimetion case:
df = \\frac{2(LCL')^{2}}{g'Ag}
```
where ``A = 2H``
where ``A = -2H``
where ``g = \\triangledown _{\\theta}(LC^{-1}L')``
Expand All @@ -41,7 +41,7 @@ function contrast(rbe::RBE, L::Matrix; numdf = 0, name = "Contrast", memopt = tr
lcl = L*rbe.C*L'
lclr = rank(lcl)
F = β'*L'*inv(lcl)*L*β/lclr
θ = theta(rbe)
θ = rvarlink(theta(rbe), rbe.vlm)

if numdf == 0 numdf = rank(L) end

Expand Down Expand Up @@ -96,7 +96,7 @@ where
C = \\sum_{i=1}^{n} X_i'V_i^{-1}X_i
```
where ``A = 2H``
where ``A = -2H``
where ``g = \\triangledown _{\\theta}(LC^{-1}L')``
Expand Down Expand Up @@ -126,7 +126,7 @@ function estimate(rbe::RBE, L::Matrix; df = :sat, name = "Estimate", memopt = tr
se = sqrt((lcl)[1])
#F = β'*L'*inv(lcl)*L*β/lclr
if df == :sat
θ = theta(rbe)
θ = rvarlink(theta(rbe), rbe.vlm)
g = ForwardDiff.gradient(x -> lclgf(L, L', rbe.Xv, rbe.Zv, varlink(x, rbe.vlm); memopt = memopt), θ)
df = 2*((lcl)[1])^2/(g'*(rbe.A)*g)
elseif df == :cont
Expand Down
Loading

2 comments on commit 1199e94

@PharmCat
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

Release notes:

  • v1.0.1 (hotfix)
    • increase stability
    • simulation stability
    • change default options: g_tol = 1e-12
    • vlm = 1.0 in rbe() & rbe!()
    • initvar() changes
    • theta() - return final estimates

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/5452

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if Julia TagBot is installed, or can be done manually through the github interface, or via:

git tag -a v1.0.1 -m "<description of version>" 1199e9415e528e30c87ac95a2b9a65076237d29f
git push origin v1.0.1

Please sign in to comment.