diff --git a/.gitignore b/.gitignore index 2c5a32f..f98edde 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,5 @@ # Validation validation/ -examples/ docs/build/ docs/img/ diff --git a/.travis.yml b/.travis.yml index 8f2a372..0ac17c6 100644 --- a/.travis.yml +++ b/.travis.yml @@ -25,7 +25,9 @@ after_success: jobs: allow_failures: - - julia: 1.4 + - julia: + - 1.4 + - 1.5 include: - stage: "Documentation" julia: 1.0 @@ -35,3 +37,7 @@ jobs: Pkg.instantiate()' - julia --project=docs/ docs/make.jl after_success: skip + +script: + - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi + - travis_wait 30 julia --project -e 'using Pkg; Pkg.build(); Pkg.test(; coverage=true)' diff --git a/Project.toml b/Project.toml index 2297e87..e4835a7 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.10" +version = "1.0.11" [deps] @@ -19,7 +19,6 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" 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] @@ -31,14 +30,13 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" test = ["CSV", "Test", "StatsBase"] [compat] -CategoricalArrays = "0.7" +CategoricalArrays = "0.7, 0.8" julia = "1.0, 1.1, 1.2, 1.3, 1.4, 1.5" StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33" -DataFrames = "0.19, 0.20" +DataFrames = "0.19, 0.20, 0.21" StatsModels = "0.6" Distributions = "0.20, 0.21, 0.22, 0.23" -PDMats = "0.9, 0.10" ForwardDiff = "0.10" -Optim = "0.19, 0.20, 0.21" +Optim = "0.19, 0.20, 0.21, 0.22, 1.0, 1.1, 1.2" LineSearches = "7.0" StaticArrays = "0.11, 0.12" diff --git a/chagelog.md b/chagelog.md index 1d3614b..1d852d8 100644 --- a/chagelog.md +++ b/chagelog.md @@ -1,3 +1,10 @@ +- v1.0.11 + * Optimizations + * Cosmetics + * PDMats exclude + * Fix experimental sigmoid function + * Add in initial variance estimate (OLS via QR decomposition) function + - v1.0.10 * Cosmetics * StatsBase.coeftable diff --git a/examples/benchmark.jl b/examples/benchmark.jl new file mode 100644 index 0000000..44b0a33 --- /dev/null +++ b/examples/benchmark.jl @@ -0,0 +1,122 @@ +using BenchmarkTools, CSV, DataFrames +using ReplicateBE + +datapath = joinpath(dirname(pathof(ReplicateBE)))*"/../test/csv/df6.csv" + +df6 = CSV.File(datapath) |> DataFrame + +#@benchmark be = ReplicateBE.rbe!(df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) +@benchmark be = ReplicateBE.rbe!($df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) seconds = 25 samples = 200 evals = 5 + +#= +julia> @benchmark be = ReplicateBE.rbe!(datadf, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) +BenchmarkTools.Trial: + memory estimate: 6.15 MiB + allocs estimate: 66530 + -------------- + minimum time: 24.307 ms (0.00% GC) + median time: 27.732 ms (0.00% GC) + mean time: 29.258 ms (3.78% GC) + maximum time: 47.647 ms (0.00% GC) + -------------- + samples: 172 + evals/sample: 1 + +julia> +=# + +#v1.5 +#= +julia> @benchmark be = ReplicateBE.rbe!(df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) +BenchmarkTools.Trial: + memory estimate: 5.40 MiB + allocs estimate: 48564 + -------------- + minimum time: 16.528 ms (0.00% GC) + median time: 18.803 ms (0.00% GC) + mean time: 20.431 ms (4.65% GC) + maximum time: 53.102 ms (50.43% GC) + -------------- + samples: 245 + evals/sample: 1 +=# + +#v1.0.11 Julia 1.4 +#= +julia> @benchmark be = ReplicateBE.rbe!(df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol += 1e-10) +BenchmarkTools.Trial: + memory estimate: 4.34 MiB + allocs estimate: 46670 + -------------- + minimum time: 10.112 ms (0.00% GC) + median time: 18.371 ms (0.00% GC) + mean time: 19.176 ms (2.40% GC) + maximum time: 83.157 ms (74.30% GC) + -------------- + samples: 261 + evals/sample: 1 +=# + +#v1.0.11d Julia 1.5.1 +#= +julia> @benchmark be = ReplicateBE.rbe!(df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol += 1e-10) +BenchmarkTools.Trial: + memory estimate: 4.48 MiB + allocs estimate: 43228 + -------------- + minimum time: 17.604 ms (0.00% GC) + median time: 19.636 ms (0.00% GC) + mean time: 20.833 ms (2.70% GC) + maximum time: 47.499 ms (0.00% GC) + -------------- + samples: 240 + evals/sample: 1 +=# + +#= +julia> @benchmark be = ReplicateBE.rbe!($df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) seconds = 25 samples = 200 evals = 5 +BenchmarkTools.Trial: + memory estimate: 4.47 MiB + allocs estimate: 43136 + -------------- + minimum time: 17.257 ms (0.00% GC) + median time: 19.049 ms (0.00% GC) + mean time: 19.897 ms (3.34% GC) + maximum time: 29.502 ms (32.83% GC) + -------------- + samples: 200 + evals/sample: 5 +=# + +#= +julia> @benchmark be = ReplicateBE.rbe!($df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) seconds = 25 samples = 200 evals = 5 +BenchmarkTools.Trial: + memory estimate: 4.44 MiB + allocs estimate: 43179 + -------------- + minimum time: 16.638 ms (0.00% GC) + median time: 19.191 ms (0.00% GC) + mean time: 20.705 ms (2.86% GC) + maximum time: 48.517 ms (14.31% GC) + -------------- + samples: 200 + evals/sample: 5 +=# + +#initvar2 +#= +julia> @benchmark be = ReplicateBE.rbe!($df6, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10) seconds = 25 samples = 200 evals = 5 +BenchmarkTools.Trial: + memory estimate: 4.97 MiB + allocs estimate: 48656 + -------------- + minimum time: 18.353 ms (0.00% GC) + median time: 21.240 ms (0.00% GC) + mean time: 22.810 ms (2.87% GC) + maximum time: 43.109 ms (10.80% GC) + -------------- + samples: 200 + evals/sample: 5 +=# diff --git a/src/ReplicateBE.jl b/src/ReplicateBE.jl index 8955758..7dd37c3 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, StaticArrays +using DataFrames, Distributions, StatsModels, StatsBase, ForwardDiff, LinearAlgebra, Random, 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 diff --git a/src/algebra.jl b/src/algebra.jl index 22d9933..d4efe0a 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -1,11 +1,11 @@ """ A' * B * A -> + θ (cache) """ -function mulαtβαinc!(θ, A::AbstractMatrix, B::Matrix, c) +function mulαtβαinc!(θ, A::AbstractMatrix, B, c) q = size(B, 1) p = size(A, 2) for i = 1:p - c .= 0 + fill!(c, zero(eltype(θ))) for n = 1:q for m = 1:q @inbounds c[n] += B[m, n] * A[m, i] @@ -25,9 +25,9 @@ A' * B * A -> + θ function mulαtβαinc!(θ, A::AbstractMatrix, B::AbstractMatrix) q = size(B, 1) p = size(A, 2) - c = zeros(eltype(B), q) + c = zeros(eltype(θ), q) for i = 1:p - c .= 0 + fill!(c, zero(eltype(θ))) for n = 1:q for m = 1:q @inbounds c[n] += B[m, n] * A[m, i] @@ -50,7 +50,7 @@ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::AbstractVe q = size(B, 1) p = size(A, 2) for i = 1:p - c .= 0 + fill!(c, zero(eltype(θ))) for n = 1:q for m = 1:q @inbounds c[n] += B[m, n] * A[m, i] @@ -69,11 +69,12 @@ end """ (y - X * β)' * V * (y - X * β) (cache) """ -function mulθ₃(y::Vector, X::AbstractMatrix, β::Vector, V::AbstractMatrix, c) +function mulθ₃(y, X, β, V, c) q = size(V, 1) p = size(X, 2) - θ = 0 - c .= 0 + θ = zero(eltype(β)) + #c .= 0 + fill!(c, zero(eltype(β))) for n = 1:q for m = 1:p @inbounds c[n] += X[n, m] * β[m] @@ -90,27 +91,27 @@ end """ A * B * A' + C """ -function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) +function mulαβαtc(A, B, C)::Symmetric q = size(B, 1) p = size(A, 1) c = zeros(eltype(B), q) mx = zeros(eltype(B), p, p) for i = 1:p - c .= 0 - for n = 1:q - for m = 1:q + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q @inbounds c[n] += A[i, m] * B[n, m] end end - for n = 1:p - for m = 1:q + @simd for n = i:p + @simd for m = 1:q @inbounds mx[i, n] += A[n, m] * c[m] end + @inbounds mx[i, n] += C[i, n] end end - mx .+= C - #Symmetric(mx) - #SMatrix{p,p}(mx) + #mx .+= C + Symmetric(mx) end """ A * B * A' + C (cache) @@ -118,26 +119,50 @@ A * B * A' + C (cache) function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c::AbstractVector) q = size(B, 1) p = size(A, 1) - #c = mem.svec[p] mx = zeros(eltype(B), p, p) for i = 1:p - c .= 0 - for n = 1:q - for m = 1:q + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q @inbounds c[n] += A[i, m] * B[n, m] end end - for n = 1:p - for m = 1:q + @simd for n = i:p + @simd for m = 1:q @inbounds mx[i, n] += A[n, m] * c[m] end + @inbounds mx[i, n] += C[i, n] end end - mx .+= C - #SMatrix{p,p,eltype(mx)}(mx) - + #mx .+= C + Symmetric(mx) +end +""" +A * B * A' + Diagonal(A*C) (cache) +""" +function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractVector, c::AbstractVector)::Symmetric + q = size(B, 1) + p = size(A, 1) + mx = zeros(eltype(B), p, p) + for i = 1:p + fill!(c, zero(eltype(c))) + @simd for n = 1:q + @simd for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + @simd for n = i:p + @simd for m = 1:q + @inbounds mx[i, n] += A[n, m] * c[m] + end + end + @simd for m = 1:length(C) + @inbounds mx[i, i] += A[i, m] * C[m] + end + end + #mx + Symmetric(mx) end - function invchol(M) q = size(M, 1) v = zeros(eltype(M), q, q) diff --git a/src/design.jl b/src/design.jl index d146165..691c012 100644 --- a/src/design.jl +++ b/src/design.jl @@ -11,10 +11,6 @@ struct RBEDataStructure mem::MemCache end -function rebuildcache(data, type) - data.mem.svec[1] = zeros(type, data.maxobs) -end - struct RBEResults{T <: AbstractFloat} reml::T # logREML β::Vector{T} # β Vector diff --git a/src/generalfunc.jl b/src/generalfunc.jl index ed1d6a4..c64bb54 100644 --- a/src/generalfunc.jl +++ b/src/generalfunc.jl @@ -5,11 +5,11 @@ """ Make X, Z matrices and vector y for each subject; """ -function sortsubjects(df::DataFrame, sbj::Symbol, X::Matrix, Z::Matrix, y::Vector) +function sortsubjects(df::DataFrame, sbj::Symbol, X::Matrix{T}, Z::Matrix{T}, y::Vector{T}) where T u = unique(df[!, sbj]) - Xa = Vector{Matrix{eltype(y)}}(undef, length(u)) - Za = Vector{Matrix{eltype(y)}}(undef, length(u)) - ya = Vector{Vector{eltype(y)}}(undef, length(u)) + Xa = Vector{SubArray{T}}(undef, length(u)) + Za = Vector{SubArray{T}}(undef, length(u)) + ya = Vector{SubArray{T}}(undef, length(u)) @simd for i = 1:length(u) @inbounds v = findall(x->x==u[i], df[!, sbj]) @inbounds Xa[i] = view(X, v, :) @@ -86,34 +86,31 @@ function mvmat(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, mem, cache)::Ma end end -function mvmatall(G::AbstractMatrix, σ::AbstractVector, Z::AbstractMatrix, mem, cache) - #h = hash(Z) - #if h in keys(cache) +function mvmatall(G::AbstractMatrix, σ::AbstractVector{T}, Z::AbstractMatrix, mem, cache) where T if Z in keys(cache) - #return cache[h] return cache[Z] else - V = mulαβαtc(Z, G, Diagonal(Z*σ), mem) + local V::Symmetric{T,Array{T,2}} + local V⁻¹::Symmetric{T,Array{T,2}} + local log│V│::T + #V = mulαβαtcupd!(mem.mvec, Z, G, σ, mem.svec) + V = mulαβαtc(Z, G, σ, mem) #V = Z * G * Z' + Diagonal(Z*σ) - V⁻¹ = nothing - try - V⁻¹ = invchol(V) - catch e - if typeof(e) <: PosDefException - V⁻¹ = inv(V) - else - throw(e) - end + + if size(V, 1) <= 14 + sV = SHermitianCompact(SMatrix{size(V, 1),size(V, 1)}(V)) + V⁻¹ = Symmetric(inv(sV)) + log│V│ = logdet(sV) + else + V⁻¹ = Symmetric(inv(V)) + log│V│ = logdet(V) end - #V⁻¹ = inv(V) - log│V│ = logdet(V) - cache[Z] = (V, V⁻¹, log│V│) - return V, V⁻¹, log│V│ + cache[Z] = (V⁻¹, log│V│) + return cache[Z] end end -#println("θ₁: ", θ1, " θ₂: ", θ2, " θ₃: ", θ3) -function minv(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, cache::Dict)::Matrix +function minv(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, cache::Dict) #h = hash(M) if Z in keys(cache) #return cache[h] @@ -121,7 +118,8 @@ function minv(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, cache::Dict)::Ma else V = mulαβαtc(Z, G, Diagonal(Z*σ)) #V = Z * G * Z' + Diagonal(Z*σ) - V⁻¹ = invchol(V) + #V⁻¹ = inv(SHermitianCompact(SMatrix{size(V, 1),size(V, 1)}(V))) + V⁻¹ = inv(V) #V⁻¹ = inv(V) #ldV = logdet(V) cache[Z] = V⁻¹ @@ -145,23 +143,23 @@ end """ -2 REML function for ForwardDiff """ -function reml2(data::RBEDataStructure, θ, β::Vector; memopt::Bool = true) +function reml2(data::RBEDataStructure, θ::Vector{T}, β::Vector; memopt::Bool = true) where T #memory optimizations to reduse allocations (cache rebuild) #empty!(data.mem.dict) - rebuildcache(data, promote_type(eltype(data.yv[1]), eltype(θ))) - cache = Dict{Matrix, Tuple{Matrix, Matrix, eltype(θ)}}() + rebuildcache(data, promote_type(eltype(data.yv[1]), T)) + cache = Dict{Matrix, Tuple{Matrix{T}, T}}() #cache = data.mem.dict #--------------------------------------------------------------------------- G = gmat(view(θ, 3:5)) - θ₁ = 0 + θ₁ = zero(T) θ₂ = zeros(promote_type(eltype(data.yv[1]), eltype(θ)), data.p, data.p) - θ₃ = 0 + θ₃ = zero(T) V⁻¹ = nothing #mVec = pmap(x -> mvmatall(G, view(θ,1:2), x, first(data.mem.svec), cache), data.Zv) @simd for i = 1:data.n if MEMOPT && memopt - V, V⁻¹, log│V│ = mvmatall(G, view(θ,1:2), data.Zv[i], first(data.mem.svec), cache) + V⁻¹, log│V│ = mvmatall(G, view(θ,1:2), data.Zv[i], first(data.mem.svec), cache) #V, V⁻¹, log│V│ =mVec[i] θ₁ += log│V│ else @@ -186,24 +184,24 @@ end function reml2bfd(data::RBEDataStructure, θ; memopt::Bool = true) return reml2b(data, θ; memopt = memopt)[1] end -function reml2b(data::RBEDataStructure, θ; memopt::Bool = true) +function reml2b(data::RBEDataStructure, θ::Vector{T}; memopt::Bool = true) where T - rebuildcache(data, promote_type(eltype(data.yv[1]), eltype(θ))) - cache = Dict() + rebuildcache(data, promote_type(eltype(data.yv[1]), T)) + cache = Dict{Matrix, Tuple{Matrix{T}, T}}() #--------------------------------------------------------------------------- G = gmat(view(θ,3:5)) - V⁻¹ = Vector{Matrix{eltype(θ)}}(undef, data.n) # Vector of V⁻¹ matrices - V = nothing - 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) + V⁻¹ = Vector{Matrix{T}}(undef, data.n) # Vector of V⁻¹ matrices + #V = nothing + local log│V│::T # Vector log determinant of V matrix + θ₁ = zero(T) + θ₂ = zeros(promote_type(eltype(first(data.yv)), T), data.p, data.p) + θ₃ = zero(T) + βm = zeros(promote_type(eltype(first(data.yv)), T), data.p) + β = zeros(promote_type(eltype(first(data.yv)), T), data.p) #mVec = map(x -> mvmatall(G, θ[1:2], x, first(data.mem.svec), cache), data.Zv) @simd for i = 1:data.n if MEMOPT && memopt - @inbounds V, V⁻¹[i], log│V│ = mvmatall(G, view(θ,1:2), data.Zv[i], first(data.mem.svec), cache) + @inbounds V⁻¹[i], log│V│ = mvmatall(G, view(θ,1:2), data.Zv[i], first(data.mem.svec), cache) #V, V⁻¹[i], log│V│ = mVec[i] θ₁ += log│V│ else @@ -235,42 +233,45 @@ end """ 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), SVector{length(θ), eltype(θ)}(θ)) - result = Vector{Matrix}(undef, 0) +function cmatgf(data, θ::Vector{T}; memopt::Bool = true) where T + p = size(data.Xv[1], 2) + jC = ForwardDiff.jacobian(x -> cmatvec(data, x; memopt = memopt), θ) + result = Vector{Matrix{T}}(undef, length(θ)) for i in 1:length(θ) - push!(result, reshape(jC[:,i], p, p)) + result[i] = reshape(view(jC, :, i), p, p) # x*x, r)/(length(r) - size(X, 2)) + var1 = zero(eltype(first(yv))) + for i = 1:length(yv) + if length(yv[i]) > 1 + var1 += var(yv[i]) + end + end + var2 = var1/length(yv) + if res - var2 > 0.0 && var2 > 0.0 + return res-var2, var2, b + elseif var2 < 1.0e-6 + return res/2.0, res/2.0, b + else + return (res+var2)/2.0, (res+var2)/2.0, b + end +end #------------------------------------------------------------------------------- function optimcallback(x) false @@ -365,30 +394,40 @@ function vlinkr(σ::T) where T <: Real log(σ) end -function rholinkpsigmoid(ρ, m) - return 1/(1 + exp(ρ * m)) +function rholinkpsigmoid(ρ::T, m) where T <: Real + return 1.0/(1.0 + exp(ρ * m)) end -function rholinkpsigmoidr(ρ, m) - return log(1/ρ - 1)/m +function rholinkpsigmoidr(ρ::T, m) where T <: Real + return log(1.0/ρ - 1.0)/m end -function rholinksigmoid(ρ, m) - return ρ/sqrt(1 + ρ^2) +function rholinksigmoid(ρ::T, m) where T <: Real + return ρ/sqrt(1.0 + ρ^2) end -function rholinksigmoidr(ρ, m) - return sign(ρ)*sqrt(ρ^2/(1 - ρ^2)) +function rholinksigmoidr(ρ::T, m) where T <: Real + return sign(ρ)*sqrt(ρ^2/(1.0 - ρ^2)) end -function rholinksigmoid2(ρ, m) - return atan(ρ) +function rholinksigmoid2(ρ::T, m) where T <: Real + return atan(ρ)/pi*2.0 end -function rholinksigmoid2r(ρ, m) - return tan(ρ) +function rholinksigmoid2r(ρ::T, m) where T <: Real + return tan(ρ*pi/2.0) end function varlinkmap(θ, r1::Union{Int, UnitRange}, r2::Union{Int, UnitRange}, f1::Function, f2::Function) - θl = similar(θ) - θl[r1] = f1.(θ[r1]) - θl[r2] = f2.(θ[r2]) - return θl + #θl = similar(θ) + @inbounds @simd for i in r1 + θ[i] = f1(θ[i]) + end + # + @inbounds @simd for i in r2 + θ[i] = f2(θ[i]) + end + return θ +end + +@inline function lvecupd!(L::AbstractVector, fac) + L .= 0. + L[fac] .= 1. end diff --git a/src/memalloc.jl b/src/memalloc.jl index 0abea34..7fa001a 100644 --- a/src/memalloc.jl +++ b/src/memalloc.jl @@ -1,8 +1,28 @@ struct MemCache svec::Vector{Vector} + mvec::Vector{Matrix} #dict::Dict function MemCache(maxobs) - svec = Vector{Vector}(undef, 1) - new(svec, #=Dict{Matrix, Tuple{Matrix, Matrix, Number}}()=#)::MemCache + #ForwardDiff.Dual{ForwardDiff.Tag{ReplicateBE.var"#23#39"{Bool,Float64,ReplicateBE.RBEDataStructure},Float64}} + svec = Vector{Vector}(undef, 1) + svec[1] = zeros(maxobs) + mvec = Vector{Matrix}(undef, 1) + mvec[1] = zeros(maxobs, maxobs) + new(svec, mvec)::MemCache + end end + +#function Base.Float64(n::T) where T <: ForwardDiff.Dual +# n.value +#end + +function rebuildcache(data, type) + if !(eltype(data.mem.svec[1]) <: type) + data.mem.svec[1] = zeros(type, data.maxobs) + data.mem.mvec[1] = zeros(type, data.maxobs, data.maxobs) + #for i = 1:data.maxobs + # data.mem.mvec[i] = Symmetric(zeros(type, i, i)) + #end + end +end diff --git a/src/randrbeds.jl b/src/randrbeds.jl index 8cccafe..0ab9adc 100644 --- a/src/randrbeds.jl +++ b/src/randrbeds.jl @@ -179,27 +179,52 @@ function randrbeds(n::Int, sequence::Vector, Mv[i] = zeros(pnum) .+ intercept .+ seqcoef[i] + periodcoef + Zv[i]*formcoef end - subjds = DataFrame(subject = Int[], formulation = String[], period = Int[], sequence = String[], var = Float64[]) - subj = 1 - subjmx = Array{Any, 2}(undef, pnum, 5) + #subjds = DataFrame(subject = Int[], formulation = String[], period = Int[], sequence = String[], var = Float64[]) + subjmx1 = Array{Int, 1}(undef, n*pnum) + subjmx2 = Array{String, 1}(undef, n*pnum) + subjmx3 = Array{Int, 1}(undef, n*pnum) + subjmx4 = Array{String, 1}(undef, n*pnum) + subjmx5 = Array{Float64, 1}(undef, n*pnum) + subj = 0 + #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] = 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, :]) - end + s = 1 + (subj - 1)*pnum + e = pnum + (subj - 1)*pnum + subjmx1[s:e] .= subj + subjmx2[s:e] = design[i,:] + subjmx3[s:e] = collect(1:pnum) + subjmx4[s:e] .= sqname[i] + #subjmx[:, 5] = rand(rng, MvNormal(PDMat(Vv[i]))) + Mv[i] + subjmx5[s:e] = rand(rng, MvNormal(Vv[i])) + Mv[i] + + #subjdsm = vcat(subjdsm, subjmx) + #for c = 1:pnum + # 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)) + #subjds = DataFrame(subject = Int.(subjdsm[:,1]), + #formulation = string.(subjdsm[:,2]), + #period = Int.(subjdsm[:,3]), + #sequence = string.(subjdsm[:,4]), + #var = Float64.(subjdsm[:,5])) + + if dropobs > 0 && dropobs < subj + dellist = sort!(sample(rng, 1:subj, dropobs, replace = false)) + deleteat!(subjmx1, dellist) + deleteat!(subjmx2, dellist) + deleteat!(subjmx3, dellist) + deleteat!(subjmx4, dellist) + deleteat!(subjmx5, dellist) + #deleterows!(subjds, sort!(dellist)) end + subjds = DataFrame(subject = subjmx1, + formulation = subjmx2, + period = subjmx3, + sequence = subjmx4, + var = subjmx5) categorical!(subjds, :subject); categorical!(subjds, :formulation); categorical!(subjds, :period); diff --git a/src/rbe.jl b/src/rbe.jl index 6311c8d..9aa8f6c 100644 --- a/src/rbe.jl +++ b/src/rbe.jl @@ -110,7 +110,12 @@ function rbe(df; dvar::Symbol, singlim = 1e-10) #Check - if any(x -> x ∉ names(df), [subject, formulation, period, sequence]) throw(ArgumentError("Names not found in DataFrame!")) end + dfn = names(df) + if eltype(dfn) <: String dfn = Symbol.(dfn) end + if any(x -> x ∉ dfn, [subject, formulation, period, sequence]) + throw(ArgumentError("Names not found in DataFrame! \n Names: $([subject, formulation, period, sequence]) \n DF names: $(names(df))")) + end + if !(eltype(df[!,dvar]) <: AbstractFloat) @warn "Responce variable ∉ AbstractFloat!" end @@ -159,13 +164,17 @@ function rbe(df; dvar::Symbol, data = RBEDataStructure([sequence, period, formulation], Xv, Zv, yv, p, N, n, (N - p) * LOG2PI, maxobs, MemCache(maxobs)) #Calculate initial variance + #iv = initvar2(df, X, yv, dvar, subject) if length(init) == 5 θvec0 = init else + intra = sum(replace!(var.(yv) .* (length.(yv) .- 1), NaN => 0))/(sum(length.(yv))-1) iv = initvar(df, dvar, formulation) iv = iv .+ eps() θvec0 = [intra, intra, iv[1], iv[2], 0.5] + + #θvec0 = [iv[2], iv[2], iv[1], iv[1], 0.5] end #Variance link function @@ -225,7 +234,7 @@ function rbe(df; dvar::Symbol, #H = Calculus.hessian(x -> reml2(data, x, β), θ) # If no varlink using can be obtained from optim results #H = O.trace[end].metadata["h(x)"] - + #print(H) A = nothing #If rho is near to 1.0 it leads to singular hessian matrix, and rho should be removed from variance-covariance matrix #It can be done another way: using varlink everywhere, but it leads to problems of calling varlink after RBE object creation with other methods @@ -239,9 +248,9 @@ function rbe(df; dvar::Symbol, #Secondary parameters calculation # if inv(H) incorrect pinv(H) used if abs(minimum(svdvals(H))) > singlim - A = 2 * inv(H) + A = 2.0 * inv(H) else - A = 2 * pinv(H) + A = 2.0 * pinv(H) end C = pinv(iC) se = Vector{eltype(C)}(undef, p) @@ -249,7 +258,7 @@ function rbe(df; dvar::Symbol, df = Vector{eltype(C)}(undef, p) t = Vector{eltype(C)}(undef, p) pval = Vector{eltype(C)}(undef, p) - gradc = cmatg(Xv, Zv, θ, C; memopt = memopt) + gradc = cmatg(data, θ, C; memopt = memopt) for i = 1:p @@ -300,7 +309,12 @@ function rbe!(df; dvar::Symbol, singlim = 1e-6) - if any(x -> x ∉ names(df), [subject, formulation, period, sequence]) throw(ArgumentError("Names not found in DataFrame!")) end + dfn = names(df) + if eltype(dfn) <: String dfn = Symbol.(dfn) end + if any(x -> x ∉ dfn, [subject, formulation, period, sequence]) + throw(ArgumentError("Names not found in DataFrame! \n Names: $([subject, formulation, period, sequence]) \n DF names: $(names(df))")) + end + if !(eltype(df[!,dvar]) <: AbstractFloat) @warn "Responce variable ∉ AbstractFloat!" df[!,dvar] = float.(df[!,dvar]) diff --git a/test/test.jl b/test/test.jl index ce903d8..beca84d 100644 --- a/test/test.jl +++ b/test/test.jl @@ -6,9 +6,9 @@ using Test, CSV, DataFrames, StatsBase path = dirname(@__FILE__) - +println("Load data...") include("testdata.jl") - +println("Start test...") @testset " Basic mixed model test " begin be = ReplicateBE.rbe!(df0, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, g_tol = 1e-10); e1 = ReplicateBE.fixed(be).est[6] @@ -66,8 +66,6 @@ include("testdata.jl") be = ReplicateBE.rbe!(df0, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence, rholink = :arctgsigmoid, singlim = 1e-4) ci2 = confint(be)[end] - - end @testset " #1 Patterson and Jones 2017 E 4.3 " begin @@ -200,6 +198,10 @@ end end @testset " # Validation with generated datasets " begin + using CSV + path = dirname(@__FILE__) + #println("File path: " , path) + #println("File rds12.csv ($(path*"/csv/rds12.csv")) exist: ", isfile(path*"/csv/rds12.csv")) #1 #TRTR/RTRT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10001) @@ -209,6 +211,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.800182 atol=1E-5 @test ci[end][2] ≈ 1.208955 atol=1E-5 + print("[.") #DF contain 46 #DF contain form 44 #2 @@ -220,7 +223,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.877858 atol=1E-5 @test ci[end][2] ≈ 1.266491 atol=1E-5 - + print(".") #3 #TTRR/RRTT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1], design = ["T" "T" "R" "R"; "R" "R" "T" "T"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10003) @@ -230,6 +233,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.838738 atol=1E-5 @test ci[end][2] ≈ 1.132936 atol=1E-5 + print(".") #4 #TRTR/RTRT/TRRT/RTTR rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1,1,1], design = ["T" "R" "T" "R"; "R" "T" "R" "T" ; "T" "R" "R" "T"; "R" "T" "T" "R"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0, 0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10004) @@ -237,10 +241,9 @@ end @test ReplicateBE.reml2(be) ≈ 178.85707596709256 atol=1E-5 @test ReplicateBE.stderror(be)[end] ≈ 0.1126566088472447 atol=1E-5 ci = confint(be, 0.1, expci = true) - @test ci[end][1] ≈ 0.901150701434849 atol=1E-5 #SPSS 0.900982 @test ci[end][2] ≈ 1.3260902090001896 atol=1E-5 # 1.326338 - + print(".") #5 #TRRT/RTTR/TTRR/RRTT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1,1,1], design = ["T" "R" "R" "T"; "R" "T" "T" "R" ; "T" "T" "R" "R"; "R" "R" "T" "T"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0, 0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10005) @@ -250,6 +253,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.8697647630450958 atol=1E-4 @test ci[end][2] ≈ 1.3438521036812334 atol=1E-4 + print(".") #6 #TRTR/RTRT/TTRR/RRTT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1,1,1], design = ["T" "R" "T" "R"; "R" "T" "R" "T" ; "T" "T" "R" "R"; "R" "R" "T" "T"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0, 0.0, 0.0], periodcoef = [0.0, 0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10006) @@ -259,6 +263,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.7280915876704905 atol=1E-5 @test ci[end][2] ≈ 1.0609077632653205 atol=1E-5 + print(".") #7 #TRT/RTR rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1], design = ["T" "R" "T"; "R" "T" "R"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10007) @@ -270,6 +275,7 @@ end @test ci[end][2] ≈ 1.269283 atol=1E-5 #DF contain 23 #DF contain form 44 + print(".") #8 #TRR/RTT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1], design = ["T" "R" "R"; "R" "T" "T"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10008) @@ -279,6 +285,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.8957103074455708 atol=1E-5 @test ci[end][2] ≈ 1.2534507598650542 atol=1E-5 + print(".") #9 #TR/RT/TT/RR rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1,1,1], design = ["T" "R"; "R" "T"; "T" "T"; "R" "R"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0, 0.0, 0.0], periodcoef = [0.0, 0.0], formcoef = [0.0, 0.0], seed = 10009) @@ -288,6 +295,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.5581809471131624 atol=1E-5 @test ci[end][2] ≈ 1.1201862625078982 atol=1E-5 + print(".") #10 #TRR/RTR/RRT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1,1], design = ["T" "R" "R"; "R" "T" "R"; "R" "R" "T"], inter=[0.5, 0.4, 0.9], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0, 0.0], periodcoef = [0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10010) @@ -297,6 +305,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.8192235933294734 atol=1E-5 @test ci[end][2] ≈ 1.1698775266746582 atol=1E-5 + print(".") #11 #TRR/RTR #SPSS REML 129.655513635398 @@ -305,17 +314,23 @@ end @test ReplicateBE.reml2(be) ≈ 129.06198795781413 atol=1E-5 @test ReplicateBE.stderror(be)[end] ≈ 0.1096846988112254 atol=1E-5 ci = confint(be, 0.1, expci = true) + #1.0.4 - @test ci[end][1] ≈ 0.8719537380504332 atol=1E-5 - @test ci[end][2] ≈ 1.2621358539962937 atol=1E-5 + #@test ci[end][1] ≈ 0.8719537380504332 atol=1E-5 + #@test ci[end][2] ≈ 1.2621358539962937 atol=1E-5 #1.0.5 #@test ci[end][1] ≈ 0.8719023753585434 atol=1E-3 #@test ci[end][2] ≈ 1.2622102048603623 atol=1E-3 - + #1.0.11 + #@test ci[end][1] ≈ 0.8719313261195931 atol=1E-5 + #@test ci[end][2] ≈ 1.2621682956584102 atol=1E-5 + #J1.4.0 + @test ci[end][1] ≈ 0.8719537380504332 atol=1E-4 + @test ci[end][2] ≈ 1.2621358539962937 atol=1E-4 #DF contain 74 #DF contain form 92 - + print(".") #Unbalanced + dropouts #12 #TRTR/RTRT @@ -328,6 +343,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.9223372675334348 atol=1E-5 @test ci[end][2] ≈ 1.183433775497828 atol=1E-5 + print(".") #13 #TRRT/RTTR rds = CSV.File(path*"/csv/rds13.csv") |> DataFrame @@ -338,6 +354,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.9109673193677009 atol=1E-5 @test ci[end][2] ≈ 1.1871890135539473 atol=1E-5 + print(".") #14 #TTRR/RRTT rds = CSV.File(path*"/csv/rds14.csv") |> DataFrame @@ -348,6 +365,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.9372828683843416 atol=1E-5 @test ci[end][2] ≈ 1.2057922163897947 atol=1E-5 + print(".") #15 #TRTR/RTRT/TRRT/RTTR rds = CSV.File(path*"/csv/rds15.csv") |> DataFrame @@ -358,6 +376,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.9431875747076479 atol=1E-5 @test ci[end][2] ≈ 1.2478691227393788 atol=1E-5 + print(".") #16 #TRRT/RTTR/TTRR/RRTT rds = CSV.File(path*"/csv/rds16.csv") |> DataFrame @@ -368,6 +387,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.9088567332123243 atol=1E-5 @test ci[end][2] ≈ 1.1476219506765315 atol=1E-5 + print(".") #17 #TRTR/RTRT/TTRR/RRTT rds = CSV.File(path*"/csv/rds17.csv") |> DataFrame @@ -379,6 +399,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.92723952053598 atol=1E-5 @test ci[end][2] ≈ 1.1951418270978003 atol=1E-5 + print(".") #18 #TRT/RTR rds = CSV.File(path*"/csv/rds18.csv") |> DataFrame @@ -390,6 +411,7 @@ end @test ci[end][1] ≈ 0.878355738329916 atol=1E-5 @test ci[end][2] ≈ 1.1975977027952331 atol=1E-5 #DF contain 34 + print(".") #19 #TRR/RTT rds = CSV.File(path*"/csv/rds19.csv") |> DataFrame @@ -400,6 +422,7 @@ end ci = confint(be, 0.1, expci = true) @test ci[end][1] ≈ 0.7755005422832147 atol=1E-5 @test ci[end][2] ≈ 1.0916738532751222 atol=1E-5 + print(".") #20 #TR/RT/TT/RR #SPSS REML 151.783195849874 @@ -414,6 +437,7 @@ end @test ci[end][1] ≈ 0.7865283528654503 atol=1E-5 @test ci[end][2] ≈ 1.6924035882963178 atol=1E-5 #DF contain 20 + print(".") #21 #TRR/RTR/RRT rds = CSV.File(path*"/csv/rds21.csv") |> DataFrame @@ -429,9 +453,10 @@ end #v1.0.5 #@test ci[end][1] ≈ 0.8293880338983111 atol=1E-3 #@test ci[end][2] ≈ 1.1480342197874598 atol=1E-3 - + #1.4.0 #DF contain 35 + print(".") #22 #TRR/RTR rds = CSV.File(path*"/csv/rds22.csv") |> DataFrame @@ -447,9 +472,10 @@ end #v1.0.5 #@test ci[end][1] ≈ 0.9064394831113237 atol=1E-3 #@test ci[end][2] ≈ 1.3055658048846954 atol=1E-3 + #1.4.0 #DF contain 32 - + print(".") # Unbalanced by sequences #23 #TRTR/RTRT @@ -460,6 +486,7 @@ end @test ReplicateBE.reml2(be) ≈ 252.0449059441563 atol=1E-5 @test ci[end][1] ≈ 0.638039 atol=1E-5 @test ci[end][2] ≈ 1.135006 atol=1E-5 + print(".") #24 #TRT/RTR rds = ReplicateBE.randrbeds(;n=36, sequence=[1,2], design = ["T" "R" "T"; "R" "T" "R"], inter=[0.4, 0.3, 0.2], intra=[0.05, 0.02], intercept = 1.0, seqcoef = [0.1, 0.0], periodcoef = [0.0, 0.0, 0.0], formcoef = [0.1, 0.0], seed = 10024) @@ -468,7 +495,7 @@ end @test ReplicateBE.reml2(be) ≈ 140.10714908682417 atol=1E-5 @test ci[end][1] ≈ 0.854985 atol=1E-5 @test ci[end][2] ≈ 1.293459 atol=1E-5 - + print(".") #25 #TRR/RTT rds = ReplicateBE.randrbeds(;n=24, sequence=[1,1], design = ["T" "R" "R"; "R" "T" "T"], inter=[0.5, 0.4, 0.5], intra=[0.1, 0.2], intercept = 1.0, seqcoef = [0.0, 0.0], periodcoef = [0.0, 0.0, 0.0], formcoef = [0.0, 0.0], seed = 10008) @@ -477,7 +504,7 @@ end @test ReplicateBE.reml2(be) ≈ 140.37671499958694 atol=1E-5 @test ci[end][1] ≈ 0.8208303872086634 atol=1E-5 @test ci[end][2] ≈ 1.3228068920153602 atol=1E-5 - + print(".") # Unbalanced by sequences #26 #TRTR/RTRT @@ -487,6 +514,7 @@ end @test ReplicateBE.reml2(be) ≈ 899.044671814185 atol=1E-5 @test ci[end][1] ≈ 0.6189149073430663 atol=1E-5 @test ci[end][2] ≈ 0.7726263704382558 atol=1E-5 + print(".") #27 #TRT/RTR rds = ReplicateBE.randrbeds(;n=128, sequence=[1,2], design = ["T" "R" "T"; "R" "T" "R"], inter=[0.4, 0.3, 0.7], intra=[0.05, 0.2], intercept = 1.0, seqcoef = [1.1, 0.0], periodcoef = [0.0, 1.0, 0.0], formcoef = [0.2, 0.0], seed = 10027) @@ -495,8 +523,9 @@ end @test ReplicateBE.reml2(be) ≈ 614.3407387528628 atol=1E-5 @test ci[end][1] ≈ 1.18696161597237 atol=1E-5 @test ci[end][2] ≈ 1.43349025480276 atol=1E-5 - + print(".") # Unbalanced by sequences + #28 #TRTR/RTRT rds = ReplicateBE.randrbeds(;n=512, sequence=[1,2], design = ["T" "R" "T" "R"; "R" "T" "R" "T"], inter=[0.5, 0.4, 0.7], intra=[0.1, 0.15], intercept = 1.0, seqcoef = [1.0, 0.0], periodcoef = [0.0, 1.0, 0.0, 0.0], formcoef = [0.0, 0.3], seed = 10028) @@ -505,6 +534,7 @@ end @test ReplicateBE.reml2(be) ≈ 3495.580039200916 atol=1E-5 @test ci[end][1] ≈ 0.6975351860484861 atol=1E-5 @test ci[end][2] ≈ 0.7698493312540238 atol=1E-5 + print(".") #29 #TRT/RTR rds = ReplicateBE.randrbeds(;n=512, sequence=[1,2], design = ["T" "R" "T"; "R" "T" "R"], inter=[0.4, 0.3, 0.7], intra=[0.05, 0.2], intercept = 1.0, seqcoef = [1.1, 0.0], periodcoef = [0.0, 0.0, 1.0], formcoef = [0.4, 0.0], seed = 10029) @@ -514,7 +544,7 @@ end @test ci[end][1] ≈ 1.361487533935354 atol=1E-5 @test ci[end][2] ≈ 1.498934137926506 atol=1E-5 - + println(".]") end @testset " # Simulation " begin @@ -537,6 +567,9 @@ end pow = ReplicateBE.simulation(task; io = io, num = 1007, seed = 1234, verbose = true) @test pow.result ≈ 0.06 atol=1E-2 + pow = ReplicateBE.simulation(task; io = io, num = 100, seed = 1234, verbose = true, rsabe = true) + #@test pow.result ≈ 0.06 atol=1E-2 + #Custom simulation task = ReplicateBE.randrbetask(;n=24,