线性混合模型系列四：矩阵求解
- 2019 年 11 月 22 日
- 筆記
混合线性模型，有两大重点，一是估算方差组分，二是矩阵求解。
估算方差组分有很多方法，最常用的是基于REML的方法。
矩阵求解有两种方法，直接法和间接法。
这篇文章通过R语言代码的形式，介绍给定方差组分的情况下，如何根据两种矩阵求解的方法分别计算BLUE值和BLUP值。
1. 混合模型矩阵求解
混合线性模型
BLUE和BLUP计算公式
2. 举个栗子
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))
2.1 数据
dat
2.2 模型介绍
模型介绍
- 固定因子：Herd
- 随机因子：Sire
- 观测值：Yield
2.3 固定因子矩阵X和随机因子Z
固定因子矩阵X
X = model.matrix(~Herd-1,data=dat) X
随机因子矩阵Z
Z = model.matrix(~Sire-1,data=dat) Z
2.4 V矩阵构建
因为这里没有系谱，关系矩阵为单位矩阵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
2.5 G矩阵
G = g*idg_mat
G
2.6 R矩阵
R = se*ide_mat
R
2.7 V矩阵
V = Z %*% G %*% t(Z) + R
V
2.8 y矩阵
y = as.vector(dat$Yield)
y
3. 固定因子b
b = solve(t(X) %*% solve(V) %*% X) %*% t(X) %*% solve(V) %*% y b
4. 随机因子u
u = G %*% t(Z) %*% solve(V) %*% (y - X %*% b) u
5. MME 混合线性方程组求解
V矩阵随着数据量的增大，对其进行求解不现实，而混合线性方程组MME，只需要对A的逆矩阵，大大降低了运算量。
5.1 等式左边计算
计算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
5.2 等号右边计算
计算MME的右边矩阵
RHS=rbind(Xpy,Zpy) #RHSRHS
5.3 求解b值和u值
sol=solve(LHS)%*%RHS #sol
对比直接矩阵形式计算的结果
# 固定因子效应值b
# 随机因子效应值u
可以看出，两种矩阵求解方法，结果一致
