估算方差组分有很多方法,最常用的是基于REML的方法。
矩阵求解有两种方法,直接法和间接法。
这篇文章通过R语言代码的形式,介绍给定方差组分的情况下,如何根据两种矩阵求解的方法分别计算BLUE值和BLUP值。
混合线性模型
BLUE和BLUP计算公式
dat = data.frame(Herd=factor(c(1,1,2,2,2,3,3,3,3)),
Sire = c("ZA","AD","BB","AD","AD","CC","CC","AD","AD"),
Yield = c(110,100,110,100,100,110,110,100,100))
dat
模型介绍
固定因子矩阵X
X = model.matrix(~Herd-1,data=dat)
X
随机因子矩阵Z
Z = model.matrix(~Sire-1,data=dat)
Z
因为这里没有系谱,关系矩阵为单位矩阵I,这里假定sire的方差组分为0.1, 残差方差组分为1. G = 0.1单位矩阵 R = 1单位矩阵
idg_mat = diag(4)
ide_mat = diag(9)
g = 0.1se = 1G = g*idg_mat
R = se*ide_mat
V = Z %*% G %*% t(Z) + R
G = g*idg_mat
G
R = se*ide_mat
R
V = Z %*% G %*% t(Z) + R
V
y = as.vector(dat$Yield)
y
b = solve(t(X) %*% solve(V) %*% X) %*% t(X) %*% solve(V) %*% y
b
u = G %*% t(Z) %*% solve(V) %*% (y - X %*% b)
u
V矩阵随着数据量的增大,对其进行求解不现实,而混合线性方程组MME,只需要对A的逆矩阵,大大降低了运算量。
计算MME方差的左边矩阵
alpha <- se/g
XpX=crossprod(X) #X’XXpZ=crossprod(X,Z) #X’ZZpX=crossprod(Z,X) #Z’XZpZ=crossprod(Z) #Z’ZXpy=crossprod(X,y) #X’yZpy=crossprod(Z,y) #Z'yLHS=rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ+diag(4)*alpha)) #LHS
LHS
计算MME的右边矩阵
RHS=rbind(Xpy,Zpy) #RHSRHS
sol=solve(LHS)%*%RHS #sol
对比直接矩阵形式计算的结果
# 固定因子效应值b
# 随机因子效应值u
可以看出,两种矩阵求解方法,结果一致