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。带有协方差的模型更适合。
发表评论