diff --git a/src/ReplicateBE.jl b/src/ReplicateBE.jl index cb7ff87..cf09e8d 100644 --- a/src/ReplicateBE.jl +++ b/src/ReplicateBE.jl @@ -17,6 +17,7 @@ LOG2PI = log(2π) struct RBE model::ModelFrame #Model frame + rmodel::ModelFrame #Random effect model factors::Array{Symbol, 1} #Factor list β::Array{Float64, 1} #β coefficients (fixed effect) θ0::Array{Float64, 1} #Initial variance paramethers @@ -80,8 +81,9 @@ function rbe(df; dvar::Symbol, Xf = @eval(@formula($dvar ~ $sequence + $period + $formulation)) Zf = @eval(@formula($dvar ~ 0 + $formulation)) MF = ModelFrame(Xf, df) + RMF = ModelFrame(Zf, df, contrasts = Dict(formulation => StatsModels.FullDummyCoding())) X = ModelMatrix(MF).m - Z = ModelMatrix(ModelFrame(Zf, df, contrasts = Dict(formulation => StatsModels.FullDummyCoding()))).m + Z = ModelMatrix(RMF).m p = rank(X) y = df[:, dvar] #Dependent variable @@ -144,7 +146,7 @@ function rbe(df; dvar::Symbol, A = 2*pinv(H) se, F, df, C = ctrst(p, Xv, Zv, iVv, θ, β, A) df2 = N / pn - sn - return RBE(MF, [sequence, period, formulation], β, θvec0, θ, remlv, se, F, df, df2, Rv, Vv, G, C, A, H, X, Z, Xv, Zv, yv, dH, pO, O) + return RBE(MF, RMF, [sequence, period, formulation], β, θvec0, θ, remlv, se, F, df, df2, Rv, Vv, G, C, A, H, X, Z, Xv, Zv, yv, dH, pO, O) end #END OF rbe() #------------------------------------------------------------------------------- #------------------------------------------------------------------------------- @@ -291,7 +293,7 @@ function ctrst(p, Xv, Zv, iVv, θ, β, A) se = Array{Float64, 1}(undef, p) F = Array{Float64, 1}(undef, p) df = Array{Float64, 1}(undef, p) - for i = 2:p + for i = 1:p L = zeros(p) L[i] = 1 L = L' diff --git a/src/show.jl b/src/show.jl index 6a636d1..61d127b 100644 --- a/src/show.jl +++ b/src/show.jl @@ -1,28 +1,22 @@ function Base.show(io::IO, obj::RBE) println(io, "Bioequivalence Linear Mixed Effect Model") println(io, "") - println(io, "REML: ", obj.reml) + println(io, "REML: ", round(obj.reml, sigdigits=6)) println(io, "") - lines = 0 - for f in obj.factors - lines += length(obj.model.contrasts[f].termnames) - end - pm = Array{Any,2}(undef, lines+2, 5) + coef = coefnames(obj.model); + rcoef = coefnames(obj.rmodel); + pm = Array{Any,2}(undef, length(coef)+1, 5); pm[1,1] = "Level"; pm[1,2] = "Value"; pm[1,3] = "SE"; pm[1,4] = "DF"; pm[1,5] = "F"; - pm[2,1] = "Intercept"; pm[2,2] = obj.β[1]; pm[2,3] = "-"; pm[2,4] = "-"; pm[2,5] = "-"; - it = 2 - pmr = 3 - for f in obj.factors - for l in obj.model.contrasts[f].termnames - pm[pmr,1] = string(f)*": "*string(l); - pm[pmr,2] = obj.β[it]; - pm[pmr,3] = obj.se[it]; - pm[pmr,4] = obj.df[it]; - pm[pmr,5] = obj.f[it]; - it += 1 - pmr += 1 - end + + for i = 1:length(coef) + pm[i+1,1] = coef[i] + pm[i+1,2] = obj.β[i]; + pm[i+1,3] = obj.se[i]; + pm[i+1,4] = obj.df[i]; + pm[i+1,5] = obj.f[i]; end + + # for r = 1:size(pm)[1] for c = 1:size(pm)[2] if isa(pm[r,c], Float64) pm[r,c] = round(pm[r,c], sigdigits=6) end @@ -32,39 +26,66 @@ function Base.show(io::IO, obj::RBE) pm = string.(pm) len = Array{Int,1}(undef, 0) vch = Array{Vector{Char},1}(undef, 0) + line = "" + #Line for c = 1:size(pm)[2] ml = maximum(length.(pm[:,c])) push!(len, ml) - ch = Vector{Char}(undef, ml) + ch = Vector{Char}(undef, ml+2) ch .= '─' - push!(vch, ch) - end - for v in vch - print(io, String(v)) - print(io, "┬") + line *= String(ch) end - println(io, "") + # + println(io, line) + #Head for c = 1:size(pm)[2] print(io, pm[1,c]) ch = Vector{Char}(undef, len[c]-length(pm[1,c])) ch .= ' ' print(io, String(ch)) - print(io, "│") + print(io, " ") end println(io, "") + # + #Line + println(io, line) + # + #Body for r = 2:size(pm)[1] - for v in vch - print(io, String(v)) - print(io, "┼") - end - println(io, "") for c = 1:size(pm)[2] print(io, pm[r,c]) ch = Vector{Char}(undef, len[c]-length(pm[r,c])) ch .= ' ' print(io, String(ch)) - print(io, "│") + print(io, " ") end println(io, "") end + # + #line + println(io, line) + # + println(io, "Intra-individual variation:") + println(io, rcoef[1], " ", round(obj.θ[1], sigdigits=6), " CVᵂ: ", round(geocv(obj.θ[1]), sigdigits=6)) + println(io, rcoef[2], " ", round(obj.θ[2], sigdigits=6), " CVᵂ: ", round(geocv(obj.θ[2]), sigdigits=6)) + println(io, "") + println(io, "Inter-individual variation:") + println(io, rcoef[1], " ", round(obj.θ[3], sigdigits=6)) + println(io, rcoef[2], " ", round(obj.θ[4], sigdigits=6)) + println(io, "Cov:", " ", round(sqrt(obj.θ[4]*obj.θ[3])*obj.θ[5], sigdigits=6)) + println(io, line) + println(io, "Confidence intervals(90%):") + println(io, "") + ci = confint(obj, 0.1, expci = true, inv = false) + println(io, rcoef[1], " / ", rcoef[2]) + println(io, round(ci[end][1]*100, digits=4), " - ", round(ci[end][2]*100, digits=4), " (%)") + println(io, "") + ci = confint(obj, 0.1, expci = true, inv = true) + println(io, rcoef[2], " / ", rcoef[1]) + println(io, round(ci[end][1]*100, digits=4), " - ", round(ci[end][2]*100, digits=4), " (%)") + +end + +function geocv(var) + return sqrt(exp(var)-1.0) end diff --git a/test/test.jl b/test/test.jl index 62058db..b662bb5 100644 --- a/test/test.jl +++ b/test/test.jl @@ -7,7 +7,7 @@ using Test, DataFrames, CSV include("testdata.jl") -@testset " Basic mixed model test " begin +@testset " Basic mixed model test " begin df = CSV.read(IOBuffer(minibe)) |> DataFrame be = ReplicateBE.rbe(df, dvar = :var, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence); @test be.β[6] ≈ -0.0791666 atol=1E-5 @@ -28,7 +28,7 @@ end end #Patterson SD, Jones B. Viewpoint: observations on scaled average bioequivalence. Pharm Stat. 2012; 11(1): 1–7. doi:10.1002/pst.498 -@testset " #5 Pub Bioequivalence Dataset " begin +@testset " #5 Pub Bioequivalence Dataset " begin #REML 321.44995530 - SAS STOP! df = CSV.read(IOBuffer(be5)) |> DataFrame be = ReplicateBE.rbe(df, dvar = :var1, subject = :subject, formulation = :formulation, period = :period, sequence = :sequence); @@ -40,7 +40,7 @@ end end #Shumaker RC, Metzler CM. The Phenytoin Trial is a Case Study of ‘Individual’ Bioequivalence. Drug Inf J. 1998; 32(4): 1063–72 -@testset " #6 Pub Bioequivalence TTRR/RRTT Dataset " begin +@testset " #6 Pub Bioequivalence TTRR/RRTT Dataset " begin #REML 329.25749378 #SE 0.04153 #DF 62