如何在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

结果是一致的

wechat