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

发布于

1, 数据

代码语言：javascript
复制
这次使用一个PPT里面的数据, 用R语言演示一下如何做BLUP值计算.

下面是生成数据的代码

代码语言：javascript
复制
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)
代码语言：javascript
复制
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, 计算亲缘关系逆矩阵

代码语言：javascript
复制
library(nadiv)

提取系谱信息

代码语言：javascript
复制
ped <- dat[,2:4]

ped

ID

Sire

Dam

1

0

0

2

0

0

3

1

0

4

1

2

5

3

2

计算亲缘关系逆矩阵

代码语言：javascript
复制
pped = prepPed(ped)


pped
代码语言：javascript
复制
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

代码语言：javascript
复制
Ainv = makeAinv(pped)$Ainv

Ainv
代码语言：javascript
复制
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</p><div class="rno-markdown-code"><div class="rno-markdown-code-toolbar"><div class="rno-markdown-code-toolbar-info"><div class="rno-markdown-code-toolbar-item is-type"><span class="is-m-hidden">代码语言：</span>javascript</div></div><div class="rno-markdown-code-toolbar-opt"><div class="rno-markdown-code-toolbar-copy"><i class="icon-copy"></i><span class="is-m-hidden">复制</span></div></div></div><div class="developer-code-block"><pre class="prism-token token line-numbers language-javascript"><code class="language-javascript" style="margin-left:0">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构建矩阵, 比较方便

代码语言：javascript
复制
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

构建单元矩阵

代码语言：javascript
复制
Z <- diag(length(unique(dat$ID)))

Z

1

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

代码语言：javascript
复制
y <- as.matrix(dat$weight)

y

140

152

135

143

160

混合线性方程组

代码语言：javascript
复制
XpZ <- crossprod(X,Z);XpZ

Chang1

1

1

1

0

0

Chang2

0

0

0

1

1

X’X

代码语言：javascript
复制
XpX <- crossprod(X) ;XpX

Chang1

Chang2

Chang1

3

0

Chang2

0

2

Z’X

代码语言：javascript
复制
ZpX <- crossprod(Z,X);ZpX

Chang1

Chang2

1

0

1

0

1

0

0

1

0

1

Z’Z

代码语言：javascript
复制
ZpZ <- crossprod(Z);ZpZ

1

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

代码语言：javascript
复制
Xpy <- crossprod(X,y);Xpy

Chang1

427

Chang2

303

Z’y

代码语言：javascript
复制
Zpy <- crossprod(Z,y);Zpy

140

152

135

143

160

K

代码语言：javascript
复制
K <- 2;K

2

代码语言：javascript
复制
LHS <- rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ+AinvK))

LHS
代码语言：javascript
复制
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左手矩阵和上图结果一致.

代码语言：javascript
复制
RHS <- rbind(Xpy,Zpy)

RHS

求解BLUP值

代码语言：javascript
复制
solve(LHS)%%RHS
代码语言：javascript
复制
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里面的结果是错误的，这该如何圆场？？？现场翻车记！！！

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

代码语言：javascript
复制
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

#na