The Animal Model
Julius Mugambe
2021 01 14
This is a page for my R publications R Pubs. You can checkout some of my few code.
Animal model example problem with R
Consider the mixed linear equation below: \[\begin{equation}y=Xb+Zu+e\end{equation}\] Variance components
sigma_e = 40
sigma_a = 20
alpha = sigma_e/sigma_a
alpha## [1] 2Set working directory
setwd("E:/R-MME")
# print working directory
getwd()## [1] "E:/R-MME"Loading the data (phenotypes)
data <- read.table("data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data## id sex wwg
## 1 4 Male 4.5
## 2 5 Female 2.9
## 3 6 Female 3.9
## 4 7 Male 3.5
## 5 8 Male 5.0y = data[,3]
y## [1] 4.5 2.9 3.9 3.5 5.0Load the pedigree file
# read ped file
ped <- read.table("ped.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
ped## calve sire dam
## 1 1 0 0
## 2 2 0 0
## 3 3 0 0
## 4 4 1 0
## 5 5 3 2
## 6 6 1 2
## 7 7 4 5
## 8 8 3 6dim(ped)## [1] 8 3head(ped)## calve sire dam
## 1 1 0 0
## 2 2 0 0
## 3 3 0 0
## 4 4 1 0
## 5 5 3 2
## 6 6 1 2tail(ped)## calve sire dam
## 3 3 0 0
## 4 4 1 0
## 5 5 3 2
## 6 6 1 2
## 7 7 4 5
## 8 8 3 6Set up the X matrix Which is a matrix of ones ( no effect) which is a 8 by 1 matrix in this case
Method
Function to make a design matrix
desgn <- function(v) {
if (is.numeric(v)) {
va = v
mrow = length(va)
mcol = max(va)
}
if (is.character(v)) {
vf = factor(v)
# Save the index value of each level
va = as.numeric(vf)
mrow = length(va)
mcol = length(levels(vf))
}
# Initializing X-matrix
X = matrix(data = c(0), nrow = mrow, ncol = mcol)
for (i in 1:mrow) {
ic = va[i]
X[i, ic] = 1
}
return(X)
}The X matrix
X = desgn(data[,2])
X## [,1] [,2]
## [1,] 0 1
## [2,] 1 0
## [3,] 1 0
## [4,] 0 1
## [5,] 0 1The Z matrix will then be:
Desgn_M <- function(vnum) {
mrow = length(vnum)
mcol = max(vnum)
W= matrix(data=c(0),nrow = mrow,ncol=mcol)
for (i in 1:mrow){ic = vnum[i]
W[i,ic]=1}
return(W)
}
anim = c(4:8) # animals with records ( the first four animals have no records)
Z = Desgn_M(anim)
Z## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 0 0 1 0 0 0 0
## [2,] 0 0 0 0 1 0 0 0
## [3,] 0 0 0 0 0 1 0 0
## [4,] 0 0 0 0 0 0 1 0
## [5,] 0 0 0 0 0 0 0 1The number of animals according to the data will be:
# number of animal
no_of_animals = nrow(ped)
no_of_animals## [1] 8Construction of the numerator relationship matrix
# First NRN Initialization
A = matrix(rep(0, no_of_animals * no_of_animals), nrow = no_of_animals)
# empty NRM
A## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0# Reading and filling the NRM one by one
for (i in 1:no_of_animals) {
animal = ped[i, 1]
sire = ped[i, 2]
dam = ped[i, 3]
if (sire == 0 & dam == 0) {
# Unknown parents
A[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
#
for (j in 1:animal - 1) {
A[j, animal] = 0.5 * (A[j, sire])
A[animal, j] = A[j, animal]
}
A[animal, animal] = 1
} else if (sire == 0 & dam != 0) {
#
for (j in 1:animal - 1) {
A[j, animal] = 0.5 * (A[j, dam])
A[animal, j] = A[j, animal]
}
A[animal, animal] = 1
} else {
#
for (j in 1:animal - 1) {
A[j, animal] = 0.5 * (A[j, sire] + A[j, dam])
A[animal, j] = A[j, animal]
}
A[animal, animal] = 1 + 0.5 * A[sire, dam]
}
}The NRM will thus be:
A## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1.00 0.00 0.00 0.500 0.000 0.50 0.25 0.250
## [2,] 0.00 1.00 0.00 0.000 0.500 0.50 0.25 0.250
## [3,] 0.00 0.00 1.00 0.000 0.500 0.00 0.25 0.500
## [4,] 0.50 0.00 0.00 1.000 0.000 0.25 0.50 0.125
## [5,] 0.00 0.50 0.50 0.000 1.000 0.25 0.50 0.375
## [6,] 0.50 0.50 0.00 0.250 0.250 1.00 0.25 0.500
## [7,] 0.25 0.25 0.25 0.500 0.500 0.25 1.00 0.250
## [8,] 0.25 0.25 0.50 0.125 0.375 0.50 0.25 1.000Now construct the inverse of the numerator relationship matrix \(A^{-1}\)
# install 'MASS' packages
# install.packages("MASS")
# activate package
require("MASS")## 필요한 패키지를 로딩중입니다: MASS\(A^{-1}\) will then be:
AINV <-as.matrix( ginv(A))
AINV## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.833333e+00 5.000000e-01 3.885781e-15 -6.666667e-01 6.661338e-16
## [2,] 5.000000e-01 2.000000e+00 5.000000e-01 4.440892e-16 -1.000000e+00
## [3,] 1.165734e-15 5.000000e-01 2.000000e+00 1.054712e-15 -1.000000e+00
## [4,] -6.666667e-01 -9.992007e-16 -1.998401e-15 1.833333e+00 5.000000e-01
## [5,] -2.220446e-16 -1.000000e+00 -1.000000e+00 5.000000e-01 2.500000e+00
## [6,] -1.000000e+00 -1.000000e+00 5.000000e-01 -7.771561e-16 -7.771561e-16
## [7,] 0.000000e+00 -1.110223e-15 -6.661338e-16 -1.000000e+00 -1.000000e+00
## [8,] -9.402201e-16 -4.378442e-15 -1.000000e+00 -3.028827e-15 1.672273e-15
## [,6] [,7] [,8]
## [1,] -1.000000e+00 -5.384582e-15 -2.681882e-15
## [2,] -1.000000e+00 -2.220446e-16 -2.185752e-16
## [3,] 5.000000e-01 -4.440892e-16 -1.000000e+00
## [4,] 4.996004e-16 -1.000000e+00 1.797174e-15
## [5,] 5.551115e-16 -1.000000e+00 -1.700029e-15
## [6,] 2.500000e+00 5.551115e-16 -1.000000e+00
## [7,] -4.440892e-16 2.000000e+00 1.002670e-15
## [8,] -1.000000e+00 3.084338e-15 2.000000e+00Setting up the Mixed Model Equations
\(\begin{bmatrix}X^\prime X & X^\prime Z\\ Z^\prime X & Z^\prime Z \end{bmatrix} \begin{bmatrix} \hat{b} \\ \hat{a} \end{bmatrix} = \begin{bmatrix} X^\prime y \\ Z^\prime y \end{bmatrix}\)
# The left hand side
LHS = rbind(
cbind(t(X) %*% X, t(X) %*% Z),
cbind(t(Z) %*% X, t(Z) %*% Z + AINV*alpha))
LHS## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 2 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 0 3 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00
## [3,] 0 0 3.666667e+00 1.000000e+00 7.771561e-15 -1.333333e+00
## [4,] 0 0 1.000000e+00 4.000000e+00 1.000000e+00 8.881784e-16
## [5,] 0 0 2.331468e-15 1.000000e+00 4.000000e+00 2.109424e-15
## [6,] 0 1 -1.333333e+00 -1.998401e-15 -3.996803e-15 4.666667e+00
## [7,] 1 0 -4.440892e-16 -2.000000e+00 -2.000000e+00 1.000000e+00
## [8,] 1 0 -2.000000e+00 -2.000000e+00 1.000000e+00 -1.554312e-15
## [9,] 0 1 0.000000e+00 -2.220446e-15 -1.332268e-15 -2.000000e+00
## [10,] 0 1 -1.880440e-15 -8.756884e-15 -2.000000e+00 -6.057654e-15
## [,7] [,8] [,9] [,10]
## [1,] 1.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00
## [2,] 0.000000e+00 0.000000e+00 1.000000e+00 1.000000e+00
## [3,] 1.332268e-15 -2.000000e+00 -1.076916e-14 -5.363765e-15
## [4,] -2.000000e+00 -2.000000e+00 -4.440892e-16 -4.371503e-16
## [5,] -2.000000e+00 1.000000e+00 -8.881784e-16 -2.000000e+00
## [6,] 1.000000e+00 9.992007e-16 -2.000000e+00 3.594347e-15
## [7,] 6.000000e+00 1.110223e-15 -2.000000e+00 -3.400058e-15
## [8,] -1.554312e-15 6.000000e+00 1.110223e-15 -2.000000e+00
## [9,] -2.000000e+00 -8.881784e-16 5.000000e+00 2.005340e-15
## [10,] 3.344547e-15 -2.000000e+00 6.168677e-15 5.000000e+00RHS = rbind(t(X) %*% y,
t(Z) %*% y)
RHS## [,1]
## [1,] 6.8
## [2,] 13.0
## [3,] 0.0
## [4,] 0.0
## [5,] 0.0
## [6,] 4.5
## [7,] 2.9
## [8,] 3.9
## [9,] 3.5
## [10,] 5.0The solution from the equations will thus be;
Sol = round(ginv(LHS) %*% RHS,3)
Sol## [,1]
## [1,] 3.404
## [2,] 4.359
## [3,] 0.098
## [4,] -0.019
## [5,] -0.041
## [6,] -0.009
## [7,] -0.186
## [8,] 0.177
## [9,] -0.249
## [10,] 0.183Plot
par(mar = c(4, 4, 0.1, 0.1))
plot(pressure, pch = 19, type = "b")
Figures
knitr::include_graphics(rep("pig.jfif", 2))
Circularize a plot