Changes in .Show()

PharmCat · Aug 27, 2019 · 68e3abb · 68e3abb
1 parent dacfc70
commit 68e3abb
Show file tree

Hide file tree

Showing 3 changed files with 62 additions and 39 deletions.
diff --git 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
@@ -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
@@ -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