Julia vs Phoenix WinNonlin (examples) [Software]

posted by Helmut Homepage – Vienna, Austria, 2022-01-13 15:18 (254 d 09:35 ago) – Posting: # 22731
Views: 1,163

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([email protected](treatment|subject), Metida.CSH),
                 repeated = Metida.VarEffect([email protected](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([email protected](Treatment|Subject), Metida.CSH),
                 repeated = Metida.VarEffect([email protected](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

Complete thread:

UA Flag
Activity
 Admin contact
22,385 posts in 4,684 threads, 1,594 registered users;
online 4 (0 registered, 4 guests [including 3 identified bots]).
Forum time: Sunday 01:53 CEST (Europe/Vienna)

You really don’t know what you don’t know until you write about it.
Then, everyone knows what you don’t know.    Rod Machado

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