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

using Metida; # hide
include(joinpath(dirname(pathof(Metida)), "..", "test", "validation_init.jl")); # hide
include(joinpath(dirname(pathof(Metida)), "..", "test", "validation_s3.jl"))

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.