diff --git a/Project.toml b/Project.toml index d73132d..e14e597 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.8" +version = "1.0.9" [deps] @@ -37,7 +37,7 @@ StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0 DataFrames = "0.19, 0.20" StatsModels = "0.6" Distributions = "0.20, 0.21, 0.22, 0.23" -PDMats = "0.9" +PDMats = "0.9, 0.10" ForwardDiff = "0.10" Optim = "0.19, 0.20, 0.21" LineSearches = "7.0" diff --git a/chagelog.md b/chagelog.md index 13c8b96..efb38a6 100644 --- a/chagelog.md +++ b/chagelog.md @@ -1,3 +1,7 @@ +- v1.0.9 + * PDMats 0.10 + * Opimizations + - v1.0.8 * Changing in RBE struct! @@ -9,7 +13,7 @@ - v1.0.7 * settings: initial - * settings: initial + * settings: initial * randrbeds "light" generation * optimizations * Distributions bump diff --git a/src/algebra.jl b/src/algebra.jl index e3543b2..e77b83d 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -43,7 +43,7 @@ function mulαtβαinc!(θ, A::AbstractMatrix, B::AbstractMatrix) end #------------------------------------------------------------------------------- """ -A' * B * A -> +θ +A' * B * A -> +θ (cache) A' * B * C -> +β """ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::Vector, c) @@ -67,26 +67,7 @@ function mulθβinc!(θ, β, A::AbstractMatrix, B::AbstractMatrix, C::Vector, c) end #------------------------------------------------------------------------------- """ -A' * B * A -""" -function mulall(A::Vector, B::AbstractMatrix) - q = size(B, 1) - θ = 0 - #c .= 0 - for n = 1:q - for m = 1:q - #@inbounds c[n] += B[m, n] * A[m] - @inbounds θ += B[m, n] * A[m] * A[n] - end - #@inbounds θ += A[n] * c[n] - end - #for n = 1:q - # @inbounds θ += A[n] * c[n] - #end - θ -end -""" -(y - X * β)' * V * (y - X * β) +(y - X * β)' * V * (y - X * β) (cache) """ function mulθ₃(y::Vector, X::AbstractMatrix, β::Vector, V::AbstractMatrix, c) q = size(V, 1) @@ -106,30 +87,6 @@ function mulθ₃(y::Vector, X::AbstractMatrix, β::Vector, V::AbstractMatrix, c return θ end - -""" -A * B * A' -""" -function mulαβαt(A::AbstractMatrix, B::AbstractMatrix) - 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 - @inbounds c[n] += A[i, m] * B[n, m] - end - end - for n = 1:p - for m = 1:q - @inbounds mx[i, n] += A[n, m] * c[m] - end - end - end - mx -end """ A * B * A' + C """ @@ -154,6 +111,9 @@ function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) mx .+= C end +""" +A * B * A' + C (cache) +""" function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c::AbstractVector) q = size(B, 1) p = size(A, 1) @@ -176,49 +136,6 @@ function mulαβαtc(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, c: #SMatrix{p,p,eltype(mx)}(mx) end -""" -A * B * A' + C -> O -""" -function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) - q = size(B, 1) - p = size(A, 1) - c = zeros(eltype(B), q) - O .= 0 - for i = 1:p - c .= 0 - for n = 1:q - for m = 1:q - @inbounds c[n] += A[i, m] * B[n, m] - end - end - for n = 1:p - for m = 1:q - @inbounds O[i, n] += A[n, m] * c[m] - end - end - end - O .+= C -end -function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, mem::MemCache) - q = size(B, 1) - p = size(A, 1) - c = mem.svec[p] - O .= 0 - for i = 1:p - c .= 0 - for n = 1:q - for m = 1:q - @inbounds c[n] += A[i, m] * B[n, m] - end - end - for n = 1:p - for m = 1:q - @inbounds O[i, n] += A[n, m] * c[m] - end - end - end - O .+= C -end function invchol(M) q = size(M, 1) diff --git a/src/deprecated.jl b/src/deprecated.jl index ad9c2e9..2243499 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -1,4 +1,220 @@ +""" +A' * B * A +""" +function mulall(A::Vector, B::AbstractMatrix) + q = size(B, 1) + θ = 0 + #c .= 0 + for n = 1:q + for m = 1:q + #@inbounds c[n] += B[m, n] * A[m] + @inbounds θ += B[m, n] * A[m] * A[n] + end + #@inbounds θ += A[n] * c[n] + end + #for n = 1:q + # @inbounds θ += A[n] * c[n] + #end + θ +end + +""" +A * B * A' +""" +function mulαβαt(A::AbstractMatrix, B::AbstractMatrix) + 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 + @inbounds c[n] += A[i, m] * B[n, m] + end + end + for n = 1:p + for m = 1:q + @inbounds mx[i, n] += A[n, m] * c[m] + end + end + end + mx +end + +""" +A * B * A' + C -> O +""" +function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) + q = size(B, 1) + p = size(A, 1) + c = zeros(eltype(B), q) + O .= 0 + for i = 1:p + c .= 0 + for n = 1:q + for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + for n = 1:p + for m = 1:q + @inbounds O[i, n] += A[n, m] * c[m] + end + end + end + O .+= C +end +function mulαβαtc!(O::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, mem::MemCache) + q = size(B, 1) + p = size(A, 1) + c = mem.svec[p] + O .= 0 + for i = 1:p + c .= 0 + for n = 1:q + for m = 1:q + @inbounds c[n] += A[i, m] * B[n, m] + end + end + for n = 1:p + for m = 1:q + @inbounds O[i, n] += A[n, m] * c[m] + end + end + end + O .+= C +end +#------------------------------------------------------------------------------- +#Return term levels count by symbol +function termmodelleveln(MF::ModelFrame, symbol::Symbol)::Int + id = findterm(MF, symbol) + return length(MF.f.rhs.terms[id].contrasts.levels) +end +#= + +""" + lsm(rbe::RBE, L::Matrix) +Deprecated. +""" +#Deprecated +function lsm(rbe::RBE, L::Matrix) + lcl = L*rbe.C*L' + return L*coef(rbe), sqrt.(lcl) +end +# +""" + emm(obj::RBE, fm::Matrix, lm::Matrix) + +Matrix mask. +""" +function emm(obj::RBE, fm::Matrix, lm::Matrix) + La = lmean(obj::RBE) + L = La .* fm + L = L .+ lm + return lsm(obj, Matrix(L)) +end +#General mean contrast L matrix 1xp +""" + lmean(obj::RBE) + +Return L-matrix for general mean. +""" +function lmean(obj::RBE) + L = zeros(1, length(fixed(obj).est)) + L[1] = 1.0 + it = 2 + for f in obj.data.factors + term = findterm(obj.model, f) + len = length(obj.model.f.rhs.terms[term].contrasts.termnames) + dev = 1/length(obj.model.f.rhs.terms[term].contrasts.levels) + for i = 1:len + L[it] = dev + it += 1 + end + end + return L +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 +=# + +#------------------------------------------------------------------------------- +#= """ Feel M with zero feeled matrices """ @@ -104,7 +320,7 @@ function lclgf(L, Lt, Xv::Vector, Zv::Vector, θ::Vector; memopt::Bool = true) end return (L * inv(C) * Lt)[1] end - +=# #= function mrmat(σ::Vector{S}, Z::Matrix{T}, cache)::Matrix where S <: Real where T <: Real h = hash(tuple(σ, Z)) diff --git a/src/generalfunc.jl b/src/generalfunc.jl index c7395ee..42c299e 100644 --- a/src/generalfunc.jl +++ b/src/generalfunc.jl @@ -6,30 +6,36 @@ Make X, Z matrices and vector y for each subject; """ function sortsubjects(df::DataFrame, sbj::Symbol, X::Matrix, Z::Matrix, y::Vector) - u = unique(df[:, sbj]) - Xa = Vector{Matrix{eltype(X)}}(undef, length(u)) - Za = Vector{Matrix{eltype(Z)}}(undef, length(u)) + 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)) for i = 1:length(u) - v = findall(x->x==u[i], df[:, sbj]) - Xs = Vector{eltype(X)}(undef, 0) - Zs = Vector{eltype(Z)}(undef, 0) - ys = Vector{eltype(y)}(undef, 0) - for r in v - append!(Xs, X[r, :]) - append!(Zs, Z[r, :]) - push!(ys, y[r]) - end - Xa[i] = Matrix(reshape(Xs, size(X)[2], :)') - Za[i] = Matrix(reshape(Zs, size(Z)[2], :)') - ya[i] = ys + v = findall(x->x==u[i], df[!, sbj]) + #Xs = Vector{eltype(X)}(undef, 0) + #Zs = Vector{eltype(Z)}(undef, 0) + #ys = Vector{eltype(y)}(undef, 0) + + #for r in v + #append!(Xs, X[r, :]) + #append!(Zs, Z[r, :]) + #push!(ys, y[r]) + #end + Xa[i] = view(X, v, :) + Za[i] = view(Z, v, :) + ya[i] = view(y, v) + #Xa[i] = Matrix(reshape(Xs, size(X)[2], :)') + #Za[i] = Matrix(reshape(Zs, size(Z)[2], :)') + #ya[i] = ys end + #= for i = 1:length(u) for c = 1:length(u) if Za[i] == Za[c] && Za[i] !== Za[c] Za[i] = Za[c] end if Xa[i] == Xa[c] && Xa[i] !== Xa[c] Xa[i] = Xa[c] end end end + =# return Xa, Za, ya end """ @@ -52,13 +58,13 @@ end """ R matrix (ForwardDiff+) """ -@inline function rmat(σ::Vector, Z::Matrix)::Matrix +@inline function rmat(σ::Vector, Z::AbstractMatrix)::Matrix return Diagonal(Z*σ) end """ R matrix (memory pre-allocation) """ -@inline function rmat!(R::AbstractMatrix{T}, σ::Vector{T}, Z::Matrix{T}) where T <: AbstractFloat +@inline function rmat!(R::AbstractMatrix{T}, σ::Vector{T}, Z::AbstractMatrix{T}) where T <: AbstractFloat copyto!(R, Diagonal(Z*σ)) return end @@ -67,17 +73,17 @@ end """ Return variance-covariance matrix V """ -@inline function vmat(G::AbstractMatrix, R::AbstractMatrix, Z::Matrix)::AbstractMatrix +@inline function vmat(G::AbstractMatrix, R::AbstractMatrix, Z::AbstractMatrix)::AbstractMatrix return mulαβαtc(Z, G, R) end -@inline function vmat!(V::Matrix{T}, G::AbstractMatrix{T}, R::AbstractMatrix{T}, Z::Matrix{T}, memc) where T <: AbstractFloat +@inline function vmat!(V::Matrix{T}, G::AbstractMatrix{T}, R::AbstractMatrix{T}, Z::AbstractMatrix{T}, memc) where T <: AbstractFloat #copyto!(V, Z*G*Z') mul!(memc[size(Z)[1]], Z, G) mul!(V, memc[size(Z)[1]], Z') V .+= R return end -function mvmat(G::AbstractMatrix, σ::Vector, Z::Matrix, cache)::Matrix +function mvmat(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, cache)::Matrix #h = hash(tuple(σ, Z)) if Z in keys(cache) return cache[Z] @@ -88,7 +94,7 @@ function mvmat(G::AbstractMatrix, σ::Vector, Z::Matrix, cache)::Matrix return V end end -function mvmat(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache)::Matrix +function mvmat(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, mem, cache)::Matrix #h = hash(tuple(σ, Z)) if Z in keys(cache) return cache[Z] @@ -100,7 +106,7 @@ function mvmat(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache)::Matrix end end -function mvmatall(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache) +function mvmatall(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, mem, cache) #h = hash(Z) #if h in keys(cache) if Z in keys(cache) @@ -127,7 +133,7 @@ function mvmatall(G::AbstractMatrix, σ::Vector, Z::Matrix, mem, cache) end #println("θ₁: ", θ1, " θ₂: ", θ2, " θ₃: ", θ3) -function minv(G::AbstractMatrix, σ::Vector, Z::Matrix, cache::Dict)::Matrix +function minv(G::AbstractMatrix, σ::Vector, Z::AbstractMatrix, cache::Dict)::Matrix #h = hash(M) if Z in keys(cache) #return cache[h] diff --git a/src/randrbeds.jl b/src/randrbeds.jl index 945c2d5..8cccafe 100644 --- a/src/randrbeds.jl +++ b/src/randrbeds.jl @@ -417,75 +417,3 @@ 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 6766b33..5a2d247 100644 --- a/src/rbe.jl +++ b/src/rbe.jl @@ -137,7 +137,7 @@ function rbe(df; dvar::Symbol, Z = ModelMatrix(RMF).m p = rank(X) zxr = rank(ModelMatrix(ModelFrame(@eval(@formula($dvar ~ $sequence + $period + $subject*$formulation)), df)).m) - y = df[:, dvar] #Dependent variable + y = df[!, dvar] #Dependent variable #Make pre located arrays with matrices for each subject Xv, Zv, yv = sortsubjects(df, subject, X, Z, y) n = length(Xv) @@ -363,8 +363,8 @@ end Returm -2logREML for rbe model with θ variance vector. """ -function reml2(rbe::RBE, θ::Vector{T}) where T <: AbstractFloat - return reml2(rbe.data, θ, coef(rbe)) +function reml2(rbe::RBE, θ::Vector{T}; memopt::Bool = true) where T <: AbstractFloat + return reml2(rbe.data, θ, coef(rbe); memopt = memopt) end """ reml2(rbe::RBE) diff --git a/src/utils.jl b/src/utils.jl index f888c32..61957bc 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -164,55 +164,7 @@ function termmodellen(MF::ModelFrame, symbol::Symbol)::Int id = findterm(MF, symbol) return length(MF.f.rhs.terms[id].contrasts.termnames) end -#Return term levels count by symbol -function termmodelleveln(MF::ModelFrame, symbol::Symbol)::Int - id = findterm(MF, symbol) - return length(MF.f.rhs.terms[id].contrasts.levels) -end # -""" - lsm(rbe::RBE, L::Matrix) - -Deprecated. -""" -#Deprecated -function lsm(rbe::RBE, L::Matrix) - lcl = L*rbe.C*L' - return L*coef(rbe), sqrt.(lcl) -end -# -""" - emm(obj::RBE, fm::Matrix, lm::Matrix) - -Matrix mask. -""" -function emm(obj::RBE, fm::Matrix, lm::Matrix) - La = lmean(obj::RBE) - L = La .* fm - L = L .+ lm - return lsm(obj, Matrix(L)) -end -#General mean contrast L matrix 1xp -""" - lmean(obj::RBE) - -Return L-matrix for general mean. -""" -function lmean(obj::RBE) - L = zeros(1, length(fixed(obj).est)) - L[1] = 1.0 - it = 2 - for f in obj.data.factors - term = findterm(obj.model, f) - len = length(obj.model.f.rhs.terms[term].contrasts.termnames) - dev = 1/length(obj.model.f.rhs.terms[term].contrasts.levels) - for i = 1:len - L[it] = dev - it += 1 - end - end - return L -end #------------------------------------------------------------------------------- function checkdata(X::Matrix, Z::Matrix, Xv::Vector, Zv::Vector, y::Vector) if size(Z)[2] != 2 error("Size random effect matrix != 2. Not implemented yet!") end diff --git a/test/test.jl b/test/test.jl index 1571091..c8c437b 100644 --- a/test/test.jl +++ b/test/test.jl @@ -15,6 +15,7 @@ include("testdata.jl") @test ReplicateBE.fixed(be).est[6] ≈ -0.0791666 atol=1E-5 @test ReplicateBE.fixed(be).se[6] ≈ 0.09037378448083119 atol=1E-5 @test ReplicateBE.reml2(be) ≈ 10.065238638105903 atol=1E-5 + @test ReplicateBE.reml2(be, ReplicateBE.theta(be), memopt = false) ≈ ReplicateBE.reml2(be) atol=1E-8 ci = confint(be, 0.1, expci = false, inv = false) @test ci[end][1] ≈ -0.25791330363201714 atol=1E-5 @test ci[end][2] ≈ 0.09957997029868393 atol=1E-5 @@ -531,5 +532,6 @@ end push!(out, ReplicateBE.theta(be)[5]) end result = ReplicateBE.simulation!(task, out, simfunc!; num = 10, seed = 10, verbose = false) + Base.show(io, result) @test mean(result) ≈ 0.8904498245891423 end