API

@covstr(ex)

Macros for random/repeated effect model.

Example

@covstr(model|subject)

@lmmformula(formula, args...)

Macro for made formula with variance-covariance structure representation. @lmmformula could be used for shorter LMM construction.

Example:

lmm = Metida.LMM(@lmmformula(var~sequence+period+formulation,
random = formulation|subject:CSH,
repeated = formulation|subject:DIAG),
df0)

equal to:

lmm = LMM(@formula(var~sequence+period+formulation), df0;
random = VarEffect(@covstr(formulation|subject), CSH),
repeated = VarEffect(@covstr(formulation|subject), DIAG),
)

@lmmformula have 3 components - 1'st is a formula for fixed effect, it's defined like in StstsModels (1st argument just provided to @formula macro). Other arguments should be defined like keywords. repeated keyword define repeated effect part, random - define random effect part. You can use several random factors as in example bellow:

lmm = LMM(@lmmformula(var~sequence+period+formulation,
random = formulation|subject:CSH,
random = 1|subject:DIAG,
repeated = formulation|subject:DIAG),
df0)

random or repeated structure made by template:

effect formula | blocking factor [/ nested factor] [: covariance structure]

| - devide effect formula form blocking factor definition (necessarily), / and : modificator are optional.

/ work like in MixedModels or in RegressionFormulae - expand factor f|a/b to f|a + f|a&b. It can't be used in repeated effect definition.

: - covariance structure defined right after : (SI, DIAG, CS, CSH, ets...), if : not used then SI used for this effect.

Terms like a+b or a*b shuould not be used as a blocking factors.

CovarianceType(cm::AbstractCovmatMethod)

Make covariance type with AbstractCovmatMethod.

LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing, wts::Union{Nothing, AbstractVector, AbstractMatrix, AbstractString, Symbol} = 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 or vector

wts: regression weights (residuals).

Weigts can be set as Symbol or String, in this case weights taken from tabular data. If weights is vector then this vector applyed to R-side part of covariance matrix (see Weights details). If weights is matrix then R-side part of covariance matrix multiplied by corresponding part of weight-matrix.

See also: @lmmformula

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)

Example

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

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

Autoregressive()

Autoregressive covariance type.

AR = Autoregressive()

\[\begin{bmatrix} 1 & \rho & \rho^2 & \rho^3 \\ \rho & 1 & \rho & \rho^2 \\ \rho^2 & \rho & 1 & \rho \\ \rho^3 & \rho^2 & \rho & 1 \end{bmatrix}\sigma^2\]

AutoregressiveMovingAverage()

Autoregressive moving average covariance type.

ARMA = AutoregressiveMovingAverage()

\[\begin{bmatrix} 1 & \gamma & \gamma\rho & \gamma\rho^2 \\ \gamma & 1 & \gamma & \gamma\rho \\ \gamma\rho & \gamma & 1 & \gamma \\ \gamma\rho^2 & \gamma\rho & \gamma & 1 \end{bmatrix}\sigma^2\]

CompoundSymmetry()

Compound symmetry covariance type.

CS = CompoundSymmetry()

\[\begin{bmatrix} 1 & \rho & \rho & \rho \\ \rho & 1 & \rho & \rho \\ \rho & \rho & 1 & \rho \\ \rho & \rho & \rho & 1 \end{bmatrix}\sigma^2\]

Diag()

Diagonal covariance type.

DIAG = Diag()

\[\begin{bmatrix} \sigma_a^2 & 0 & 0 \\ 0 & \sigma_b^2 & 0 \\ 0 & 0 & \sigma_c^2 \end{bmatrix}\]

HeterogeneousAutoregressive()

Heterogeneous autoregressive covariance type.

ARH = HeterogeneousAutoregressive()

\[\begin{bmatrix} \sigma_a^2 & \rho\sigma_a\sigma_b & \rho^2\sigma_a\sigma_c & \rho^3\sigma_a\sigma_d \\ \rho\sigma_b\sigma_a & \sigma_b^2 & \rho\sigma_b\sigma_c & \rho^2\sigma_b\sigma_d \\ \rho^2\sigma_c\sigma_a & \rho\sigma_c\sigma_b & \sigma_c^2 & \rho\sigma_c\sigma_d \\ \rho^3\sigma_d\sigma_a & \rho^2\sigma_d\sigma_b & \rho\sigma_d\sigma_c & \sigma_d^2 \end{bmatrix}\]

HeterogeneousCompoundSymmetry()

Heterogeneous compound symmetry covariance type.

CSH = HeterogeneousCompoundSymmetry()

\[\begin{bmatrix} \sigma_a^2 & \rho\sigma_a\sigma_b & \rho\sigma_a\sigma_c & \rho\sigma_a\sigma_d \\ \rho\sigma_b\sigma_a & \sigma_b^2 & \rho\sigma_b\sigma_c & \rho\sigma_b\sigma_d \\ \rho\sigma_c\sigma_a & \rho\sigma_c\sigma_b & \sigma_c^2 & \rho\sigma_c\sigma_d \\ \rho\sigma_d\sigma_a & \rho\sigma_d\sigma_b & \rho\sigma_d\sigma_c & \sigma_d^2 \end{bmatrix}\]

HeterogeneousToeplitz()

Heterogeneous toeplitz covariance type. Only for G matrix.

TOEPH = HeterogeneousToeplitz()

\[\begin{bmatrix} \sigma_a^2 & \rho_1 \sigma_a \sigma_b & \rho_2 \sigma_a \sigma_c & \rho_3 \sigma_a \sigma_d \\ \rho_1 \sigma_b \sigma_a & \sigma_b^2 & \rho_1 \sigma_b \sigma_c & \rho_2 \sigma_b \sigma_d \\ \rho_2 \sigma_c \sigma_a & \rho_1 \sigma_c \sigma_b & \sigma_c^2 & \rho_1 \sigma_c \sigma_d \\ \rho_3 \sigma_d \sigma_a & \rho_2 \sigma_d \sigma_b & \rho_1 \sigma_d \sigma_c & \sigma_d^2 \end{bmatrix}\]

HeterogeneousToeplitzParameterized(p::Int)

Heterogeneous toeplitz covariance type with parameter p, (number of bands = p - 1, if p = 1 it's equal DIAG structure).

TOEPHP(p) = HeterogeneousToeplitzParameterized(p)

ScaledIdentity()

Scaled identity covariance type.

SI = ScaledIdentity()

\[\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\sigma^{2}\]

SpatialExponential()

Spatian Exponential covariance structure. Used only for repeated effect.

\[R_{i,j} = \sigma^{2} * exp(-dist(i,j)/\theta)\]

where dist - Euclidean distance between row-vectors of repeated effect matrix for subject i and j, θ > 0.

SPEXP = SpatialExponential()

SpatialGaussian()

Spatian Gaussian covariance structure. Used only for repeated effect.

\[R_{i,j} = \sigma^{2} * exp(- dist(i,j)^2 / \theta^2)\]

where dist - Euclidean distance between row-vectors of repeated effect matrix for subject i and j, θ ≠ 0.

SPGAU = SpatialGaussian()

SpatialPower()

Spatian Power covariance structure. Used only for repeated effect.

\[R_{i,j} = \sigma^{2} * \rho^{dist(i,j)}\]

where dist - Euclidean distance between row-vectors of repeated effect matrix for subject i and j, 1 > ρ > -1.

SPPOW = SpatialPower()

Toeplitz()

Toeplitz covariance type. Only for G matrix.

TOEP = Toeplitz()

\[\begin{bmatrix} 1 & \rho_1 & \rho_2 & \rho_3 \\ \rho_1 & 1 & \rho_1 & \rho_2 \\ \rho_2 & \rho_1 & 1 & \rho_1 \\ \rho_3 & \rho_2 & \rho_1 & 1 \end{bmatrix}\sigma^2\]

ToeplitzParameterized(p::Int)

Toeplitz covariance type with parameter p, (number of bands = p - 1, if p = 1 it's equal SI structure).

TOEPP(p) = ToeplitzParameterized(p)

Unstructured()

Unstructured covariance structure with t*(t+1)/2-t paremeters where t - number of factor levels, t*(t+1)/2-2t of them is covariance (ρ) patemeters. All levels for repeated effect should be unique within each subject.

UN = Unstructured()

caic(lmm::LMM)

Conditional Akaike Information Criterion.

coefn(lmm)

Coef number.

dof_satter(lmm::LMM{T}, l) where T

Return Satterthwaite approximation for the denominator degrees of freedom, where l is a contrast vector (estimable linear combination of fixed effect coefficients vector (β).

\[df = \frac{2(LCL')^{2}}{g'Ag}\]

Where: $A = 2H^{-1}$, $g = \triangledown_{\theta}(LC^{-1}_{\theta}L')$

dof_satter(lmm::LMM{T}, i::Int) where T

Return Satterthwaite approximation for the denominator degrees of freedom, where n - coefficient number.

dof_satter(lmm::LMM{T}) where T

Return Satterthwaite approximation for the denominator degrees of freedom for all coefficients.

dof_satter(lmm::LMM{T}, l::Matrix) where T

Return Satterthwaite approximation for the denominator degrees of freedom for conrast matrix l.

For size(l, 1) > 1:

\[df = \frac{2E}{E - rank(LCL')}\]

where:

• let $LCL' = QΛQ^{-1}$, where $QΛQ^{-1}$ - spectral decomposition of $LCL'$
• $Lq_i$ is the i-th row of $Q^{-1}L$
• $A = 2H^{-1}$, $g = \triangledown_{\theta}(Lq_i C^{-1}_{\theta} Lq_i')$
• $v_i = \frac{2*Λ_{i,i}^2}{g' * A * g}$
• $E = \sum_{i=1}^n {\frac{v_i}(v_i - 2)}$ for $v_i > 2$

estimate(lmm, l::AbstractVector; level = 0.95, name = "Estimate")

Estimate table for l vector. Satter DF used.

estimate(lmm; level = 0.95)

Estimates table. Satter DF used.

getlog(lmm::LMM)

Return fitting log.

gmatrix(lmm::LMM{T}, r::Int) where T

hessian(lmm, theta)

Calculate Hessian matrix of REML for theta.

lcontrast(lmm::LMM, i::Int)

L-contrast matrix for i fixed effect.

Number of blocks

rand(rng::AbstractRNG, lmm::LMM{T}) where T

Generate random responce vector for fitted 'lmm' model.

rand(rng::AbstractRNG, lmm::LMM{T}; theta) where T

Generate random responce vector 'lmm' model, theta covariance vector, and zero means.

rand(rng::AbstractRNG, lmm::LMM{T}; theta, beta) where T

Generate random responce vector 'lmm' model, theta covariance vector and mean's vector.

rand!(v::AbstractVector, lmm::LMM) = rand!(default_rng(), v, lmm, lmm.result.theta, lmm.result.beta)

Generate random responce vector for fitted 'lmm' model, store results in v.

rand!(rng::AbstractRNG, lmm::LMM{T}; theta) where T

Generate random responce vector 'lmm' model, theta covariance vector, and zero means, store results in v.

raneff(lmm::LMM{T}, i)

Vector of random effect coefficients for block i.

raneff(lmm::LMM{T})

Vector of random effect coefficients for all subjects by each random effect.

raneffn(lmm)

Retuen number of random effects.

rankx(lmm::LMM)

Return rank of X matrix.

rmatrix(lmm::LMM{T}, i::Int) where T

theta(lmm::LMM)

Return theta vector.

thetalength(lmm::LMM)

Length of theta vector.

vmatrix(lmm::LMM, i::Int)

Return variance-covariance matrix V for i bolock.

vmatrix!(V, θ, lmm, i)

Update variance-covariance matrix V for i bolock. Upper triangular updated.

m2logreml(lmm::LMM, θ = theta(lmm); maxthreads::Int = num_cores())

-2 logREML

logreml(lmm::LMM, θ = theta(lmm); maxthreads::Int = num_cores())

logREML

StatsBase.aic(lmm::LMM)

Akaike Information Criterion.

StatsBase.aicc(lmm::LMM)

Corrected Akaike Information Criterion.

StatsBase.bic(lmm::LMM)

Bayesian information criterion.

StatsBase.coef(lmm::LMM) = copy(lmm.result.beta)

Model coefficients (β).

StatsBase.coefnames(lmm::LMM) = StatsBase.coefnames(lmm.mf)

Coefficients names.

coeftable(lmm::LMM)

Return coefficients table.

StatsBase.confint(lmm::LMM{T}; level::Real=0.95, ddf::Symbol = :satter) where T

Confidece interval for coefficients.

ddf = :satter/:residual

\[CI_{U/L} = β ± SE * t_{ddf, 1-α/2}\]

See also: dof_satter, dof_residual

StatsBase.confint(lmm::LMM{T}, i::Int; level::Real=0.95, ddf::Symbol = :satter) where T

Confidece interval for coefficient i.

StatsBase.confint(br::BootstrapResult, n::Int; level::Float64=0.95, method=:bp, metric = :coef, delrml = false)

Confidence interval for bootstrap result.

*method:

• :bp - bootstrap percentile;
• :rbp - reverse bootstrap percentile;
• :norm - Normal distribution;
• :bcnorm - Bias corrected Normal distribution;
• :jn - bias corrected (jackknife resampling).

crossmodelmatrix(lmm::LMM)

Return X'X.

StatsBase.dof(lmm::LMM)

DOF.

StatsBase.dof_residual(lmm::LMM)

DOF residuals: N - rank(X), where N - total number of observations.

fit(::Type{T}, f::FormulaTerm, data;
contrasts=Dict{Symbol,Any}(),  
random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, 
repeated::Union{Nothing, VarEffect} = nothing,
kwargs...)

Fit LMM model with @formula.

Keywords see fit!

fit(::Type{T}, f::LMMformula, data;
contrasts=Dict{Symbol,Any}(),  
kwargs...) where T <: LMM

Fit LMM model with @lmmformula.

Keywords see fit!

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 - - 1 - only log, 2 - log and print, 3 - print only errors, other log, 0 (or any other value) - no logging
• varlinkf - :exp / :sq / :identity ref
• rholinkf - :sigm / :atan / :sqsigm / :psigm
• aifirst - first iteration with AI-like method - :default / :ai / :score
• aifmax - maximum pre-optimization steps
• 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 (false by default)
• optmethod - Optimization method. Look at Optim.jl documentation. (Newton by default)
• singtol - singular tolerance = 1e-8
• maxthreads - maximum threads = min(num_cores(), Threads.nthreads())

StatsBase.islinear(model::LMM)

StatsBase.isfitted(lmm::LMM)

StatsBase.loglikelihood(lmm::LMM)

Return loglikelihood value.

StatsBase.modelmatrix(lmm::LMM)

Fixed effects matrix.

StatsBase.nobs(lmm::MetiaModel)

Number of observations.

StatsBase.response(lmm::LMM)

Response vector.

responsename(lmm::LMM)

Responce varible name.

StatsBase.stderror(lmm::LMM)

Standard error

StatsBase.vcov(lmm::LMM)

Variance-covariance matrix of β.

SpatialExponentialD()

Same as SpatialExponential, but add D to all diagonal elements.

SPEXPD = SpatialExponentialD()

SpatialGaussianD()

Same as SpatialGaussianD, but add D to all diagonal elements.

SPGAUD = SpatialGaussianD()

SpatialPowerD()

Same as SpatialPower, but add D to all diagonal elements.

SPPOWD = SpatialPowerD()

ScaledWeightedCov(wtsm::AbstractMatrix{T})

Scaled weighted covariance matrix, where wtsm - NxN within block correlation matrix (N - total number of observations). Used only for repeated effect.

SWC = ScaledWeightedCov

\[R = Corr(W) * \sigma_c^2\]

where $Corr(W)$ - diagonal correlation matrix.

example:

matwts = Symmetric(UnitUpperTriangular(rand(size(df0,1), size(df0,1))))
lmm = LMM(@formula(var~sequence+period+formulation), df0;
    repeated = VarEffect(@covstr(1|subject), SWC(matwts)))
fit!(lmm)

There is no wtsm checks for symmetricity or values.

ACOV(c, action = 0)

Augmented (adjusted) covariance. Add additional correlations to existed R-part of variance covariance matrix. Can be used with AR or CS types.

action if existed covariance not equal sero:

• 0 - add
• 1 - replace
• 2 - do nothing (use existed value)
• 3 - warning and add
• 4 - warning and replace
• 5 - warning and use existed value
• other - error

dof_contain(lmm, i)

typeiii(lmm::LMM{T}; ddf::Symbol = :satter) where T

Type III table.

MILMM(lmm::LMM, data)

Multiple imputation model.

mutable struct RawCoding <: AbstractContrasts

Contrast for CategoricalTerm to get column "as it is" for model matrix.

contrast(lmm, l::AbstractMatrix; name::String = "Contrast", ddf = :satter)

User contrast table. ddf = :satter or :residual or any number for direct ddf setting.

fvalue(lmm::LMM, l::AbstractMatrix)

F value for contrast matrix l.

\[F = \frac{\beta'L'(LCL')^{-1}L\beta}{rank(LCL')}\]

mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

θ + A * B * A'

Change θ (only upper triangle). B is symmetric.

mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, alpha)

θ + A * B * A' * alpha

Change θ (only upper triangle). B is symmetric.

mulαβαtinc!(θ::AbstractVector{T}, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha) where T

θ + A * B * (a - b) * alpha

Change θ (only upper triangle). B is symmetric.

mulθ₃(y, X, β, V::AbstractMatrix{T})::T where T

(y - X * β)' * (-V) * (y - X * β)

use only upper triangle of V

mulαtβinc!(θ::AbstractVector{T}, A::AbstractMatrix, b::AbstractVector) where T

θ + A' * b

Change θ.

Term name.

Return number of subject foe each random effet in current block.

edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T

Distance between vector mx[i, :] and mx[j, :].

m2logml(lmm::LMM, β = coef(lmm), θ = theta(lmm); maxthreads::Int = num_cores())

-2 logML

logml(lmm::LMM, beta = coef(lmm), θ = theta(lmm); maxthreads::Int = num_cores())

logML