N3 | 动物模型计算BLUP值
1, 数据
这次使用一个PPT里面的数据, 用R语言演示一下如何做BLUP值计算.下面是生成数据的代码
Chang <- c(1,1,1,2,2)
ID <- c(1,2,3,4,5)
Sire <- c(0,0,1,1,3)
Dam <- c(0,0,0,2,2)
weight <- c(140,152,135,143,160)
dat <- data.frame(Chang,ID,Sire,Dam,weight)datChang | ID | Sire | Dam | weight |
|---|---|---|---|---|
1 | 1 | 0 | 0 | 140 |
1 | 2 | 0 | 0 | 152 |
1 | 3 | 1 | 0 | 135 |
2 | 4 | 1 | 2 | 143 |
2 | 5 | 3 | 2 | 160 |
2, 计算亲缘关系逆矩阵
library(nadiv)提取系谱信息
ped <- dat[,2:4]
pedID | Sire | Dam |
|---|---|---|
1 | 0 | 0 |
2 | 0 | 0 |
3 | 1 | 0 |
4 | 1 | 2 |
5 | 3 | 2 |
计算亲缘关系逆矩阵
pped = prepPed(ped)
ppedWarning message in prepPed(ped):
"Zero in the dam column interpreted as a missing parent"Warning message in prepPed(ped):
"Zero in the sire column interpreted as a missing parent"首先, 将系谱进行一下转换, 使用nadiv的prepPed函数, 预处理. 它会自动不齐没有亲本的个体, 变为NA.
ID | Sire | Dam |
|---|---|---|
1 | NA | NA |
2 | NA | NA |
3 | 1 | NA |
4 | 1 | 2 |
5 | 3 | 2 |
如果是计算逆矩阵的矩阵形式, 可以使用makeAinv(pped)$Ainv
Ainv = makeAinv(pped)$Ainv
Ainv5 x 5 sparse Matrix of class "dgCMatrix"
1 1.8333333 0.5 -0.6666667 -1 .
2 0.5000000 2.0 0.5000000 -1 -1
3 -0.6666667 0.5 1.8333333 . -1
4 -1.0000000 -1.0 . 2 .
5 . -1.0 -1.0000000 . 2如果是计算逆矩阵的行列形式, 可以使用makeAinv(pped)$listAinv
makeAinv(pped)$listAinvrow | column | Ainv | |
|---|---|---|---|
1 | 1 | 1 | 1.8333333 |
5 | 2 | 1 | 0.5000000 |
6 | 2 | 2 | 2.0000000 |
10 | 3 | 1 | -0.6666667 |
11 | 3 | 2 | 0.5000000 |
12 | 3 | 3 | 1.8333333 |
14 | 4 | 1 | -1.0000000 |
15 | 4 | 2 | -1.0000000 |
16 | 4 | 4 | 2.0000000 |
17 | 5 | 2 | -1.0000000 |
18 | 5 | 3 | -1.0000000 |
19 | 5 | 5 | 2.0000000 |
教科书的结果, 两者一样
3, 构建模型
$$ y = Xb + Zu + e $$
构建固定因子矩阵
这里使用函数model.matrix构建矩阵, 比较方便
for(i in 1:4) dat[,i] <- as.factor(dat[,i])
X <- model.matrix(~Chang-1,dat)
XChang1 | Chang2 | |
|---|---|---|
1 | 1 | 0 |
2 | 1 | 0 |
3 | 1 | 0 |
4 | 0 | 1 |
5 | 0 | 1 |
构建单元矩阵
Z <- diag(length(unique(dat$ID)))
Z1 | 0 | 0 | 0 | 0 |
|---|---|---|---|---|
0 | 1 | 0 | 0 | 0 |
0 | 0 | 1 | 0 | 0 |
0 | 0 | 0 | 1 | 0 |
0 | 0 | 0 | 0 | 1 |
构建y的矩阵
y <- as.matrix(dat$weight)
y140 |
|---|
152 |
135 |
143 |
160 |
混合线性方程组
XpZ <- crossprod(X,Z);XpZChang1 | 1 | 1 | 1 | 0 | 0 |
|---|---|---|---|---|---|
Chang2 | 0 | 0 | 0 | 1 | 1 |
X’X
XpX <- crossprod(X) ;XpXChang1 | Chang2 | |
|---|---|---|
Chang1 | 3 | 0 |
Chang2 | 0 | 2 |
Z’X
ZpX <- crossprod(Z,X);ZpXChang1 | Chang2 |
|---|---|
1 | 0 |
1 | 0 |
1 | 0 |
0 | 1 |
0 | 1 |
Z’Z
ZpZ <- crossprod(Z);ZpZ1 | 0 | 0 | 0 | 0 |
|---|---|---|---|---|
0 | 1 | 0 | 0 | 0 |
0 | 0 | 1 | 0 | 0 |
0 | 0 | 0 | 1 | 0 |
0 | 0 | 0 | 0 | 1 |
X’y
Xpy <- crossprod(X,y);XpyChang1 | 427 |
|---|---|
Chang2 | 303 |
Z’y
Zpy <- crossprod(Z,y);Zpy140 |
|---|
152 |
135 |
143 |
160 |
K
K <- 2;K2
LHS <- rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ+Ainv*K))
LHS7 x 7 sparse Matrix of class "dgCMatrix"
Chang1 Chang2
Chang1 3 . 1.000000 1 1.000000 . .
Chang2 . 2 . . . 1 1
1 1 . 4.666667 1 -1.333333 -2 .
2 1 . 1.000000 5 1.000000 -2 -2
3 1 . -1.333333 1 4.666667 . -2
4 . 1 -2.000000 -2 . 5 .
5 . 1 . -2 -2.000000 . 5可以看到, 里面的LHS左手矩阵和上图结果一致.
RHS <- rbind(Xpy,Zpy)
RHS
求解BLUP值
solve(LHS)%*%RHS7 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 142.842105
[2,] 151.118421
[3,] -2.462551
[4,] 3.052632
[5,] -2.116397
[6,] -1.387652
[7,] 2.150810最后，为了方便大家重演，我将相关的生产数据代码，运行代码汇总如下：
Chang <- c(1,1,1,2,2)ID <- c(1,2,3,4,5)Sire <- c(0,0,1,1,3)Dam <- c(0,0,0,2,2)weight <- c(140,152,135,143,160)dat <- data.frame(Chang,ID,Sire,Dam,weight)dat
library(nadiv)ped <- dat[,2:4]ped
pped = prepPed(ped)pped
Ainv = makeAinv(pped)$AinvAinv
makeAinv(pped)$listAinv
for(i in 1:4) dat[,i] <- as.factor(dat[,i])X <- model.matrix(~Chang-1,dat)X
Z <- diag(length(unique(dat$ID)))Z
y <- as.matrix(dat$weight)y
XpZ <- crossprod(X,Z);XpZ
XpX <- crossprod(X) ;XpX
ZpX <- crossprod(Z,X);ZpX
ZpZ <- crossprod(Z);ZpZ
Xpy <- crossprod(X,y);Xpy
Zpy <- crossprod(Z,y);Zpy
K <- 2;K
LHS <- rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ+Ainv*K))LHS
RHS <- rbind(Xpy,Zpy)RHS
solve(LHS)%*%RHS