Installation
import Pkg; Pkg.add("Metida")
Simple example
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;
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
— TypeLMM(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)
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 specifiedVarEffect(@covstr(1|1), SI)
used. If no repeated effects specified vector of ones used.
Random/repeated model
Metida.@covstr
— Macro@covstr(ex)
Macros for random/repeated effect model.
Example
@covstr(factor|subject)
Random/repeated effect construction
Metida.VarEffect
— TypeVarEffect(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)
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()))
Fitting
StatsBase.fit!
— Functionfit!(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 / :cudaverbose
- :auto / 1 / 2 / 3varlinkf
- not implementedrholinkf
- :sigm / :atan / :sqsigm / :psigmaifirst
- first iteration with AI-like method - :default / :ai / :scoreg_tol
- absolute tolerance in the gradientx_tol
- absolute tolerance of theta vectorf_tol
- absolute tolerance in changes of the REMLhes
- calculate REML Hessianinit
- initial theta valuesio
- uotput IO
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
— TypeCovmatMethod(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
Step 2: make custom CovarianceType
CovarianceType(cm::AbstractCovmatMethod)
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