专栏首页育种数据分析之放飞自我线性混合模型系列四:矩阵求解

线性混合模型系列四:矩阵求解

混合线性模型,有两大重点,一是估算方差组分,二是矩阵求解。

估算方差组分有很多方法,最常用的是基于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

可以看出,两种矩阵求解方法,结果一致

本文分享自微信公众号 - 育种数据分析之放飞自我(R-breeding),作者:邓飞2013

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-11-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 一步法中混合线性模型方程组构建和控制--blupf90

    http://nce.ads.uga.edu/wiki/lib/exe/fetch.php?media=singlestepblupf90.pdf

    邓飞
  • 软件介绍: BLUPF90的无敌和寂寞

    前同桌说, 要带粉丝了, 一个换一个, 我就把题目改了一下, 原来的题目是《遗传育种软件三剑客之一:BLUPF90》,但我还是喜欢现在这个题目《BLUPF90的...

    邓飞
  • R语言中矩阵常用的操作(笔记)

    发现好久没有更新微信文了, 所谓才思枯竭, 黔驴技穷就是我现在的状态. 记得看过这样一句话: "如果你不知道写什么东西, 那就写不知道写什么事情这件事吧". 深...

    邓飞
  • 数学实验(预习)

    也可以用初等变换求逆矩阵,构造一个n行2n列的矩阵(A E),并进行初等变换,A编程单位矩阵的时候,E就变成了A的逆矩阵.

    云深无际
  • 数据分析与数据挖掘 - 06线性代数

    导数是高等数学中非常重要的知识点,也是人工智能的算法应用中比较常用的一个知识,这一章我们的重点就是讲解一下导数和其求导法则。首先我们来看一下导数的基本概念:函数...

    马一特
  • 吹弹牛皮之Unity 引擎基础 - 矩阵(一)

    沉迷于硬笔的练习偷懒了很长时间。过去的7月份仅仅更新了一篇文章,实在是深表遗憾。接着之前的向量篇小菜继续向下探索。谢谢大家长久来的鼓励和支持。

    用户7698595
  • 一起来学matlab-matlab学习笔记10 10_1一般运算符

    本文为matlab自学笔记的一部分,之所以学习matlab是因为其真的是人工智能无论是神经网络还是智能计算中日常使用的,非常重要的软件。也许最近其带来的一...

    DrawSky
  • 吴恩达机器学习笔记18-逆矩阵、矩阵转置

    “Linear Algebra review(optional)——Inverse and transpose”

    讲编程的高老师
  • 吴恩达机器学习笔记16-矩阵与矩阵的乘法

    “Linear Algebra review(optional)——Matrix-matrix multiplication”

    讲编程的高老师
  • Matlab系列之矩阵秀

    上次讲完了数组的基本操作,不知道是否熟悉使用了,本篇将要对矩阵部分的操作再进行介绍,这部分的内容我觉得蛮有意思的,不过你们觉不觉得我就不知了,但还是想让你们可以...

    狂人V

扫码关注云+社区

领取腾讯云代金券