Skip to content

Commit

Permalink
v1.0.5
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Dec 2, 2019
1 parent 9a6b6ff commit 1b446fe
Show file tree
Hide file tree
Showing 10 changed files with 132 additions and 109 deletions.
3 changes: 2 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.4"
version = "1.0.5"

[deps]
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
Expand Down
20 changes: 10 additions & 10 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,18 +51,18 @@ Bioequivalence Linear Mixed Effect Model (status: converged)
-2REML: 329.257 REML: -164.629
Fixed effect:
──────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────
Effect Value SE F DF t P|t|
──────────────────────────────────────────────────────────────────────────────────────────────────────
(Intercept) 4.42158 0.119232 1375.21 68.6064 37.0838 4.0203900000000003e-47*
sequence: 2 0.360591 0.161776 4.96821 62.0 2.22895 0.0294511*
period: 2 0.027051 0.0533388 0.257206 122.73 0.507155 0.612956
period: 3 -0.00625777 0.0561037 0.012441 153.634 -0.111539 0.911334
period: 4 0.036742 0.0561037 0.428886 153.634 0.654894 0.513515
formulation: 2 0.0643404 0.0415345 2.39966 62.0 1.54908 0.126451
──────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────
(Intercept) 4.42158 0.119232 1375.21 68.6064 37.0838 4.02039E-47*
sequence: 2 0.360591 0.161776 4.96821 62.0 2.22895 0.0294511*
period: 2 0.027051 0.0533388 0.257206 122.73 0.507155 0.612956
period: 3 -0.00625777 0.0561037 0.012441 153.634 -0.111539 0.911334
period: 4 0.036742 0.0561037 0.428886 153.634 0.654894 0.513515
formulation: 2 0.0643404 0.0415345 2.39966 62.0 1.54908 0.126451
───────────────────────────────────────────────────────────────────────────────────────────
Intra-individual variance:
formulation: 1 0.108629 CVᵂ: 33.87 %
formulation: 1 0.108629 CVᵂ: 33.87 %
formulation: 2 0.0783544 CVᵂ: 28.55 %
Inter-individual variance:
Expand Down
2 changes: 2 additions & 0 deletions chagelog.md
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
- v1.0.5
* inv(H) used where possible
* DF for typeiii checks
* Type optimization
* Output fix

- v1.0.4
* CSV test dataset
Expand Down
22 changes: 11 additions & 11 deletions docs/src/details.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,24 +143,24 @@ Bioequivalence Linear Mixed Effect Model (status: converged)
-2REML: 329.257 REML: -164.629
Fixed effect:
──────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────
Effect Value SE F DF t P|t|
──────────────────────────────────────────────────────────────────────────────────────────────────────
(Intercept) 4.42158 0.119232 1375.21 68.6064 37.0838 4.0203900000000003e-47*
sequence: 2 0.360591 0.161776 4.96821 62.0 2.22895 0.0294511*
period: 2 0.027051 0.0533388 0.257206 122.73 0.507155 0.612956
period: 3 -0.00625777 0.0561037 0.012441 153.634 -0.111539 0.911334
period: 4 0.036742 0.0561037 0.428886 153.634 0.654894 0.513515
formulation: 2 0.0643404 0.0415345 2.39966 62.0 1.54908 0.126451
──────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────────────────
(Intercept) 4.42158 0.119232 1375.21 68.6064 37.0838 4.02039E-47*
sequence: 2 0.360591 0.161776 4.96821 62.0 2.22895 0.0294511*
period: 2 0.027051 0.0533388 0.257206 122.73 0.507155 0.612956
period: 3 -0.00625777 0.0561037 0.012441 153.634 -0.111539 0.911334
period: 4 0.036742 0.0561037 0.428886 153.634 0.654894 0.513515
formulation: 2 0.0643404 0.0415345 2.39966 62.0 1.54908 0.126451
───────────────────────────────────────────────────────────────────────────────────────────
Intra-individual variance:
formulation: 1 0.108629 CVᵂ: 33.87 %
formulation: 1 0.108629 CVᵂ: 33.87 %
formulation: 2 0.0783544 CVᵂ: 28.55 %
Inter-individual variance:
formulation: 1 0.377846
formulation: 2 0.421356
ρ: 0.980288 Cov: 0.391143
ρ: 0.980288 Cov: 0.391143
Confidence intervals(90%):
formulation: 1 / formulation: 2
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, CategoricalArrays
using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, PDMats, Optim, LineSearches, CategoricalArrays, Printf

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
12 changes: 6 additions & 6 deletions src/memalloc.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
struct MemAlloc{T <: AbstractFloat}
mem1::Array{Array{T, 2}, 1}
mem2::Array{Array{T, 2}, 1}
mem3::Array{Array{T, 1}, 1}
mem1::Vector{Matrix{T}}
mem2::Vector{Matrix{T}}
mem3::Vector{Vector{T}}
#mem4::Array{Float64, 2}
function MemAlloc(p, zs, yv::Vector{Vector{T}}) where T <: AbstractFloat
maxobs = maximum(length.(yv))
yvtype = eltype(yv[1])
memc1 = Array{Array{yvtype, 2}, 1}(undef, maxobs)
memc2 = Array{Array{yvtype, 2}, 1}(undef, maxobs)
memc3 = Array{Array{yvtype, 1}, 1}(undef, maxobs)
memc1 = Vector{Matrix{yvtype}}(undef, maxobs)
memc2 = Vector{Matrix{yvtype}}(undef, maxobs)
memc3 = Vector{Vector{yvtype}}(undef, maxobs)
for i = 1:maxobs
memc1[i] = zeros(i, zs)
memc2[i] = zeros(p, i)
Expand Down
4 changes: 2 additions & 2 deletions src/randrbeds.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ mutable struct RandRBEDS
seqcoef::Vector
periodcoef::Vector
formcoef::Vector
dropsubj::Float64 #Deprecated
dropsubj::Real #Deprecated
dropobs::Int
seed
function RandRBEDS(;n=24, sequence=[1,1],
Expand Down Expand Up @@ -127,7 +127,7 @@ function randrbeds(n::Int, sequence::Vector,
design::Matrix,
θinter::Vector, θintra::Vector,
intercept::Real, seqcoef::Vector, periodcoef::Vector, formcoef::Vector,
dropsubj::Float64, dropobs::Int, seed)
dropsubj::Real, dropobs::Int, seed)
if seed != 0
rng = MersenneTwister(seed)
else
Expand Down
155 changes: 81 additions & 74 deletions src/rbe.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,62 +3,62 @@
#
"""
```julia
struct RBE
model::ModelFrame #Model frame
rmodel::ModelFrame #Random effect model
design::Design #Design description
factors::Array{Symbol, 1} #Factor list
θ0::Array{Float64, 1} #Initial variance paramethers
vlm::Real
θ::Tuple{Vararg{Float64}} #Final variance paramethers
reml::Float64 #-2REML
fixed::EffectTable #Fixed Effect table
typeiii::ContrastTable #Type III table
R::Array{Matrix{Float64},1} #R matrices for each subject
V::Array{Matrix{Float64},1} #V matrices for each subject
G::Matrix{Float64} #G matrix
C::Matrix{Float64} #C var(β) p×p variance-covariance matrix
A::Matrix{Float64} #asymptotic variance-covariance matrix of b θ
H::Matrix{Float64} #Hessian matrix
X::Matrix #Matrix for fixed effects
Z::Matrix #Matrix for random effects
Xv::Array{Matrix{Float64},1} #X matrices for each subject
Zv::Array{Matrix{Float64},1} #Z matrices for each subject
yv::Array{Array{Float64, 1},1} #responce vectors for each subject
detH::Float64 #Hessian determinant
preoptim::Union{Optim.MultivariateOptimizationResults, Nothing} #Pre-optimization result object
optim::Optim.MultivariateOptimizationResults #Optimization result object
end
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
R::Vector{Matrix{T}} # R matrices for each subject
V::Vector{Matrix{T}} # V matrices for each subject
G::Matrix{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
Xv::Vector{Matrix{T}} # X matrices for each subject
Zv::Vector{Matrix{T}} # Z matrices for each subject
yv::Vector{Vector{T}} # responce vectors for each subject
detH::T # Hessian determinant
preoptim::Union{Optim.MultivariateOptimizationResults, Nothing} # Pre-optimization result object
optim::Optim.MultivariateOptimizationResults # Optimization result object
end
```
Replicate bioequivalence structure.
"""
struct RBE
model::ModelFrame #Model frame
rmodel::ModelFrame #Random effect model
struct RBE{T <: AbstractFloat}
model::ModelFrame # Model frame
rmodel::ModelFrame # Random effect model
design::Design
factors::Array{Symbol, 1} #Factor list
θ0::Array{Float64, 1} #Initial variance paramethers
vlm::Real
θ::Tuple{Vararg{Float64}} #Final variance paramethers
reml::Float64 #-2REML
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
R::Array{Matrix{Float64},1} #R matrices for each subject
V::Array{Matrix{Float64},1} #V matrices for each subject
G::Matrix{Float64} #G matrix
C::Matrix{Float64} #C var(β) p×p variance-covariance matrix
A::Matrix{Float64} #asymptotic variance-covariance matrix ofb θ
H::Matrix{Float64} #Hessian matrix
X::Matrix #Matrix for fixed effects
Z::Matrix #Matrix for random effects
Xv::Array{Matrix{Float64},1} #X matrices for each subject
Zv::Array{Matrix{Float64},1} #Z matrices for each subject
yv::Array{Array{Float64, 1},1} #responce vectors for each subject
detH::Float64 #Hessian determinant
preoptim::Union{Optim.MultivariateOptimizationResults, Nothing} #Pre-optimization result object
optim::Optim.MultivariateOptimizationResults #Optimization result object
R::Vector{Matrix{T}} # R matrices for each subject
V::Vector{Matrix{T}} # V matrices for each subject
G::Matrix{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
Xv::Vector{Matrix{T}} # X matrices for each subject
Zv::Vector{Matrix{T}} # Z matrices for each subject
yv::Vector{Vector{T}} # responce vectors for each subject
detH::T # Hessian determinant
preoptim::Union{Optim.MultivariateOptimizationResults, Nothing} # Pre-optimization result object
optim::Optim.MultivariateOptimizationResults # Optimization result object
end

"""
Expand Down Expand Up @@ -111,7 +111,6 @@ V_{i} = Z_{i}GZ_i'+R_{i}
* memopt = true - memory optimization (function cache)
* init = [] - initial variance paremeters
* postopt = false - post optimization
* maxopttry = 100 - maximum attempts to optimize
"""
function rbe(df; dvar::Symbol,
Expand All @@ -132,7 +131,6 @@ function rbe(df; dvar::Symbol,
if !(eltype(df[!,dvar]) <: AbstractFloat)
@warn "Responce variable ∉ Float!"
end
vartype = eltype(df[!,dvar])
if !(typeof(df[!,subject]) <: CategoricalArray)
@warn "Subject variable not Categorical, use rbe!()!"
end
Expand Down Expand Up @@ -202,14 +200,14 @@ function rbe(df; dvar::Symbol,
end


θvec0 = rvarlink(θvec0, vlm)
θvec0 = rvarlink(θvec0, vlm)
#Prelocatiom for G, R, V, V⁻¹ matrices
G = zeros(2, 2)
Rv = Array{Array{vartype,2}, 1}(undef, n)
Vv = Array{Array{vartype,2}, 1}(undef, n)
iVv = Array{Array{vartype,2}, 1}(undef, n)
matvecz!(Rv, Zv)
matvecz!(Vv, Zv)
G = zeros(2, 2)
Rv = Vector{Matrix{eltype(df[!,dvar])}}(undef, n)
Vv = Vector{Matrix{eltype(df[!,dvar])}}(undef, n)
iVv = Vector{Matrix{eltype(df[!,dvar])}}(undef, n)
matvecz!(Rv, Zv)
matvecz!(Vv, Zv)
matvecz!(iVv, Zv)
#Optimization
pO = nothing
Expand Down Expand Up @@ -275,11 +273,11 @@ function rbe(df; dvar::Symbol,

#A = 2 * pinv(H)
C = cmat(Xv, Zv, iVv, θ)
se = Array{vartype, 1}(undef, p)
F = Array{vartype, 1}(undef, p)
df = Array{vartype, 1}(undef, p)
t = Array{vartype, 1}(undef, p)
pval = Array{vartype, 1}(undef, p)
se = Vector{eltype(C)}(undef, p)
F = Vector{eltype(C)}(undef, p)
df = Vector{eltype(C)}(undef, p)
t = Vector{eltype(C)}(undef, p)
pval = Vector{eltype(C)}(undef, p)
for i = 1:p
L = zeros(1, p)
L[i] = 1
Expand All @@ -296,23 +294,32 @@ function rbe(df; dvar::Symbol,
end
fixed = EffectTable(coefnames(MF), β, se, F, df, t, pval)
fac = [sequence, period, formulation]
F = Array{vartype, 1}(undef, length(fac))
df = Array{vartype, 1}(undef, length(fac))
ndf = Array{vartype, 1}(undef, length(fac))
pval = Array{vartype, 1}(undef, length(fac))
F = Vector{eltype(C)}(undef, length(fac))
df = Vector{eltype(C)}(undef, length(fac))
ndf = Vector{eltype(C)}(undef, length(fac))
pval = Vector{eltype(C)}(undef, length(fac))
for i = 1:length(fac)
L = lmatrix(MF, fac[i])
lcl = L*C*L'
lclr = rank(lcl)
F[i] = β'*L'*inv(lcl)*L*β/lclr
if lclr ≥ 2
vm = Array{vartype, 1}(undef, lclr)
vm = Vector{eltype(C)}(undef, lclr)
for i = 1:lclr
g = ForwardDiff.gradient(x -> lclgf(L[i:i,:], L[i:i,:]', Xv, Zv, x; memopt = memopt), θ)
dfi = 2*((L[i:i,:]*C*L[i:i,:]')[1])^2/(g'*A*g)
vm[i] = dfi/(dfi-2)
g = ForwardDiff.gradient(x -> lclgf(L[i:i,:], L[i:i,:]', Xv, Zv, x; memopt = memopt), θ)
dfi = 2*((L[i:i,:]*C*L[i:i,:]')[1])^2/(g'*A*g)
if dfi > 2
vm[i] = dfi/(dfi-2)
else
vm[i] = 0
end
end
E = sum(vm)
if E > lclr
dfi = 2 * E / (E - lclr)
else
dfi = 0
end
dfi = 2*sum(vm)/(sum(vm)-lclr)
else
g = ForwardDiff.gradient(x -> lclgf(L, L', Xv, Zv, x; memopt = memopt), θ)
dfi = 2*((lcl)[1])^2/(g'*(A)*g)
Expand Down Expand Up @@ -379,11 +386,11 @@ end
#-------------------------------------------------------------------------------
#returm -2REML
"""
reml2(rbe::RBE, θ::Array{T, 1}) where T <: AbstractFloat
reml2(rbe::RBE, θ::Vector{T}) where T <: AbstractFloat
Returm -2logREML for rbe model with θ variance vector.
"""
function reml2(rbe::RBE, θ::Array{T, 1}) where T <: AbstractFloat
function reml2(rbe::RBE, θ::Vector{T}) where T <: AbstractFloat
return -2*reml(rbe.yv, rbe.Zv, rank(ModelMatrix(rbe.model).m), rbe.Xv, θ, coef(rbe))
end
"""
Expand Down Expand Up @@ -504,7 +511,7 @@ function StatsBase.confint(obj::RBE, alpha::Real; expci::Bool = false, inv::Bool
df = obj.fixed.df
end
end
a = Array{Tuple{Float64, Float64},1}(undef, length(obj.fixed.est))
a = Array{Tuple{eltype(obj.fixed.est), eltype(obj.fixed.est)},1}(undef, length(obj.fixed.est))
for i = 1:length(obj.fixed.est)
a[i] = calcci(obj.fixed.est[i]*ifv, obj.fixed.se[i], df[i], alpha, expci)
end
Expand Down
Loading

2 comments on commit 1b446fe

@PharmCat
Copy link
Owner

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.5
    • inv(H) used where possible
    • DF for typeiii checks
    • Type optimization
    • Output fix

@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/6138

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.5 -m "<description of version>" 1b446fe081db06debbbf5e6e67117d1aff84faca
git push origin v1.0.5

Please sign in to comment.