原因之一是，动物的建模群体规模通常超过1000，而在植物中，1000个个体组成的建模群体在实际育种过程中可遇而不可求。另一个重要的原因是，植物比动物有更显著的基因型与环境互作效应，会对预测造成极大的干扰。
One of the reasons is that the training population size of animals usually exceeds 1000, while in plants, a training population composed of 1000 individuals can be difficult in the actual breeding process. Another important reason is that plants have a more significant genotype environment interaction (GEI) effect than animals, which can cause significant interference in predictions.
原文：http://potatobreeding.cals.wisc.edu/wp-content/uploads/sites/161/2014/01/GS_tutorial.pdf
作者：Jeffrey Endelman
版本：4
更新：20130615
翻译：张敖
翻译更新：20221026 12:44:13
本文介绍如何使用第四版中新加入rrBLUP的特性(Endelman 2011)。
本包的基本核心仍然是mixed.solve，它求解了除残差之外的一个方差分量的混合模型。该函数估计两个方差分量，通过ML或REML模型估计，估计的执行使用Kang等(2008)描述的光谱分析算法。在这个过程中，很容易创建表型协方差矩阵的逆矩阵，逆矩阵随后用于固定效应和随机效应的BLUE和BLUP求解计算(Searle et al. 1992)。在(Endelman 2011)中，作者展示了mixed.solve如何用于全基因组预测，要么建模标记作为随机效应，要么家系随机效应。
第4版是围绕A.mat和kin.blup函数设计的。
A.mat用标记估计真实的亲缘关系矩阵（A），对于没有缺失数据的高密度标记，A.mat的VanRaden(2008)建议的第一个公式：
这里，矩阵W中等位基因的表示数值通过群体均值进行中心化。Endelman and Jannink (2012)证明了该公式的平均对角元素为1+f，f是近交系数。
缺失标记数据
当基因型有缺失时，A.mat有两个补缺失的选项。一个是缺失值用标记的群体均值代替，这对SNP芯片类，只有少量缺失值是足够的。对GBS标记，缺失值的水平可能过高。这种情况下，A.mat可以使用基于多元正态分布的EM算法估计亲缘关系矩阵(Poland et al. 2012)。
为了证明EM算法，我下载了Poland et al. (2012) 的GBS数据https://www.crops.org/publications/tpg/supplements/5/tpg12-06-0006-dataset-s2.gz（链接已失效）。
下列代码会读取GBS数据，并转换为rrBLUP需要的{-1，0，1}格式。
GBS<-read.csv("tpg12-06-0006-dataset-s2",header=T,as.is=TRUE,row.names=1)
alleles <- setdiff(unique.x,union("H","N")){
unique.x <- unique(x)
y <- rep(0,length(x))
y[which(x==alleles[1])] <- -1
y[which(x==alleles[2])] <- 1
y[which(x=="N")] <- NA
return(y)
}
X <- apply(GBS[,-c(1:3)],1,parse.GBS)
dim(X) #lines by markers
frac.missing <- apply(X,2,function(z){length(which(is.na(z)))/length(z)})
length(which(frac.missing<0.5))
hist(frac.missing)
有1.6万的标记缺失率小于50%，足够估计出254个家系的关系矩阵。因为作者是在一个具有多个处理器的unix兼容的系统上运行这段代码的，所以我可以使用12个核来加速EM算法：
library(rrBLUP)
system.time(A1 <- A.mat(X,impute.method="EM",n.core=12,max.missing=0.5)) EM算法在其进行过程中显示收敛序列，它表示亲缘关系系数的根均方误差。默认停止标准为0.02，但是可以通过tol参数改变（输入?A.mat获得更多信息）。在本例中，当它达到0.0151时，计算停止，这仅仅需要1分钟的时间。如果只有一个核，可以简单的省略n.core选项，因为默认为1核。译者注：windows下指定核数量无效。
用均值进行补缺失，可以使用下列语法：
system.time(A2 <- A.mat(X, max.missing=0.5))用均值补缺失肯定会更快，在很多情况下，在GEBV的预测精度上，均值的表现与EM算法和其他更先进的方法一样好。然而，与EM方法相比，均值的方法的育种值往往更有偏向性(Poland et al. 2012)。对两个矩阵的平均对角线元素比较表明，EM算法的结果更加接近给定1%杂合率的期望，表达式1+f≈2。提出
round(mean(diag(A2)),2) # imputed with mean
round(mean(dia(A1)),2) # imputed with EMA矩阵的收缩估计
A.mat的另一个特性是收缩估计，它可以用于低密度的标记，例如来自384个SNP的芯片。当家系的数量与标记的数量相当或更多时，上面的方程可能是最优化的亲缘关系矩阵估计。Endelman and Jannink (2012)建议将估算值降低到(1+f)I，用收缩强度选择最小化均方误差。
为了说明这一点，我将使用BLR包中的一个数据集，该数据集由1279个DArT标记对599个小麦家系进行了基因分型。
library(BLR)
data(wheat)
M <- 2*X-1 #convert markers to {-1,1}
dim(M) #=[1] 599 1279
A1 <- A.mat(M,shrink=TRUE) #= [1] "Shrinkage intensity: 0.03"
A2 <- A.mat(M[,sample(1:1279,384)],shrink=TRUE) #= "Shrinkage intensity: 0.1"
A3 <- A.mat(M[,sample(1:1279,192)],shrink=TRUE) #= "Shrinkage intensity: 0.17"如上例所示，收缩强度从0（无收缩）到1（完全收缩），标记密度减少，缩减强度增加。所有1279个标记，使用了非常小的缩减（3%），而384和192个标记的随机集合，缩减强度是10%和17%。
Endelman(2001)中，我介绍了一个mixed.solve的“包装器”，叫做kinship.BLUP，这是为家系的预测而设计的，使用加性遗传模型或高斯核。然而，“包装器”没有那么方便，所以我又设计了一个新的函数代替它，这个函数叫做kin.blup。该函数不需要用户创建设计矩阵，新函数会自动从数据框中完成创建部分。另一个不同是，用户用kin.blup传递的是亲缘关系矩阵，而不是标记，这样允许我们更好的计算A矩阵（参见A.mat）。
为了说明基础功能，我将继续以上述的小麦数据为例，其中，A1矩阵用所有的标记估计。这599个自交系已经在BLR包中分成了10个集合用于交叉验证。为了预测集合1的育种值，使用集合2至10的自交系表型，首先，表型和相应的基因型标识符必须组装成数据框。
test <- which(sets==1)
yNA <- Y[,1] # grain yield in environment 1
yNA[test] <- NA # mask yields for validation set
data1 <- data.frame(y=yNA, gid=1:599)验证集合需要遮盖表型数据，需要将他们的值设成NA，如上面的例子。最小的数据集有两列，一列是表型值，一列是基因型的标签（注：基因型名称）。确保数据框中基因型标签（名称）与关系矩阵的行名对应，我们可以做全基因组预测：
rownames(A1) <- 1:599
ans1 <- kin.blup(data1, K=A1, geno="gid", pheno="y")
str(ans1) 正如上面的例子，亲缘关系矩阵传入参数“K”（给kinship）到kin.blup函数，以及数据框中的表型和基因型变量名。该函数返回方差组分（$Vg, $Ve）的REML估计，以及遗传值（$g）的BLUP，即本例中的育种值。同时，还将返回来自混合模型的残差。如上例所示，BLUP值在亲缘关系矩阵中返回599个条目（观测值），尽管这些系中10%的系没有表型（由于遮盖）。BLUP的顺序由K矩阵的顺序决定（array也如此命名）。
为了评估GEBV的预测精度，计算验证数据集的预测值与遮盖的表型值的相关性：
round(cor(ans1$g[test),Y[test,1],2)对于具有高斯核的预测，而不是传递关系矩阵到K参数，使用欧式距离，并将GAUSS参数者只为TRUE。
D <- as.matrix(dist(M)) # Euclidean distance
system.time(ans2 <- kin.blup(data1, K=D, GAUSS=TRUE, geno="gid", pheno="y"))
system.time(ans2 <- kin.blup(data1, K=D, GAUSS=TRUE, geno="gid", pheno="y", n.core=10))
round(cor(ans$g[test],Y[test,1],2) 正如所见，多核可以提高GAUSS核的运算速度。高斯核的预测精度是0.63，高于使用亲缘关系矩阵。该结果与Endelman（2011）的结果不完全相同，用于确定最优尺度网络点（grid point）参数的并不相同；输入?kin.blup获得设置网络点的更多信息。
多环境试验
kin.blup函数能够处理不同环境中的重复测量。来自BLR包的小麦数据集平均在4个环境下测定599个自交系。环境2-4是相关的，本例中，将考虑为育种计划的一个目标环境。下面的代码将创建一个不平衡的数据集，使用跨环境的部分重复。
y <- c(Y[1:400,2],Y[101,400,3],Y[201:500,4])
env <- c(rep(2,400),rep(3,300),rep(4,300))
gid <- c(1:400,101:400,201:500)
data2 <- data.frame(y=y,env=env,gid=gid)
nrow(data2)为了在预测模型中将环境效应作为固定效应，数据框的列被传递给函数：
system.time(ans <- kin.blup(data2,K=A1,geno="gid",pheno="y",fixed="env"))
round(cor(ans$g[501:599],rowMeans(Y[501:599,2:4])),3) 当有多个固定效应用于建模时，例如年份和地点，只需要传递一个列名数组，例如fixed=c(“year”,”location”)。
上面的例子是一步预测，这比首先计算线平均值的两步方法的计算要求更高。对于不平衡数据，kin.blup可以用两步法的速度做预测，同时保留关于不同重复水平的信息。这通过reduce=TRUE实现，转换混合模型的维度等于自交系的数量（参看使用手册获得更多细节）：
system.time(ans2<-kin.blup(data2,K=A1,geno="gid",pheno="y",fixed="env",reduce=TRUE))
round(cor(ans2$g[501:599],ans$g[501:599]),3)
round(cor(ans2$g[501:599],rowMeans(Y[501:599,2:4])),3)在上面的例子中，采用约简方法的计算时间减少了近5倍。两种方法的预测非常相似（r = 0.96），但在本例中没有降低。
Endelman, J.B. 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4:250–255. doi:10.3835/plantgenome2011.08.0024
Endelman, J.B., and J.-L. Jannink. 2012. Shrinkage estimation of the realized relationship matrix. G3:Genes, Genomes, Genetics. 2:1405-1413. doi:10.1534/g3.112.004259
Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709–1723.
Pérez et al. 2010. Genomic-enabled prediction based on molecular markers and pedigree using the Bayesian Linear Regression package in R. Plant Genome 3:106–116.
Poland, J., J. Endelman et al. 2012. Genomic selection in wheat breeding using genotyping-by-sequencing. Plant Genome 5:103–113. doi: 10.3835/plantgenome2012.06.0006．
Searle et al. 1992. Variance Components. John Wiley & Sons, Hoboken.
VanRaden, P.M. 2008. Efficient methods to compute genomic predictions. J. Dairy Science
91:4414–4423.
原版更新：2018.06.08 翻译更新：2022.05.11
翻译 刘翼宁、张敖
题目 使用深度学习从表型预测基因型
版本 1.2（注：最新版本为1.2）
作者 Chuang Ma, Zhixu Qiu, Qian Cheng and Wenlong Ma
维护 Zhixu Qiu zhixuqiu2015@gmail.com, Chuang Ma
描述 名为‘DeepGS’的R包可以执行全基因组选择（Genomic selection，GS），GS在动植物育种中前途光明。 DeepGS使用整个基因组的全部分子标记结合先进的机器学习技术（深度学习）预测表型。DeepGS用国际玉米小麦改良中心（CIMMYT）基因库的2000个伊朗面包小麦的8个表型预测完成了有效性验证。此外，基于粒子群优化（ELBPSO）的集成学习，可以用于不同GS模型的预测的线性整合。
依赖 R(>=3.0.0)
建议 mxnet
许可 GPL-2 | GPL-3
编码 UTF-8
集成数据 true
RoxygenNote 6.0.1
需要编译 no
‘DeepGS’包的中文说明R主题文档cvSampleIndexELBPSOmeanNDCGpredict_GSModeltrain_deepGSModelwheat_example
生成建模集和测试集的索引。
描述
该函数为建模和测试集中的样本生成索引，以执行N倍交叉验证实验。
使用
cvSampleIndex(sampleNum, cross = 5, seed = 1, randomSeed = FALSE)
参数
| 参数 | 描述 |
|---|---|
| sampleNum | 需要划分为建模集和测试集的样本数。 |
| cross | 交叉验证的倍数（折叠数）。 |
| seed | 一个整数种子，用于重复结果，默认为1。 |
| randomSeed | 逻辑变量。默认值为FALSE。 |
值
每个元素的列表，包括$trainIdx、$testIdx和$cvIdx。
$trainIdx 建模样本的索引。
$testIdx 测试样本的索引。
$cvIdx 交叉验证的索引。
| 值 | 描述 |
|---|---|
| $trainIdx | 建模样本的索引。 |
| $testIdx | 测试样本的索引。 |
| $cvIdx | 交叉验证的索引。 |
作者
Chuang Ma, Zhixu Qiu, Qian Cheng and Wenlong Ma
例子
#' ## 读取数据 ##
data(wheat_example)
## 5倍交叉验证
b <- cvSampleIndex(sampleNum = 2000, cross = 5, seed = 1)
## 提取其中一次验证的数据集
b$cv1
基于粒子群优化的集成学习（Ensemble Learning Based on Particle Swarm Optimization）。
描述
该算法结合了不同GS模型的预测表型值，并返回这些模型的最佳权重。
使用
ELBPSO(rep_times = 100, interation_times, weight_dimension, weight_min, weight_max, rate_min, rate_max, paticle_number, pred_matrix, IW = 1, AF1 = 2, AF2 = 2)
参数
| 参数 | 描述 |
|---|---|
| rep_times | ELBPSO重复次数，默认是100。 |
| interation_times | 每次ELBPSO优化迭代次数。 |
| weight_dimension | 权重维度。 |
| weight_min | 最小权重。 |
| weight_max | 最大权重。 |
| rate_min | 最小更新率。 |
| rate_max | 最大更新率。 |
| paticle_number | 帕特里克数量。 |
| pred_matrix | 第一列代表真实值，其他列代表其他GS模型的预测值。 |
| IW | 惯性（惰性）权重。默认为1。 |
| AF1 | 加速因子1。默认是2。 |
| AF2 | 加速因子2。默认是2。 |
例子1.2
# # Not run
# library(DeepGS)
# library(rrBLUP)
# data("wheat_example")
# Markers <- wheat_example$Markers
# y <- wheat_example$y
# cvSampleList <- cvSampleIndex(length(y),10,1)
# # select one fold
# cvIdx <- 1
# trainIdx <- cvSampleList[[cvIdx]]$trainIdx
# testIdx <- cvSampleList[[cvIdx]]$testIdx
# trainMat = Markers[trainIdx,]
# trainPheno = y[trainIdx]
# validIdx <- sample(1:length(trainIdx),floor(length(trainIdx)*0.1))
# validMat <- trainMat[validIdx,]
# validPheno <- trainPheno[validIdx]
# testMat = Markers[testIdx,]
# testPheno = y[testIdx]
# # design DeepGS architecture
# conv_kernel <- c("1*18") # convolution kernels (fileter shape)
# conv_stride <- c("1*1")
# conv_num_filter <- c(8) # number of filters
# pool_act_type <- c("relu") # active function for next pool
# pool_type <- c("max") # max pooling shape
# pool_kernel <- c("1*4") # pooling shape
# pool_stride <- c("1*4") # number of pool kernerls
# fullayer_num_hidden <- c(32,1)
# fullayer_act_type <- c("sigmoid")
# drop_float <- c(0.2,0.1,0.05)
# cnnFrame <- list(conv_kernel =conv_kernel,conv_num_filter = conv_num_filter,
# conv_stride = conv_stride,pool_act_type = pool_act_type,
# pool_type = pool_type,pool_kernel =pool_kernel,
# pool_stride = pool_stride,fullayer_num_hidden= fullayer_num_hidden,
# fullayer_act_type = fullayer_act_type,drop_float = drop_float)
# markerImage = paste0("1*",ncol(trainMat))
# # train DeepGS model
# DeepGS_obj <- train_deepGSModel(trainMat = trainMat,trainPheno = trainPheno,
# validMat = validMat,validPheno = validPheno, markerImage = markerImage,
# cnnFrame = cnnFrame,device_type = "cpu",gpuNum = 1, eval_metric = "mae",
# num_round = 6000,array_batch_size= 30,learning_rate = 0.01,
# momentum = 0.5,wd = 0.00001, randomseeds = 0,initializer_idx = 0.01,
# verbose =TRUE)
# # make predictions based on the trained model
# DeepGS_pred <- predict_GSModel(GSModel = DeepGS_obj,testMat = Markers[testIdx,],
# markerImage = markerImage )
# # train RR-BLUP model#'
# rrBLUP_obj <-mixed.solve(trainPheno, Z=trainMat, K=NULL, SE = FALSE, return.Hinv=FALSE)
# # make predictions based on the trained model
# rrBLUP_pred <- testMat %*% rrBLUP_obj$u + as.numeric(rrBLUP_obj$beta )
# # prepare the prediction matrix
# test_predMat <- cbind(t(DeepGS_pred), rrBLUP_pred)
# train_predMat <- cbind(testPheno, t(DeepGS_pred), rrBLUP_pred)
# colnames(train_predMat) <- c("real", "DeepGS", "RR-BLUP")
## End not run
例子1.1
## 不运行
# rrBLUP例子
# library(DeepGS)
# library(rrBLUP)
# data("wheat_example")
# Markers <- wheat_example$Markers
# y <- wheat_example$y
# cvSampleList <- cvSampleIndex(length(y),10,1)
## 选择第N倍（折）
# cvIdx <- 1
# trainIdx <- cvSampleList[[cvIdx]]$trainIdx
# testIdx <- cvSampleList[[cvIdx]]$testIdx
# trainMat = Markers[trainIdx,]
# trainPheno = y[trainIdx]
# testMat = Markers[testIdx,]
# testPheno = y[testIdx]
# rrBLUP_obj <-mixed.solve(trainPheno, Z=trainMat, K=NULL, SE = FALSE, return.Hinv=FALSE)
# rrBLUP_pred <- testMat %*% rrBLUP_obj$u + as.numeric(rrBLUP_obj$beta)
## 不运行结束
# 使用他们的预测得分计算不同模型的权重
test_datapath <- system.file("exdata", "test_ELBPSO.RData", package = "DeepGS")
load(test_datapath)
weight <- ELBPSO(rep_times = 100,interation_times = 25,weight_dimension = 2, weight_min = 0,weight_max=1,rate_min = -0.01,rate_max = 0.01, paticle_number = 10, pred_matrix = train_predMat,IW = 1, AF1 = 2, AF2 = 2)
new_pre <- (test_predMat %*% weight)/sum(weight)
# calculating the weight of different training model by using their predict socres
test_datapath <- system.file("exdata", "test_ELBPSO.RData",
package = "DeepGS")
load(test_datapath)
weight <- ELBPSO(rep_times = 100,interation_times = 25,weight_dimension = 2, weight_min = 0,weight_max=1,rate_min = -0.01,rate_max = 0.01, paticle_number = 10, pred_matrix = train_predMat,IW = 1, AF1 = 2, AF2 = 2)
ensemble_pred <- (test_predMat %*% weight)/sum(weight)
predMat <- cbind(testPheno, t(DeepGS_pred), rrBLUP_pred, ensemble_pred)
colnames(predMat) <- c("real", "DeepGS", "RR-BLUP", "ensemble")
cor(predMat)
计算平均归一化折扣累积增益。
描述
该函数是在选择育种值高的前k个个体时，用平均归一化折扣累积增益评估全基因组选择预测模型表现的。
使用
meanNDCG(realScores, predScores, topAlpha = c(10))
参数
| 参数 | 描述 |
|---|---|
| realScores | 一个数值化向量，一个性状验证个体的育种值。 |
| predScores | 一个数值化向量或矩阵，用全基因组选择预测模型预测的个体的预测育种值。 |
| topAlpha | 一个数值化的向量（一个或多个参数），优良个体百分比，默认是10。 |
作者
Chuang Ma, Zhixu Qiu, Qian Cheng and Wenlong Ma
例子
## Not run
refer_value <- runif(100)
pred_value <- sin(refer_value) + cos(refer_value)
meanNDCG(realScores = refer_value,predScores = pred_value, topAlpha = c(10,20,30))
## End not run
预测GS模型。
描述
使用经过训练的深度学习全基因组选择预测模型预测性状值。
使用
predict_GSModel(GSModel, testMat, markerImage)
参数
GSModel 从DeepGSModel函数获得的经过训练的预测模型。
testMat 基因型矩阵（T×M；T是个体，M是标记）。
markerImage （字符串）给出“i×j”的图像格式，每个个体的标记信息（M×1）将被编码。
作者
Chuang Ma, Qian Cheng, Zhixu Qiu and Wenlong Ma
利用深度学习技术建立全基因组选择预测模型。
描述
该函数应用深度卷积神经网络建立全基因组选择的预测模型。
使用
train_deepGSModel(trainMat, trainPheno, validMat, validPheno, markerImage, cnnFrame, device_type = "cpu", gpuNum = "max", eval_metric = "mae", num_round = 6000, array_batch_size = 30, learning_rate = 0.01, momentum = 0.5, wd = 1e-05, randomseeds = NULL, initializer_idx = 0.01, verbose = TRUE...)
参数
trainMat 建模群体的基因型矩阵（N×M，N是个体，M是标记）。
trainPheno 建模群体表型的向量（N×1）。
validMat 验证群体的基因型矩阵。
validPheno 验证群体的表型向量（N*1）。
markerImage （字符串）这个给出“i×j”的图像格式，每个个体用标记信息编码（M×1）。如果图像大小超过SNP标记数，用0修改缺失的部分，如果图像大小小于SNP标记数，多余的标记被去除。
cnnFrame 包含卷积神经网络（CNN）框架的以下元素的列表：
| trainMat | 建模群体的基因型矩阵（N×M，N是个体，M是标记）。 |
|---|---|
| trainPheno | 建模群体表型的向量（N×1）。 |
| validMat | 验证群体的基因型矩阵。 |
| validPheno | 验证群体的表型向量（N*1）。 |
| markerImage | （字符串）这个给出“i×j”的图像格式，每个个体用标记信息编码（M×1）。如果图像大小超过SNP标记数，用0修改缺失的部分，如果图像大小小于SNP标记数，多余的标记被去除。 |
| cnnFrame | 包含卷积神经网络（CNN）框架的以下元素的列表： – conv_kernel：向量（K×1）给出卷积核大小（宽度x高度），分别用于过滤K个卷积层的图像矩阵。 – conv_num_filter：向量（K×1）分别给出K个卷积层的卷积核数。 – pool_act_type：给出活动函数的类型的向量（K×1），将定义K个卷积层的输出，这将分别作为相应池层的输入。它包括“relu”、“sigmoid”、“softrelu”和“tanh”。 – conv_stride：K个卷积核的一个特征（K×1）步长。 – pool_type：分别从“avg”、“max”、“sum”中选择的一个特征（K×1）类型的K个池层。 – pool_kernel：一个特征（K×1）K个池层的K池核大小（宽×高）。 – pool_stride：一个特征（K×1）K池核的步长。 – fullayer_number-hidden：一个数值化（H×1），H完全连接层分别隐藏神经元的数量，最后一个完全连接层的隐藏神经元数量必须是1。 – fullayer_act_type：一个数值型（（H-1）×1），从所有连接层的“relu”，“sigmoid”，“softrelu”和“tanh”的选择活动函数的类别。 – drop_float：数值型。 |
| device_type | 选择”cpu“或”gpu“设备去构建预测模型。 |
| gpuNum | （数值型）GPU设备的数量，如果使用多GPU（gpuNum>1），momentum参数必须大于0。 |
| eval_metric | （字符串）一种估计建模过程表现的方法，包括”mae“、”rmse“和”accuracy“，默认为”mae“。 |
| num_round | （数值型）超过建模数据去建模的迭代次数，默认=10。 |
| array_batch_size | （数值型）它为每个更新权重定义将通过网络传播的样本数，默认为128。 |
| learning_rate | 建模过程的学习率。 |
| momentum | （浮点，0~1）移动平均线动量，平均为0.9。 |
| wd | （浮点，0~1）权重衰减，默认为0。 |
| randomseeds | 设置mxnet设备特定随机数使用的种子。 |
| initializer_idx | 参数的初始化主题。 |
| verbose | 逻辑值（默认true）指定是否在建模期间的迭代时输出信息。 |
| … | 用于构造包“mxnet”中使用的神经网络的参数（http://mxnet.io/）。 |
作者
Chuang Ma , Zhixu Qiu, Qian Cheng and Wenlong Ma
例子
data(wheat_example)
Markers <- wheat_example$Markers
y <- wheat_example$y
cvSampleList <- cvSampleIndex(length(y),10,1)
# cross validation set
cvIdx <- 1
trainIdx <- cvSampleList[[cvIdx]]$trainIdx
testIdx <- cvSampleList[[cvIdx]]$testIdx
trainMat <- Markers[trainIdx,]
trainPheno <- y[trainIdx]
validIdx <- sample(1:length(trainIdx),floor(length(trainIdx)*0.1))
validMat <- trainMat[validIdx,]
validPheno <- trainPheno[validIdx]
trainMat <- trainMat[-validIdx,]
trainPheno <- trainPheno[-validIdx]
conv_kernel <- c("1*18") ## convolution kernels (fileter shape)
conv_stride <- c("1*1")
conv_num_filter <- c(8) ## number of filters
pool_act_type <- c("relu") ## active function for next pool
pool_type <- c("max") ## max pooling shape
pool_kernel <- c("1*4") ## pooling shape
pool_stride <- c("1*4") ## number of pool kernerls
fullayer_num_hidden <- c(32,1)
fullayer_act_type <- c("sigmoid")
drop_float <- c(0.2,0.1,0.05)
cnnFrame <- list(conv_kernel =conv_kernel,conv_num_filter = conv_num_filter,
conv_stride = conv_stride,pool_act_type = pool_act_type,
pool_type = pool_type,pool_kernel =pool_kernel,
pool_stride = pool_stride,fullayer_num_hidden= fullayer_num_hidden,
fullayer_act_type = fullayer_act_type,drop_float = drop_float)
markerImage = paste0("1*",ncol(trainMat))
trainGSmodel <- train_deepGSModel(trainMat = trainMat,trainPheno = trainPheno, validMat = validMat,validPheno = validPheno, markerImage = markerImage, cnnFrame = cnnFrame,device_type = "cpu",gpuNum = 1, eval_metric = "mae", num_round = 6000,array_batch_size= 30,learning_rate = 0.01,momentum = 0.5,wd = 0.00001, randomseeds = 0,initializer_idx = 0.01, verbose =TRUE)
predscores <- predict_GSModel(GSModel = trainGSmodel,testMat = Markers[testIdx,], markerImage = markerImage )
在开发功能中运行的例子。
描述
列表：
- 标记：矩阵（599×1225），每行代表一个个体的1225个标记。
- y：每个个体的真实表型值。
使用
data(wheat_example)
模型：
y = Xβ + Zu + ε代码：
# sommer的GS预测
library(sommer) # 加载sommer包，第一次安装请先运行install.packages("sommer")
data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
colnames(DT) <- paste0("X",1:ncol(DT)) # 列名为X1~X4
DT <- as.data.frame(DT) # 转换为数据框
DT$id <- as.factor(rownames(DT)) # 把行号作为列并作为因子
# select environment 1
rownames(GT) <- rownames(DT) # 给基因型赋值相同的行号（名称）
K <- A.mat(GT) # additive relationship matrix
colnames(K) <- rownames(K) <- rownames(DT) # K的行列相同，都赋值名称
# GBLUP pedigree-based approach
set.seed(12345)
y.trn <- DT
vv <- sample(rownames(DT),round(nrow(DT)/5)) # 抽取20%作为验证群体（被预测）
y.trn[vv,"X1"] <- NA
head(y.trn)
## GBLUP
ans <- mmer(X1~1,
random=~vs(id,Gu=K),
rcov=~units,
data=y.trn, verbose = FALSE) # kinship based
ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")
## rrBLUP
ans2 <- mmer(X1~1,
random=~vs(list(GT)),
rcov=~units,
data=y.trn, verbose = FALSE) # kinship based
u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
rownames(u) <- rownames(GT)
cor(u[vv,],DT[vv,"X1"]) # same correlation
数据使用的是CIMMYT的全球小麦项目组提供的数据，共有599个个体。y.trn表示籽粒产量，每一列表示一个环境，共有4个环境，例子中只计算了80%预测20%的结果。
杂交种的预测使用GCA效应和SCA效应。数据来自1254年测定的玉米单交试验。123个dent和86个flint，用测定的1254材料预测未测定的9324个材料。例子中只计算了第一个性状GY。
模型：
y = Xβ + Zu1 + Zu2 + ZuS + ε代码：
library(sommer)
data(DT_technow)
DT <- DT_technow # 表型数据
Md <- Md_technow # dent马齿基因型矩阵
Mf <- Mf_technow # flint硬粒基因型矩阵
Ad <- Ad_technow # dent的加性效应矩阵
Af <- Af_technow # flint的加性效应矩阵
y.trn <- DT
vv1 <- which(!is.na(DT$GY)) # 获得有数据的材料行号
vv2 <- sample(vv1, 100) # 在有数据的材料行号中随机选择100个材料
y.trn[vv2,"GY"] <- NA # 将抽取的100个材料设为NA
anss2 <- mmer(GY~1, # 固定效应
random=~vs(dent,Gu=Ad) + vs(flint,Gu=Af), # 随机效应，并指定方差协方差矩阵
rcov=~units,
data=y.trn, verbose = FALSE)
summary(anss2)$varcomp
zu1 <- model.matrix(~dent-1,y.trn) %*% anss2$U$`u:dent`$GY
zu2 <- model.matrix(~flint-1,y.trn) %*% anss2$U$`u:flint`$GY
u <- zu1+zu2+anss2$Beta[1,"Estimate"]
cor(u[vv2,], DT$GY[vv2])
VarComp VarCompSE Zratio Constraint
u:dent.GY-GY 16.06423 2.5737578 6.241548 Positive
u:flint.GY-GY 11.42070 2.1591718 5.289390 Positive
units.GY-GY 16.81801 0.7689509 21.871368 Positive
> cor(u[vv2,], DT$GY[vv2])
[1] 0.922328预测时，所有有数据的个体作为建模群体，这里随机减去100个个体的测定值用于验证，即建模群体是有数据的个体-100。zu1是模型中的Zu1，Z是model.matrix(~dent-1,y.trn)，即dent的设计矩阵；u1是anss2$U$u:dent$GY是GY性状dent的BLUP值。相应的，zu2是flint的设计矩阵×GY性状flint的BLUP。
anns2$Beta[1,”Estimate”]是固定效应的BLUP。
最后，100个个体的测定值和这100个个体的BLUP值的相关系数作为预测结果。
上面的例子只使用了GCA进行预测，因为SCA在这里面没有对预测精度起到多大的影响，相反，由于使用了很大的矩阵，导致计算时间明显增加。
sommer是Solving Mixed Model Equations in R的首字母缩写，译为在R中解混合线性方程。
sommer的诞生是为了给R语言用户提供一个高效、稳定的多变量混合模型求解工具。该包主要关注p>n（估计的效应多于观测值），其核心算法由C++编写。该包可以拟合包含随机效应的方差-协方差结构的混合模型，指定非齐次方差（一个变量X对另一个变量Y的影响可能因个体而异），还能获得固定/随机效应的BLUP、BLUE、残差、拟合值、方差等参数。
包的核心是解混合模型，利用一个接口调用the NR Direct-Inversion Newton-Raphson or Average Information algorithms（直接反演牛顿-拉夫逊或平均信息算法）。支持多变量模型。多变量或拓展的单变量形式为：
这里，yi是性状的表型向量；βi是固定效应向量；ui是个体的随机效应向量；ei是第i个性状的残差（i=1…t）。假设随机效应ui和ei服从正态分布且均值为0。X和Z分别是固定效应和随机效应的发生率矩阵。多元变量的响应分布和表型的方差协方差（V）是：
这里，K是第K个随机效应（u=1…k）的关系矩阵或协方差矩阵；H=I是残差项的单位矩阵或部分单位矩阵；σgi2和σei2分别表示第i个性状的遗传（或第k个任意随机项）和残差的方差；σgij和σeij分别是遗传（或第k个任意随机项）和残差的协方差（i=1…t, j=1…t）。算法用对数似然数优化。
这里||是矩阵的行列式，用牛顿优化算法更新REML的估计。
这里，θ是性状间随机效应方差分量和协方差分量的向量。H-1是第k个循环二阶导数Hessian矩阵的逆。dL/dσi2是关于方差协方差组分的似然一阶倒数向量。Newton-Raphson算法包含了Lee和Van Der Werf (2016)提出的关系矩阵的特征分析，提高了时间效率。方差分量线性组合的流行的标准差估计函数（例如遗传力和遗传相关）也被纳入该包。
线性模型的优势在于灵活的指定方差、协方差结构。一般来说，混合模型的方差结构可以看作是多重方差协方差结构的张量（kronecker）积。例如，一个多响应模型（例如2个性状），包含g个个体（例如100个）在e个处理（例如3个环境），个体作为随机效应的方差协方差可以表示为下列乘法模型：
T⨂G⨂A
T是性状间个体的协方差结构。
G是环境中个体的协方差结构。
A是一个方形矩阵，表示个体水平之间的协方差（任何已知的关系矩阵）。
上述T和G矩阵需要估计，而A矩阵已知。T和G叫做非结构化（unstructured, US）矩阵。
参考，未必是这个意思（结构化矩阵（二维表）是大小固定，取值的类型固定的矩阵，实际指的是（可能和原定义不完全符合）每一列（变量）的取值是固定类型，如年龄，只能取整数数字。非结构化矩阵，指每一列（性状）可能取不确定的类型，即值的类型不同，如长度可能取小数，病害可能取级数，开花期可能取时间等。）
协方差结构除了上述例子外，还有其他形式：
对角线（Diagonal，DIAG）协方差结构：实际除了方差外，其他都是0。
复合对称性（Compound symmetry，CS）协方差结构：貌似加入了互作。
一阶自回归（First order autoregressive，AR1）协方差结构：
前面提到过的非结构化的方差结构：
sommer可以拟合上述协方差结构。
（1）单变量齐次方差模型
表型数据的形式为：
这类模型指的是单个响应变量模型，其中兴趣变量（即基因型）需要分析与第二个随机效应（即环境）的相互作用，但假设基因型跨环境具有相同的方差分量。这就是所谓的复合对称（CS）模型。
library(sommer)
data("DT_example")
DT <- DT_example
ans1 <- mmer(Yield~Env,
random=~Name+Env:Name,
rcov=~units,
data=DT,
verbose=FALSE)
summary(ans1)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value -20.14538 46.29075 55.95182 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
Name.Yield-Yield 3.682 1.691 2.177 Positive
Env:Name.Yield-Yield 5.173 1.495 3.460 Positive
units.Yield-Yield 4.366 0.647 6.748 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.496 0.6855 24.065
2 Yield EnvCA.2012 -5.777 0.7558 -7.643
3 Yield EnvCA.2013 -6.380 0.7960 -8.015
============================================================
Groups and observations:
Yield
Name 41
Env:Name 123
============================================================
Use the '$' sign to access results and parametersmmer函数是sommer中拟合单变量或多变量混合线性模型的函数。
logLik（Log-likelihood），说明书中没有写，应该是似然函数取对数的最大值，即最大似然对数函数。
AIC（Akaike information criterion）， 赤池信息量标准。是评估统计模型的复杂度和衡量统计模型“拟合”资料之优良性（英语：Goodness of Fit）的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则建立在信息熵的概念基础上。
一般而言，模型越复杂，模型能解释的点的误差就越小，当然，这是在拟合建模群体时的情况。但是，用另一个测试集来测试，此时，越复杂的点可能误差越大，因为在拟合的时候，模型会把真实情况当作真是情况调整模型，而另外一组真实值，其误差还是误差，此时，复杂的模型很难拟合的很好。从下图看，似然值越小，说明拟合越好。在train的角度，似然值是一直下降的。在test的角度上，有下降、上升的趋势。在AIC角度上也有类似的趋势，因此，AIC作为选择模型的指标，最小值可以得到最优的模型。AIC=2k-2ln(L)，这里k是参数的个数，L是似然数。解读过来就是参数不能太少，也不能太多，找到一个平衡点。
BIC（Bayesian information criterion），贝叶斯信息量准则。是在有限集合中进行模型选择的准则：BIC最低的模型是最好的。该准则部分基于似然函数并与赤池信息量准则（AIC）紧密相关。
BIC=2k×ln(n)-2ln(L)，n是样本数。除了考虑模型参数外，还考虑的样本数。
拟合模型时，增加参数可提高似然，但如此下去可能导致过拟合。BIC与AIC都致力于向模型中引入关于参数数量的惩罚项；其中，BIC中的惩罚项会大于AIC中的惩罚项。
Method，方差估计的方法或算法，包括NR（Direct-inversion Newton-Raphson）和AI（Average Information）。
units是随机误差项。
zratio貌似是z检验比例？未确定。
（2）单变量非齐次方差模型
在多环境实验中，假设跨地点的遗传方差和残差齐次不太合理。因此，指定一个一般的遗传组分和地点特异性遗传方差是一种解决途径。这需要CS+DIAG模型（非齐次CS模型）。
library(sommer)
data("DT_example")
DT <- DT_example
ans2 <- mmer(Yield~Env,
random=~Name+vs(ds(Env),Name),
rcov=~ vs(ds(Env),units),
data=DT,verbose=FALSE)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value -15.42983 36.85965 46.52072 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
Name.Yield-Yield 2.963 1.496 1.980 Positive
CA.2011:Name.Yield-Yield 10.146 4.507 2.251 Positive
CA.2012:Name.Yield-Yield 1.878 1.870 1.004 Positive
CA.2013:Name.Yield-Yield 6.629 2.503 2.649 Positive
CA.2011:units.Yield-Yield 4.942 1.525 3.242 Positive
CA.2012:units.Yield-Yield 5.725 1.312 4.363 Positive
CA.2013:units.Yield-Yield 2.560 0.640 4.000 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.508 0.8268 19.965
2 Yield EnvCA.2012 -5.817 0.8575 -6.783
3 Yield EnvCA.2013 -6.412 0.9356 -6.854
============================================================
Groups and observations:
Yield
Name 41
CA.2011:Name 41
CA.2012:Name 41
CA.2013:Name 41
============================================================
Use the '$' sign to access results and parameters如你所见，at和diag函数都能用于表示每个环境的基因型有不同的方差。残差也一样。不同的是，at可以指定用来指定特别的水平或不同方差的环境。
（3）非结构化模型
除了上面的 CS+DIAG模型假设外，还假设一个特定的随机效应（如环境）与另一个随机效应（如基因型）有协方差结构。可以用us()函数获得。us是非结构化（unstructured model）的缩写。
library(sommer)
data("DT_example")
DT <- DT_example
ans3 <- mmer(Yield~Env,
random=~vs(us(Env),Name),
rcov=~ vs(us(Env),units),
data=DT,verbose=FALSE)
summary(ans3)
===================================================================
Multivariate Linear Mixed Model fit by REML
************************** sommer 4.1 **************************
===================================================================
logLik AIC BIC Method Converge
Value -11.49971 28.99943 38.66049 NR TRUE
===================================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
CA.2011:Name.Yield-Yield 15.665 5.4207 2.890 Positive
CA.2012:CA.2011:Name.Yield-Yield 6.110 2.4851 2.459 Unconstr
CA.2012:Name.Yield-Yield 4.530 1.8208 2.488 Positive
CA.2013:CA.2011:Name.Yield-Yield 6.384 3.0659 2.082 Unconstr
CA.2013:CA.2012:Name.Yield-Yield 0.393 1.5234 0.258 Unconstr
CA.2013:Name.Yield-Yield 8.597 2.4838 3.461 Positive
CA.2011:units.Yield-Yield 4.970 1.5323 3.243 Positive
CA.2012:CA.2011:units.Yield-Yield 4.087 0.0000 Inf Unconstr
CA.2012:units.Yield-Yield 5.673 1.3008 4.361 Positive
CA.2013:CA.2011:units.Yield-Yield 4.087 0.0000 Inf Unconstr
CA.2013:CA.2012:units.Yield-Yield 4.087 0.0000 Inf Unconstr
CA.2013:units.Yield-Yield 2.557 0.6393 4.000 Positive
===================================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.331 0.8137 20.070
2 Yield EnvCA.2012 -5.696 0.7404 -7.693
3 Yield EnvCA.2013 -6.271 0.8191 -7.656
===================================================================
Groups and observations:
Yield
CA.2011:Name 41
CA.2012:CA.2011:Name 82
CA.2012:Name 41
CA.2013:CA.2011:Name 82
CA.2013:CA.2012:Name 82
CA.2013:Name 41
===================================================================
Use the '$' sign to access results and parameters正如我们所见，us(Env)表明，基因型(name)在环境（Env）中有协方差结构。
（4）多变量齐次方差模型
多响应变量模型由某些变量隐藏的相关性决定，从预测角度看有作用。在sommer中，指定多变量模型，响应变量需要在响应中使用cbind()函数；而模型的随机部分使用us(trait)，diag(trait)或at(trait)。
library(sommer)
data("DT_example")
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
ans4 <- mmer(cbind(Yield,Weight)~Env,
random=~vs(Name,Gtc=unsm(2))+vs(EnvName,Gtc=unsm(2)),
rcov=~ vs(units,Gtc=unsm(2)),
data=DT,verbose=FALSE)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value 167.0252 -322.0505 -298.5695 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:Name.Yield-Yield 3.7089 1.68117 2.206 Positive
u:Name.Yield-Weight 0.9071 0.37944 2.391 Unconstr
u:Name.Weight-Weight 0.2243 0.08775 2.557 Positive
u:EnvName.Yield-Yield 5.0921 1.47879 3.443 Positive
u:EnvName.Yield-Weight 1.0269 0.30767 3.338 Unconstr
u:EnvName.Weight-Weight 0.2101 0.06661 3.154 Positive
u:units.Yield-Yield 4.3837 0.64941 6.750 Positive
u:units.Yield-Weight 0.9077 0.14145 6.417 Unconstr
u:units.Weight-Weight 0.2280 0.03377 6.751 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.4093 0.6783 24.191
2 Weight (Intercept) 0.9806 0.1497 6.550
3 Yield EnvCA.2012 -5.6844 0.7474 -7.606
4 Weight EnvCA.2012 -1.1846 0.1593 -7.439
5 Yield EnvCA.2013 -6.2952 0.7850 -8.019
6 Weight EnvCA.2013 -1.3559 0.1681 -8.065
============================================================
Groups and observations:
Yield Weight
u:Name 41 41
u:EnvName 94 94
============================================================
Use the '$' sign to access results and parameters可能注意到了，我们在随机效应后面增加了us(trait)，这表明应该在多变量模型中假设结构。随机效应（如name）后面的diag(trait)表明性状模型（Yield and Weight）没有协方差组分，不应该被估计，而us(trait)假设随机效应有协方差组分（如随机效应Name的产量和重量的协方差）。残差部分相同。
（5）多变量非齐次方差模型
该模型拓展了单变量非齐次方差模型。这是一个CS+DIAG多变量模型。
library(sommer)
data("DT_example")
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
ans5 <- mmer(cbind(Yield,Weight)~Env,
random=~vs(Name,Gtc=unsm(2))+vs(ds(Env),Name,Gtc=unsm(2)),
rcov=~ vs(ds(Env),units,Gtc=unsm(2)),
data=DT,verbose=FALSE)
summary(ans5)
=============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
=============================================================
logLik AIC BIC Method Converge
Value 177.8154 -343.6308 -320.1497 NR TRUE
=============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:Name.Yield-Yield 3.31936 1.45269 2.2850 Positive
u:Name.Yield-Weight 0.79393 0.32621 2.4338 Unconstr
u:Name.Weight-Weight 0.19085 0.07503 2.5438 Positive
CA.2011:Name.Yield-Yield 8.70657 4.01470 2.1687 Positive
CA.2011:Name.Yield-Weight 1.77892 0.83926 2.1196 Unconstr
CA.2011:Name.Weight-Weight 0.35966 0.17903 2.0089 Positive
CA.2012:Name.Yield-Yield 2.57109 1.94951 1.3188 Positive
CA.2012:Name.Yield-Weight 0.33245 0.39840 0.8345 Unconstr
CA.2012:Name.Weight-Weight 0.03842 0.08595 0.4470 Positive
CA.2013:Name.Yield-Yield 5.46908 2.16307 2.5284 Positive
CA.2013:Name.Yield-Weight 1.34713 0.50479 2.6687 Unconstr
CA.2013:Name.Weight-Weight 0.32902 0.12208 2.6952 Positive
CA.2011:units.Yield-Yield 4.93852 1.52318 3.2422 Positive
CA.2011:units.Yield-Weight 0.99447 0.32150 3.0932 Unconstr
CA.2011:units.Weight-Weight 0.23982 0.07394 3.2433 Positive
CA.2012:units.Yield-Yield 5.73887 1.31533 4.3631 Positive
CA.2012:units.Yield-Weight 1.28009 0.30157 4.2448 Unconstr
CA.2012:units.Weight-Weight 0.31806 0.07286 4.3652 Positive
CA.2013:units.Yield-Yield 2.56127 0.63993 4.0024 Positive
CA.2013:units.Yield-Weight 0.44569 0.12645 3.5246 Unconstr
CA.2013:units.Weight-Weight 0.12232 0.03057 4.0009 Positive
=============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.4243 0.7891 20.815
2 Weight (Intercept) 0.9866 0.1683 5.863
3 Yield EnvCA.2012 -5.7339 0.8266 -6.937
4 Weight EnvCA.2012 -1.1998 0.1698 -7.066
5 Yield EnvCA.2013 -6.3128 0.8757 -7.209
6 Weight EnvCA.2013 -1.3621 0.1915 -7.114
=============================================================
Groups and observations:
Yield Weight
u:Name 41 41
CA.2011:Name 41 41
CA.2012:Name 41 41
CA.2013:Name 41 41
=============================================================
Use the '$' sign to access results and parameters（6）多变量非结构化方差模型
模型是单变量非结构化方差模型的拓展，是一个us多变量模型。
library(sommer)
data("DT_example")
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
ans6 <- mmer(cbind(Yield,Weight)~Env,
random=~vs(us(Env),Name,Gtc=unsm(2)),
rcov=~ vs(ds(Env),units,Gtc=unsm(2)),
data=DT,verbose=FALSE)
====================================================================
Multivariate Linear Mixed Model fit by REML
************************** sommer 4.1 **************************
====================================================================
logLik AIC BIC Method Converge
Value 181.7945 -351.5889 -328.1079 NR TRUE
====================================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
CA.2011:Name.Yield-Yield 15.6451 5.35692 2.921 Positive
CA.2011:Name.Yield-Weight 3.3586 1.14633 2.930 Unconstr
CA.2011:Name.Weight-Weight 0.7182 0.24871 2.888 Positive
CA.2012:CA.2011:Name.Yield-Yield 6.5289 2.48615 2.626 Positive
CA.2012:CA.2011:Name.Yield-Weight 1.3505 0.52388 2.578 Unconstr
CA.2012:CA.2011:Name.Weight-Weight 0.2842 0.11259 2.524 Positive
CA.2012:Name.Yield-Yield 4.7893 1.86183 2.572 Positive
CA.2012:Name.Yield-Weight 0.8640 0.38377 2.251 Unconstr
CA.2012:Name.Weight-Weight 0.1693 0.08354 2.027 Positive
CA.2013:CA.2011:Name.Yield-Yield 5.9934 2.93830 2.040 Positive
CA.2013:CA.2011:Name.Yield-Weight 1.4232 0.64973 2.190 Unconstr
CA.2013:CA.2011:Name.Weight-Weight 0.3379 0.14680 2.302 Positive
CA.2013:CA.2012:Name.Yield-Yield 2.0987 1.44034 1.457 Positive
CA.2013:CA.2012:Name.Yield-Weight 0.5240 0.32356 1.619 Unconstr
CA.2013:CA.2012:Name.Weight-Weight 0.1342 0.07572 1.772 Positive
CA.2013:Name.Yield-Yield 8.6257 2.47811 3.481 Positive
CA.2013:Name.Yield-Weight 2.1048 0.58748 3.583 Unconstr
CA.2013:Name.Weight-Weight 0.5125 0.14285 3.588 Positive
CA.2011:units.Yield-Yield 4.9516 1.52694 3.243 Positive
CA.2011:units.Yield-Weight 0.9993 0.32286 3.095 Unconstr
CA.2011:units.Weight-Weight 0.2411 0.07432 3.244 Positive
CA.2012:units.Yield-Yield 5.7790 1.32423 4.364 Positive
CA.2012:units.Yield-Weight 1.2914 0.30408 4.247 Unconstr
CA.2012:units.Weight-Weight 0.3212 0.07356 4.366 Positive
CA.2013:units.Yield-Yield 2.5567 0.63883 4.002 Positive
CA.2013:units.Yield-Weight 0.4452 0.12631 3.524 Unconstr
CA.2013:units.Weight-Weight 0.1223 0.03056 4.001 Positive
====================================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.3342 0.8254 19.790
2 Weight (Intercept) 0.9677 0.1770 5.466
3 Yield EnvCA.2012 -5.6637 0.7449 -7.604
4 Weight EnvCA.2012 -1.1855 0.1604 -7.390
5 Yield EnvCA.2013 -6.2153 0.8340 -7.453
6 Weight EnvCA.2013 -1.3406 0.1806 -7.425
====================================================================
Groups and observations:
Yield Weight
CA.2011:Name 41 41
CA.2012:CA.2011:Name 82 82
CA.2012:Name 41 41
CA.2013:CA.2011:Name 82 82
CA.2013:CA.2012:Name 82 82
CA.2013:Name 41 41
====================================================================
Use the '$' sign to access results and parameters任何随机效应都可以指定不同的结构。
（7）方差模型特殊函数的细节
方差模型主函数vs()和辅助函数
sommer的函数vs()允许构建复杂的方差模型并传递给mmer()函数。vs()函数的格式如下：
random=~vs(..., Gu, Gti, Gtc)
这表示，vs()函数反映了矩阵中每个随机效应可能有特别的方差结构。
var(u)=T⊗E⊗…⊗A
这里，…表示vs()函数的参数，用来指定随机效应方差矩阵的keronecker products乘积。辅助函数ds()、us()、cs()、at()可用于定义这些方差结构。因为，随机变量x（如个体）的方差模型可能需要更加灵活，且参数不应只有一个：
random=~x~被读作预测，~x表示用x预测。
例如，如果个体在不同的时间点或环境下测试，我们可以假设在不同的环境-时间点组合中，个体间有不同的方差和协方差组分。例如：
var(u)=T⊗E⊗S⊗A
可以在vs()函数中指定：
random=~vs(us(e),us(s),x,Gu=A,Gtc=T)
这里，e是数据框中环境的列向量，s是数据框中时间点的列向量，x是数据框中个体的向量，A是个体间已知的方差协方差矩阵（通常是一个单位矩阵；默认不指定），T是协方差矩阵且行和列的数量等于性状的数量。
为随机效应构建协方差模型的辅助函数是+ds()对角线协方差结构+us()非结构化协方差+at()是指定水平的非齐次方差结构+cs()自定义协方差结构。
ds()指定对角线（DIAG）协方差结构
对角线协方差结构类似下面的形式：
考虑一个随机效应的例子，g为个体，在e个环境中测定，模型则为：
random=~vs(ds(e),g)
其意义是用向量g与对角线方差矩阵e的张量积做预测。
us()指定非结构化协方差
非结构化协方差表示如下：
考虑相同的例子，个体在环境中的模型为：
random=~vs(us(e),g)
at()用来指定特殊水平的非齐次方差
第二随机效应（方差和协方差）的个对角线协方差结构是这样个样子的：
还是之前的例子，g是个体，有ABC三个处理（环境），则模型可能是：
random=~vs(at(e,c("A","B")),g)
这里方差组分只为g拟合水平A和B。
cs()指定水平的方差协方差结构
第二个随机变量的自定义协方差结构可能形式如下：
依然用之前的例子，g是个体，e是三个处理（环境）A、B、C，模型可能是：
random=~vs(cs(e,mm),g)
这里mm表明为g估计的方差和协方差组分。
（8）方差组分的约束规则（Gtc参数）
sommer的一大优势是在多性状框架下能够非常灵活的指定方差协方差结构。sommer 3.7版后可以非常容易的通过vs()和Gtc参数实现。Gtc通过以下填充数字，用于随机效应方差协方差组分的约束矩阵的预测。
0：不顾及参数
1：有约束的估计
2：无约束的估计
3：不估计，但是在Gti中提供固定值
unsm()用来快速指定非结构化的约束矩阵，fixm()用于固定值约束，fcm()用于固定效应约束。
考虑到4个性状的（y1，y2，y3，y4）多性状模型，一个随机效应（u）和1个固定效应x
fixed=cbind(y1,y2,y3,y4)~x
random=~vs(u,Gtc=?)
4×4的方差协方差组分估计的约束可以被估计为：
（a）非结构化（方差组分必须是正值，协方差可正可负）
random=~vs(u,Gtc=unsm(4))
unsm(4)
[,1] [,2] [,3] [,4]
[1,] 1 2 2 2
[2,] 2 1 2 2
[3,] 2 2 1 2
[4,] 2 2 2 1（b）非协方差（任意组分方差或协方差可以正或负）
random=~vs(u,Gtc=uncm(4))
uncm(4)
[,1] [,2] [,3] [,4]
[1,] 2 2 2 2
[2,] 2 2 2 2
[3,] 2 2 2 2
[4,] 2 2 2 2（c）固定值（方差或协方差组分表明3被作为固定值）
random=~vs(u,Gtc=fixm(4),Gti=mm)
fixm(4)
[,1] [,2] [,3] [,4]
[1,] 3 3 3 3
[2,] 3 3 3 3
[3,] 3 3 3 3
[4,] 3 3 3 3这里mm是4×4的矩阵，该矩阵具有方差分量的初始值。
（d）考虑固定效应
fixed=cbind(y1,y2,y3,y4)~vs(x,Gtc=fcm(c(1,0,1,0)))
fcm(c(1,0,1,0))
[,1] [,2]
[1,] 1 0
[2,] 0 0
[3,] 0 1
[4,] 0 0这里1和0表示性状，固定效应1被估计0不估计。
（9）特殊模型的特殊函数
随机回归模型
拟合随机回归模型，用户可以使用leg()函数拟合Legendre polynomials多项式。这可以用其他协方差结构，如ds()和us()等。
library(orthopolynom)
data(DT_legendre)
DT <- DT_legendre
mRR2<-mmer(Y~ 1 + Xf,
random=~ vs(us(leg(X,1)),SUBJECT),
rcov=~vs(units),
data=DT, verbose = FALSE)
summary(mRR2)$varcomp
VarComp VarCompSE Zratio Constraint
leg0:SUBJECT.Y-Y 2.5782969 0.6717242 3.838326 Positive
leg1:leg0:SUBJECT.Y-Y 0.4765431 0.2394975 1.989763 Unconstr
leg1:SUBJECT.Y-Y 0.3497299 0.2183229 1.601893 Positive
u:units.Y-Y 2.6912226 0.3825197 7.035513 Positive这里，协方差X用于解释主体的轨迹，结合了一个非结构化的协方差矩阵。细节可以查看理论。
GWAS模型
尽管基因组相关分析可以考虑通过各种方法，混合模型仍然是最流行的，我们采用这个方法。最流行和经典的两个方法，第一个通过标记混合建模，获得标记效应，提供p值-log10。第二个假设所有标记的遗传方差组分相似，因此，方差组分只被估计一次，在一般线性模型中使用表型方差矩阵的倒数（V-inverse）测试所有标记，b=(X’V-X)-XV-y。此时的GWAS更快更高效。sommer提供直接的拓展，GWAS函数可以拟合以上两种方法。库中已经有很多类似的方法。sommer中获得标记效应的一般线性模型的形式是：
b=(X’V-X)X’V-y
X=ZMi
这里，b是标记效应（维度1×mt），y是响应变量（单或多变量）（维度1×nt），V是表型方差矩阵的逆（维度nt×nt），Z是随机效应选择的关系矩阵，用来执行GWAS（维度nt×ut），Mi是标记矩阵（M参数）第i列（维度u×m）。
t是性状数量，n是观测值数量，m是标记数量，u是随机效应水平数量。如果P3D=TRUE，计算一次V矩阵并用于所有的标记测试；FALSE，每个标记使用REML估计。
这里我们展示一个简单的GWAS模型，单个性状：
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
#### create the variance-covariance matrix
A <- A.mat(GT) # additive relationship matrix
#### look at the data and fit the model
head(DT,3)
id Row Col Year color Yield FruitAver Firmness Rowf Colf
P003 P003 3 1 2014 0.10075269 154.67 41.93 588.917 3 1
P004 P004 4 1 2014 0.13891940 186.77 58.79 640.031 4 1
P005 P005 5 1 2014 0.08681502 80.21 48.16 671.523 5 1
head(MP,3)
Locus Position Chrom
1 scaffold_77830_839 0 1
2 scaffold_39187_895 0 1
3 scaffold_50439_2379 0 1
GT[1:3,1:4]
scaffold_50439_2381 scaffold_39344_153 uneak_3436043 uneak_2632033
P003 0 0 0 1
P004 0 0 0 1
P005 0 -1 0 1
mix1 <- GWAS(color~1,
random=~vs(id,Gu=A)
+ Rowf + Colf,
rcov=~units,
data=DT,
M=GT, gTerm = "u:id",
verbose = FALSE)
## Performing GWAS evaluation
ms <- as.data.frame(mix1$scores)
ms$Locus <- rownames(ms)
MP2 <- merge(MP,ms,by="Locus",all.x = TRUE);
manhattan(MP2, pch=20,cex=.5, PVCN = "color")
需要注意，标记矩阵M必须补缺（不允许缺失值）。确保M矩阵中的行数等于g项，如id有300个个体，M矩阵为有300×m，m是标记数量。
叠加模型[overlay()函数]
overlay()叠加不同随机效应的矩阵，估计叠加项的单个方差组分。
data("DT_halfdiallel")
DT <- DT_halfdiallel
head(DT)
rep geno male female sugar
1 1 12 1 2 13.950509
2 2 12 1 2 9.756918
3 1 13 1 3 13.906355
4 2 13 1 3 9.119455
5 1 14 1 4 5.174483
6 2 14 1 4 8.452221
DT$femalef <- as.factor(DT$female)
DT$malef <- as.factor(DT$male)
DT$genof <- as.factor(DT$geno)
#### model using overlay
modh <- mmer(sugar~1,
random=~vs(overlay(femalef,malef))
+ genof,
data=DT,verbose = FALSE)
这里，famalef和malef随机效应被叠加，变成了单个随机效应，具有相同的方差组分。
空间模型（使用二维条线）
在田间实验设计中，使用CPdata展示3维条线的使用，用于适应田间设计的空间效应。早期的各种试验，可用的种子很少，因此，需要使用非重复的设计。实验设计例如增强设计和部分重复设计（p-rep）现在非常流行。
适应空间趋势，田间空间协方差矩阵被提出（即自回归残差；arl）。不幸的是，这些协方差矩阵使得模型非常不稳定。最近，其他研究组提出了使用2维空间去克服上述问题，空间项建模更强大。
下面的例子，我们假设没有重复的群体，row和range信息可用，允许拟合2维条线模型。
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
### mimic two fields
A <- A.mat(GT)
mix <- mmer(Yield~1,
random=~vs(id, Gu=A) +
vs(Rowf) +
vs(Colf) +
vs(spl2D(Row,Col)),
rcov=~vs(units),
data=DT,verbose = FALSE)
summary(mix)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value -151.2011 304.4021 308.2938 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:id.Yield-Yield 783.4 319.3 2.4536 Positive
u:Rowf.Yield-Yield 814.7 390.5 2.0863 Positive
u:Colf.Yield-Yield 182.2 129.7 1.4053 Positive
u:Row.Yield-Yield 513.6 694.7 0.7393 Positive
u:units.Yield-Yield 2922.6 294.1 9.9368 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 132.1 8.791 15.03
============================================================
Groups and observations:
Yield
u:id 363
u:Rowf 13
u:Colf 36
u:Row 168
============================================================
Use the '$' sign to access results and parametersGT矩阵作为随机效应被封装在一个列表中，在vs()中使用。
（10）全基因组选择
你可以用rrBLUP和GBLUP模型拟合。本例子在群体中指定个体。基本形式为：
利用已知的基因型效应方差协方差矩阵作为加性关系矩阵（A），使用A.met函数建立所有个体的联系，并预测未测定个体的BLUP。本例中，标记矩阵是随机效应，标记的关系也可以被指定，但是在这里假设是一个对角线矩阵。
data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
colnames(DT) <- paste0("X",1:ncol(DT))
DT <- as.data.frame(DT);DT$id <- as.factor(rownames(DT))
# select environment 1
rownames(GT) <- rownames(DT)
K <- A.mat(GT) # additive relationship matrix
colnames(K) <- rownames(K) <- rownames(DT)
# GBLUP pedigree-based approach
set.seed(12345)
y.trn <- DT
vv <- sample(rownames(DT),round(nrow(DT)/5))
y.trn[vv,"X1"] <- NA
## GBLUP
ans <- mmer(X1~1,
random=~vs(id,Gu=K),
rcov=~units,
data=y.trn,verbose = FALSE) # kinship based
ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")
[1] 0.5737594
## rrBLUP
ans2 <- mmer(X1~1,
random=~vs(list(GT)),
rcov=~units,
data=y.trn,verbose = FALSE) # kinship based
u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
rownames(u) <- rownames(GT)
cor(u[vv,],DT[vv,"X1"]) # same correlation
# the same can be applied in multi-response models in GBLUP or rrBLUP
[1] 0.5737594~1的1表示截距，基本公式为y=C+ε，即常数。
（11）似然比测试
似然比测试（LRT）是调查随机效应或特定方差协方差组分显著性的方式。
（11.1）方差组分显著性测试
例如，想象研究的模型在加入了空间内核效应后改善了，那么模型可能是：
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
### mimic two fields
A <- A.mat(GT)
mix1 <- mmer(Yield~1,
random=~vs(id, Gu=A) +
vs(Rowf) +
vs(Colf),
rcov=~vs(units),
data=DT, verbose = FALSE)
带有空间内核的模型为：
mix2 <- mmer(Yield~1,
random=~vs(id, Gu=A) +
vs(Rowf) +
vs(Colf) +
vs(spl2D(Row,Col)),
rcov=~vs(units),
data=DT,verbose = FALSE)
似然比测试，检测第二个模型：
lrt <- anova(mix1, mix2)
Likelihood ratio test for mixed models
==============================================================
Df AIC BIC loLik Chisq ChiDf PrChisq
mod2 8 304.4021 308.2938 -151.2011
mod1 7 305.0477 308.9393 -151.5238 0.64554 1 0.42171
==============================================================
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1似然值增加了，不显著。
（11.2）测试协方差组分的显著性
有些时候，研究人员对协方差结构是否感相关感兴趣。假设我们有两个多性状模型，1）拟合非协方差（独立），2）拟合群体中产量和颜色的遗传协方差：
data(DT_example)
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
modelBase <- mmer(cbind(Yield, Weight) ~ Env,
random= ~ vs(Name, Gtc=diag(2)), # here is diag()
rcov= ~ vs(units, Gtc=unsm(2)),
data=DT,verbose = FALSE)
modelCov <- mmer(cbind(Yield, Weight) ~ Env,
random= ~ vs(us(Env),Name, Gtc=unsm(2)), # here is unsm()
rcov= ~ vs(ds(Env),units, Gtc=unsm(2)),
data=DT,verbose = FALSE)
lrt <- anova(modelBase, modelCov)
Likelihood ratio test for mixed models
==============================================================
Df AIC BIC loLik Chisq ChiDf PrChisq
mod2 45 -351.5889 -328.1079 181.7945
mod1 23 -253.4383 -229.9573 132.7192 98.15058 22 1.3552747104066e-11 ***
==============================================================
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1基因型的协方差显著改善了拟合，卡方分布小于0.05。带有协方差的模型更适合。