First step · Metida.jl

Installation

import Pkg; Pkg.add("Metida")

Simple example

LMM
Metida.@covstr
Metida.VarEffect
fit!

Step 1: Load data

Load provided data with CSV and DataFrames:

using Metida, CSV, DataFrames

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

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

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    -2 logREML: 10.0652    BIC: 23.9282

    Fixed-effects parameters:
Name             Estimate     SE          z           Pr(>|z|)     
(Intercept)      1.57749      0.334543    4.71536     2.41283e-6   
sequence: 2      -0.170833    0.384381    -0.444437   0.656726     
period: 2        0.195984     0.117228    1.67182     0.0945606    
period: 3        0.145014     0.109171    1.32832     0.184073     
period: 4        0.157363     0.117228    1.34236     0.179478     
formulation: 2   -0.0791667   0.0903709   -0.876019   0.38102      

    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

5-element Array{Metida.LMMLogMsg,1}:
   INFO  : Initial θ: [0.2108585059318441, 0.2108585059318441, 0.0001, 0.2108585059318441, 0.2108585059318441]

   ERROR : DomainError (-9.735556609752802e32, log will only return a complex result if called with a complex argument. Try log(Complex(x)).) during REML calculation.

   INFO  : Resulting θ: [0.45558429062258543, 0.3676561978014236, 0.9999999997954805, 0.14368242865088818, 0.20565739028506555]; 22 iterations.

   INFO  : Model fitted.

   WARN  : det(G) of random effect 1 is less 1e-08.

Model construction

Metida.LMM — Type

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)

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. 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 effect construction

Metida.VarEffect — Type

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

StatsBase.fit! — Function

fit!(lmm::LMM{T};
solver::Symbol = :default,
verbose = :auto,
varlinkf = :exp,
rholinkf = :sigm,
aifirst = :default,
g_tol::Float64 = 1e-12,
x_tol::Float64 = 1e-12,
f_tol::Float64 = 1e-12,
hes::Bool   = true,
init = nothing,
io::IO = stdout,
) where T

Fit LMM model.

solver - :default / :nlopt / :cuda
verbose - :auto / 1 / 2 / 3
varlinkf - not implemented
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 - uotput IO

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.

Custom structure

Step 1: make custom CovarianceMethod for R and G matrix

Metida.CovmatMethod — Type

CovmatMethod(nparamf::Function, xmat!::Function)

nparamf - function type (t, q) -> (a, b)

where: t - size(z, 2) - number of levels for effect (number of columns of individual z matriz); q - number of factors in the effect model; a - number of variance parameters; b - number of ρ parameters;

Example: (t, q) -> (t, 1) for CSH structure; (t, q) -> (1, 1) for AR, ets.

Tuple{Int, Int} should be returned.

xmat! - construction function

G matrix function should update mx, where mx is zero matrix, p - parameter of CovarianceType structure.

Exapple

function gmat_diag!(mx, θ::Vector{T}, p) where T
    for i = 1:size(mx, 1)
            mx[i, i] = θ[i] ^ 2
        end
    nothing
end

Exapple

R matrix function should add R part to mx, rz - is repeated effect matrix.

Example

function rmatp_diag!(mx, θ::Vector{T}, rz, p) where T
    for i = 1:size(mx, 1)
        for c = 1:length(θ)
            mx[i, i] += rz[i, c] * θ[c] * θ[c]
        end
    end
    nothing
end

source

Step 2: make custom CovarianceType

  CovarianceType(cm::AbstractCovmatMethod)

See: Metida.CovarianceType

Step 3 Fit your model

#Make methods for G and R matrix and CovarianceType struct
CCTG = CovarianceType(CovmatMethod((q,p) -> (q, 1), Metida.gmat_csh!))
CCTR = CovarianceType(CovmatMethod((q,p) -> (q, 0), Metida.rmatp_diag!))

#Make model
lmm = LMM(@formula(var~sequence+period+formulation), df;
random = VarEffect(@covstr(formulation|subject), CCTG),
repeated = VarEffect(@covstr(formulation|subject), CCTR),
)

#Fit model
fit!(lmm)

Linear Mixed Model: var ~ sequence + period + formulation
Random 1: 
    Model: formulation|subject
    Type: FUNC (3), Subjects: 5
Repeated: 
    Model: formulation|subject
    Type: FUNC (2)
    Blocks: 5, Maximum block size: 4

Status: converged    -2 logREML: 10.0652    BIC: 23.9282

    Fixed-effects parameters:
Name             Estimate     SE          z           Pr(>|z|)     
(Intercept)      1.57749      0.334543    4.71536     2.41283e-6   
sequence: 2      -0.170833    0.384381    -0.444437   0.656726     
period: 2        0.195984     0.117228    1.67182     0.0945606    
period: 3        0.145014     0.109171    1.32832     0.184073     
period: 4        0.157363     0.117228    1.34236     0.179478     
formulation: 2   -0.0791667   0.0903709   -0.876019   0.38102      

    Variance components:
    θ vector: [0.455584, 0.367656, 1.0, 0.143682, 0.205657]
Random 1   NA   var   0.207557    
Random 1   NA   var   0.135171    
Random 1   NA   rho   1.0         
Residual   NA   var   0.0206445   
Residual   NA   var   0.0422948