Bioequivalence and Bioavailability Forum • Julia vs Phoenix WinNonlin (examples)

Julia vs Phoenix WinNonlin (examples) [Software]

posted by Helmut – Vienna, Austria, 2022-01-13 16:18 (1005 d 16:10 ago) – Posting: # 22731
Views: 3,282

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 🖖🏼 Довге життя Україна!
Helmut Schütz

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

Complete thread:

Metida.jl: experimental Julia LinearMixedModel software PharmCat 2020-10-03 19:59 [Software]
- Metida.jl: validation in progress… PharmCat 2021-02-08 21:33
  - user structure PharmCat 2021-02-18 18:44
    - Validation results PharmCat 2021-03-05 14:17
      - Validation results Weidson 2022-01-11 23:57
        
        Julia ≠ 🇷 Helmut 2022-01-12 01:02
        
        Julia ≠ 🇷 Weidson 2022-04-04 14:50
        
        Julia vs Phoenix WinNonlin (examples)Helmut 2022-01-13 15:18
        
        Julia vs Phoenix WinNonlin (examples) PharmCat 2022-05-20 17:54