線性混合模型系列四：矩陣求解
- 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
可以看出，兩種矩陣求解方法，結果一致
