分类： R语言

library(sommer)
data(DT_example)
DT <- DT_example  # 41个土豆材料，1个地点，3个年份=3个环境
A <- A_example   # 系谱矩阵
# 1. 单地点
# BLUP
ansSingle <- mmer(Yield~1,
                  random= ~ vs(Name, Gu=A),
                  rcov= ~ units,
                  data=DT, verbose = FALSE)
summary(ansSingle)

# 2. 主效应，假设没有基因型与环境互作
# 环境作为固定效应+BLUP
ansMain <- mmer(Yield~Env,
                random= ~ vs(Name, Gu=A),
                rcov= ~ units,
                data=DT, verbose = FALSE)
summary(ansMain)

# 3. 对角线模型
# 假设存在基因型与环境互作
# 假设每个环境基因型的方差不同
# 每个环境单独拟合基因型效应
# 缺点：没有环境间的协方差（每个环境为不同个体）
# 环境作为固定效应+单环境的BLUP（环境的对角线矩阵）
ansDG <- mmer(Yield~Env,
              random= ~ vs(ds(Env),Name, Gu=A),
              rcov= ~ units,
              data=DT, verbose = FALSE)
summary(ansDG)

# 4. 复合对称模型
# 假设存在基因型与环境互作
# 主要的基因型方差协方差分量跨所有环境
# 假设主要的基因型效应被跨环境效应影响
# 缺点：各环境之间有相同的方差-协方差
# 环境作为固定效应+BLUP+GEI
E <- diag(length(unique(DT$Env)))   # 环境的对角线矩阵
rownames(E) <- colnames(E) <- unique(DT$Env)   #  添加环境名称
# 克罗内克乘积，矩阵1的每一项×矩阵2的每一项
# 若矩阵1是m×n，矩阵2是p×q，则克罗内克乘积是mp×nq的方形矩阵
# make.dimnames表示添加行名列名
EA <- kronecker(E,A, make.dimnames = TRUE)
ansCS <- mmer(Yield~Env,
              random= ~ vs(Name, Gu=A) + vs(Env:Name, Gu=EA),
              rcov= ~ units,
              data=DT, verbose = FALSE)
summary(ansCS)

# 5. 复合对称+对角线模型
# 假设存在基因型与环境互作
# 假设每个环境都有不同的基因型与环境互作方差
# 缺点：各环境有相同的协方差
# 环境作为固定效应+BLUP+每个地点的基因型与环境互作效应
ansMain <- mmer(Yield~Env,
                random= ~ vs(Name, Gu=A) + vs(ds(Env),Name, Gu=A),
                rcov= ~ units,
                data=DT, verbose = FALSE)
summary(ansMain)

# 6. 非结构化模型
# 假设存在基因型与环境互作
# 非结构化模型可以得到尽可能多的环境到环境组合的协方差
# 缺点：方差分量很大，方差组分很难收敛，一些方差协方差组分是零，且很难有良好的初始值
# 环境作为固定效应+单环境BLUP（协方差调整）
ansUS <- mmer(Yield~Env,
              random= ~ vs(us(Env),Name, Gu=A),
              rcov= ~ units,
              data=DT, verbose = FALSE)
summary(ansUS)

# 7. 随机回归模型
# 假设环境可以视为连续变量，因此，截距和斜率的方差组分可以被拟合
# 方差分量的数量将依赖于Legendre多项式拟合的阶数(leg函数)
GBIT.library("orthopolynom")
DT$EnvN <- as.numeric(as.factor(DT$Env))
ansRR <- mmer(Yield~Env,
              random= ~ vs(leg(EnvN,1),Name),
              rcov= ~ units,
              data=DT, verbose = FALSE)
summary(ansRR)
# 此外，其他结构的方差协方差矩阵可以应用于这个多项式
GBIT.library("orthopolynom")
DT$EnvN <- as.numeric(as.factor(DT$Env))
ansRR <- mmer(Yield~Env,
              random= ~ vs(us(leg(EnvN,1)),Name),
              rcov= ~ units,
              data=DT, verbose = FALSE)
summary(ansRR)

# 8. 基因型环境互作协方差结构
# 不常用，但1阶自回归和其他协方差结构可以应用于GEI模型
# 
E <- AR1(DT$Env) # can be AR1() or CS(), etc.
rownames(E) <- colnames(E) <- unique(DT$Env)
EA <- kronecker(E,A, make.dimnames = TRUE)
ansCS <- mmer(Yield~Env,
              random= ~ vs(Name, Gu=A) + vs(Env:Name, Gu=EA),
              rcov= ~ units,
              data=DT, verbose = FALSE)
summary(ansCS)