Installation

import Pkg; Pkg.add("Metida")

Simple example

Step 1: Load data

Load provided data with CSV and DataFrames:

using Metida, CSV, DataFrames, CategoricalArrays

df = CSV.File(joinpath(dirname(pathof(Metida)), "..", "test", "csv", "df0.csv")) |> DataFrame;
Note

Check that all categorical variables are categorical.

transform!(df, :subject => categorical, renamecols=false)
transform!(df, :period => categorical, renamecols=false)
transform!(df, :sequence => categorical, renamecols=false)
transform!(df, :formulation => categorical, renamecols=false)

Step 2: Make model

Make model with @formula macro from StatsModels. Define random and repreated effects with Metida.VarEffect using Metida.@covstr macros. Left side of @covstr is model of effect and right side is a effect itself. Metida.HeterogeneousCompoundSymmetry and Metida.Diagonal in example bellow is a model of variance-covariance structure. See also Metida.@lmmformula macro.

Note

In some cases levels of repeated effect should not be equal inside each level of subject or model will not have any sense. For example, it is assumed that usually CSH or UN (Unstructured) using with levels of repeated effect is different inside each level of subject. Metida does not check this!

lmm = LMM(@formula(var~sequence+period+formulation), df;
random = VarEffect(@covstr(formulation|subject), CSH),
repeated = VarEffect(@covstr(formulation|subject), DIAG));
Linear Mixed Model: var ~ sequence + period + formulation
Random 1: 
    Model: formulation|subject
    Type: CSH (3), Subjects: 5
Repeated: 
    Model: formulation|subject
    Type: DIAG (2)
Blocks: 5, Maximum block size: 4
Not fitted.

Step 3: Fit

Just fit the model.

fit!(lmm)
Linear Mixed Model: var ~ sequence + period + formulation
Random 1: 
    Model: formulation|subject
    Type: CSH (3), Subjects: 5
Repeated: 
    Model: formulation|subject
    Type: DIAG (2)
Blocks: 5, Maximum block size: 4
Status: converged (No Errors)
    -2 logREML: 10.0652    BIC: 23.9282

    Fixed-effects parameters:
───────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)      1.57749     0.334543    4.72    <1e-05
sequence: 2     -0.170833    0.384381   -0.44    0.6567
period: 2        0.195984    0.117228    1.67    0.0946
period: 3        0.145014    0.109171    1.33    0.1841
period: 4        0.157363    0.117228    1.34    0.1795
formulation: 2  -0.0791667   0.0903709  -0.88    0.3810
───────────────────────────────────────────────────────
    Variance components:
    θ vector: [0.455584, 0.367656, 1.0, 0.143682, 0.205657]
  Random 1   σ² formulation: 1   var   0.207557
  Random 1   σ² formulation: 2   var   0.135171
  Random 1   ρ                   rho   1.0
  Residual   σ² formulation: 1   var   0.0206445
  Residual   σ² formulation: 2   var   0.0422948
Check warnings and errors in log.
lmm.log
3-element Vector{Metida.LMMLogMsg}:
   INFO  : Initial θ: [0.21085850593184408, 0.21085850593184408, 0.0001, 0.21085850593184408, 0.21085850593184408]

   INFO  : Resulting θ: [0.4555842906225857, 0.3676561978014238, 0.9999999999972058, 0.1436824286508881, 0.2056573902850657]; 26 iterations.

   INFO  : Model fitted.
Type III Tests of Fixed Effects
Warning

Experimental

typeiii(lmm)
 ------------- ---------- ----- --------- -----------
         Name          F   ndf       ddf        pval
 ------------- ---------- ----- --------- -----------
  (Intercept)    22.2346   1.0   3.16597   0.0159885
     sequence   0.197525   1.0   3.00055    0.686828
       period     1.0789   3.0   8.66853    0.407683
  formulation   0.767409   1.0   5.46311    0.417867
 ------------- ---------- ----- --------- -----------

Model construction

Metida.LMMType
LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing)

Make Linear-Mixed Model object.

model: is a fixed-effect model (@formula)

data: tabular data

contrasts: contrasts for fixed factors

random: vector of random effects or single random effect

repeated: is a repeated effect (only one)

See also: @lmmformula

source
  • model - example: @formula(var ~ sequence + period + formulation)

  • random - effects can be specified like this: VarEffect(@covstr(formulation|subject), CSH). @covstr is a effect model: @covstr(formulation|subject). CSH is a CovarianceType structure. Premade constants: SI, DIAG, AR, ARH, CS, CSH, ARMA, TOEP, TOEPH, UN. If not specified only repeated used.

  • repeated - can be specified like random effect. If not specified VarEffect(@covstr(1|1), SI) used. If no repeated effects specified vector of ones used.

Random/repeated model

Metida.@covstrMacro
@covstr(ex)

Macros for random/repeated effect model.

Example

@covstr(factor|subject)
source

Random/repeated effect construction

Metida.VarEffectType
VarEffect(formula, covtype::T, coding) where T <: AbstractCovarianceType

VarEffect(formula, covtype::T; coding = nothing) where T <: AbstractCovarianceType

VarEffect(formula; coding = nothing)

Random/repeated effect.

  • formula from @covstr(ex) macros.

  • covtype - covariance type (SI, DIAG, CS, CSH, AR, ARH, ARMA, TOEP, TOEPH, TOEPP, TOEPHP)

Note

Categorical factors are coded with FullDummyCoding() by default, use coding for other contrast codeing.

Example

VarEffect(@covstr(1+factor|subject), CSH)

VarEffect(@covstr(1 + formulation|subject), CSH; coding = Dict(:formulation => StatsModels.DummyCoding()))
source

Fitting

StatsAPI.fit!Function
fit!(lmm::LMM{T}; kwargs...
) where T

Fit LMM model.

Keywords:

  • solver - :default / :nlopt for using with MetidaNLopt.jl/ :cuda for using with MetidaCu.jl
  • verbose - :auto / 1 / 2 / 3
  • varlinkf - :exp / :sq / :identity ref
  • rholinkf - :sigm / :atan / :sqsigm / :psigm
  • aifirst - first iteration with AI-like method - :default / :ai / :score
  • g_tol - absolute tolerance in the gradient
  • x_tol - absolute tolerance of theta vector
  • f_tol - absolute tolerance in changes of the REML
  • hes - calculate REML Hessian
  • init - initial theta values
  • io - output IO
  • time_limit - time limit = 120 sec
  • iterations - maximum iterations = 300
  • refitinit - true/false - if true - use last values for initial condition
  • optmethod - Optimization method. Look at Optim.jl documentation. (Newton by default)
  • singtol - singular tolerance = 1e-8
source
  • solver - :default solving with Optim.jl, for :nlopt and :cuda MetidaNLopt.jl and MetidaCu.jl should be installed.

  • verbose - 1 - only log, 2 - log and print, 3 - print only errors, other log, 0 (or any other value) - no logging.