\[ \newcommand{\E}{\text{E}} \newcommand{\mbf}{\mathbf} \newcommand{\bs}{\boldsymbol} \newcommand{\Cov}{\text{Cov}} \newcommand{\Var}{\text{Var}} \newcommand{\Cov}{\text{Cov}} \newcommand{\A}[1]{\mathbf{A}_{#1}} \newcommand{\Ai}[1]{\mathbf{A}^{#1}} \]

18. Extension to Multi-trait Models#

18.1. Full-conditional posteriors for \(\mathbf{G}_0\) and \(\mathbf{R}_0\)#

Suppose the covariance matrix of \(\mbf{u}\) can be written as:

\[\begin{split} \begin{align} \Var(\mbf{u}) &=\mbf{G}_0 \otimes \mbf{A} \\ &= \begin{bmatrix} g_{011}\mbf{A} & \cdots & g_{01m}\mbf{A} \\ \vdots & \ddots & \vdots \\ g_{0m}\mbf{A} & \cdots & g_{0mm}\mbf{A} \end{bmatrix}, \end{align} \end{split}\]

where \(\mbf{G}_0\) and \(\mbf{A}\) are positive definite matricies, and \(g_{0ij}\) denotes the element from row \(i\) and column \(j\) of \(\mbf{G}_0\).

Then, the mutivariate normal prior for \(\mbf{u}\) is:

\[\begin{split} \begin{align} f(\mbf{u}|\mbf{G}_0) &\propto |\mbf{G}_0|^{-k/2} \exp\left\{-\frac{1}{2}\mbf{u}'\mbf{G}^{-1}\mbf{u} \right\} \\ &\propto |\mbf{G}_0|^{-k/2} \exp\left\{-\frac{1}{2}\mbf{u}'(\mbf{G}_0\otimes\mbf{A})^{-1}\mbf{u}\right\} \\ &\propto |\mbf{G}_0|^{-k/2} \exp\left\{-\frac{1}{2}\mbf{u}'(\mbf{G}_0^{-1}\otimes\mbf{A}^{-1})\mbf{u}\right\} \\ &\propto |\mbf{G}_0|^{-k/2} \exp\left\{-\frac{1}{2} \begin{bmatrix} \mbf{u}_1 & \mbf{u}_2 &\cdots &\mbf{u}_m \end{bmatrix} \begin{bmatrix} g_0^{11}\mbf{A}^{-1} & \cdots & g_0^{1m}\mbf{A}^{-1} \\ \vdots & \ddots & \vdots \\ g_0^{1m}\mbf{A}^{-1} & \cdots & g_0^{mm}\mbf{A}^{-1} \end{bmatrix} \begin{bmatrix} \mbf{u}_1\\ \mbf{u}_2\\ \vdots\\ \mbf{u}_m \end{bmatrix} \right\}\\ &\propto |\mbf{G}_0|^{-k/2} \exp\left\{ -\frac{1}{2}\sum_{ij} g_0^{ij}\mbf{u}_i'\mbf{A}^{-1}\mbf{u}_j \right\}\\ &\propto |\mbf{G}_0|^{-k/2} \exp\left\{ -\frac{1}{2}\sum_{ij} g_0^{ij}s_{ij} \right\}\\ &\propto |\mbf{G}_0|^{-k/2} \exp\left\{ -\frac{1}{2}\text{tr}(\mbf{G}_0^{-1}\mbf{S}) \right\}, \tag{1} \end{align} \end{split}\]

where the \(\mbf{S} = \{s_{ij} = \mbf{u}_i'\mbf{A}^{-1}\mbf{u}_j\}\). In this prior, the matrix \(\mbf{A}\) is assumed to be known, and the matrix \(\mbf{G}_0\) may be considered unknown with a scaled inverted Wishart prior:

\[ f(\mbf{G}_0|\nu_u,\mbf{V_u}) \propto |\mbf{G}_0|^{-(\nu_u+m+2)/2} \exp\left\{ -\frac{1}{2}\text{tr}(\mbf{G}_0^{-1}\mbf{V}_u) \right\}, \]

where \(\nu_u\) and \(\mbf{V}_u\) are the degrees of freedom and scale matrix of this prior. The scale matrix \(\mbf{V}_u\) in this prior is the inverse of the scale matrix in Sorensen and Gianola (2002). The parameterization presented here is that used in the Julia Distributions package. Combining this prior with (1) gives full-conditional posterior for \(\mbf{G}_0\):

\[\begin{split} \begin{align} f(\mbf{G}_0|\text{ELSE}) &\propto |\mbf{G}_0|^{-k/2} \exp\left\{ -\frac{1}{2}\text{tr}(\mbf{G}_0^{-1}\mbf{S}) \right\} \\ &\times |\mbf{G}_0|^{-(\nu_u+m+2)/2} \exp\left\{ -\frac{1}{2}\text{tr}(\mbf{G}_0^{-1}\mbf{V}_u) \right\}\\ &\propto |\mbf{G}_0|^{-(\nu_u+m + k +2)/2} \exp\left\{ -\frac{1}{2}\text{tr}[\mbf{G}_0^{-1}(\mbf{V}_u + \mbf{S})] \right\}, \end{align} \end{split}\]

which can be recognized as an inverted Wishart distribution with \(\tilde{\nu}_u = (\nu_u+k)\) degrees of freemom and scale matrix \(\tilde{\mbf{V}}_u = (\mbf{V}_u + \mbf{S})\).

Similarly, it can be shown that the full-conditional posterior for \(\mbf{R}_0\) is a scaled inverted Wishart distribution with \(\tilde{\nu}_e = (\nu_e + n)\) degrees of freedom and scale matrix \(\tilde{\mbf{V}}_e = (\mbf{V}_e + \mbf{S}_e)\), when the residual covariance matrix has the form:

\[ \begin{align} \Var(\mbf{e}) &= \mbf{R}_0 \otimes \mbf{I}, \end{align} \]

where \(\mbf{R}_0\) is an \(m \times m\) positive definite matrix, \(\mbf{I}\) is an \(n\times n\) identity matrix, \(\mbf{S}_e = \{\mbf{e}_i'\mbf{e}_j \}\), and \(\nu_e\) and \(\mbf{V}_e\) are the hyperparameters of the scaled inverted Wishart prior for \(\mbf{R}_0\).

Consider the mixed linear model

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Zu} + \mathbf{e}, \]

where \(\boldsymbol{\beta}\) is a vector of fixed effects, \(\mathbf{X}\) is an observed matrix that relates \(\boldsymbol{\beta}\) to \(\mathbf{y}\), \(\mathbf{u}\) is a vector of multivariate normal random effects with null means, and it is assumed that its covariance matrix can be written as

\[\begin{split} \begin{align} \Var(\mbf{u}) &=\mbf{G}_0 \otimes \mbf{A} \\ &= \begin{bmatrix} g_{011}\mbf{A} & \cdots & g_{01m}\mbf{A} \\ \vdots & \ddots & \vdots \\ g_{0m}\mbf{A} & \cdots & g_{0mm}\mbf{A} \end{bmatrix}, \end{align} \end{split}\]

where \(\mbf{G}_0\) and \(\mbf{A}\) are positive definite matricies, and \(g_{0ij}\) denotes the element from row \(i\) and column \(j\) of \(\mbf{G}_0\). The matrix \(\mbf{A}\) is assumed to be known, and the matrix \(\mbf{G}_0\) may be considered unknown with a scaled inverted Wishart prior:

\[ f(\mbf{G}_0|\nu_u,\mbf{V_u}) \propto |\mbf{G}_0|^{-(\nu_u+m+2)/2} \exp\left\{ -\frac{1}{2}\text{tr}(\mbf{G}_0^{-1}\mbf{V}_u) \right\}, \]

where \(\nu_u\) and \(\mbf{V}_u\) are the degrees of freedom and scale matrix of this prior. The scale matrix \(\mbf{V}_u\) in this prior is the inverse of the scale matrix in Sorensen and Gianola (2002). The parameterization presented here is that used in the Julia Distributions package.

18.2. Full-conditional Posteriors for Location Parameters#

If we denote the location parameters (\(\boldsymbol{\beta}, \mbf{u}\)) of the model by \(\boldsymbol{\theta}\), as we have seen before, the posteior distribution of \(\boldsymbol{\theta}\) is multivariate normal with mean given by the solution to

\[ \mbf{C}\hat{\bs{\theta}} = \mbf{r} \]

and the variance is \(\mbf{C}^{-1}\sigma^2_e\), where \(\mbf{C}\) and \(\mbf{r}\) are the coefficient matrix and right-hand-side of the MME. So, the quadratic in the exponent of this normal distribution is:

\[ Q = \frac{-1}{2\sigma^2_e}(\bs{\theta} - \hat{\boldsymbol{\theta}})'\mathbf{C}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}). \tag{2} \]

In order to use the single-site Gibbs sampler to construct a Markov chain for \(\bs{\theta}\), we need to derive the full-conditional posterior distribution for \((\theta_i|\text{ELSE})\), where \(\text{ELSE}\) stands for all other variables in the model other than \(\theta_i\). To simplify this derivation, \(\bs{\theta}\) is partitioned as:

\[ \bs{\theta} = [\bs{\theta}_1, \bs{\theta}_2], \]

where \(\bs{\theta}_1 = \theta_i\) and \(\bs{\theta}_2\) contains the remaining elements of \(\bs{\theta}\). Now, the quadratic (2) can be written as:

\[\begin{split} Q = \frac{-1}{2\sigma^2_e} \begin{bmatrix} (\bs{\theta}_1 - \hat{\bs{\theta}}_1) & (\bs{\theta}_2 - \hat{\bs{\theta}}_2) \end{bmatrix} \begin{bmatrix} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \end{bmatrix} \begin{bmatrix} (\bs{\theta}_1 - \hat{\bs{\theta}}_1) \\ (\bs{\theta}_2 - \hat{\bs{\theta}}_2) \end{bmatrix}, \tag{3} \end{split}\]

where the coefficient matrix \(\mbf{C}\) has also been partitioned to correspond to the partitioning of \(\bs{\theta}\).

To derive the distribution of \((\bs{\theta}_1|\bs{\theta}_2)\), we need to only consider the terms in (3) that involve \(\bs{\theta}_1\). These terms are (after dropping \(\frac{-1}{2\sigma^2_e}\)):

\[\begin{split} \begin{align} Q_1 &= \bs{\theta}_1'\mbf{C}_{11}\bs{\theta}_1 - 2\bs{\theta}_1'\mbf{C}_{11}\hat{\bs{\theta}}_1 + 2\bs{\theta}_1'\mbf{C}_{12}\bs{\theta}_2 - 2\bs{\theta}_1'\mbf{C}_{12}\hat{\bs{\theta}}_2 \tag{4} \\ &= (\bs{\theta}_1 - \tilde{\bs{\theta}}_1)'\mbf{C}_{11}(\bs{\theta}_1 - \tilde{\bs{\theta}}_1) - \tilde{\bs{\theta}}_1'\mbf{C}_{11}\tilde{\bs{\theta}}_1, \tag{5} \end{align} \end{split}\]

where for (4) to equal (5), \(\tilde{\bs{\theta}}_1\) must satisfy the equality:

\[ \bs{\theta}_1'\mbf{C}_{11}\tilde{\bs{\theta}}_1 = \bs{\theta}_1'\mbf{C}_{11}\hat{\bs{\theta}}_1 + \bs{\theta}_1'\mbf{C}_{12}\hat{\bs{\theta}}_2 - \bs{\theta}_1'\mbf{C}_{12}\bs{\theta}_2. \]

This equlity must hold for all \(\bs{\theta}_1\), and so

\[\begin{split} \begin{align} \mbf{C}_{11}\tilde{\bs{\theta}}_1 &= \mbf{C}_{11}\hat{\bs{\theta}}_1 + \mbf{C}_{12}\hat{\bs{\theta}}_2 - \mbf{C}_{12}\bs{\theta}_2 \\ \mbf{C}_{11}\tilde{\bs{\theta}}_1 &= \mbf{r}_1 - \mbf{C}_{12}\bs{\theta}_2, \label{eq:singleSiteGibbs} \end{align} \end{split}\]

because the right-hand-side of the MME corresponding to \(\hat{\bs{\theta}}_1\) is equal to

\[ \mbf{C}_{11}\hat{\bs{\theta}}_1 + \mbf{C}_{12}\hat{\bs{\theta}}_2 = \mbf{r}_1. \]

Thus, \((\bs{\theta}_1|\bs{\theta}_2) \sim \text{N}(\tilde{\bs{\theta}},\mbf{C}_{11}^{-1}\sigma^2_e)\). An alternative derivation is given in Sorensen and Gianola (2002).

18.3. Gibbs Sampler for Location Parameters#

function Gibbs!(A,x,b,varRes::Float64,nIter::Int64;outFreq=100)
    n = size(x,1)
    xMean = zeros(n)
    for iter = 1:nIter
        if iter%outFreq==0
            println("at sample: ",iter)
        end
        for i=1:n
            cVarInv = 1.0/A[i,i]
            cMean   = cVarInv*(b[i] - A[:,i]'x) + x[i]
            x[i]    = randn()*sqrt(cVarInv) + cMean
        end
        xMean += (x - xMean)/iter
    end
    return xMean
end

Gibbs! (generic function with 1 method)

using NBInclude
@nbinclude("../MME/3.7.3.BuildMME.ipynb"; regex=r"#\s*EXECUTE")

updateLhsRhs! (generic function with 1 method)

18.4. Structures for Building MME#

#EXECUTE 
using DataFrames, SparseArrays, LinearAlgebra, Random, JWAS, CSV
using Printf

mutable struct ModelTerm
    trmString::AbstractString
    iModel::Int64                           # the model (trait) this term belongs to
    startPos::Int64
    endPos::Int64
    randomType::String
end

mutable struct RandomEffect
    modelTermVec::Array{ModelTerm,1}      
    νpr::Float64
    Spr::Array{Float64,2}
    Vi::Array{Float64,2}
    meanV::Array{Float64,2}
    Ai
    estimate::Bool
    outSamples::Bool
    outFreq::Int64
    outStream
end

mutable struct MME
    modelEquations::Array{AbstractString,1} #"[y1 = A + B + A*B ;y2 = A + B"]
    covVec::Array{String}                   #["age"]
    modelTermVec::Array{ModelTerm,1}        #[modelTerm("A") , modelTerm("A*B")]
    modelTermDict::Dict{AbstractString,ModelTerm}
    randomEffectsVec::Array{RandomEffect,1} #[RandomEffect("Ind"), RandomEffect("Mat")]
    ped                                     # false or a pedigee structure defined in JWAS.PedModule
    depVarVec::Array{AbstractString,1}      #["y1","y2"]
    X::SparseMatrixCSC{Float64,Int64}       #design matrix
    mmeLhs::SparseMatrixCSC{Float64,Int64}  #left-hand side for MME
    mmeRhs::Array{Float64,1}                #right-hand side for MME
    mmeSpl::Array{Float64,1}                #sample of location effects
    varRes::Array{Float64,2}                #residual variance 
    νRes::Float64                           #prior degrees of freedom for varRes
    SRes::Array{Float64,2}                  #prior scale parameter for resVar
    y::Array{Float64,1}                     #matrix of values for dependent variable
    meanVarRes::Array{Float64,2}                     
end 

function ModelTerm(str::AbstractString,iModel)
    ModelTerm(str,iModel,0,0,"fixed")
end  

ModelTerm

18.5. Functions for Building MME#

#EXECUTE 
function initMME(modelEquations::AbstractString,varRes;νRes=5.0)  # "y1 = A + A*B; y2 = A + B"
    if modelEquations==""
        error("modelEquations is empty\n")
    end
    modelVec   = [strip(i) for i in split(modelEquations,[';','\n'],keepempty=false)]
    nModels    = size(modelVec,1)
    depVarVec  = Array{AbstractString,1}()
    modelTermVec  = Array{ModelTerm,1}()                 # modelTerms across all the models
    modelTermDict = Dict{AbstractString,ModelTerm}()  
    for (m,model) = enumerate(modelVec)
        leftRight  = split(model,"=")                 # "y", "A+A*B"
        depVarVec  = [depVarVec; strip(leftRight[1])] # ["y1","y2"]
        modelParts = strip(leftRight[2])              # "A+A*B" 
        termsVec   = split(modelParts,"+")            # "A","A*B"
        modelTermVeci = [ModelTerm(string(m)*":"*strip(trmStr),m) for trmStr in termsVec]
        modelTermVec  = [modelTermVec;modelTermVeci]    
    end
    for i in modelTermVec
        modelTermDict[i.trmString] = i
    end
    randomEffectsVec = []
    ped = false
    covVec = []
    X = spzeros(0,0)
    mmeLhs = spzeros(0,0)
    mmeRhs = []
    mmeSol = []
    y      = Array{Float64,1}(undef,0)
    m = size(varRes,1)
    varRes = m==1 ? fill(varRes,1,1) : varRes
    SRes   = varRes*(νRes - m - 1)
    meanVarRes = zero(varRes)
    return MME(modelVec,covVec,modelTermVec,modelTermDict,randomEffectsVec,ped,
               depVarVec,X,mmeLhs,mmeRhs,mmeSol,varRes,νRes,SRes,y,meanVarRes)
end 

initMME (generic function with 1 method)

varRes = 1.0
mme = initMME("y = sex + breed + age",varRes)
mme.covVec = ["age"];

mme.modelTermVec[1].trmString

"1:sex"

mme.modelTermVec[1].trmString[3:end]

"sex"

#EXECUTE 
# This function returns returns a dictionary with the names in the pedigree as the keys and their 
# sequential numbers as the associated values
# It also returns vector with the keys in sequential order. 

function mkDict(ped::JWAS.PedModule.Pedigree)
    d = Dict()
    names = Array{String}(undef,length(ped.idMap))
    for i in ped.idMap    
        d[i.first] = i.second.seqID
        names[i.second.seqID] = i.first
    end
    return d,names
end

# This function returns 
# a dictionary with the unique values in the vector "a" as the keys and their 
# sequential numbers as the associated values
# It also returns vector with the keys in sequential order. 
function mkDict(a)
    d = Dict()
    aUnique = unique(a)
    names = Array{String}(undef,size(aUnique,1))
    for (i,s) in enumerate(aUnique)
    names[i] = s
    d[s] = i
    end
    return d,names
end

function getX(mme,modelTerm::ModelTerm,covariables,df)
    n = size(df,1)
    nModels = size(mme.modelEquations,1)
    trmString = split(modelTerm.trmString,':')[2]
    if trmString == "intercept"
        ii = (modelTerm.iModel-1)*n .+ (1:n) # row numbers
        jj = ones(n)
        val = 1.0
        p = 1
        X  = sparse(ii,jj,val,n*nModels,p)
        colNames = modelTerm.trmString
        return X,colNames
    end
    factors = strip.(split(trmString,"*"))
    covs = [i in covariables for i in factors]
    
    if covs[1] == false
        str = string.(df[:,Symbol(factors[1])])
        val = 1.0
    else
        str = fill(factors[1],n) 
        val = df[:,Symbol(factors[1])]    
    end       

    for i in 2:length(factors)
        if covs[i] == false
            str = str .*" x ".*string.(df[:,Symbol(factors[i])])
            val = val .* 1.0 
        else
            str = str .*" x ".*fill(factors[i],n) 
            val = val .* df[:,Symbol(factors[i])]    
        end 
    end
    dict,colNames   = modelTerm.randomType != "A" ? mkDict(str) : mkDict(mme.ped)
    str = modelTerm.randomType != "A" ? str : string.(df[:,Symbol(factors[1])])        
    ii = (modelTerm.iModel-1)*n .+ (1:n)      # row numbers 
    jj = [dict[i] for i in str]  # column numbers
    p = length(colNames)        
    X  = sparse(ii,jj,val,n*nModels,p)
    return X, strip(modelTerm.trmString)*": ".*colNames   
end

function setRandom!(mme::MME,trmStrings::String,V,ped=false;
                estimate::Bool=false,
                νPrior::Float64=4.1,
                outSamples::Bool=false,
                outFreq::Int64=100
            )
    m = size(V,1)        
    if outSamples==true
        fileName = "V($trmStrings).samples"
        outStream = open(fileName,"w")
        for i=1:m,j=i:m
            if i==j==1
                @printf(outStream,"%6s", "V[$i,$j]")
            else
                @printf(outStream," %6s", "V[$i,$j]")    
            end        
        end
        @printf(outStream,"\n")        
    else
        outStream = nothing
    end        
    var = m==1 ? fill(V,1,1) : V        
    trmVec = strip.(split(trmStrings,[',';' '],keepempty=false))
    modelTerms = []
    [begin
        trmString = split(modelTerm.trmString,':')[2]
        if trmString in trmVec
            push!(modelTerms,modelTerm)
        end
     end
     for modelTerm in mme.modelTermVec
    ] 
    if length(modelTerms) != size(var,1)
        println("size of var: $(size(var,1)) does not match the number of random effects: $(length(modelTerms))")
        return
    end 
                
    if ped==false
        [trm.randomType = "I" for trm in modelTerms]
        Ai = I
    else
        [trm.randomType = "A" for trm in modelTerms] 
        mme.ped = pedigree
        Ai = JWAS.PedModule.AInverse(pedigree)    
    end
    Spr = (νPrior - m - 1)*var    
    randomEffect = RandomEffect(modelTerms,νPrior,Spr,inv(var),zero(var),Ai,
                        estimate,outSamples,outFreq,outStream)
    push!(mme.randomEffectsVec,randomEffect)   
end

function addGiMats!(mme::MME)
    for randomEffect in mme.randomEffectsVec
        for (i,modelTermi) in enumerate(randomEffect.modelTermVec), (j,modelTermj) in enumerate(randomEffect.modelTermVec)
            starti = modelTermi.startPos
            startj = modelTermj.startPos
            endi = modelTermi.endPos
            endj = modelTermj.endPos
            mme.mmeLhs[starti:endi,startj:endj] = mme.mmeLhs[starti:endi,startj:endj] + randomEffect.Ai*randomEffect.Vi[i,j]
        end
    end
end

addGiMats! (generic function with 1 method)

#EXECUTE
function getLhsRhs!(mme,df)
    n = size(df,1)
    X,colNames = getX(mme,mme.modelTermVec[1],mme.covVec,df)
    mme.modelTermVec[1].startPos = 1
    mme.modelTermVec[1].endPos  = mme.modelTermVec[1].startPos + size(X,2) - 1
    for i = 2:size(mme.modelTermVec,1)
        Xi,namesi = getX(mme,mme.modelTermVec[i],mme.covVec,df)
        X = [X Xi]
        mme.modelTermVec[i].startPos = mme.modelTermVec[i-1].endPos + 1
        mme.modelTermVec[i].endPos   = mme.modelTermVec[i].startPos + size(Xi,2) - 1
        colNames = [colNames; namesi]
    end
    y = Array{Float64,1}()
    for v in mme.depVarVec
        y = [y;df[:,Symbol(v)]]
    end
    Ri = kron(inv(mme.varRes),sparse(I, n, n))
    mme.X = X
    mme.mmeLhs = X'Ri*X
    mme.mmeRhs = X'Ri*y
    mme.mmeSpl = zero(mme.mmeRhs)
    mme.y      = y
    addGiMats!(mme)
    return mme.mmeLhs,mme.mmeRhs,colNames
end
function updateLhs!(mme)
    m = size(mme.varRes,1)
    n = Int(size(mme.y,1)/m)
    Ri = kron(inv(mme.varRes),sparse(I, n, n))
    mme.mmeLhs = mme.X'Ri*mme.X
    mme.mmeRhs = mme.X'Ri*mme.y
    addGiMats!(mme)
end

updateLhs! (generic function with 1 method)

Random.seed!(31415)
animal= ["animal1","animal2","animal3","animal4","animal5","animal6"]
sex   = ["m","f","f","m","f","f"]
breed = ["Angus","Angus","Hereford","Hereford","Angus","Angus"]
age   = [40,36,38,42,40,36]
df    = DataFrame(animal=animal,sex=sex,breed=breed,age=age,y=round.(randn(6),digits=3))

6 rows × 5 columns

animalsexbreedagey
StringStringStringInt64Float64
1animal1mAngus40-0.859
2animal2fAngus360.631
3animal3fHereford381.775
4animal4mHereford42-0.17
5animal5fAngus40-0.561
6animal6fAngus360.214
mme = initMME("y = intercept + sex + breed + age",varRes)
mme.covVec = ["age"];
lhs,rhs,names = getLhsRhs!(mme,df)
[names Matrix(lhs) rhs]

6×8 Array{Any,2}:
 "1:intercept"          6.0   2.0    4.0    4.0   2.0   232.0   1.03 
 "1:sex: m"             2.0   2.0    0.0    1.0   1.0    82.0  -1.029
 "1:sex: f"             4.0   0.0    4.0    3.0   1.0   150.0   2.059
 "1:breed: Angus"       4.0   1.0    3.0    4.0   0.0   152.0  -0.575
 "1:breed: Hereford"    2.0   1.0    1.0    0.0   2.0    80.0   1.605
 "1:age: age"         232.0  82.0  150.0  152.0  80.0  9000.0  33.93 

mme.modelTermVec

4-element Array{ModelTerm,1}:
 ModelTerm("1:intercept", 1, 1, 1, "fixed")
 ModelTerm("1:sex", 1, 2, 3, "fixed")      
 ModelTerm("1:breed", 1, 4, 5, "fixed")    
 ModelTerm("1:age", 1, 6, 6, "fixed")      

pedigree   = get_pedigree("pedFile",separator=",",header=false);

The delimiter in pedFile is ','.

coding pedigree... 100%|████████████████████████████████| Time: 0:00:00

Finished!

data = CSV.read("data.phen")

3 rows × 5 columns

IndMaty1y2x
Int64Int64Float64Float64Float64
1328.99.211.9
2429.75.710.8
3548.88.511.7
varRes = 1.0
mme = initMME("y1 = intercept + Ind",varRes);
setRandom!(mme,"Ind",1.0,pedigree);
#setRandom!(mme,"Ind",1.0);

lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs) rhs]

6×8 Array{Any,2}:
 "1:intercept"  3.0   0.0   0.0   1.0   1.0   1.0  27.4
 "1:Ind: 1"     0.0   2.0   1.0  -1.0  -1.0   0.0   0.0
 "1:Ind: 2"     0.0   1.0   2.0  -1.0  -1.0   0.0   0.0
 "1:Ind: 4"     1.0  -1.0  -1.0   3.5   0.5  -1.0   9.7
 "1:Ind: 3"     1.0  -1.0  -1.0   0.5   3.5  -1.0   8.9
 "1:Ind: 5"     1.0   0.0   0.0  -1.0  -1.0   3.0   8.8

Matrix(mme.randomEffectsVec[1].Ai)

5×5 Array{Float64,2}:
  2.0   1.0  -1.0  -1.0   0.0
  1.0   2.0  -1.0  -1.0   0.0
 -1.0  -1.0   2.5   0.5  -1.0
 -1.0  -1.0   0.5   2.5  -1.0
  0.0   0.0  -1.0  -1.0   2.0

size(mme.mmeLhs)

(6, 6)

Matrix(mme.X'mme.X)

6×6 Array{Float64,2}:
 3.0  0.0  0.0  1.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  0.0  0.0  1.0  0.0
 1.0  0.0  0.0  0.0  0.0  1.0

mme.modelTermVec

2-element Array{ModelTerm,1}:
 ModelTerm("intercept", 1, 1, "fixed")
 ModelTerm("Ind", 2, 6, "A")          

varRes = 1.0
mme = initMME("y1 = intercept + Ind + Mat",varRes);
setRandom!(mme,"Ind",1.0,pedigree)
setRandom!(mme,"Mat",1.0,pedigree);

lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs)]

11×12 Array{Any,2}:
 "intercept"  3.0   0.0   0.0   1.0   1.0   1.0   0.0   2.0   1.0   0.0   0.0
 "Ind: 1"     0.0   2.0   1.0  -1.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0
 "Ind: 2"     0.0   1.0   2.0  -1.0  -1.0   0.0   0.0   0.0   0.0   0.0   0.0
 "Ind: 4"     1.0  -1.0  -1.0   3.5   0.5  -1.0   0.0   1.0   0.0   0.0   0.0
 "Ind: 3"     1.0  -1.0  -1.0   0.5   3.5  -1.0   0.0   1.0   0.0   0.0   0.0
 "Ind: 5"     1.0   0.0   0.0  -1.0  -1.0   3.0   0.0   0.0   1.0   0.0   0.0
 "Mat: 1"     0.0   0.0   0.0   0.0   0.0   0.0   2.0   1.0  -1.0  -1.0   0.0
 "Mat: 2"     2.0   0.0   0.0   1.0   1.0   0.0   1.0   4.0  -1.0  -1.0   0.0
 "Mat: 4"     1.0   0.0   0.0   0.0   0.0   1.0  -1.0  -1.0   3.5   0.5  -1.0
 "Mat: 3"     0.0   0.0   0.0   0.0   0.0   0.0  -1.0  -1.0   0.5   2.5  -1.0
 "Mat: 5"     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -1.0  -1.0   2.0

Matrix(mme.X'mme.X)

11×11 Array{Float64,2}:
 3.0  0.0  0.0  1.0  1.0  1.0  0.0  2.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 2.0  0.0  0.0  1.0  1.0  0.0  0.0  2.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

V = [1   0.5
     0.5 2.0]

2×2 Array{Float64,2}:
 1.0  0.5
 0.5  2.0

varRes = 1.0
mme = initMME("y1 = intercept + Ind + Mat",varRes);
setRandom!(mme,"Ind, Mat",1.0,pedigree);

size of var: 1 does not match the number of random effects: 2

varRes = 1.0
mme = initMME("y1 = intercept + Ind + Mat",varRes);
setRandom!(mme,"Ind, Mat",V,pedigree);

lhs,rhs,names = getLhsRhs!(mme,data)
Matrix(lhs)

11×11 Array{Float64,2}:
 3.0   0.0        0.0        1.0       …   1.0        0.0        0.0     
 0.0   2.28571    1.14286   -1.14286       0.285714   0.285714   0.0     
 0.0   1.14286    2.28571   -1.14286       0.285714   0.285714   0.0     
 1.0  -1.14286   -1.14286    3.85714      -0.714286  -0.142857   0.285714
 1.0  -1.14286   -1.14286    0.571429     -0.142857  -0.714286   0.285714
 1.0   0.0        0.0       -1.14286   …   1.28571    0.285714  -0.571429
 0.0  -0.571429  -0.285714   0.285714     -0.571429  -0.571429   0.0     
 2.0  -0.285714  -0.571429   1.28571      -0.571429  -0.571429   0.0     
 1.0   0.285714   0.285714  -0.714286      2.42857    0.285714  -0.571429
 0.0   0.285714   0.285714  -0.142857      0.285714   1.42857   -0.571429
 0.0   0.0        0.0        0.285714  …  -0.571429  -0.571429   1.14286 

varRes = 1.0
mme = initMME("y1 = intercept + Ind + Ind*x",varRes);
mme.covVec = ["x"]
setRandom!(mme,"Ind",1.0,pedigree);

size(mme.randomEffectsVec,1)

1

setRandom!(mme,"Ind*x",1.0,pedigree);

size(mme.randomEffectsVec,1)

2

lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs)]

11×12 Array{Any,2}:
 "1:intercept"   3.0   0.0   0.0   1.0  …   0.0   0.0   10.8    11.9    11.7 
 "1:Ind: 1"      0.0   2.0   1.0  -1.0      0.0   0.0    0.0     0.0     0.0 
 "1:Ind: 2"      0.0   1.0   2.0  -1.0      0.0   0.0    0.0     0.0     0.0 
 "1:Ind: 4"      1.0  -1.0  -1.0   3.5      0.0   0.0   10.8     0.0     0.0 
 "1:Ind: 3"      1.0  -1.0  -1.0   0.5      0.0   0.0    0.0    11.9     0.0 
 "1:Ind: 5"      1.0   0.0   0.0  -1.0  …   0.0   0.0    0.0     0.0    11.7 
 "1:Ind*x: 1"    0.0   0.0   0.0   0.0      2.0   1.0   -1.0    -1.0     0.0 
 "1:Ind*x: 2"    0.0   0.0   0.0   0.0      1.0   2.0   -1.0    -1.0     0.0 
 "1:Ind*x: 4"   10.8   0.0   0.0  10.8     -1.0  -1.0  119.14    0.5    -1.0 
 "1:Ind*x: 3"   11.9   0.0   0.0   0.0     -1.0  -1.0    0.5   144.11   -1.0 
 "1:Ind*x: 5"   11.7   0.0   0.0   0.0  …   0.0   0.0   -1.0    -1.0   138.89

Matrix(mme.X'mme.X)

11×11 Array{Float64,2}:
  3.0  0.0  0.0   1.0   1.0   1.0  0.0  0.0   10.8    11.9    11.7 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
  1.0  0.0  0.0   1.0   0.0   0.0  0.0  0.0   10.8     0.0     0.0 
  1.0  0.0  0.0   0.0   1.0   0.0  0.0  0.0    0.0    11.9     0.0 
  1.0  0.0  0.0   0.0   0.0   1.0  0.0  0.0    0.0     0.0    11.7 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
 10.8  0.0  0.0  10.8   0.0   0.0  0.0  0.0  116.64    0.0     0.0 
 11.9  0.0  0.0   0.0  11.9   0.0  0.0  0.0    0.0   141.61    0.0 
 11.7  0.0  0.0   0.0   0.0  11.7  0.0  0.0    0.0     0.0   136.89

varRes = 1.0
mme = initMME("y1 = intercept + Ind + Ind*x",varRes);
mme.covVec = ["x"]
setRandom!(mme,"Ind, Ind*x",V,pedigree);

size(mme.randomEffectsVec,1)

1

lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs) rhs]

11×13 Array{Any,2}:
 "1:intercept"   3.0   0.0        0.0       …   11.9        11.7        27.4 
 "1:Ind: 1"      0.0   2.28571    1.14286        0.285714    0.0         0.0 
 "1:Ind: 2"      0.0   1.14286    2.28571        0.285714    0.0         0.0 
 "1:Ind: 4"      1.0  -1.14286   -1.14286       -0.142857    0.285714    9.7 
 "1:Ind: 3"      1.0  -1.14286   -1.14286       11.1857      0.285714    8.9 
 "1:Ind: 5"      1.0   0.0        0.0       …    0.285714   11.1286      8.8 
 "1:Ind*x: 1"    0.0  -0.571429  -0.285714      -0.571429    0.0         0.0 
 "1:Ind*x: 2"    0.0  -0.285714  -0.571429      -0.571429    0.0         0.0 
 "1:Ind*x: 4"   10.8   0.285714   0.285714       0.285714   -0.571429  104.76
 "1:Ind*x: 3"   11.9   0.285714   0.285714     143.039      -0.571429  105.91
 "1:Ind*x: 5"   11.7   0.0        0.0       …   -0.571429  138.033     102.96

Matrix(mme.X)

3×11 Array{Float64,2}:
 1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0   0.0  11.9   0.0
 1.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  10.8   0.0   0.0
 1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0   0.0   0.0  11.7

Matrix(mme.X'mme.X)

11×11 Array{Float64,2}:
  3.0  0.0  0.0   1.0   1.0   1.0  0.0  0.0   10.8    11.9    11.7 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
  1.0  0.0  0.0   1.0   0.0   0.0  0.0  0.0   10.8     0.0     0.0 
  1.0  0.0  0.0   0.0   1.0   0.0  0.0  0.0    0.0    11.9     0.0 
  1.0  0.0  0.0   0.0   0.0   1.0  0.0  0.0    0.0     0.0    11.7 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
  0.0  0.0  0.0   0.0   0.0   0.0  0.0  0.0    0.0     0.0     0.0 
 10.8  0.0  0.0  10.8   0.0   0.0  0.0  0.0  116.64    0.0     0.0 
 11.9  0.0  0.0   0.0  11.9   0.0  0.0  0.0    0.0   141.61    0.0 
 11.7  0.0  0.0   0.0   0.0  11.7  0.0  0.0    0.0     0.0   136.89

varRes = 1.0
mme = initMME("y1 = intercept + Ind + Ind*x",varRes);
mme.covVec = ["x"]
setRandom!(mme,"Ind",1.0,pedigree)
setRandom!(mme,"Ind*x",1.0);

lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs) rhs]

9×11 Array{Any,2}:
 "1:intercept"      3.0   0.0   0.0   1.0  …   11.9    10.8    11.7    27.4 
 "1:Ind: 1"         0.0   2.0   1.0  -1.0       0.0     0.0     0.0     0.0 
 "1:Ind: 2"         0.0   1.0   2.0  -1.0       0.0     0.0     0.0     0.0 
 "1:Ind: 4"         1.0  -1.0  -1.0   3.5       0.0    10.8     0.0     9.7 
 "1:Ind: 3"         1.0  -1.0  -1.0   0.5      11.9     0.0     0.0     8.9 
 "1:Ind: 5"         1.0   0.0   0.0  -1.0  …    0.0     0.0    11.7     8.8 
 "1:Ind*x: 3 x x"  11.9   0.0   0.0   0.0     142.61    0.0     0.0   105.91
 "1:Ind*x: 4 x x"  10.8   0.0   0.0  10.8       0.0   117.64    0.0   104.76
 "1:Ind*x: 5 x x"  11.7   0.0   0.0   0.0       0.0     0.0   137.89  102.96

varRes =[1   0.5
         0.5 2.0]
V = [1.0 0.0 0.0
     0.0 2.0 0.0
     0.0 0.0 3.0]
mme = initMME("y1 = intercept + Ind + Ind*x;
               y2 = intercept + Ind",varRes);
mme.covVec = ["x"]
#setRandom!(mme,"Ind Ind*x",V,pedigree);

lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs)]

11×12 Array{Any,2}:
 "1:intercept"      3.42857    1.14286   …  -0.285714  -0.285714  -0.285714
 "1:Ind: 3"         1.14286    1.14286      -0.285714   0.0        0.0     
 "1:Ind: 4"         1.14286    0.0           0.0       -0.285714   0.0     
 "1:Ind: 5"         1.14286    0.0           0.0        0.0       -0.285714
 "1:Ind*x: 3 x x"  13.6       13.6          -3.4        0.0        0.0     
 "1:Ind*x: 4 x x"  12.3429     0.0       …   0.0       -3.08571    0.0     
 "1:Ind*x: 5 x x"  13.3714     0.0           0.0        0.0       -3.34286 
 "2:intercept"     -0.857143  -0.285714      0.571429   0.571429   0.571429
 "2:Ind: 3"        -0.285714  -0.285714      0.571429   0.0        0.0     
 "2:Ind: 4"        -0.285714   0.0           0.0        0.571429   0.0     
 "2:Ind: 5"        -0.285714   0.0       …   0.0        0.0        0.571429

mme.randomEffectsVec[1].modelTermVec

3-element Array{ModelTerm,1}:
 ModelTerm("1:Ind", 1, 2, 6, "A")   
 ModelTerm("1:Ind*x", 1, 7, 11, "A")
 ModelTerm("2:Ind", 2, 13, 17, "A") 

inv(V)

3×3 Array{Float64,2}:
 1.0  0.0  0.0     
 0.0  0.5  0.0     
 0.0  0.0  0.333333

varRes =[1   0.5
         0.5 2.0]
V = [1.0 0.0 0.0
     0.0 2.0 0.0
     0.0 0.0 3.0]
mme = initMME("y1 = intercept + Ind + Ind*x;
               y2 = intercept + Ind",varRes);
mme.covVec = ["x"]
setRandom!(mme,"Ind Ind*x",V,pedigree);
lhs,rhs,names = getLhsRhs!(mme,data)
[names Matrix(lhs)]

17×18 Array{Any,2}:
 "1:intercept"   3.42857    0.0   0.0  …  -0.285714  -0.285714  -0.285714
 "1:Ind: 1"      0.0        2.0   1.0      0.0        0.0        0.0     
 "1:Ind: 2"      0.0        1.0   2.0      0.0        0.0        0.0     
 "1:Ind: 4"      1.14286   -1.0  -1.0     -0.285714   0.0        0.0     
 "1:Ind: 3"      1.14286   -1.0  -1.0      0.0       -0.285714   0.0     
 "1:Ind: 5"      1.14286    0.0   0.0  …   0.0        0.0       -0.285714
 "1:Ind*x: 1"    0.0        0.0   0.0      0.0        0.0        0.0     
 "1:Ind*x: 2"    0.0        0.0   0.0      0.0        0.0        0.0     
 "1:Ind*x: 4"   12.3429     0.0   0.0     -3.08571    0.0        0.0     
 "1:Ind*x: 3"   13.6        0.0   0.0      0.0       -3.4        0.0     
 "1:Ind*x: 5"   13.3714     0.0   0.0  …   0.0        0.0       -3.34286 
 "2:intercept"  -0.857143   0.0   0.0      0.571429   0.571429   0.571429
 "2:Ind: 1"      0.0        0.0   0.0     -0.333333  -0.333333   0.0     
 "2:Ind: 2"      0.0        0.0   0.0     -0.333333  -0.333333   0.0     
 "2:Ind: 4"     -0.285714   0.0   0.0      1.40476    0.166667  -0.333333
 "2:Ind: 3"     -0.285714   0.0   0.0  …   0.166667   1.40476   -0.333333
 "2:Ind: 5"     -0.285714   0.0   0.0     -0.333333  -0.333333   1.2381  

pedigree = get_pedigree("simData.ped",separator=",",header=false);

The delimiter in simData.ped is ','.
Pedigree information:
#individuals: 1248
#sires:       50
#dams:        198
#founders:    248

using Statistics

function sampleLoc!(mme)
    A = mme.mmeLhs
    b = mme.mmeRhs
    x = mme.mmeSpl
    n = size(x,1)
    for i=1:n
        cVarInv = 1.0/A[i,i]
        cMean   = cVarInv*(b[i] - A[:,i]'x) + x[i]
        x[i]    = randn()*sqrt(cVarInv) + cMean
    end
end

sampleLoc! (generic function with 1 method)

18.6. Full-conditional Posteriors for Covariance Matrices in Multi-trait Models#

Full-conditional posteriors for covariance matrices in multi-trait models are given here

18.7. Gibbs Sampler for Covariance Matrices#

function sampleVar!(mme,iter=0)
    for randomEffect in mme.randomEffectsVec
        if randomEffect.estimate == false continue end
        modelTerm1 = randomEffect.modelTermVec[1]      
        k = modelTerm1.endPos - modelTerm1.startPos
        m = size(randomEffect.modelTermVec,1)
        
        S = zeros(m,m)
        for i=1:m
            modelTermi = randomEffect.modelTermVec[i]
            starti = modelTermi.startPos
            endi = modelTermi.endPos
            ui = mme.mmeSpl[starti:endi]
            for j=i:m
                modelTermj = randomEffect.modelTermVec[j]
                startj = modelTermj.startPos
                endj = modelTermj.endPos
                uj = mme.mmeSpl[startj:endj]
                S[i,j] = ui'randomEffect.Ai*uj
                S[j,i] = S[i,j]
            end
        end
        Vpo = randomEffect.Spr + S
        νpo = randomEffect.νpr + k 
        V = rand(InverseWishart(νpo,Vpo))
        randomEffect.Vi = inv(V)
        if iter>0
            randomEffect.meanV += (V - randomEffect.meanV)/iter
        end
        if randomEffect.outSamples==true && iter%randomEffect.outFreq==0
            for i=1:m,j=i:m
                if i==j==1
                    @printf(randomEffect.outStream,"%10.5e", V[i,j])
                else
                    @printf(randomEffect.outStream," %10.5e", V[i,j])    
                end
            end
            @printf(randomEffect.outStream,"\n")
        end
            
    end
end

sampleVar! (generic function with 2 methods)

18.8. Gibbs Sampler for Residual Covariance Matrix#

function sampleVarRes!(mme,iter=0)
    m = size(mme.varRes,1)
    n = Int(size(mme.y,1)/m)
    S = zeros(m,m)
    e = mme.y - mme.X*mme.mmeSpl
    for i=1:m
        starti = (i - 1)*n + 1
        endi = starti + n - 1
        ei = e[starti:endi]
        for j=i:m
            startj = (j - 1)*n + 1
            endj = startj + n - 1
            ej = e[startj:endj]
            S[i,j] = ei'ej
            S[j,i] = S[i,j]
        end
    end
    νpo = n + mme.νRes
    Spo = S + mme.SRes
    mme.varRes = rand(InverseWishart(νpo,Spo))
    if iter>0
        mme.meanVarRes += (mme.varRes - mme.meanVarRes)/iter
    end  
end

sampleVarRes! (generic function with 2 methods)

using Distributions

using Plots

using ProgressMeter
function runMCMC!(mme,nIter,burnIn)
    @showprogress "MCMC sampling" for iter = 1:nIter
        sampleLoc!(mme)
        sampleVar!(mme,iter-burnIn)
        sampleVarRes!(mme,iter-burnIn)
        updateLhsRhs!(mme)
    end
    for randomEffect in mme.randomEffectsVec
        if randomEffect.outStream != nothing close(randomEffect.outStream) end
    end
end

runMCMC! (generic function with 1 method)

function runMean(x)
    meanx = 0
    [meanx += (v-meanx)/i for (i,v) in enumerate(x)]
end        

runMean (generic function with 1 method)

18.9. Numerical Example#

dataTwoTrait = CSV.read("dataTwoTrait.phen",DataFrame)
first(dataTwoTrait,5)

5 rows × 3 columns

Indy1y2
Int64Float64Float64
12514.9750.312
22522.069-4.109
32533.413-0.052
42548.806-1.066
52555.759-1.057
varGen = [10.0 2.0; 2.0 3.0]
varRes = [10.0 0.0; 0.0 3.0]
mme = initMME("y1 = intercept + Ind
               y2 = intercept + Ind",varRes);
setRandom!(mme,"Ind",varGen,pedigree,estimate=true,outSamples=true);
lhs,rhs,names = getLhsRhs!(mme,dataTwoTrait);
nIter  = 10_000
burnIn = 1_000

1000

mme.randomEffectsVec[1].modelTermVec

2-element Vector{ModelTerm}:
 ModelTerm("1:Ind", 1, 2, 1249, "A")
 ModelTerm("2:Ind", 2, 1251, 2498, "A")

runMCMC!(mme,nIter,burnIn)

MCMC sampling 100%|██████████████████████████████████████| Time: 0:00:15

mme.randomEffectsVec[1].meanV

2×2 Matrix{Float64}:
 11.3768   2.11806
  2.11806  2.9284

mme.meanVarRes

2×2 Matrix{Float64}:
 9.19382  1.80055
 1.80055  3.05443