\[ \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}} \]

19. Extensions for Models with SNP Covariates#

We have already seen how the MME can be constructed from the Data in a DataFrame and a higher-level description of the model. The model description included a string representation of the model equation, identification of the quantitative factors in the model, and specification of the random terms in the model. Here, we will see how substitution effects of SNPs can be included in the model as random effects in the model. In multiple-trait analyses, a single command will be used to include SNP effects in all the models.

19.1. Breeding value models and SNP effects models#

When genomic data is not available, consider the pedigree-based breeding value model,

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

where \(\mbf{u}\) is a vector of breeding values. Then, given a pedigree \(\mbf{P}\) for the individuals, the covariance matrix of \(\mbf{u}\) is

\[ \Var(\mbf{u}|\mbf{P}) = \mbf{A}\sigma^2_a, \]

where \(\mbf{A}\) is the additive genetic relationship matrix and \(\sigma^2_a\) is the additive genetic variance.

Suppose SNP genotypes are available on each of the \(n\) individuals at \(k\) loci. Let \(\mbf{M}\) denote the \(n\times k\) matrix of SNP genotype covariates and \(\bs{\alpha}\) denote their effects. Although from a frequentist point of view \(\bs{\alpha}\) is a fixed effect, for Bayesian analyses, \(\bs{\alpha}\) is often assigned a normal prior with null mean and covariance matrix: \(\Var(\bs{\alpha}) = \mbf{I}\sigma^2_{\alpha}\). Then, the breeding values can be modeled as:

\[ \mbf{u} = \mbf{M}\bs{\alpha}. \]

Assuming that the \(\mbf{Z}\) matrix in (1) is an identity matrix, the model for the phenotypes can be wrtten as either:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{M}\bs{\alpha} + \mathbf{e}, \label{eq.MEM} \]

or as

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mbf{u} + \mathbf{e}. \tag{2} \]

In (2), conditional on the SNP genotypes, the covariance matrix of \(\mbf{u}\) is:

\[\begin{split} \begin{align} \Var(\mbf{u}|\mbf{M}) &= \mbf{MM}'\sigma^2_{\alpha}\\ &= \frac{\mbf{MM}'}{k\overline{2pq}}k\overline{2pq}\sigma^2_{\alpha}\\ &= \mbf{G}\sigma^2_u, \end{align} \end{split}\]

where \(\mbf{G} = \frac{\mbf{MM}'}{k\overline{2pq}}\) and \(\sigma^2_u = k\overline{2pq}\sigma^2_{\alpha}\). The model given by (2) is a genomic breeding value model.

Assuming, for simplicity that \(\mbf{R} = \mbf{I}\), mixed model equations (MME) that correspend to model (1) is:

\[\begin{split} \begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{X}'\mathbf{M} \\ \mathbf{M}'\mathbf{X} & \mathbf{M}'\mathbf{M} + \frac{\sigma_{e}^{2}}{\sigma_{\alpha}^2}\mathbf{I} \end{bmatrix} \begin{bmatrix} \hat{\bs{\beta}}\\ \hat{\bs{\alpha}} \end{bmatrix} = \begin{bmatrix} \mathbf{X}'\mathbf{y}\\ \mathbf{M}'\mathbf{y} \end{bmatrix}, \label{eq.MMEMEM} \end{split}\]

and the MME corresponding to (2) is:

\[\begin{split} \begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{X}' \\ \mathbf{X} & \mathbf{I} + \frac{\sigma_{e}^{2}}{\sigma_u^2}\mathbf{G}^{-1} \end{bmatrix} \begin{bmatrix} \hat{\bs{\beta}}\\ \hat{\mbf{u}} \end{bmatrix} = \begin{bmatrix} \mathbf{X}'\mathbf{y}\\ \mathbf{y} \end{bmatrix}, \label{MMEBVM} \end{split}\]

where \(\mbf{G} = \frac{\mbf{MM}'}{k\overline{2pq}}\) and \(\sigma^2_u = k\overline{2pq}\sigma^2_{\alpha}\).

19.2. Extension for SNP Effects#

19.2.1. Single-trait model#

We will initially consider a univariate mixed linear model of the form:

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Zu} + \mathbf{W}\boldsymbol{\alpha} + \mathbf{e}, \tag{model} \]

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 covariance matrix \(\mathbf{G}\sigma^2_u\), \(\mathbf{Z}\) is an observed matrix that relates \(\mathbf{u}\) to \(\mathbf{y}\), \(\boldsymbol{\alpha}\) is a vector of random SNP effects with null means and covariance matrix \(\mathbf{I}\sigma^2_{\alpha}\), \(\mathbf{W}\) is an observed matrix of SNP covariates, and \(\mathbf{e}\) is a vector of multivariate normal residuals with null means and covariance matrix \(\mathbf{R}\sigma^2_e\). The fixed effects are usually assigned a flat prior distribution, and the matrices \(\mathbf{G}\) and \(\mathbf{R}\) are assumed to be known. It is assumed that vectors \(\mbf{u}\), \(\bs{\alpha}\) and \(\mbf{e}\) are mutually independent. The variance components \(\sigma^2_u\), \(\sigma^2_{\alpha}\) and \(\sigma^2_e\) are assigned scaled inverted chi-square prior distributions.

The MME for this model are:

\[\begin{split} \begin{bmatrix} \mathbf{X}'\mathbf{R}^{-1}\mathbf{X} & \mathbf{X}'\mathbf{R}^{-1}\mathbf{Z} & \mathbf{X}'\mathbf{R}^{-1}\mathbf{W} \\ \mathbf{Z}'\mathbf{R}^{-1}\mathbf{X} & \mathbf{Z}'\mathbf{R}^{-1}\mathbf{Z} + \frac{\sigma_{e}^{2}}{\sigma_u^2}\mathbf{G}^{-1}& \mathbf{Z}'\mathbf{R}^{-1}\mathbf{W} \\ \mathbf{W}'\mathbf{R}^{-1}\mathbf{X} & \mathbf{W}'\mathbf{R}^{-1}\mathbf{Z} & \mathbf{W}'\mathbf{R}^{-1}\mathbf{W} + \frac{\sigma_{e}^{2}}{\sigma_{\alpha}^2}\mathbf{I} \end{bmatrix} \begin{bmatrix} \hat{\bs{\beta}}\\ \hat{\mbf{u}} \\ \hat{\bs{\alpha}} \end{bmatrix} = \begin{bmatrix} \mathbf{X}'\mathbf{R}^{-1}\mathbf{y}\\ \mathbf{Z}'\mathbf{R}^{-1}\mathbf{y}\\ \mathbf{W}'\mathbf{R}^{-1}\mathbf{y} \end{bmatrix}. (MME) \end{split}\]

The rows and columns in Equation (MME) corresponding to \(\hat{\bs{\alpha}}\) are dense. Thus, an alternative sampling scheme is described below for sampling the location parameters in Equation (model) without constrcting the MME in Equation (MME).

In iteration \(i\) of the Gibbs sampler, conditional on the sampled values of \(\bs{\alpha} = \bs{\alpha}^{(i-1)}\) in the previous iteration, the model given by Equation (model) can be written as:

\[ \begin{align} \mathbf{y} - \mathbf{W}\boldsymbol{\alpha}^{(i-1)} &= \mathbf{X}\boldsymbol{\beta} + \mathbf{Zu} + \mathbf{e}, \tag{model2} \end{align} \]

where distributions of \(\mbf{u}\), and \(\mbf{e}\) in Equation (model2) are identical to their distibutions in Equation (model), because \(\bs{\alpha}\) is independent of \(\mbf{u}\), and \(\mbf{e}\). The MME for Equation (model2) will be sparse when \(\mathbf{X}\), \(\mathbf{Z}\) and \(\mathbf{G}^{-1}\) are sparse as in the usual non-genomic models. Thus, the function sampleLoc! can be applied to the sparse MME for Equation (model2) to obtain single-site Gibbs samples of the elements in \(\bs{\beta}\) and \(\mbf{u}\). Recall that this function can be used for single-trait and multi-trait models.

Similarly, for a single-trait model, to apply single-site Gibbs sampling to element \(j\) of \(\bs{\alpha}\), the model in Equation (model) can be written as

\[\begin{split} \begin{align} \mbf{y} - \mbf{X}\bs{\beta}^{(i)} - \mbf{Z}\mbf{u}^{(i)} - \sum_{l < j}\mbf{w}_l\alpha_l^{(i)} - \sum_{l >j}\mbf{w}_l\alpha_l^{(i-1)} &= \mbf{w}_j\alpha_j+ \mathbf{e}\\ \mbf{y}_{\text{adj-j}}^{(j)} &= \mbf{w}_j\alpha_j + \mathbf{e}, \tag{model3} \end{align} \end{split}\]

where \(\mbf{w}_j\) is column \(j\) of \(\mbf{W}\). In this model, \(\alpha_j\) is the only unknown, and the MME corresponding to Equation (model3) is:

\[ (\mbf{w}_j'\mbf{w}_j + \frac{\sigma^2_e}{\sigma^2_{\alpha}})\hat{\alpha}_j = \mbf{w}_j'\mbf{y}_{\text{adj-j}}^{(j)}. \tag{mmeSNP} \]

Thus, as we have seen previously, the full-conditional distribution of \(\alpha_j\) can be shown to have a normal distribution with mean \(\hat{\alpha}_j\) and variance \(\frac{\sigma^2_e}{\mbf{w}_j'\mbf{w}_j + \frac{\sigma^2_e}{\sigma^2_{\alpha}}}\).

19.2.1.1. An Efficient Approach#

To compute the right-hand-side (RHS) of (mmeSNP) efficiently, before sampling \(\alpha_j\), let

\[ \mbf{y}_{\text{adj}}^{(j-1)} = \mbf{y} - \mbf{X}\bs{\beta}^{(i)} - \mbf{Z}\mbf{u}^{(i)} - \sum_{l < j}\mbf{w}_l\alpha_l^{(i)} - \mbf{w}_j\alpha_j^{(i-1)} - \sum_{l > j}\mbf{w}_l\alpha_l^{(i-1)}, \]

where \(\mbf{y}\) has been adjusted for all the effects in the model with their sampled values, including \(\alpha_j = \alpha_j^{(i-1)}\). Then, the RHS of Equation (mmeSNP) can be written as

\[ \mbf{w}_j'\mbf{y}_{\text{adj-j}}^{(j)} = \mbf{w}_j'\mbf{y}_{\text{adj}}^{(j)} + \mbf{w}_j'\mbf{w}_j\alpha_j^{(i-1)}. \]

After sampling \(\alpha_j\), \(\mbf{y}_{\text{adj}}^{{j}}\) is obtained efficiently as

\[\begin{split} \begin{align} \mbf{y}_{\text{adj}}^{(j)} &= \mbf{y} - \mbf{X}\bs{\beta}^{(i)} - \mbf{Z}\mbf{u}^{(i)} - \sum_{l < j}\mbf{w}_l\alpha_l^{(i)} - \mbf{w}_j\alpha_j^{(i)} - \sum_{l > j}\mbf{w}_l\alpha_l^{(i-1)}\\ &= \mbf{y}_{\text{adj}}^{(j-1)} + \mbf{w}_j\alpha_j^{(i-1)} - \mbf{w}_j\alpha_j^{(i)} \\ &= \mbf{y}_{\text{adj}}^{(j-1)} + \mbf{w}_j(\alpha_j^{(i-1)} - \alpha_j^{(i)}). \end{align} \end{split}\]

After all elements of \(\bs{\alpha}\) have been sampled, adding \(\mbf{X}\bs{\beta}^{(i)} + \mbf{Z}\mbf{u}^{(i)}\) to \(\mbf{y}_{\text{adj}}^{{k}}\) gives:

\[ \mathbf{y} - \mathbf{W}\boldsymbol{\alpha}^{(i)} = \mbf{y}_{\text{adj}}^{{k}} + \mbf{X}\bs{\beta}^{(i)} + \mbf{Z}\mbf{u}^{(i)}, \]

which will be used in the next round of the Gibbs sampler for obtaining samples of \(\bs{\beta}\) and \(\mbf{u}\).

19.2.2. Multi-trait model#

Consider now a two-trait model, where, to simplify the notation, we will denote \(\mbf{y}_{\text{adj-j}}^{{j}}\) by \(\mbf{y}_1\) for the first trait and by \(\mbf{y}_2\) for the second trait. Then, the model for sampling the effects for locus \(j\) can be written as:

\[\begin{split} \begin{align} \begin{bmatrix} \mbf{y}_1\\ \mbf{y}_2 \end{bmatrix} &= \begin{bmatrix} \mbf{w}_j & \mbf{0}\\ \mbf{0} & \mbf{w}_j \end{bmatrix} \begin{bmatrix} \alpha_{j1}\\ \alpha_{j2} \end{bmatrix} + \begin{bmatrix} \mbf{e}_1\\ \mbf{e}_2 \end{bmatrix}\\ &= (\mbf{I}_2\otimes \mbf{w}_j)\bs{\alpha}_j + \mbf{e}, \tag{21} \end{align} \end{split}\]

where \(\bs{\alpha}_j \sim \text{N}(\mbf{0},\bs{\Sigma})\), and \(\mbf{e} \sim \text{N}(\mbf{0},\mbf{R}_0\otimes\mbf{I}_n)\). The MME for this model is:

\[\begin{split} \begin{align} [(\mbf{I}_2\otimes \mbf{w}'_j)(\mbf{R}_0^{-1}\otimes\mbf{I}_n)(\mbf{I}_2\otimes \mbf{w}_j) + \bs{\Sigma}^{-1}]\hat{\bs{\alpha}}_j &= (\mbf{I}_2\otimes \mbf{w}'_j)(\mbf{R}_0^{-1}\otimes\mbf{I}_n)\mbf{y} \\ (\begin{bmatrix} \mbf{w}'_j & \mbf{0}\\ \mbf{0} & \mbf{w}'_j \end{bmatrix} \begin{bmatrix} r_0^{11}\mbf{I} & r_0^{12}\mbf{I}\\ r_0^{21}\mbf{I} & r_0^{22}\mbf{I} \end{bmatrix} \begin{bmatrix} \mbf{w}_j & \mbf{0}\\ \mbf{0} & \mbf{w}_j \end{bmatrix} + \bs{\Sigma}^{-1}) \hat{\bs{\alpha}}_j &= (\begin{bmatrix} \mbf{w}'_j & \mbf{0}\\ \mbf{0} & \mbf{w}'_j \end{bmatrix} \begin{bmatrix} r_0^{11}\mbf{I} & r_0^{12}\mbf{I}\\ r_0^{21}\mbf{I} & r_0^{22}\mbf{I} \end{bmatrix} \begin{bmatrix} \mbf{y}_1\\ \mbf{y}_2 \end{bmatrix} \\ (\mbf{w}'_j\mbf{w}_j\mbf{R}_0^{-1} + \bs{\Sigma}^{-1})\hat{\bs{\alpha}}_j &= \begin{bmatrix} r_0^{11}\mbf{w}' & r_0^{12}\mbf{w}'\\ r_0^{21}\mbf{w}' & r_0^{22}\mbf{w}' \end{bmatrix} \begin{bmatrix} \mbf{y}_1\\ \mbf{y}_2 \end{bmatrix} \\ (\mbf{w}'_j\mbf{w}_j\mbf{R}_0^{-1} + \bs{\Sigma}^{-1})\hat{\bs{\alpha}}_j &= \mbf{R}_0^{-1} \begin{bmatrix} \mbf{w}'_j\mbf{y}_1\\ \mbf{w}'_j\mbf{y}_2 \end{bmatrix} %\label{eq:MME2Trait} \end{align} \end{split}\]

When the multi-trait model is for \(m\) traits, the coefficient matrix of the MME in Equation (21) expands to an \(m\times m\) matrix and the RHS to an \(m\times 1\) vector. Then, the sampleLoc! function can be applied to this MME to obtain single-site Gibbs samples of the elements in \(\bs{\alpha}_j\).

19.2.3. Structures for Building MME#

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

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

# Composite type for including SNP effects in the model
mutable struct Genotypes                    # ZM where Z is a sparse incidence matrix and M is a dense matrix of covariates           
    trmString::AbstractString               # "Ztrm"
    Z::SparseMatrixCSC{Float64,Int64}       # incidence matrix Z to map rows of M to phenotypes according to "Ind" 
    M::Array{Float64,2}                     # dense matrix of covariates
    MArray::Array{Any,1}                    # column views of M
    MPMArray::Array{Float64,1}              # diagonals of M'M
    idGeno::Array{String}                   # row names of M (Ind must be a subset of idGeno)
    α::Array{Float64,2}                     # nCovariates x nTraits matrix of effects
    αMean::Array{Float64,2}                 # mean of sampled values  
    yAdjArray::Array{Any,1}                 # Array of Arrays that will be used to access adjusted y-vector
    νpr::Float64
    Spr::Array{Float64,2}
    Vi::Array{Float64,2}
    meanV::Array{Float64,2}
    estimate::Bool
    outSamples::Bool
    outFreq::Int64
    outStream   
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}
    genotypes                               # dense component of model
    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 parameters
    meanEffects::Array{Float64,1}           #mean of location parameters
    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}                     #dependent variable
    yAdj::Array{Float64,1}
    meanVarRes::Array{Float64,2}                     
end 

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

┌ Info: Precompiling DataFrames [a93c6f00-e57d-5684-b7b6-d8193f3e46c0]
└ @ Base loading.jl:1273
┌ Info: Precompiling JWAS [c9a035f4-d403-5e6b-8649-6be755bc4798]
└ @ Base loading.jl:1273

ModelTerm

Ri = [1.0 0.5;0.5 2.0]

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

19.2.4. Functions for Building MME#

#EXECUTE 
function initMME(modelEquations::AbstractString,varRes;νRes=5.0)  # "y1 = A + A*B; y2 = A + B + !(A)M
    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
    genotypes = false
    randomEffectsVec = []
    ped = false
    covVec = []
    X = spzeros(0,0)
    mmeLhs = spzeros(0,0)
    mmeRhs = []
    mmeSpl = []
    meanEffects = []
    y      = Array{Float64,1}(undef,0)
    yAdj   = 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,genotypes,randomEffectsVec,ped,
               depVarVec,X,mmeLhs,mmeRhs,mmeSpl,meanEffects,varRes,νRes,SRes,y,yAdj,meanVarRes)
end 

initMME (generic function with 1 method)

19.2.5. Function for Adding SNP Effects to Model#

#EXECUTE
function addGenotypes!(mme::MME,df::DataFrame,ZTrm::AbstractString,M::Array{Float64,2},
            idGeno::Array{String,1},V;
            estimate::Bool=false,
            νPrior::Float64=4.1,
            outSamples::Bool=false,
            outFreq::Int64=100
        )
        m = size(V,1)        
        if outSamples==true
            fileName = "G.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
        d = Dict()
        for (i,s) = enumerate(idGeno)
            d[s] = i
        end
        n = size(df,1) 
        ii = collect(1:n) 
        jj = Array{Int64,1}(undef,n)
        idPheno = string.(df[!,Symbol(ZTrm)])    
        for (i,v) = enumerate(idPheno)
            jj[i] = d[v]
        end
        Z = sparse(ii,jj,1.0,n,size(idGeno,1))
        α = zeros(0,0)
        αMean = zeros(0,0)    
        Spr = (νPrior - m - 1)*var 
        yAdjArray = Array{Array{Float64,1},1}(undef,0)
        MArray = Array{Any,1}(undef,0)
        MPMArray = Array{Any,1}(undef,0)  
        mme.genotypes = Genotypes(ZTrm,Z,M,MArray,MPMArray,idGeno,α,αMean,yAdjArray,νPrior,Spr,
            inv(var),zero(var),estimate,outSamples,outFreq,outStream)
        return
end

addGenotypes! (generic function with 1 method)

#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])])
    #str = modelTerm.randomType != "A" ? str : string.(df[:, Symbol(factors[factors .!= mme.covVec][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 = []
    for modelTerm in mme.modelTermVec
        trmString = split(modelTerm.trmString,':')[2]
        if trmString in trmVec
            push!(modelTerms,modelTerm)
        end                
    end
    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.meanEffects = zero(mme.mmeRhs)
    mme.y      = y
    mme.yAdj = copy(mme.y)
    addGiMats!(mme)
    if mme.genotypes != false
        mme.genotypes.M = mme.genotypes.M .- mean(mme.genotypes.M,dims=1)
        mme.genotypes.M = mme.genotypes.Z * mme.genotypes.M
        nModels  = size(mme.varRes,1)
        nMarkers = size(mme.genotypes.M,2)
        yAdjArray = Array{Any,1}(undef,nModels)
        for i=1:nModels
            startPos = (i-1)*n + 1
            endPos = startPos + n - 1
            yAdjArray[i] = @view mme.yAdj[startPos:endPos]
        end
        MArray   = Array{Any,1}(undef,nMarkers)
        MPMArray = Array{Float64,1}(undef,nMarkers)
        for i=1:nMarkers
            m = @view mme.genotypes.M[:,i]
            MArray[i] = m
            MPMArray[i] = m'm
        end
        mme.genotypes.yAdjArray = yAdjArray
        mme.genotypes.MArray = MArray
        mme.genotypes.MPMArray = MPMArray
        mme.genotypes.α = zeros(nMarkers,nModels)
        mme.genotypes.αMean = zeros(nMarkers,nModels)
    end
    return mme.mmeLhs,mme.mmeRhs,colNames
end
function updateLhsRhs!(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.yAdj
    addGiMats!(mme)
end

updateLhsRhs! (generic function with 1 method)

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

The delimiter in pedFile is ','.

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

Finished!

typeof(pedigree)

JWAS.PedModule.Pedigree

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);

M = float.(rand(Binomial(2,0.5),5,10))

5×10 Array{Float64,2}:
 1.0  2.0  1.0  0.0  2.0  1.0  2.0  1.0  1.0  1.0
 1.0  2.0  2.0  1.0  2.0  0.0  0.0  0.0  0.0  1.0
 2.0  2.0  1.0  1.0  1.0  2.0  1.0  2.0  2.0  1.0
 2.0  0.0  2.0  1.0  1.0  2.0  1.0  0.0  2.0  1.0
 1.0  0.0  1.0  2.0  0.0  1.0  1.0  0.0  2.0  2.0

idGeno = string.(1:5);
addGenotypes!(mme,data,"Ind",M,idGeno,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

mme.genotypes.M

3×10 Array{Float64,2}:
  0.6   0.8  -0.4  0.0  -0.2   0.8  0.0   1.4  0.6  -0.2
  0.6  -1.2   0.6  0.0  -0.2   0.8  0.0  -0.6  0.6  -0.2
 -0.4  -1.2  -0.4  1.0  -1.2  -0.2  0.0  -0.6  0.6   0.8

mean(mme.genotypes.M,dims=1)

1×10 Array{Float64,2}:
 0.266667  -0.533333  -0.0666667  0.333333  …  0.0  0.0666667  0.6  0.133333

mme.genotypes.MArray

10-element Array{Any,1}:
 [0.6, 0.6, -0.4]  
 [0.8, -1.2, -1.2] 
 [-0.4, 0.6, -0.4] 
 [0.0, 0.0, 1.0]   
 [-0.2, -0.2, -1.2]
 [0.8, 0.8, -0.2]  
 [0.0, 0.0, 0.0]   
 [1.4, -0.6, -0.6] 
 [0.6, 0.6, 0.6]   
 [-0.2, -0.2, 0.8] 

mme.genotypes.MPMArray

10-element Array{Float64,1}:
 0.8799999999999998
 3.5200000000000005
 0.68              
 1.0               
 1.5200000000000005
 1.3199999999999994
 0.0               
 2.6799999999999997
 1.0799999999999994
 0.7199999999999999

Matrix(mme.genotypes.Z)

3×5 Array{Float64,2}:
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

mme.genotypes.yAdjArray[1]'

1×3 Adjoint{Float64,SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}}:
 8.9  9.7  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

mme.genotypes.yAdjArray

1-element Array{Array{Float64,1},1}:
 [8.9, 9.7, 8.8]

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);

mme.randomEffectsVec[1].modelTermVec

3-element Array{ModelTerm,1}:
 ModelTerm("1:Ind", 1, 0, 0, "A")  
 ModelTerm("1:Ind*x", 1, 0, 0, "A")
 ModelTerm("2:Ind", 2, 0, 0, "A")  

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  

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

mme.randomEffectsVec[1].modelTermVec

2-element Array{ModelTerm,1}:
 ModelTerm("1:Ind", 1, 0, 0, "A")
 ModelTerm("2:Ind", 2, 0, 0, "A")

y = float.(1:10)

10-element Array{Float64,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0

yArray = []
push!(yArray,unsafe_wrap(Array,pointer(y,1),5))
push!(yArray,unsafe_wrap(Array,pointer(y,6),5))

2-element Array{Any,1}:
 [1.0, 2.0, 3.0, 4.0, 5.0] 
 [6.0, 7.0, 8.0, 9.0, 10.0]

yArray[2][2] = 7.7
y

10-element Array{Float64,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.7
  8.0
  9.0
 10.0

v = []
push!(v,@view y[1:5])

1-element Array{Any,1}:
 [1.0, 2.0, 3.0, 4.0, 5.0]

@view y[1:5]

5-element view(::Array{Float64,1}, 1:5) with eltype Float64:
 1.0
 2.0
 3.0
 4.0
 5.0

push!(v,@view y[6:10])

2-element Array{Any,1}:
 [1.0, 2.0, 3.0, 4.0, 5.0] 
 [6.0, 7.0, 8.0, 9.0, 10.0]

typeof(v[1])

SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}

v[1][2] = 2.3

2.3

y

10-element Array{Float64,1}:
  1.0
  2.3
  3.0
  4.0
  5.0
  6.0
  7.7
  8.0
  9.0
 10.0

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

The delimiter in pedFile is ','.

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

Finished!

#JWAS.PedModule.getinfo(pedigree)

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

3 rows × 5 columns

IndMaty1y2x
Int64Int64Float64Float64Float64
1328.99.211.9
2429.75.710.8
3548.88.511.7
using Statistics, Distributions

M = float.(rand(Binomial(2,0.5),5,10))

5×10 Array{Float64,2}:
 1.0  1.0  1.0  0.0  2.0  2.0  1.0  2.0  1.0  1.0
 0.0  0.0  2.0  0.0  2.0  0.0  1.0  0.0  2.0  2.0
 0.0  0.0  1.0  0.0  2.0  1.0  2.0  1.0  1.0  0.0
 0.0  1.0  0.0  0.0  1.0  2.0  1.0  0.0  1.0  1.0
 2.0  1.0  0.0  0.0  2.0  1.0  0.0  2.0  2.0  1.0

var(data[:,3])

0.24333333333333262

function sampleLoc!(mme,iter=0)
    A = mme.mmeLhs
    b = mme.mmeRhs
    x = mme.mmeSpl
    n = size(x,1)
    sampleLoc!(A,b,x,n)
    if iter > 0
        mme.meanEffects += (x - mme.meanEffects)/iter
    end
end

sampleLoc! (generic function with 2 methods)

function sampleLoc!(A,b,x,n)
    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 3 methods)

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 = @view mme.mmeSpl[starti:endi]
            for j=i:m
                modelTermj = randomEffect.modelTermVec[j]
                startj = modelTermj.startPos
                endj = modelTermj.endPos
                uj = @view 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)

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 = @view e[starti:endi]
        for j=i:m
            startj = (j - 1)*n + 1
            endj = startj + n - 1
            ej = @view 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)

function sampleGenEffects!(mme,iter=0)
    # y has already been corrected for SNP effects from previous iteration
    mme.yAdj[:] -= mme.X*mme.mmeSpl # adjust for effects in sparse MME
    nMarkers = size(mme.genotypes.M,2)
    Ri = inv(mme.varRes)
    nModels = size(mme.varRes,1)
    rhs = zeros(nModels)
    mmeGen = mme.genotypes
    for i=1:nMarkers
        lhs = mmeGen.MPMArray[i]*Ri + mmeGen.Vi
        for j=1:nModels
            rhs[j] = mmeGen.MArray[i]'mmeGen.yAdjArray[j] + mmeGen.MPMArray[i]*mmeGen.α[i,j]
        end
        rhs = Ri*rhs
        x = @view mmeGen.α[i,:]
        old = copy(x)
        sampleLoc!(lhs,rhs,x,nModels)
        for j=1:nModels
            mmeGen.yAdjArray[j][:] += mmeGen.MArray[i]*(old[j] - x[j])
        end
    end
    mme.yAdj[:] += mme.X*mme.mmeSpl # unadjust for effects in sparse MME
    if iter > 0
        mmeGen.αMean += (mmeGen.α - mmeGen.αMean)/iter
    end
    nothing
end

sampleGenEffects! (generic function with 2 methods)

using ProgressMeter
function runMCMC!(mme,nIter,burnIn)
    @showprogress "MCMC sampling" for iter = 1:nIter
        sampleLoc!(mme,iter-burnIn)
        if mme.genotypes !=false sampleGenEffects!(mme,iter-burnIn) end
        #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
    if mme.genotypes !=false && mme.genotypes.outStream != nothing close(mme.genotypes.outStream) end
end

runMCMC! (generic function with 1 method)

varGen = 5.0
varRes = 10.0
mme = initMME("y1 = intercept",varRes);
#setRandom!(mme,"Ind",varGen,pedigree,estimate=false);
#setRandom!(mme,"Ind",varGen,estimate=false);

idGeno = string.(1:5);
varSNPEffect = varGen/size(M,2)
addGenotypes!(mme,data,"Ind",M,idGeno,varSNPEffect,estimate=false,outSamples=false)

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

nIter  = 100_000
burnIn = 1_0000
runMCMC!(mme,nIter,burnIn)

MCMC sampling100%|██████████████████████████████████████| Time: 0:00:13

for i=1:size(M,2)
    colNum  = 5+i
    colName = Symbol("m$i")
    insertcols!(data,colNum,colName => mme.genotypes.M[:,i])
end

data

3 rows × 15 columns (omitted printing of 5 columns)

IndMaty1y2xm1m2m3m4m5
Int64Int64Float64Float64Float64Float64Float64Float64Float64Float64
1328.99.211.9-0.6-0.60.20.00.2
2429.75.710.8-0.60.4-0.80.0-0.8
3548.88.511.71.40.4-0.80.00.2
modelEqn = "y1 = intercept"
for i=1:size(M,2)
    modelEqn = modelEqn *  " + m$i"
end

modelEqn

"y1 = intercept + m1 + m2 + m3 + m4 + m5 + m6 + m7 + m8 + m9 + m10"

mme1 = initMME(modelEqn,varRes);
#setRandom!(mme1,"Ind",varGen,pedigree,estimate=false);
#setRandom!(mme1,"Ind",varGen,estimate=false);

varSNPEffect

0.5

mme1.covVec = []
for i=1:size(M,2)
    randomTrm = "m$i"
    push!(mme1.covVec,randomTrm)
    setRandom!(mme1,randomTrm,varSNPEffect)
end

nIter  = 100_000
burnIn = 1_000
lhs,rhs,names = getLhsRhs!(mme1,data)
runMCMC!(mme1,nIter,burnIn)

MCMC sampling100%|██████████████████████████████████████| Time: 0:00:05

names

11-element Array{String,1}:
 "1:intercept"
 "1:m1: m1"   
 "1:m2: m2"   
 "1:m3: m3"   
 "1:m4: m4"   
 "1:m5: m5"   
 "1:m6: m6"   
 "1:m7: m7"   
 "1:m8: m8"   
 "1:m9: m9"   
 "1:m10: m10" 

[names[1] mme1.meanEffects[1] mme.meanEffects]

1×3 Array{Any,2}:
 "1:intercept"  9.13158  9.12796

[mme1.meanEffects[2:11] mme.genotypes.αMean]

10×2 Array{Float64,2}:
 -0.0274914   -0.0246814  
  0.00936227   0.0102446  
 -0.00985962  -0.0128209  
  0.00381277  -0.00354914 
 -0.0224857   -0.0261506  
  0.0223116    0.0234557  
  0.00152241   0.000699012
 -0.0375149   -0.0357298  
 -0.0118666   -0.0147853  
  0.0131752    0.0123046