单性状动物模型矩阵形式计算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)

dat

Chang	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]

ped

ID	Sire	Dam
1	0	0
2	0	0
3	1	0
4	1	2
5	3	2

计算亲缘关系逆矩阵

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.

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

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

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

X

	Chang1	Chang2
1	1	0
2	1	0
3	1	0
4	0	1
5	0	1

构建单元矩阵

Z <- diag(length(unique(dat$ID)))
Z

0	0	0	0
1	0	0	0
0	1	0	0
0	0	1	0
0	0	0	1

构建y的矩阵

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

140
152
135
143
160

混合线性方程组

XpZ <- crossprod(X,Z);XpZ

Chang1	1	1	1	0	0
Chang2	0	0	0	1	1

X’X

XpX <- crossprod(X) ;XpX

	Chang1	Chang2
Chang1	3	0
Chang2	0	2

Z’X

ZpX <- crossprod(Z,X);ZpX

Chang1	Chang2
1	0
1	0
1	0
0	1
0	1

Z’Z

ZpZ <- crossprod(Z);ZpZ

0	0	0	0
1	0	0	0
0	1	0	0
0	0	1	0
0	0	0	1

X’y

Xpy <- crossprod(X,y);Xpy

Chang1	427
Chang2	303

Z’y

Zpy <- crossprod(Z,y);Zpy

140
152
135
143
160

K <- 2;K

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)$Ainv
Ainv

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