如何在R语言中构建动物模型混合线性模型方程并计算BLUP值
数据
参考书籍：王金玉, 陈国宏. 《数量遗传与动物育种》. 东南大学出版社, 2004，P200，第四节：单性状BLUP育种值估计
dat <- data.frame(id=c(4,5,6),sire = c(1,3,3),dam=c(2,2,4),y=c(200,170,180))
dat
dat
A的逆矩阵
library(asreml)
ainv <- asreml.Ainverse(dat[,1:3])$ginv
ainv_mat <- asreml.sparse2mat(ainv)
row.names(ainv_mat) = colnames(ainv_mat) <- attr(ainv,"rowNames")
ainv_mat
Pedigree insert: Individual "1" inserted at record 1
Pedigree insert: Individual "3" inserted at record 2
Pedigree insert: Individual "2" inserted at record 3
A
Y 矩阵形式
Y
X 矩阵形式
X
Z矩阵
Z
混合线性模型矩阵形式
HMM
其中，左边的矩阵为LHS
LHS
右边的矩阵为RHS
RHS
R 语言实现代码
X <- matrix(1,3,1)
X
X
Z = matrix(c(0,0,0,1,0,0,
0,0,0,0,1,0,
0,0,0,0,0,1),3,byrow = T);Z
Z
H2
LHS的计算方法
tXX = t(X)%*%X
tXZ = t(X)%*%Z
tZX = t(Z)%*%X
tZZk = t(Z)%*%Z + ainv_mat*2
LHS = rbind(cbind(tXX,tXZ),cbind(tZX,tZZk))
LHS
LHS
RHS的计算方法
Y <- matrix(c(200,170,180),3,1)
tXY = t(X)%*%Y
tZY = t(Z)%*%Y
RHS = rbind(tXY,tZY)
RHS
RHS
参数的计算方法
f
ab = (solve(LHS))%*%RHS
ab
blup
用ASreml运算作为对照，结果一样
dat$id <- as.factor(dat$id)
moda <- asreml(y ~ 1, random = ~ped(id),ginverse=list(id=ainv),data=dat,
start.values = T)
tt <- moda$gammas.table
tt
varcomp
# 遗传力为1/3，这里加性假定为1，残差假定为2
tt$Value <- c(1,2)
tt$Constraint <- "F"
modb <- asreml(y ~ 1, random = ~ped(id),ginverse=list(id=ainv),data=dat,
R.param = tt,G.param = tt)
ASReml: Tue Dec 19 16:34:54 2017
LogLik S2 DF wall cpu
-87.6387 2.0000 2 16:34:54 0.0
-87.6387 2.0000 2 16:34:54 0.0
-87.6387 2.0000 2 16:34:54 0.0
-87.6387 2.0000 2 16:34:54 0.0
Finished on: Tue Dec 19 16:34:54 2017
LogLikelihood Converged
coef(modb)$fixed
coef(modb)$random
asreml-blup
ab
r-blup
结果是一致的
r-breeding.png
