完结篇 | GWAS计算BLUE值4--联合方差分析演示

GWAS计算BLUE值4--联合方差分析演示 #2021.12.14

本篇，用书籍中的数据和结论，用R语言的一般线性模型和混合线性模型，做一下一年多点的联合方差分析的演示。

1. 参考书籍《农业试验设计与统计分析》

这本书是压箱底的好书，里面是用最原始的方式（公式推导和分解）介绍农业试验设计与统计分析相关的知识点，这次使用第九章 多点试验结果的联合分析，第一节，一年多点随机区组试验结果的分析的数据，作为演示。

2. 什么是联合方差分析

联合方差分析是一个比较内行的、比较封闭的、也可以说是比较专业的称谓，它主要应用于农业多点试验，用来评价品种的平均表现、对环境变化的反应特点以及其适应区域。

它首先是方差分析，但由于它是一年多点的试验，包括品种和环境的互作。问题来了，方差分析的假定是数据正态性、方差齐次性、数据独立性，由于是多点的试验，你如何在没有证明方差齐次性的前提下就进行多点数据的分析呢？

所以联合的前提是地点间误差方差的齐次性，只有齐次了，才可以进行联合方差分析，才叫做联合方差分析。如果不齐次，就要考虑其它方法了，比如去掉误差方差过大的点，比如对数据进行转化，比如用加权联合方差分析，比如用混合线性模型来分析（齐次和非齐次）。

3. 联合方差分析的误差同质性检验

使用Barrlett检验方法，下面两个截图是这种统计方法的介绍。

4. 联合方差分析的手动实现

步骤：

• 1，将数据分割为单地点数据
• 2，对单地点数据分别进行随机区组方差分析，提取残差方差均方
• 3，对各个地点的误差方差进行同质性检验
• 4，如果满足同质，则合并进行联合方差分析。如果不满足，则不能进行联合方差分析。

「读取数据：」

library(data.table)
library(tidyverse)
library(openxlsx)
library(asreml)

dat = read.xlsx("wang9.1-lianhefangchafenxi.xlsx")
head(dat)
str(dat)

col = 1:3
dat[,col] = dat %>% select(all_of(col)) %>% map_df(as.factor)
str(dat)

re1 = dat %>% split(.$loc) %>% map(.,~aov(yield ~ cul + block,.) %>% summary()) 
mse = re1 %>% map(~.[[1]][3,3]) %>% unlist()
mse

地点的残差如上图所示，和书中的结果一致：

计算方差齐性的Bartlett检验的显著性：

s = mean(mse)
sj = sum(log(mse))
X_2 = (9-1)*(4-1)*(10*log(s) - sj)
X_2
1 - pchisq(q = X_2,df = 9)

可以看出，p值为0.07，不显著，说明个地点残差方差齐次，可以进行联合方差分析。

5. 用LMM模型进行方差齐性检验

用混合线性模型，检验方差是否齐性，比较简单，构建两个模型：

• 模型1：方差齐性的模型
• 模型2：方差不齐性的模型
• 比较两个模型的显著性

「模型1：」

mod1 = asreml(yield ~ cul*loc + loc:block, data = dat)
summary(mod1)$varcomp
summary(mod1)$bic

「模型2：」

dd = arrange(dat,loc,cul)
mod2 = asreml(yield ~ cul*loc + loc:block, residual = ~ dsum(~units|loc), data = dd)
summary(mod2)$varcomp

「注意，这里的方差组分，和使用单地点方差分析的残差方差组分是一致的。」

比较两个模型，使用似然比检验：

# 比较两个模型
summary(mod1)$bic
summary(mod2)$bic
lrt.asreml(mod1,mod2,boundary = F)

可以看到，这里的P值为0.075，和我们手动计算的Bartlet检验的p值一致。

6. 计算BLUE值

# 计算BLUE值
re = predict(mod1,"cul")$pval %>% as.data.frame()
head(re)

保存到本地excel文件：

write.xlsx(re,"re_blue.xlsx")

混合线性模型进行一年多点数据分析的优越性介绍

这里看一下专家是怎么说的：

❝[1]王春平, 胡希远, 沈琨仑. 玉米区域试验中误差方差的异质性及其对品种评价的影响[J]. 作物学报, 2013, 39(3):6.❞

如果误差同质模型，就认为e~N（0，σ²），分析的方法采用的是普通最小二乘法（OLSE），这也是传统的方差分析方法。

如果误差是异质的，即认为试验误差随着试验环境的不同而变化，即e~N（0，σ²），e1，e2……不完全相等，这时应用的统计方法为加权最小二乘估计（WLSE）。

「无论误差方差同质还是误差方差异质模型，都是线性混合模型误差的方差协方差的特例，因此两者都可以用混合线性模型的框架进行分析，这时试验的固定效应的估计是最佳线性无偏估计（BLUE）」，当方差同质时，BLUE相当于OLSE，当误差方差异质时，BLUE相当于WLSE。

两模型统一在线性模型原理框架下分析的优点在于，不仅能得到试验效应的最佳线性无偏估计，而且能对试验误差特性与模型的适用性进行评价，以及方便的取得不同误差特性时试验效应差异显著性测验的结果。

「在之前计算条件有限的情况下，联合方差分析不失为作物区试试验分析最佳的选择，但是在误差异质性广泛存在的情况下，之前所做的努力比如对误差点较大的试验点去除、对数据进行转化以及加权联合方差分析都是比较麻烦和不灵活的，而如果对异质性不考虑直接进行方差分析的话，又会对品种和环境互作影响较大，因此推荐用混合线性模型来代替联合方差分析，而且它还可以计算BLUE值，有很大的优越性。」