基于R语言进行AMMI分析2

news2024/9/22 11:37:26

接续上文【基于R语言进行AMMI分析1】

1、AMMI()函数的结果解读

# 加载agricolae包
library(agricolae)
# 加载数据
data(plrv)
# 查看数据
head(plrv)
model<-with(plrv,AMMI(Locality,Genotype,Rep,Yield,PC=TRUE))
# 查看方差分析结果
model$ANOVA
# 查看主成分的方差分析结果
model$analysis
# 查看基因型和环境的均值
model$means
# 查看双标图
model$biplot

（1）model$ANOVA：主效应和互作效应的方差分析

结果显示：ENV（环境）间差异显著；GEN（品种）间差异显著；ENV:GEN（环境与品种互作）显著。

（2）model$analysis：环境与品种交互作用的主成分方差分析

结果显示：环境与品种交互作用可分解为5个主成分，其中通过显著检验的主成分有3个（PC1、PC2、PC3），这三个主成分对交互作用总变异（交互作用平方和）的累计贡献率为92.7%，对数据具有很大的代表性。用于做AMMI2双标图的PC1和PC2的累计贡献率为83.3%，此双标图的反映的品种稳产性和环境辨别力结果的可信度较高。

【注：这里的PC1实际指的是IPCA1，依此类推】

（3）model$means：基因型和环境的均值

（4）model$biplot：用于制作双标图的数据

下图中可以看出：

AMMI1双标图，plot(model,0,1)：0对应Yield列，1对应PC1列。

AMMI2双标图，plot(model,1,2)：1对应的PC1列，2对应PC2列。

（5）示例1

data(CIC)
data1<-CIC$comas[,c(1,6,7,17,18)]
data2<-CIC$oxapampa[,c(1,6,7,19,20)]
head(data1)
head(data2)
cic<-rbind(data1,data2)
model<-with(cic,
            AMMI(Locality,Genotype,Rep,relative))
model$ANOVA

（6）示例2

# 加载agricolae包
library(agricolae)
# 加载数据集
data("sinRepAmmi")
# 展示数据集
head(sinRepAmmi)
# 设置重复数
REP<-3
# 设置MSE
MSerror<-93.24224
# AMMI模型分析
model<-with(sinRepAmmi,AMMI(ENV,GEN,REP,YLD,MSE=MSerror,PC=TRUE))
# 查看方差分析结果
print(model$ANOVA,na.print="")
# 查看AMMI1双标图
plot(model,0,1)
# 查看G×E交互作用的主成分分析结果
pc<-model$PC 
pc$loadings
summary(pc)
biplot(pc)

2、index.AMMI()函数

index.AMMI()函数的使用方法为：index.AMMI(model)，其中model为AMMI()函数的运算结果。

index.AMMI()函数的输出结果为数据框，包括的指标如下：

①ASV：AMMI stability value

②YSI：Yield stability index

③rASV：Rank of AMMI stability value

④rYSI：Rank of yield stability index

⑤means：基因型的均值

计算公式可参考：

结果解释：

ASV是衡量基因型（品种）稳定性的指标，其值越小，基因型（品种）的稳定性越好（【个人评价】此时对IPCA1和IPCA2对G×E交互作用的方差解释率有较高要求，解释率越高，ASV的评价就越准确）；而YSI则是在稳定性评价（ASV从小到大排序的Rank值）的基础上增加了产量数据（means从大到小排序的Rank值）（【个人评价】YSI对品种丰产性和稳产性的评价过于简单，其结果不及GGE双标图）。而根据案例给出的计算结果rYSI实际是YSI计算公式中的RY。

# 加载agricolae包
library(agricolae)
# 加载数据集
data(plrv)
# 查看数据
head(plrv)
# AMMI分析
model<-with(plrv,AMMI(Locality,Genotype,Rep,Yield,console = FALSE))
# 产量稳定性分析
Idx<-index.AMMI(model)
# 查看稳定性分析指标
names(Idx)
# 查看稳定性分析结果
Idx

3、稳定性参数D值

根据本人查看的文献，国内对品种稳定性的评价更多的是使用【基于AMMI模型的品种稳定性分析_张泽】的稳定性参数。解释如下：

当统计检验显著的IPCA轴多于2个以上时，用一般的方法不可能做出像双标图那样直观的坐标图来分析，并且为了更全面、定量地反映GEI（基因型与环境互作）信息，需要计算出一个能度量稳定性的综合指标Di(j)。Di(j)就是在IPCA多维空间中品种(或试点)离原点的距离(即欧氏距离)，计算公式如下：

$D_{i}=\sqrt{\sum^{N}_{r=1}S^2_{ir}}$

$D_{j}=\sqrt{\sum^{N}_{r=1}S^2_{jr}}$

其中，N为显著的IPCA个数， $S^2_{i(j)r}$ 为第i个基因型（或第j个试点）在r个IPCA上的得分。由于Di(j)实际上是在N个IPCA上得分的累加，因此Di值越小、则品种越稳定。对应地，Dj值越大、则表示试点对品种差异的判别力越强。（【个人评价】此方法对G×E互作的主成分进行更多的综合考虑，即所有统计显著的主成分都要求计算在内，主成分对G×E互作的方差解释率更高，评价也更准确）

# 加载agricolae包
library(agricolae)
# 加载数据集
data(plrv)
# 查看数据
head(plrv)
# AMMI分析
# 将console设置为TRUE，可以直接查看主成分分析结果
model<-with(plrv,AMMI(Locality,Genotype,Rep,Yield,console =TRUE))
# 由结果可知PC1、PC2、PC3结果显著
# 查看基因型和环境在各主成分上的得分
model$biplot
# 产量稳定性分析
df<-model$biplot
df<-df[,-c(6,7)]
df$D<-with(df,sqrt(PC1^2+PC2^2+PC3^2))
library(dplyr)
df<-df%>%
  group_by(type) %>%
  arrange(D) %>%
  mutate(rank_D=row_number())
# 结果展示
print(df,n=37)