GitHub

Validation

Validation provided with 3 sections:

  • REML validation for public datasets with Metida & SPSS
  • Parameters validation for public datasets Metida & SPSS & MixedModels
  • Validation with bioequivalence datasets with Metida & SPSS

To run validation:

using Metida; include(joinpath(dirname(pathof(Metida)), "..", "test", "validation.jl"))

Section 1: REML validation for public datasets Metida & SPSS

REML result table

ModelDataSetUsed cov. typesREML MetidaREML SPSS
1sleepstudy.csvSI/SI1729.49256023670251729.492560
2sleepstudy.csvCS/SI1904.32651707221321904.327
3sleepstudy.csvCSH/SI1772.09532519970461772.095
4sleepstudy.csvARH/SI1730.18954273983221730.189543
5Pastes.csvSI,SI/SI246.99074585348623246.990746
6Pastes.csvARMA/SI246.81895071012508246.818951
7Penicillin.csvSI,SI/SI330.86058899109184330.860589
8RepeatedPulse.csvSI/AR453.3395435627574453.339544
9RepeatedPulse.csv0/AR471.85107712169827471.851077
10RepeatedPulse.csvAR/SI453.3395560121246453.339555

sleepstudy.csv

Model 1
lmm = LMM(@formula(Reaction~Days), df;
  random = VarEffect(@covstr(1|Subject), SI),
  )
  fit!(lmm)

SPSS:

MIXED Reaction BY Days
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=Days | SSTYPE(3)
  /METHOD=REML
  /RANDOM=INTERCEPT | SUBJECT(Subject) COVTYPE(ID).
Model 2
lmm = LMM(@formula(Reaction~1), df;
  random = VarEffect(Metida.@covstr(Days|Subject), CS),
  )
  fit!(lmm)

SPSS:

MIXED Reaction BY Days
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=INTERCEPT | SSTYPE(3)
  /METHOD=REML
  /RANDOM=Days | SUBJECT(Subject) COVTYPE(CS).
Model 3
lmm = LMM(@formula(Reaction~1), df;
  random = VarEffect(@covstr(Days|Subject), CSH)
  )
  fit!(lmm)

SPSS:

MIXED Reaction BY Days
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=INTERCEPT | SSTYPE(3)
  /METHOD=REML
  /RANDOM=Days | SUBJECT(Subject) COVTYPE(CSH).
Model 4
lmm = LMM(@formula(Reaction~1), df;
  random = VarEffect(@covstr(Days|Subject), ARH)
  )
   fit!(lmm)

SPSS:

MIXED Reaction BY Days
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=INTERCEPT | SSTYPE(3)
  /METHOD=REML
  /RANDOM=Days | SUBJECT(Subject) COVTYPE(ARH1).

pastes.csv

Model 5
lmm =  LMM(@formula(strength~1), df;
random = [VarEffect(@covstr(1|batch), SI),  VarEffect(@covstr(1|batch & cask), SI)]
)
fit!(lmm)

SPSS:

MIXED strength
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /RANDOM=INTERCEPT | SUBJECT(batch) COVTYPE(ID)
  /RANDOM=INTERCEPT | SUBJECT(cask * batch) COVTYPE(ID).
Model 6
lmm =  LMM(@formula(strength~1), df;
random = VarEffect(Metida.@covstr(cask|batch),  ARMA),
)
fit!(lmm)

SPSS:

MIXED strength by cask
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /RANDOM=cask | SUBJECT(batch) COVTYPE(ARMA11).

penicillin.csv

Model 7
lmm =  LMM(@formula(diameter ~ 1), df;
random = [VarEffect(@covstr(1|plate), SI), VarEffect(@covstr(1|sample), SI)]
)
fit!(lmm)

SPSS:

MIXED diameter
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /RANDOM=INTERCEPT | SUBJECT(plate) COVTYPE(ID)
  /RANDOM=INTERCEPT | SUBJECT(sample) COVTYPE(ID).

RepeatedPulse.csv

Model 8
sort!(df, :Day)
lmm =  LMM(@formula(Pulse~1), df;
random =  VarEffect(Metida.@covstr(Time|Time),  SI),
repeated =  VarEffect(Metida.@covstr(Day|Time),  AR),
)
fit!(lmm)

SPSS:

MIXED Pulse BY Day Time
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /RANDOM=Time | SUBJECT(Time) COVTYPE(ID)
  /REPEATED = Day | SUBJCET(Time) COVTYPE(AR1).
Model 9
sort!(df, :Day)
lmm =  LMM(@formula(Pulse~1), df;
repeated = VarEffect(Metida.@covstr(Day|Time),  AR),
)
 fit!(lmm)

SPSS:

MIXED Pulse BY Day Time
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /REPEATED = Day | SUBJCET(Time) COVTYPE(AR1).
Model 10
sort!(df, :Day)
lmm =  LMM(@formula(Pulse~1), df;
random =  VarEffect(Metida.@covstr(Day|Time),  AR),
)
 fit!(lmm)

SPSS:

MIXED Pulse BY Day Time
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /RANDOM=Day | SUBJECT(Time) COVTYPE(AR1).

Section 2: Parameters validation for public datasets Metida & SPSS & MixedModels

Model 7

Metida:

lmm = LMM(@formula(diameter ~ 1), df;
random = [VarEffect(@covstr(1|plate), SI), VarEffect(@covstr(1|sample), SI)]
)
fit!(lmm)

MixedModels:

fm = @formula(diameter ~ 1 + (1|plate) + (1|sample))
mm  = fit(MixedModel, fm, df, REML=true)

SPSS:

MIXED diameter
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=| SSTYPE(3)
  /METHOD=REML
  /EMMEANS=TABLES(OVERALL)
  /RANDOM=INTERCEPT | SUBJECT(plate) COVTYPE(ID)
  /RANDOM=INTERCEPT | SUBJECT(sample) COVTYPE(ID).
ModelParameterValue MetidaValue MMValue SPSS
7(Intercept) estimate22.972222.972222.972
7(Intercept) SE0.8085730.8085720.809
7plate σ²0.7169080.7169080.716908
7sample σ²3.730923.7309013.730918
7Residual σ²0.3024150.3024160.302415

Section 3: Validation with bioequivalence datasets with Metida & SPSS

Model BE-B

lmm =  LMM(@formula(lnpk~sequence+period+treatment), dfrds;
random =  VarEffect(Metida.@covstr(1|subject),  SI),
)
fit!(lmm)

Model BE-C

lmm =  LMM(@formula(lnpk~sequence+period+treatment), dfrds;
    random =  VarEffect(Metida.@covstr(treatment|subject),  CSH),
    repeated =  VarEffect(Metida.@covstr(treatment|subject),  DIAG),
    )
     fit!(lmm)

Typical SPSS code

MIXED lnpk BY period sequence treatment subject
  /CRITERIA=DFMETHOD(SATTERTHWAITE) CIN(95) MXITER(100) MXSTEP(10) SCORING(1)
    SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=period sequence treatment | SSTYPE(3)
  /METHOD=REML
  /RANDOM= subject(sequence) | COVTYPE(ID)
  /EMMEANS=TABLES(treatment) COMPARE REFCAT(FIRST) ADJ(LSD).

MIXED lnpk BY period treatment sequence subject
  /CRITERIA=CIN(90) MXITER(200) MXSTEP(20) SCORING(2) SINGULAR(0.000000000001) HCONVERGE(0,
    RELATIVE) LCONVERGE(0.0000000000001, RELATIVE) PCONVERGE(0, RELATIVE)
  /FIXED=period treatment sequence | SSTYPE(3)
  /METHOD=REML
  /RANDOM=treatment | SUBJECT(subject) COVTYPE(CSH)
  /REPEATED=treatment | SUBJECT(subject*period) COVTYPE(DIAG)
  /EMMEANS=TABLES(treatment) COMPARE REFCAT(FIRST) ADJ(LSD).

Results

include(joinpath(dirname(pathof(Metida)), "..", "test", "validation_s3.jl"))
Test Summary:                                               | Pass  Total  Time
  RDS Test                                                  |  120    120  7.7s

Bioequivalence Reference Datasets - REML - type B
.-----.----------.----------.-------------.
| RDS |   Metida |    SPSS  |        DIFF |
|  N  |   REML B |   REML B |             |
:-----+----------+----------+-------------:
|   1 |  536.201 |  536.201 |          OK |
|   2 | -28.5882 | -28.5882 |          OK |
|   3 |  436.512 |  436.512 |          OK |
|   4 |  316.051 |  316.051 |          OK |
|   5 | -74.5417 | -74.5417 |          OK |
|   6 |  536.201 |  536.201 |          OK |
|   7 |  1388.54 |  1388.54 |          OK |
|   8 |  2345.96 |  2345.96 |          OK |
|   9 |  2985.72 |  2985.72 |          OK |
|  10 | -15.8938 | -15.8938 |          OK |
|  11 |  253.051 |  253.051 |          OK |
|  12 |  1143.36 |  1143.36 |          OK |
|  13 |   2090.1 |   2090.1 |          OK |
|  14 |  1050.31 |  1050.31 |          OK |
|  15 |   2090.1 |   2090.1 |          OK |
|  16 |  326.536 |  326.536 |          OK |
|  17 |  78.8447 |  78.8447 |          OK |
|  18 |  937.986 |  937.986 |          OK |
|  19 |  815.225 |  815.225 |          OK |
|  20 |  827.985 |  827.985 |          OK |
|  21 |  474.088 |  474.088 |          OK |
|  22 |  249.254 |  249.254 |          OK |
|  23 |  129.163 |  129.163 |          OK |
|  24 |  283.595 |  283.595 |          OK |
|  25 |  678.745 |  678.745 |          OK |
|  26 |  435.445 |  435.445 |          OK |
|  27 |  1125.08 |  1125.09 | -0.00895071 |
|  28 |  331.808 |  331.808 |          OK |
|  29 |  45.3049 |  45.3049 |          OK |
|  30 |   32.418 |   32.418 |          OK |
'-----'----------'----------'-------------'

Bioequivalence Reference Datasets - REML - type C
.-----.----------.----------.-------------.-------.
| RDS |   Metida |    SPSS  |        DIFF | Comm. |
|  N  |   REML C |   REML C |             |       |
:-----+----------+----------+-------------+-------:
|   1 |  530.145 |  530.145 |          OK |     * |
|   2 | -30.6746 | -30.6746 |  6.15692e-6 |    ** |
|   3 |  425.447 |  425.447 |          OK |     * |
|   4 |  314.222 |  314.222 |          OK |       |
|   5 |   -74.88 |   -74.88 |          OK |    ** |
|   6 |  530.145 |  530.145 |          OK |     * |
|   7 |  1387.09 |  1387.09 |          OK |    ** |
|   8 |   2342.6 |   2342.6 |          OK |       |
|   9 |  2983.26 |  2983.26 |          OK |       |
|  10 | -16.4173 | -16.4173 |          OK |    ** |
|  11 |  250.945 |  250.945 |          OK |       |
|  12 |  1140.38 |  1140.38 |          OK |       |
|  13 |  2087.48 |  2087.48 |          OK |     * |
|  14 |  1012.35 |  1012.35 |          OK |       |
|  15 |  2087.48 |  2087.48 |          OK |     * |
|  16 |  323.998 |  323.998 |          OK |     * |
|  17 |   77.569 |   77.569 |          OK |       |
|  18 |  904.874 |  904.874 |          OK |       |
|  19 |   782.94 |   782.94 |          OK |       |
|  20 |  796.312 |  796.312 |          OK |       |
|  21 |  470.591 |  470.591 |          OK |       |
|  22 |   248.99 |   248.99 |          OK |    ** |
|  23 |  119.806 |  119.806 |          OK |       |
|  24 |  274.306 |  274.306 |          OK |     * |
|  25 |  660.047 |  660.047 |          OK |     * |
|  26 |  433.841 |  433.841 |          OK |       |
|  27 |  1123.66 |  1123.66 | -0.00870748 |       |
|  28 |  329.257 |  329.257 |          OK |       |
|  29 |  26.9661 |  26.9661 |          OK |       |
|  30 |  26.3165 |   14.944 |     11.3725 |    ** |
'-----'----------'----------'-------------'-------'

*  - The final Hessian matrix is not positive definite although all convergence criteria
are satisfied. The MIXED procedure continues despite this warning. Validity of subsequent
results cannot be ascertained.

** - Iteration was terminated but convergence has not been achieved. The MIXED procedure
continues despite this warning. Subsequent results produced are based on the last iteration.
Validity of the model fit is uncertain.

DataSet 27 - MixedModels.jl result:

Linear mixed model fit by REML
 lnpk ~ 1 + sequence + period + treatment + (1 | subject)
 REML criterion at convergence: 1125.083432363973

Variance components:
            Column   Variance Std.Dev.
subject  (Intercept)  0.489831 0.699879
Residual              0.110231 0.332011
 Number of obs: 623; levels of grouping factors: 312

  Fixed-effects parameters:
─────────────────────────────────────────────────────
                   Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)    4.58171     0.0846385  54.13    <1e-99
sequence: RT   0.135734    0.119701    1.13    0.2568
sequence: TR  -0.0165249   0.119701   -0.14    0.8902
sequence: TT   0.103687    0.124075    0.84    0.4033
period: 2     -0.010685    0.0266212  -0.40    0.6881
treatment: T  -0.175322    0.0377032  -4.65    <1e-05
─────────────────────────────────────────────────────

Bioequivalence Reference Datasets - 90% Confidence intervals
.-----.----------.----------.----------.----------.--------.--------.
| RDS | Metida B | Metida B | Metida C | Metida C |   REF* |   REF* |
|  N  |      LCI |      UCI |      LCI |     UCI  |    LCI |    UCI |
:-----+----------+----------+----------+----------+--------+--------:
|   1 |   107.17 |   124.97 |    107.1 |   124.89 | 107.17 | 124.97 |
|   2 |    97.32 |   107.46 |    97.02 |   107.79 |  97.32 | 107.46 |
|   3 |   113.31 |   136.73 |   113.17 |   136.49 | 113.31 | 136.73 |
|   4 |    117.9 |   159.69 |   118.67 |   158.65 |  117.9 | 159.69 |
|   5 |   103.82 |   112.04 |    103.8 |   112.06 | 103.82 | 112.04 |
|   6 |    80.02 |    93.31 |    80.07 |    93.37 |  80.02 |  93.31 |
|   7 |    86.46 |    92.81 |    86.44 |    92.83 |  86.46 |  92.81 |
|   8 |    75.69 |     87.6 |    75.57 |    87.74 |  75.69 |   87.6 |
|   9 |    75.69 |     87.6 |    75.57 |    87.74 |  75.69 |   87.6 |
|  10 |    96.27 |   107.59 |    95.96 |   107.94 |  96.27 | 107.59 |
|  11 |    80.64 |   100.38 |    79.55 |   101.75 |  80.64 | 100.38 |
|  12 |    90.35 |   157.88 |    89.13 |   159.16 |  90.35 | 157.88 |
|  13 |    72.87 |    85.51 |    72.94 |     85.6 |  72.87 |  85.51 |
|  14 |    69.21 |   121.27 |    65.23 |   137.67 |  69.21 | 121.27 |
|  15 |    72.87 |    85.51 |    72.94 |     85.6 |  72.87 |  85.51 |
|  16 |    69.54 |    89.37 |    69.32 |    89.66 |  69.54 |  89.37 |
|  17 |   115.97 |   155.09 |   113.64 |   158.87 | 115.97 | 155.09 |
|  18 |    59.13 |    107.2 |    59.84 |   130.49 |  59.13 |  107.2 |
|  19 |    53.85 |    98.77 |    49.86 |   112.63 |  53.85 |  98.77 |
|  20 |    50.92 |    95.62 |    47.26 |   109.53 |  50.92 |  95.62 |
|  21 |   111.72 |   127.73 |   111.35 |   128.69 | 111.72 | 127.73 |
|  22 |    77.98 |   106.09 |    77.98 |   106.09 |  77.98 | 106.09 |
|  23 |    97.13 |   128.41 |    95.76 |   130.01 |  97.13 | 128.41 |
|  24 |    87.24 |   109.85 |    86.64 |   110.54 |  87.24 | 109.85 |
|  25 |    77.93 |     98.1 |    77.93 |     98.1 |  77.93 |   98.1 |
|  26 |   133.51 |   171.42 |   132.26 |   172.37 | 133.51 | 171.42 |
|  27 |    78.86 |     89.3 |    78.75 |    89.49 |  78.86 |   89.3 |
|  28 |    87.86 |   100.07 |    87.49 |    100.5 |  87.86 | 100.07 |
|  29 |    88.43 |   121.59 |    81.12 |   128.18 |  88.43 | 121.59 |
|  30 |    79.58 |   108.07 |    79.56 |   107.76 |  79.58 | 108.07 |
'-----'----------'----------'----------'----------'--------'--------'

* Reference: Schütz H, Labes D, Tomashevskiy M, la Parra MG, Shitova A,
Fuglsang A. Reference Datasets for Studies in a Replicate Design Intended for
Average Bioequivalence with Expanding Limits. AAPS J. 2020 Feb 7;22(2):44.
doi: 10.1208/s12248-020-0427-6. PMID: 32034551.
Table II, Method B, SPSS section.

Full SPSS code provided in validation folder (here).

SPSS output in DOCX format.

Validation dataset available here, 122482020427MOESM2ESM.xls.

Validation report

Validation and report can be done on local machine with Weave.jl and Pandoc.

using Weave, Metida
weave(joinpath(dirname(pathof(Metida)), "..", "test", "validation_report.jmd");
doctype = "pandoc2pdf",
out_path = :pwd,
pandoc_options=["--toc", "-V colorlinks=true" , "-V linkcolor=blue", "-V urlcolor=red", "-V toccolor=gray"])

Report will be saved in Julia working directory. For your own purpose you can edit validation_report.jmd template.