Mixed Model Equations (MME) from Data and Model (known variance components)
Contents
14. Mixed Model Equations (MME) from Data and Model (known variance components)#
Algorithms to Construct and Solve Henderson’s Mixed Model Equations (MME).
Concepts to build and solve Henderson’s MME are introduced here.
Here, we will see how the MME can be constructed from the Data in a DataFrame and a higher-level description of the model. The model description will include a string representation of the model equation, identification of the quantitative factors in the model, and specification of the random terms in the model. In the previous notebook, we constructed the left-hand-side (LHS) and the right-hand-side (RHS) of the normal equations for a fixed-effects model. In order to construct the MME for a mixed model, we need to identify the random effects in the model and also determine the row and column positions for these random effects in the LHS of the MME. In order to store this information needed to construct the MME, we will use the following stuctures (composite types).
14.1. Structures for Building MME#
#EXECUTE
using DataFrames, SparseArrays, LinearAlgebra, Random, JWAS, CSV
mutable struct ModelTerm
trmString::AbstractString
startPos::Int64
endPos::Int64
randomType::String
end
mutable struct RandomEffect
modelTermVec::Array{ModelTerm,1}
Vi
Ai
end
mutable struct MME
modelEquation::AbstractString #"y = A + B + 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
depVar::Symbol #:y
X::SparseMatrixCSC{Float64,Int64} #design matrix
mmeLhs::SparseMatrixCSC{Float64,Int64} #left-hand side for MME
mmeRhs::Array{Float64,1} #right-hand side for MME
mmeSol::Array{Float64,1} #a solution of the MME
varRes::Float64
end
function ModelTerm(str::AbstractString)
ModelTerm(str,0,0,"fixed")
end
ModelTerm
14.2. Functions for Building MME#
#EXECUTE
function initMME(modelEquation::AbstractString,varRes) # "y = A + A*B"
if modelEquation==""
error("modelEquation is empty\n")
end
leftRight = split(modelEquation,"=") # "y", "A+A*B"
depVar = Symbol(strip(leftRight[1])) # :y
modelParts = strip(leftRight[2]) #"A+A*B"
termsVec = split(modelParts,"+") #"A","A*B"
modelTermVec = [ModelTerm(strip(trmStr)) for trmStr in termsVec]
modelTermDict = Dict{AbstractString,ModelTerm}()
for i in modelTermVec
modelTermDict[i.trmString] = i
end
randomEffectsVec = []
ped = false
covVec = []
X = spzeros(0,0)
mmeLhs = spzeros(0,0)
mmeRhs = []
mmeSol = []
return MME(modelEquation,covVec,modelTermVec,modelTermDict,randomEffectsVec,ped,depVar,X,mmeLhs,mmeRhs,mmeSol,varRes)
end
initMME (generic function with 1 method)
varRes = 1.0
mme = initMME("y = sex + breed + age",varRes)
mme.covVec = ["age"];
mme.modelEquation
"y = sex + breed + age"
mme.modelTermVec
3-element Array{ModelTerm,1}:
ModelTerm("sex", 0, 0, "fixed")
ModelTerm("breed", 0, 0, "fixed")
ModelTerm("age", 0, 0, "fixed")
mme.depVar
:y
a = mme.modelTermDict["sex"]
ModelTerm("sex", 0, 0, "fixed")
a.startPos
0
#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)
trmString = modelTerm.trmString
if trmString == "intercept"
X = ones(n,1)
colNames = ["intercept"]
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 = 1:n # row numbers
jj = [dict[i] for i in str] # column numbers
p = length(colNames)
X = sparse(ii,jj,val,n,p)
return X, strip(modelTerm.trmString)*": ".*colNames
end
function setRandom(mme::MME,trmStrings::String,V,ped=false)
m = size(V,1)
var = m==1 ? fill(V,1,1) : V
trmVec = strip.(split(trmStrings,[',';' '],keepempty=false))
if length(trmVec) != size(var,1)
println("size of var: $(size(var,1)) does not match the number of random effects: $(length(trmVec))")
return
end
modelTerms = [mme.modelTermDict[trm] for trm in trmVec]
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
randomEffect = RandomEffect(modelTerms,inv(var),Ai)
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]*mme.varRes
end
end
end
addGiMats (generic function with 1 method)
#EXECUTE
function getLhsRhs(mme,df)
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 = df[:,mme.depVar]
mme.X = X
mme.mmeLhs = X'X
mme.mmeRhs = X'y
addGiMats(mme)
return mme.mmeLhs,mme.mmeRhs,colNames
end
getLhsRhs (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))
| animal | sex | breed | age | y | |
|---|---|---|---|---|---|
| String | String | String | Int64 | Float64 | |
| 1 | animal1 | m | Angus | 40 | -0.859 |
| 2 | animal2 | f | Angus | 36 | 0.631 |
| 3 | animal3 | f | Hereford | 38 | 1.775 |
| 4 | animal4 | m | Hereford | 42 | -0.17 |
| 5 | animal5 | f | Angus | 40 | -0.561 |
| 6 | animal6 | f | Angus | 36 | 0.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}:
"intercept" 6.0 2.0 4.0 4.0 2.0 232.0 1.03
"sex: m" 2.0 2.0 0.0 1.0 1.0 82.0 -1.029
"sex: f" 4.0 0.0 4.0 3.0 1.0 150.0 2.059
"breed: Angus" 4.0 1.0 3.0 4.0 0.0 152.0 -0.575
"breed: Hereford" 2.0 1.0 1.0 0.0 2.0 80.0 1.605
"age: age" 232.0 82.0 150.0 152.0 80.0 9000.0 33.93
mme.modelTermVec
4-element Array{ModelTerm,1}:
ModelTerm("intercept", 1, 1, "fixed")
ModelTerm("sex", 2, 3, "fixed")
ModelTerm("breed", 4, 5, "fixed")
ModelTerm("age", 6, 6, "fixed")
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!
data = CSV.read("data.phen")
| Ind | Mat | y1 | y2 | x | |
|---|---|---|---|---|---|
| Int64 | Int64 | Float64 | Float64 | Float64 | |
| 1 | 3 | 2 | 8.9 | 9.2 | 11.9 |
| 2 | 4 | 2 | 9.7 | 5.7 | 10.8 |
| 3 | 5 | 4 | 8.8 | 8.5 | 11.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}:
"intercept" 3.0 0.0 0.0 1.0 1.0 1.0 27.4
"Ind: 1" 0.0 2.0 1.0 -1.0 -1.0 0.0 0.0
"Ind: 2" 0.0 1.0 2.0 -1.0 -1.0 0.0 0.0
"Ind: 4" 1.0 -1.0 -1.0 3.5 0.5 -1.0 9.7
"Ind: 3" 1.0 -1.0 -1.0 0.5 3.5 -1.0 8.9
"Ind: 5" 1.0 0.0 0.0 -1.0 -1.0 3.0 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}:
"intercept" 3.0 0.0 0.0 1.0 … 0.0 0.0 10.8 11.9 11.7
"Ind: 1" 0.0 2.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0
"Ind: 2" 0.0 1.0 2.0 -1.0 0.0 0.0 0.0 0.0 0.0
"Ind: 4" 1.0 -1.0 -1.0 3.5 0.0 0.0 10.8 0.0 0.0
"Ind: 3" 1.0 -1.0 -1.0 0.5 0.0 0.0 0.0 11.9 0.0
"Ind: 5" 1.0 0.0 0.0 -1.0 … 0.0 0.0 0.0 0.0 11.7
"Ind*x: 1" 0.0 0.0 0.0 0.0 2.0 1.0 -1.0 -1.0 0.0
"Ind*x: 2" 0.0 0.0 0.0 0.0 1.0 2.0 -1.0 -1.0 0.0
"Ind*x: 4" 10.8 0.0 0.0 10.8 -1.0 -1.0 119.14 0.5 -1.0
"Ind*x: 3" 11.9 0.0 0.0 0.0 -1.0 -1.0 0.5 144.11 -1.0
"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}:
"intercept" 3.0 0.0 0.0 … 11.9 11.7 27.4
"Ind: 1" 0.0 2.28571 1.14286 0.285714 0.0 0.0
"Ind: 2" 0.0 1.14286 2.28571 0.285714 0.0 0.0
"Ind: 4" 1.0 -1.14286 -1.14286 -0.142857 0.285714 9.7
"Ind: 3" 1.0 -1.14286 -1.14286 11.1857 0.285714 8.9
"Ind: 5" 1.0 0.0 0.0 … 0.285714 11.1286 8.8
"Ind*x: 1" 0.0 -0.571429 -0.285714 -0.571429 0.0 0.0
"Ind*x: 2" 0.0 -0.285714 -0.571429 -0.571429 0.0 0.0
"Ind*x: 4" 10.8 0.285714 0.285714 0.285714 -0.571429 104.76
"Ind*x: 3" 11.9 0.285714 0.285714 143.039 -0.571429 105.91
"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}:
"intercept" 3.0 0.0 0.0 1.0 … 11.9 10.8 11.7 27.4
"Ind: 1" 0.0 2.0 1.0 -1.0 0.0 0.0 0.0 0.0
"Ind: 2" 0.0 1.0 2.0 -1.0 0.0 0.0 0.0 0.0
"Ind: 4" 1.0 -1.0 -1.0 3.5 0.0 10.8 0.0 9.7
"Ind: 3" 1.0 -1.0 -1.0 0.5 11.9 0.0 0.0 8.9
"Ind: 5" 1.0 0.0 0.0 -1.0 … 0.0 0.0 11.7 8.8
"Ind*x: 3 x x" 11.9 0.0 0.0 0.0 142.61 0.0 0.0 105.91
"Ind*x: 4 x x" 10.8 0.0 0.0 10.8 0.0 117.64 0.0 104.76
"Ind*x: 5 x x" 11.7 0.0 0.0 0.0 0.0 0.0 137.89 102.96