API

MRVCModels stands for Multivariate Response Variance Components linear mixed Models. MRVCModels.jl permits maximum likelihood (ML) or residual maximum likelihood (REML) estimation and inference.

MRTVCModel(
    Y::AbstractVecOrMat,
    X::Union{Nothing, AbstractVecOrMat},
    V::Vector{<:AbstractMatrix}
    )

Create a new MRTVCModel instance from response matrix Y, predictor matrix X, and kernel matrices V. The number of variance components must be two.

Keyword arguments

se::Bool        calculate standard errors; default true
reml::Bool      pursue REML estimation instead of ML estimation; default false

MRVCModel(
    Y::AbstractVecOrMat,
    X::Union{Nothing, AbstractVecOrMat},
    V::Union{AbstractMatrix, Vector{<:AbstractMatrix}}
    )

Create a new MRVCModel instance from response matrix Y, predictor matrix X, and kernel matrices V.

Keyword arguments

se::Bool        calculate standard errors; default true
reml::Bool      pursue REML estimation instead of ML estimation; default false

Extended help

When there are two variance components, computation can be reduced by avoiding large matrix inversion per iteration, which is achieved with MRTVCModel instance. MRTVCModels stands for Multivariate Response Two Variance Components linear mixed Models. MRVCModel is more general and is not limited to two variance components case. For MRTVCModel, the number of variance components must be two.

fisher_B!(model::MRTVCModel)

Compute the sampling variance-covariance model.Bcov of regression coefficients model.B, assuming model.Λ and model.Φ are updated.

fisher_B!(model::MRVCModel)

Compute the sampling variance-covariance model.Bcov of regression coefficients model.B, assuming inverse of covariance matrix model.Ω is available at model.storage_nd_nd.

fisher_B_reml!(model::MRTVCModel)

Compute the sampling variance-covariance model.Bcov_reml of regression coefficients model.B_reml, assuming model.Λ and model.Φ are updated.

fisher_Σ!(model::MRTVCModel)

Compute the sampling variance-covariance model.Σcov of variance component estimates model.Σ, assuming model.Λ and model.Φ are updated.

fisher_Σ!(model::MRVCModel)

Compute the sampling variance-covariance model.Σcov of variance component estimates model.Σ, assuming inverse of covariance matrix model.Ω is available at model.storage_nd_nd.

fit!(model::MRVCModel)
fit!(model::MRTVCModel)

Fit a multivariate response variance components model by MM or EM algorithm.

Keyword arguments

maxiter::Int        maximum number of iterations; default 1000
reltol::Real        relative tolerance for convergence; default 1e-6
verbose::Bool       display algorithmic information; default true
init::Symbol        initialization strategy; :default initializes by least squares, while
    :user uses user-supplied values at model.B and model.Σ
algo::Symbol        optimization algorithm; :MM (default) or :EM (for MRVCModel)
log::Bool           record iterate history or not; default false

Extended help

MM algorithm is provably faster than EM algorithm in this setting, so recommend trying MM algorithm first, which is by default, and switching to EM algorithm if there are convergence issues.

h2(model::MRVCModel)

Calculate heritability estimates and their standard errors, assuming that all variance components capture genetic effects except the last term. Also return total heritability from sum of individual contributions and its standard error.

kron_axpy!(A, X, Y)

Overwrite Y with A ⊗ X + Y. Same as Y += kron(A, X), but more memory efficient.

kron_reduction!(A, B, C; sym = false)

Overwrite C with the derivative of tr(A' (X ⊗ B)) wrt X. C[i, j] = dot(Aij, B), where Aij is the (i , j) block of A. sym = true indicates A and B are symmetric.

loglikelihood!(model::MRVCModel)

Overwrite model.storage_nd_nd by inverse of the covariance matrix model.Ω, overwrite model.storage_nd by U' \ vec(model.R), and return the log-likelihood. Assume model.Ω and model.R are already updated according to model.Σ and model.B.

loglikelihood(model::MRTVCModel)

Return the residual log-likelihood, assuming model.Λ, model.Φ, model.logdetΣ2, and model.R̃Φ are updated.

loglikelihood_miss!(model::MRVCModel)

Overwrite model.storage_nd_nd by inverse of the covariance matrix model.Ω, overwrite model.storage_nd by U' \ vec(model.R), overwrite model.storage_n_miss_n_miss_1 by conditional variance, precompute model.storage_n_miss_n_obs_1 for conditional mean, and return the value of surrogate Q-function of log-likelihood. Assume model.Ω and model.R are already updated according to model.Σ and model.B.

loglikelihood_reml(model::MRTVCModel)

Return the log-likelihood, assuming model.Λ, model.Φ, model.logdetΣ2, and model.R̃Φ_reml are updated.

lrt(model1::MRVCModel, model0::MRVCModel)

Perform a variation of the likelihood ratio test for univariate variance components models as in Molenberghs and Verbeke 2007 with model1 and model0 being the full and nested models, respectively.

permute(Y::AbstractVecOrMat)

Return permutation P such that vec(Y)[P] rearranges vec(Y) with missing values spliced after non-missing values. Also return inverse permutation invP such that vec(Y)[P][invP] = vec(Y).

rg(model::MRVCModel)

Calculate genetic/residual correlation estimates and their standard errors.

update_B!(model::MRTVCModel)

Update the regression coefficients model.B, assuming model.Λ and model.Φ are updated.

update_B!(model::MRVCModel)

Update the regression coefficients model.B, assuming inverse of covariance matrix model.Ω is available at model.storage_nd_nd.

update_B_miss!(model::MRVCModel)

Update the regression coefficients model.B, assuming inverse of covariance matrix model.Ω is available at model.storage_nd_nd and model.storage_n_miss_n_obs_1 for conditional mean is precomputed.

update_B_reml!(model::MRTVCModel)

Update the regression coefficients model.B_reml, assuming model.Λ and model.Φ are updated.

update_Γk!(model, k)

Update the parameters model.Γ[k] and model.Ψ[:,k] assuming a diagonal plus low rank structure for model.Σ[k]. Assumes covariance matrix model.Ω is available at model.storage_nd_nd and model.Ω⁻¹R precomputed.

update_Σ!(model::MRVCModel)

Update the variance component parameters model.Σ, assuming inverse of covariance matrix model.Ω is available at model.storage_nd_nd. If missing response, assume conditional variance model.storage_n_miss_n_miss_1 is precomputed.

update_Σk!(model, k, rk)

Update the model.Σ[k] assuming it has rank rk < d, assuming inverse of covariance matrix model.Ω is available at model.storage_nd_nd and model.Ω⁻¹R precomputed.

update_Σk!(model::MRVCModel, k, Val(:EM))

EM update the model.Σ[k] assuming it has full rank d, inverse of covariance matrix model.Ω is available at model.storage_nd_nd, and model.Ω⁻¹R is precomputed.

update_Σk!(model::MRVCModel, k, Val(:MM))

MM update the model.Σ[k] assuming it has full rank d, inverse of covariance matrix model.Ω is available at model.storage_nd_nd, and model.Ω⁻¹R is precomputed.

update_Σk_miss!(model::MRVCModel, k, Val(:MM))

MM update the model.Σ[k] assuming it has full rank d, inverse of covariance matrix model.Ω is available at model.storage_nd_nd, model.Ω⁻¹R is precomputed, and Ω⁻¹P'CPΩ⁻¹ is available at model.storage_nd_nd_miss.

vech(A::AbstractVecOrMat)

Return lower triangular part of A.

◺(n::Int)

Triangular number n * (n + 1) / 2.