PharmCat
★    

Russia,
2020-10-03 21:59
(1531 d 09:57 ago)

Posting: # 21968
Views: 5,701
 

 Metida.jl: experimental Julia LinearMixedModel software [Software]

Post reply
Hello, dear colleagues!

Some time ago I released the mixed model solution to bioequivalence for Julia (ReplicateBE.jl) and it works fine but this package is not flexible and can calculate only one kind of task. Now I try to make a general instrument for linear mixed models with some approaches that were used in ReplicateBE.jl.

Metida.jl is an experimental Julia package for fitting linear mixed model with flexible covariance structure. Main goal - to implement the ability to develop custom covariance structures and easy way for contributors to do it. Top goal - to realize the most used PROC MIXED abilities from SAS.

Now Metida.jl in the early development stage and checked on one small dataset (SAS code). Covariance type implemented: Scaled identity, Variance components, Heterogeneous Compound Symmetry.

If Julia is already installed Metida.jl can be installed very simple way:

import Pkg; Pkg.add("Metida")

Version 0.1.3 on this moment available.

To check you can use this example:

using Metida, StatsBase, StatsModels, CSV, DataFrames
df = CSV.File(dirname(pathof(Metida))*"\\..\\test\\csv\\df0.csv") |> DataFrame
categorical!(df, :subject);
categorical!(df, :period);
categorical!(df, :sequence);
categorical!(df, :formulation);

lmm = LMM(@formula(var ~ sequence + period + formulation), df;
random = VarEffect(@covstr(formulation), CSH),
repeated = VarEffect(@covstr(formulation), VC),
subject = :subject)

fit!(lmm)

There are many work to do: check other validation datasets, add other covariance types, performance things and other... So I will be very appreciate for any issues on github or here, advices ets.
It will be very cool if you guided me for any validation aspects.

Note: One of the most notable problems is the long initial start-up time. The second run of the same command is much faster. Worst of all, at the moment, changing the covariance structure also requires recompilation and takes a certain amount of time.

Some comparations with PROC MIXED SAS:

using Metida, StatsBase, StatsModels, CSV, DataFrames
df = CSV.File(dirname(pathof(Metida))*"\\..\\test\\csv\\df0.csv") |> DataFrame

# EXAMPLE 1
################################################################################
################################################################################
################################################################################
#=
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL  var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM  formulation/TYPE=CSH SUB=subject G V;
REPEATED/GRP=formulation SUB=subject R;
RUN;

REML: 10.06523862
=#

lmm = LMM(@formula(var ~ sequence + period + formulation), df;
random   = VarEffect(@covstr(formulation), CSH),
repeated = VarEffect(@covstr(formulation), VC),
subject  = :subject)

fit!(lmm)
#=
Linear Mixed Model: var ~ sequence + period + formulation
Random 1:
   Model: formulation
   Type: HeterogeneousCompoundSymmetry (3)
   Coefnames: ["formulation: 1", "formulation: 2"] Repeated:
   Model: formulation
   Type: VarianceComponents (2)
   Coefnames: ["formulation: 1", "formulation: 2"]
Status: converged

   -2 logREML: 10.0652

   Fixed effects:

Name             Estimate     SE
(Intercept)      1.57749      0.334543
sequence: 2      -0.170833    0.384381
period: 2        0.195984     0.117228
period: 3        0.145014     0.109171
period: 4        0.157363     0.117228
formulation: 2   -0.0791667   0.0903709

Random effects:

   θ vector: [0.455584, 0.367656, 1.0, 0.143682, 0.205657]
Random 1   formulation: 1   var   0.207557
Random 1   formulation: 2   var   0.135171
Random 1   Rho              rho   1.0
Repeated   formulation: 1   var   0.0206445
Repeated   formulation: 2   var   0.0422948
=#

# EXAMPLE 2
################################################################################
################################################################################
################################################################################
#=
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL  var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM  formulation/TYPE=VC SUB=subject G V;
REPEATED/GRP=formulation SUB=subject R;
RUN;

REML: 16.06148160
=#

lmm = LMM(
    @formula(var ~ sequence + period + formulation), df;
    random   = VarEffect(@covstr(formulation), SI),
    repeated = VarEffect(@covstr(formulation), VC),
    subject  = :subject,
)
fit!(lmm)
#=
Linear Mixed Model: var ~ sequence + period + formulation
Random 1:
   Model: formulation
   Type: ScaledIdentity (1)
   Coefnames: ["formulation: 1", "formulation: 2"] Repeated:
   Model: formulation
   Type: VarianceComponents (2)
   Coefnames: ["formulation: 1", "formulation: 2"]
Status: converged

   -2 logREML: 16.0615

   Fixed effects:

Name             Estimate     SE
(Intercept)      1.57212      0.305807
sequence: 2      -0.170833    0.279555
period: 2        0.204087     0.289957
period: 3        0.155769     0.11308
period: 4        0.160015     0.289957
formulation: 2   -0.0791667   0.279555

Random effects:

   θ vector: [0.412436, 0.145184, 0.220819]
Random 1   Var              var   0.170103
Repeated   formulation: 1   var   0.0210784
Repeated   formulation: 2   var   0.048761
=#

#EXAMPLE 3
################################################################################
################################################################################
################################################################################
#=
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL  var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM  subject/TYPE=VC G V;
RUN;

REML: 10.86212458
=#

lmm = LMM(@formula(var ~ sequence + period + formulation), df;
    random = VarEffect(@covstr(subject), SI)
    )
fit!(lmm)
#=
Linear Mixed Model: var ~ sequence + period + formulation
Random 1:
   Model: subject
   Type: ScaledIdentity (1)
   Coefnames: ["subject: 1", "subject: 2", "subject: 3", "subject: 4", "subject: 5"] Repeated:
   Model: nothing
   Type: ScaledIdentity (1)
   Coefnames: -

Status: converged

   -2 logREML: 10.8621

   Fixed effects:

Name             Estimate     SE
(Intercept)      1.61         0.309774
sequence: 2      -0.170833    0.383959
period: 2        0.144167     0.116706
period: 3        0.08         0.115509
period: 4        0.144167     0.116706
formulation: 2   -0.0791667   0.0833617

Random effects:

   θ vector: [0.410574, 0.182636]
Random 1   Var   var   0.168571
Repeated   Var   var   0.0333559
=#

#EXAMPLE 4
################################################################################
################################################################################
################################################################################
#=
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL  var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM  period/TYPE=VC G V;
RANDOM  formulation/TYPE=VC G V;
RUN;

REML: 25.12948063
=#
lmm = LMM(
    @formula(var ~ sequence + period + formulation), df;
    random   = [VarEffect(@covstr(period), VC), VarEffect(@covstr(formulation), VC)] )
fit!(lmm)

#=
Linear Mixed Model: var ~ sequence + period + formulation
Random 1:
   Model: period
   Type: VarianceComponents (4)
   Coefnames: ["period: 1", "period: 2", "period: 3", "period: 4"] Random 2:
   Model: formulation
   Type: VarianceComponents (2)
   Coefnames: ["formulation: 1", "formulation: 2"] Repeated:
   Model: nothing
   Type: ScaledIdentity (1)
   Coefnames: -

Status: converged

   -2 logREML: 25.1295

   Fixed effects:

Name             Estimate     SE
(Intercept)      1.61         0.249491
sequence: 2      -0.170833    0.192487
period: 2        0.144167     0.269481
period: 3        0.08         0.266717
period: 4        0.144167     0.269481
formulation: 2   -0.0791667   0.192487

Random effects:

   θ vector: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.421717]
Random 1   period: 1        var   0.0
Random 1   period: 2        var   0.0
Random 1   period: 3        var   0.0
Random 1   period: 4        var   0.0
Random 2   formulation: 1   var   0.0
Random 2   formulation: 2   var   0.0
Repeated   Var              var   0.177845
=#

#EXAMPLE 5
################################################################################
################################################################################
################################################################################
#=
PROC MIXED data=df0;
CLASSES subject sequence period formulation;
MODEL  var = sequence period formulation/ DDFM=SATTERTH s;
RANDOM  formulation/TYPE=UN(1) SUB=subject G V;
RUN;

REML: 16.24111264
=#

lmm = LMM(@formula(var ~ sequence + period + formulation), df;
    random  = VarEffect(@covstr(formulation), VC),
    subject = :subject)
fit!(lmm)

#=
Linear Mixed Model: var ~ sequence + period + formulation
Random 1:
   Model: formulation
   Type: VarianceComponents (2)
   Coefnames: ["formulation: 1", "formulation: 2"] Repeated:
   Model: nothing
   Type: ScaledIdentity (1)
   Coefnames: -

Status: converged

   -2 logREML: 16.2411

   Fixed effects:

Name             Estimate     SE
(Intercept)      1.61         0.334718
sequence: 2      -0.170833    0.277378
period: 2        0.144167     0.289463
period: 3        0.08         0.117047
period: 4        0.144167     0.289463
formulation: 2   -0.0791667   0.277378

Random effects:

   θ vector: [0.447322, 0.367367, 0.185068]
Random 1   formulation: 1   var   0.200097
Random 1   formulation: 2   var   0.134959
Repeated   Var              var   0.0342502
=#
################################################################################

Best regards,
 
PharmCat
★    

Russia,
2021-02-08 22:33
(1403 d 08:24 ago)

@ PharmCat
Posting: # 22205
Views: 4,570
 

 Metida.jl: validation in progress…

Post reply
Version 0.2.4 released 10 minutes ago :)

The package is documented here.

Now covariance types available:
  • Scaled Identity (SI)
  • Diagonal (DIAG)
  • Autoregressive (AR)
  • Heterogeneous Autoregressive (ARH)
  • Compound Symmetry (CS)
  • Heterogeneous Compound Symmetry (CSH)
  • Autoregressive Moving Average (ARMA)
Validation in progress. Interim results say that REML function fully reproducible with SPSS for all datasets in Helmut paper.
 
PharmCat
★    

Russia,
2021-02-18 19:44
(1393 d 11:12 ago)

@ PharmCat
Posting: # 22213
Views: 4,495
 

 user structure

Post reply
Version 0.4.0 released ...

New covariance types available:
  • Toeplitz
  • ToeplitzParameterized
  • CustomCovarianceType
For bioequivalence CSH+DIAG avialible since 0.1.1 ;-)
 
PharmCat
★    

Russia,
2021-03-05 15:17
(1378 d 15:39 ago)

@ PharmCat
Posting: # 22246
Views: 4,465
 

 Validation results

Post reply
Version 0.6.0 released ...

Add cov types:
  • Heterogeneous Toeplitz (TOEPH)
  • Heterogeneous Toeplitz Parameterized (TOEPHP)
Validation results with reference datasets.*

  • Reference: Schütz, H., Labes, D., Tomashevskiy, M. et al.
    Reference Datasets for Studies in a Replicate Design Intended for Average Bioequivalence with Expanding Limits.
    Table II, Method B, SPSS section.
 
Weidson
☆    

Brazil,
2022-01-12 00:57
(1066 d 06:00 ago)

@ PharmCat
Posting: # 22729
Views: 3,945
 

 Validation results

Post reply
Dear Friends,

I'm not able to install Metida.jl package and I would like to test your functiolities. I downloaded the archive zip (Metida.jl-master.zip) from page: https://github.com/PharmCat/Metida.jl#readme on the botton Code and save it on my diretory (my documents). When I use the option Install from Achive File (.zip; targ.gz) of Rstudio Panel and choose the file at save directory the following message appears on RStudio's console :

Installing package into ‘C:/Users/weidson.carlo/Documents/R/win-library/3.6’
(as ‘lib’ is unspecified)
Warning in install.packages :
  cannot open compressed file 'Metida.jl-master/DESCRIPTION', probable reason 'No such file or directory'
Error in install.packages : cannot open the connection

Could someone help me with problem?:-)
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2022-01-12 02:02
(1066 d 04:55 ago)

@ Weidson
Posting: # 22730
Views: 3,998
 

 Julia ≠ 🇷

Post reply
Hi Weidson,

❝ When I use the option Install from Achive File (.zip; targ.gz) of Rstudio Panel and choose the file at save directory the following message appears on RStudio's console:

❝

❝ Installing package into ‘C:/Users/weidson.carlo/Documents/R/win-library/3.6’

❝ (as ‘lib’ is unspecified)

❝ Warning in install.packages :

❝   cannot open compressed file 'Metida.jl-master/DESCRIPTION', probable reason 'No such file or directory'

❝ Error in install.packages : cannot open the connection


As expected (see the subject line). It’s like presenting a Fortran-source to a C-compiler. Can’t work.
First you have to download and install Julia. Then start julia.exe (in the installation’s /bin-folder) and in its terminal:

julia> import Pkg; Pkg.add("Metida")

Looks on my machine like this:

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.1 (2021-12-22)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> import Pkg; Pkg.add("Metida")
    Updating registry at `C:\Users\HS\.julia\registries\General`
    Updating git-repo `https://github.com/JuliaRegistries/General.git`
  Downloaded artifact: Hwloc
    Updating `C:\Users\HS\.julia\environments\v1.7\Project.toml`
  [a1dec852] + Metida v0.12.4
    Updating `C:\Users\HS\.julia\environments\v1.7\Manifest.toml`
  [79e6a3ab] + Adapt v3.3.3
  [4fba245c] + ArrayInterface v3.2.2
  [2a0fbf3d] + CPUSummary v0.1.6
  [bbf7d656] + CommonSubexpressions v0.3.0
  [b429d917] + DensityInterface v0.4.0
  [85a47980] + Dictionaries v0.3.17
  [163ba53b] + DiffResults v1.0.3
  [b552c78f] + DiffRules v1.9.0
  [31c24e10] + Distributions v0.25.37
  [1a297f60] + FillArrays v0.12.7
  [6a86dc24] + FiniteDiff v2.9.0
  [f6369f11] + ForwardDiff v0.10.24
  [0e44f5e4] + Hwloc v2.0.0
  [615f187c] + IfElse v0.1.1
  [313cdc1a] + Indexing v1.1.1
  [d3d80556] + LineSearches v7.1.1
  [1914dd2f] + MacroTools v0.5.9
  [a1dec852] + Metida v0.12.4
  [075456b7] + MetidaBase v0.5.1
  [d41bc354] + NLSolversBase v7.8.2
  [77ba4419] + NaNMath v0.3.6
  [429524aa] + Optim v1.6.0
  [90014a1f] + PDMats v0.11.5
  [d96e819e] + Parameters v0.12.3
  [85a6dd25] + PositiveFactorizations v0.2.4
  [1fd47b50] + QuadGK v2.4.2
  [03a91e81] + SplitApplyCombine v1.2.0
  [aedffcd0] + Static v0.4.1
  [90137ffa] + StaticArrays v1.3.1
  [9d95f2ec] + TypedTables v1.4.0
  [3a884ed6] + UnPack v1.0.2
  [e33a78d0] + Hwloc_jll v2.7.0+0
  [4607b0f0] + SuiteSparse
Precompiling project...
  32 dependencies successfully precompiled in 26 seconds (51 already precompiled)

Before you can run PharmCat’s example in GitHub you have to install three other libraries.

julia> import Pkg; Pkg.add("CSV")
julia> import Pkg; Pkg.add("DataFrames")
julia> import Pkg; Pkg.add("CategoricalArrays")

Then with

julia> using Metida, CSV, DataFrames, CategoricalArrays
julia> df = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","df0.csv")) |> DataFrame
julia> transform!(df, :subject => categorical, renamecols=false)
julia> transform!(df, :period => categorical, renamecols=false)
julia> transform!(df, :sequence => categorical, renamecols=false)
julia> transform!(df, :formulation => categorical, renamecols=false)
julia> lmm = LMM(@formula(var~sequence+period+formulation), df;
       random = VarEffect(@covstr(formulation|subject), CSH),
       repeated = VarEffect(@covstr(formulation|subject), DIAG),
       )

You should get:

Linear Mixed Model: var ~ sequence
Random 1:
    Model: formulation|subject
    Type: CSH (3), Subjects: 5
Repeated:
    Model: formulation|subject
    Type: DIAG (2)
Blocks: 5, Maximum block size: 4
Not fitted.

Then:

julia> fit!(lmm)

You should get:

Linear Mixed Model: var ~ sequence + period + formulation
Random 1:
    Model: formulation|subject
    Type: CSH (3), Subjects: 5
Repeated:
    Model: formulation|subject
    Type: DIAG (2)
Blocks: 5, Maximum block size: 4
Status: converged See warnings in log.
    -2 logREML: 10.0652    BIC: 23.9282

    Fixed-effects parameters:
───────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)      1.57749     0.334543    4.72    <1e-05
sequence: 2     -0.170833    0.384381   -0.44    0.6567
period: 2        0.195984    0.117228    1.67    0.0946
period: 3        0.145014    0.109171    1.33    0.1841
period: 4        0.157363    0.117228    1.34    0.1795
formulation: 2  -0.0791667   0.0903709  -0.88    0.3810
───────────────────────────────────────────────────────
    Variance components:
    θ vector: [0.455584, 0.367656, 1.0, 0.143682, 0.205657]
  Random 1   σ² formulation: 1   var   0.207557
  Random 1   σ² formulation: 2   var   0.135171
  Random 1   ρ                   rho   1.0
  Residual   σ² formulation: 1   var   0.0206445
  Residual   σ² formulation: 2   var   0.0422948


The terminal REPL (Read-Eval-Print-Loop) of Julia is pain in the back. I suggest Atom/Juno or Visual Studio Code.

Didn’t try package JuliaCall, except in very simple examples. According to the documentation one could install Julia with the function install_julia() and packages with julia_install_package(). Never tried that.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
Weidson
☆    

Brazil,
2022-04-04 16:50
(983 d 15:06 ago)

@ Helmut
Posting: # 22908
Views: 3,384
 

 Julia ≠ 🇷

Post reply
Hi Helmut,

During instalation I had problem with this two comands lines:

❝ julia> import Pkg; Pkg.add("DataFrames")

❝ julia> import Pkg; Pkg.add("CategoricalArrays")


Not executed :confused:. You can see at this figure:
[image]

How can I solver this? I am a beginner user of this application and I need your help (My OS is Windows 10 (today :clap:))

best regards.


Weidson
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2022-01-13 16:18
(1064 d 14:39 ago)

@ Weidson
Posting: # 22731
Views: 3,845
 

 Julia vs Phoenix WinNonlin (examples)

Post reply
Hi Weidson,

examples in Julia. Note that I suppressed most of the intermediate output with a semicolon ; and the end of commands.

# The (in)famous EMA reference 'data set I'
using Metida, CSV, DataFrames, CategoricalArrays;
df = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","berds","rds1.csv"), delim=";") |> DataFrame;
transform!(df, :subject => categorical, renamecols=false);
transform!(df, :period => categorical, renamecols=false);
transform!(df, :sequence => categorical, renamecols=false);
transform!(df, :treatment => categorical, renamecols=false);
lmm = Metida.LMM(@formula(log(PK)~sequence+period+treatment), df;
                 random = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.CSH),
                 repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.DIAG),);
fit!(lmm);
pe = 100 * exp(coef(lmm)[end]);
ci_satt = 100 * [exp(Metida.confint(lmm, level = 0.9)[end][1]),
                 exp(Metida.confint(lmm, level = 0.9)[end][2])];
ci_res  = 100 * [exp(Metida.confint(lmm; ddf = :residual, level = 0.9)[end][1]),
                exp(Metida.confint(lmm; ddf = :residual, level = 0.9)[end][2])];
print(lmm);
println("Point estimate           =  ", pe); println("90% CI, Satterthwaite df = ", ci_satt); print("90% CI, Residual df      = ", ci_res);

Gives
Linear Mixed Model: :(log(PK)) ~ sequence + period + treatment
Random 1:
    Model: treatment|subject
    Type: CSH (3), Subjects: 77
Repeated:
    Model: treatment|subject
    Type: DIAG (2)
Blocks: 77, Maximum block size: 4
Status: converged See warnings in log.
    -2 logREML: 530.145    BIC: 558.562

    Fixed-effects parameters:
────────────────────────────────────────────────────────
                      Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────────────
(Intercept)      7.64968      0.151214   50.59    <1e-99
sequence: TRTR  -0.0204706    0.197223   -0.10    0.9173
period: 2        0.0446376    0.0636371   0.70    0.4830
period: 3        0.00659015   0.0644219   0.10    0.9185
period: 4        0.0727201    0.0639714   1.14    0.2556
treatment: T     0.145464     0.0465012   3.13    0.0018
────────────────────────────────────────────────────────
    Variance components:
    θ vector: [0.852995, 0.828407, 1.0, 0.449575, 0.342628]
  Random 1   σ² treatment: R   var   0.7276
  Random 1   σ² treatment: T   var   0.686258
  Random 1   ρ                 rho   1.0
  Residual   σ² treatment: R   var   0.202118
  Residual   σ² treatment: T   var   0.117394

Point estimate           =  115.65764182681824
90% CI, Satterthwaite df = [107.10446639731114, 124.89385889211053]
90% CI, Residual df      = [107.11502945421051, 124.88154259117144]

In Phoenix WinNonlin I got with Heterogenous Compound Symmetry:
cshSD(1)_11                      0.85274463
cshSD(2)_11                      0.82831758
cshCorr_11                       1.0005616
Var(period*treatment*subject)_21 0.20252422
Var(period*treatment*subject)_22 0.11754154
PE                               115.656778312678
90% CL (Satterthwaite df)        107.051689268515, 124.953566459994
Close but not identical.
Residual df did not work: *** WARNING 11088: Residual DF is 0. Satterthwaite DF method will be used.
However, with Banded No-Diagonal Factor Analytic (f=2) like in the FDA’s SAS-code FA0(2):
Var(period*treatment*subject)_21 0.20211806
Var(period*treatment*subject)_2  0.11739425
PE                               115.657641801223
90% CL (Satterthwaite df)        107.104404788005, 124.893930679106
That’s much closer to PharmCat’s Metida.jl.

Now your data:
df = CSV.read("E:/Public/Documents/BEBAC/R/Weidson/IPC_Drug.Cmax.csv", delim="\t", DataFrame);
transform!(df, :Subject => categorical, renamecols=false);
transform!(df, :Period => categorical, renamecols=false);
transform!(df, :Sequence => categorical, renamecols=false);
transform!(df, :Treatment => categorical, renamecols=false);
lmm = Metida.LMM(@formula(log(PK)~Sequence+Period+Treatment), df;
                 random = Metida.VarEffect(Metida.@covstr(Treatment|Subject), Metida.CSH),
                 repeated = Metida.VarEffect(Metida.@covstr(Treatment|Subject), Metida.DIAG),);
fit!(lmm);
pe = 100 * exp(coef(lmm)[end]);
ci_satt = 100 * [exp(Metida.confint(lmm, level = 0.9)[end][1]),
                 exp(Metida.confint(lmm, level = 0.9)[end][2])];
ci_res  = 100 * [exp(Metida.confint(lmm; ddf = :residual, level = 0.9)[end][1]),
                exp(Metida.confint(lmm; ddf = :residual, level = 0.9)[end][2])];
print(lmm);
println("Point estimate           =  ", pe); println("90% CI, Satterthwaite df = ", ci_satt); print("90% CI, Residual df      = ", ci_res);

Linear Mixed Model: :(log(PK)) ~ Sequence + Period + Treatment
Random 1:
    Model: Treatment|Subject
    Type: CSH (3), Subjects: 50
Repeated:
    Model: Treatment|Subject
    Type: DIAG (2)
Blocks: 50, Maximum block size: 4
Status: converged (No Errors)
    -2 logREML: 558.71    BIC: 585.101

    Fixed-effects parameters:
───────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)      4.20291      0.249867  16.82    <1e-62
Sequence: TRTR   0.163992     0.314302   0.52    0.6018
Period: 2        0.310384     0.150476   2.06    0.0391
Period: 3        0.0440086    0.140033   0.31    0.7533
Period: 4        0.183732     0.150476   1.22    0.2221
Treatment: T    -0.238037     0.11275   -2.11    0.0348
───────────────────────────────────────────────────────
    Variance components:
    θ vector: [0.99404, 1.12122, 0.951444, 0.74508, 0.670061]
  Random 1   σ² Treatment: R   var   0.988116
  Random 1   σ² Treatment: T   var   1.25713
  Random 1   ρ                 rho   0.951444
  Residual   σ² Treatment: R   var   0.555144
  Residual   σ² Treatment: T   var   0.448982

Point estimate           =  78.81736099725553
90% CI, Satterthwaite df = [65.23698724552828, 95.22475909550096]
90% CI, Residual df      = [65.4172844390996, 94.96230924038038]

In Phoenix WinNonlin with Heterogenous Compound Symmetry:
cshSD(1)_11                      0.99404029
cshSD(2)_11                      1.1212233
cshCorr_11                       0.95144445
Var(period*treatment*subject)_21 0.55514404
Var(period*treatment*subject)_22 0.44898200
PE                               78.8173609972554
90% CL (Satterthwaite df)        65.2370226163479, 95.2247074656494
Again, close but not identical; residual df did not work.
This time with FA0(2)
Var(Period*Treatment*Subject)_21 0.55514404
Var(Period*Treatment*Subject)_22 0.44898200
PE                               78.8173609972554
90% CL (Satterthwaite df)        65.2369873029156, 95.2247590117338

An alternative is PharmCat’s ReplicateBE.jl (great documentation):
using Pkg; Pkg.add("ReplicateBE") # download/install first
using CSV, DataFrames, Metida, ReplicateBE, StatsBase;
df = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","berds","rds1.csv"), delim=";") |> DataFrame;
be = rbe!(df, dvar = :logPK, subject = :subject, formulation = :treatment, period = :period, sequence = :sequence);
println(design(be)); print(be)

Trial design
  Observations:          298
  Subjects:              77
  Sequences:             2
  Periods:               4
  Formulations:          2
  Subjects(Formulation): 77/77
  Rank X:                6
  Rank XZ:               156
  DF (robust):           75
  DF (contain):          142
  DF (Subj(Form)):       151
Bioequivalence Linear Mixed Effect Model (status:converged)
Hessian not positive defined!
Rho ~ 1.0!

-2REML: 530.145    REML: -265.072

Fixed effect:
────────────────────────────────────────────────────────────────────────────────────────
Effect           Value        SE          F           DF        t           P|t|
────────────────────────────────────────────────────────────────────────────────────────
(Intercept)      7.64968      0.151214    2559.2      81.3295   50.5885     3.08793E-63*
sequence: TRTR   -0.0204707   0.197223    0.0107733   74.735    -0.103795   0.91761
period: 2        0.0446377    0.0636371   0.492019    214.16    0.701441    0.483789
period: 3        0.00659021   0.0644219   0.0104648   189.589   0.102298    0.918629
period: 4        0.0727201    0.0639714   1.29222     214.259   1.13676     0.256909
treatment: T     0.145464     0.0465012   9.78552     208.081   3.12818     0.00201093*
────────────────────────────────────────────────────────────────────────────────────────
Intra-individual variance:
treatment: R   0.202118   CVᵂ:   47.33   %
treatment: T   0.117394   CVᵂ:   35.29   %

Inter-individual variance:
treatment: R   0.7276
treatment: T   0.686258
ρ:             1.0        Cov: 0.706627

Confidence intervals(90%):
treatment: R / treatment: T
Ratio: 86.46, CI: 80.07 - 93.37 (%)
treatment: T / treatment: R
Ratio: 115.66, CI: 107.1 - 124.89 (%)

using CSV, DataFrames, ReplicateBE, StatsBase;
df = CSV.read("E:/Public/Documents/BEBAC/R/Weidson/IPC_Drug.Cmax.csv", delim="\t", DataFrame);
df[!, :logPK] .= ( v <= 0. ? 0. : log(v) for v in df.PK);
be = rbe!(df, dvar = :logPK, subject = :Subject, formulation = :Treatment, period = :Period, sequence = :Sequence);
println(design(be)); print(be)

Trial design
  Observations:          200
  Subjects:              50
  Sequences:             2
  Periods:               4
  Formulations:          2
  Subjects(Formulation): 50/50
  Rank X:                6
  Rank XZ:               102
  DF (robust):           48
  DF (contain):          98
  DF (Subj(Form)):       97
Bioequivalence Linear Mixed Effect Model (status:converged)

-2REML: 558.71    REML: -279.355

Fixed effect:
─────────────────────────────────────────────────────────────────────────────────────
Effect           Value       SE         F           DF        t          P|t|
─────────────────────────────────────────────────────────────────────────────────────
(Intercept)      4.20291     0.249867   282.932     56.3186   16.8206    1.26492E-23*
Sequence: TRTR   0.163992    0.314302   0.27224     48.0      0.521766   0.604233
Period: 2        0.310384    0.150476   4.25464     116.951   2.06268    0.0413577*
Period: 3        0.0440086   0.140033   0.0987677   93.3646   0.314273   0.754015
Period: 4        0.183732    0.150476   1.49084     116.951   1.221      0.224542
Treatment: T     -0.238037   0.11275    4.45715     48.0      -2.1112    0.039989*
─────────────────────────────────────────────────────────────────────────────────────
Intra-individual variance:
Treatment: R   0.555144   CVᵂ:   86.15   %
Treatment: T   0.448982   CVᵂ:   75.28   %

Inter-individual variance:
Treatment: R   0.988116
Treatment: T   1.25714
ρ:             0.951444   Cov: 1.06042

Confidence intervals(90%):
Treatment: R / Treatment: T
Ratio: 126.88, CI: 105.01 - 153.29 (%)
Treatment: T / Treatment: R
Ratio: 78.82, CI: 65.24 - 95.22 (%)


[image]

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
PharmCat
★    

Russia,
2022-05-20 19:54
(937 d 12:03 ago)

@ Helmut
Posting: # 23009
Views: 3,221
 

 Julia vs Phoenix WinNonlin (examples)

Post reply
Hi Helmut!

Thank you for this great explanation and comparison! I didn't test Metida vs Phoenix, only vs SPSS and SAS. And as I know CSH and FA(2) are completely identical for data with 2 levels in random model (in this case G have 2x2 dim). And I think the difference in results with Phoenix for CSH and FA(2) is the issue (unfortunately it can't be submitted on GitHub). Also if we look at SPSS/SAS documentation we can find that rho optimized in '-1<rho<1' range. And it can't be more 1.0, but we have 'cshCorr_11 1.0005616', rly? It is difficult to make good optimization with Newton methods if parameters are limited, in Metida sigmoid function is used to change the optimization region for rho from -1:1 to -inf:+inf. In SAS and SPSS this feature works fine and I never get rho more than 1 and less -1. But maybe in Phoenix, this limiting work not so good if we see 1.0005616. From one side we can just look at REML values and get a better fitting, from the other - only the fact that we see covariance coefficient = 1.0005616 in the variance-covariance matrix leads that other computations and REML estimation have no sense.
In this situation maybe it is a good choice - never use Heterogenous Compound Symmetry in Phoenix.

If we see at point estimate and variance values for Metida and Phoenix with FA(0) - they are identical. But the minimal difference may be in DF estimation. I observe differences for DF in all software - between SPSS, SAS and Phoenix, sometimes it depends on data (especially if the covariance matrix is ill-conditioned). Maybe the cause is different approaches for derivating REML function, but it is only an assumption.

So also I can say that spatial covariance was added to Metida, and interface to any custom structure. Now I try to work on bootstrapping and multiple imputations for Metida. Maybe some other packages can be interesting:

ODMXMLTools.jl - experimental package for ODM-XML
MetidaNCA.jl - Julia NCA package
 
UA Flag
Activity
 Admin contact
23,336 posts in 4,902 threads, 1,698 registered users;
55 visitors (0 registered, 55 guests [including 12 identified bots]).
Forum time: 06:57 CET (Europe/Vienna)

Only dead fish go with the current.    Scuba divers' proverb

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5