Julia vs Phoenix WinNonlin (examples) [Software]

posted by Helmut Homepage – Vienna, Austria, 2022-01-13 16:18 (1218 d 14:22 ago) – Posting: # 22731
Views: 4,449

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

Complete thread:

UA Flag
Activity
 Admin contact
23,424 posts in 4,927 threads, 1,672 registered users;
159 visitors (0 registered, 159 guests [including 15 identified bots]).
Forum time: 07:40 CEST (Europe/Vienna)

A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions.    Micheal J. Moroney

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