单性状动物模型矩阵形式计算BLUP值

公众号叫：“育种数据分析之放飞自我”，但是竟然没有一个相关的专辑，因此，决定将之前写过的关于育种数据分析相关的博文，重新整理一下，包括软件包有些更新后语法变化了，有些软件包不在维护有了更好的替代品。经常有人给我发过来很久以前的博客，告诉我运行报错，然后我发现有些排版时错误的，有些R包已经不再维护退出CRAN了，因此重新整理一下，确保发出去的东西能够运行，整理一下之前的博文。

专业性的博文增加生活的段子，其乐无穷。

段子:

同事去大骨头, 吃过之后进行了抽奖, 然后向我炫耀, 看我抽到了一个小猪佩奇. 照片如下:

我看过之后, 冷静的说, 首先, 我们去吃饭, 奖品都是送的, 不需要抽奖的.

其次, 这是乔治, 不是佩奇.

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)

dat

2, 计算亲缘关系逆矩阵

library(nadiv)

提取系谱信息

ped <- dat[,2:4]

ped

计算亲缘关系逆矩阵

pped = prepPed(ped)

pped

Warning 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.

如果是计算逆矩阵的矩阵形式, 可以使用makeAinv(pped)$Ainv

Ainv = makeAinv(pped)$Ainv

Ainv

5 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)$listAinv

教科书的结果, 两者一样

3, 构建模型

$$ y = Xb + Zu + e $$

构建固定因子矩阵

这里使用函数model.matrix构建矩阵, 比较方便

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的矩阵

y <- as.matrix(dat$weight)
y

混合线性方程组

XpZ <- crossprod(X,Z);XpZ

X’X

XpX <- crossprod(X) ;XpX

Z’X

ZpX <- crossprod(Z,X);ZpX

Z’Z

ZpZ <- crossprod(Z);ZpZ

X’y

Xpy <- crossprod(X,y);Xpy

Z’y

Zpy <- crossprod(Z,y);Zpy

K

K <- 2;K

2

LHS <- rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ+Ainv*K))
LHS

7 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)%*%RHS

7 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

可以看到, 结果虽然结果不一致, 但是PPT里面的结果是错误的…

所以说，PPT里面的内容也不一定是正确的，现场演示之后，发现PPT里面的结果是错误的，这该如何圆场？？？现场翻车记！！！

最后，为了方便大家重演，我将相关的生产数据代码，运行代码汇总如下：

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