GWAS计算BLUE值3--LMM考虑残差异质计算BLUE值

本节，介绍如何使用R语言的asreml包拟合混合线性模型，定义残差异质，计算最佳线性无偏估计（blue）

1. 试验数据

数据来源： Isik F , Holland J , Maltecca C . Genetic Data Analysis for Plant and Animal Breeding. Springer International Publishing, 2017.

数据及代码下载，请关注公众号：育种数据分析之放飞自我，进入知识星球进行相关下载和学习

该数据有62个重组自交系（RIL），在4个地点进行试验，随机区组，每个地点2个重复，每个小区种植20株，随机选择5株的表型平均值作为观测值。

2. 读取数据及转换为因子

library(lme4)
library(emmeans)
library(data.table)
library(tidyverse)
library(asreml)dat = fread("MaizeRILs.csv",data.table = F)
head(dat)
str(dat)col = 1:5
dat[,col] = dat %>% select(all_of(col)) %>% map_df(as.factor)
str(dat)

3. 使用asreml计算BLUE值（定义残差同质）

library(asreml)
m1 = asreml(height ~ RIL, random = ~ location + location:RIL + location:rep,data=dat)
summary(m1)$varcomp
re1 = predict(m1,classify = "RIL")$pval %>% as.data.frame()

4. 使用asreml计算BLUE值（定义残差异质）

m2 = asreml(height ~ RIL, random = ~ location + location:RIL + location:rep,residual = ~ dsum(~units|location),data=dat)
summary(m2)$varcomp

从方差组分可以看到，四个地点的方差组分分别是：

ARC： 45.13
CLR：114.70
PPAC：56
TPAC：54

差别还是比较大的。那这两个模型有没有显著性差异呢，哪个模型最优呢？

5. 比较BIC和似然比检验（LRT）

summary(m1)$bic
summary(m2)$bic
lrt.asreml(m1,m2)

结果可以看出：

定义地点内残差同质的BIC为：2531.222
定义地点内残差异质的BIC为：2530.491
两个模型的LRT的P值<0.001，达到极显著

BIC越小越好。两个模型达到极显著，所以定义残差异质的模型是更好的。

所以，该数据，应该选择地点异质的模型作为计算BLUE值的模型。

6. 计算最优模型的BLUE值

re2 = predict(m2,classify = "RIL")$pval %>% as.data.frame()
head(re2)

7. 更复杂的模型：定义品种与地点互作异质

m3 = asreml(height ~ RIL, random = ~ location + at(location):RIL + location:rep,residual = ~ dsum(~units|location),data=dat)
summary(m3)$varcomp

它和模型2，哪个模型更优呢？

我们可以比较BIC和LRT：

summary(m2)$bic
summary(m3)$bic
lrt.asreml(m2,m3)

结果可以看出：

模型2（只考虑地点残差异质）的BIC为：2530.491
模型3（同时考虑互作的残差异质和地点的残差异质）的BIC为2541.703
两模型达到极显著。

这里模型2更优，并且和模型3达到极显著。所以，我们选择模型2为最优模型。

8. 选择模型不是越复杂越好，而是越合适越好

选择模型不是越复杂越好，而是越合适越好，怎么看合适不合适呢？看一下模型的BIC值。

下一节，我们用教科书的示例，介绍一下联合方差分析的计算方法。其实，从统计角度，很多区试多地点的数据进行一年多点的方差分析，这之前没有进行地点残差一致性检验，是不严谨的。

下一节，我们演示一下，手动计算各个地点的残差和LMM模型定义地点异质，两者是等价的。

数据及代码下载，请关注公众号：育种数据分析之放飞自我，进入知识星球进行相关下载和学习