N3 | 动物模型计算BLUP值

这次使用一个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

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"

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

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

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

RHS <- rbind(Xpy,Zpy)
RHS

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

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