随机模拟R语言 r语言随机效应

R中sommer包中mmer()函数

R中sommer包中的mmer函数

• R中sommer包中mmer()函数

• 1、函数random()详解

• 1.1vsr()
• 1.2 dsr(x), usr(x), csr(x) and atr(x,levs)
• 1.3 unsm(x), fixm(x) and diag(x)
• 1.4 gvsr(...,Gu,Guc,Gti,Gtc)

• 2、特殊的方差结构

• 2.1 vsr(atr(x,levels),y)
• 2.2 vsr(dsr(x),y)
• 2.3 vsr(usr(x),y)
• 2.4 vsr(overlay(...,rlist=NULL,prefix=NULL))

mmer()是sommer包中一个函数。mmer函数中的算法用于求解p >N(即要估计的系数比数据集中的记录多)或密集矩阵。

语法如下：

mmer(fixed, random, rcov, data, weights, W, nIters=20, tolParConvLL = 1e-03, 
     tolParInv = 1e-06, init=NULL, constraints=NULL,method="NR", getPEV=TRUE,
     naMethodX="exclude", naMethodY="exclude",returnParam=FALSE, 
     dateWarning=TRUE,date.warning=TRUE,verbose=TRUE, reshapeOutput=TRUE, stepWeight=NULL,
     emWeight=NULL)

简单释义：

参数	含义
fixed	指定响应变量和固定变量的公式，形如：response~covariate
random	指定随机效应名称的公式，形如：random= ~ genotype + year
rcov	误差项
data	指定数据集，其中包含响应效应、固定效应和随机效应
weights	权重的协变量名称。用于乘积R = Wsi RWsi，其中*为矩阵乘积，Wsi为W逆的平方根，R为残差矩阵
nIters	允许的最大迭代次数
tolParConvLL	对数似然变化的收敛准则
tolParInv	在估计过程中遇到奇异点时，矩阵逆的公差参数
init	方差分量的初始值。默认情况下，这是NULL，方差组件的初始值由算法提供，但如果用户想为所有var-cov组件提供初始值，这个参数是有用的。它必须作为一个列表提供，其中每个列表元素对应一个随机效应(1x1矩阵)，如果采用多性状模型，列表中的每个元素都是这种随机效应的性状间方差协方差分量的矩阵。初始值也可以在vsr函数的Gti参数中提供
constraints	当提供初始值时，这些初始值必须伴随着它们的约束条件。有关约束的更多细节，请参阅vsr函数。强烈建议使用vsr函数的Gti和Gtc参数来代替这个参数，但这些参数可以用来一次性提供所有约束
method	用于估计方差分量的算法，包括牛顿-拉弗森方法和平均信息法
getPEV	一个TRUE/FALSE值，指示程序是否应该返回预测的误差方差和随机效应的方差。之所以提供这个选项，是因为在某些模型中p在很大程度上为> n时，这可能需要很长时间
naMethodX	有两个参数值;“include”或“exclude”。如果选择“include”，则函数将用中位数对固定效应的X矩阵进行插值。如果“exclude”被选中，它将去掉所有缺少X(固定)协变量值的行。默认值为“exclude”。“include”选项应该谨慎使用
naMethodY	有三个可能值;“include”、“include2"或"exclude”(默认值)用于处理响应变量中的观察值，以便用于方差分量的估计。第一个选项“include”将为所有具有中值的行输入响应变量，而“include2”仅为至少有一个响应(仅在多响应模型中可用)存在观察值的行输入响应。如果选择“exclude”(默认)，它将删除响应中至少一个响应缺少值的行。
returnParam	一个TRUE/FALSE值，指示程序是否应该返回用于拟合模型而不是拟合模型的参数
dateWarning	TRUE/FALSE值，指示程序是否应在更新sommer包时发出警告
verbose	TRUE/FALSE值，指示程序是否应返回迭代算法的进度
stepWeight	一个值的向量(长度等于迭代的次数)，表示在每次迭代中用于乘以方差组件的更新(增量)的权重。如果为NULL，第一次迭代将乘以0.5，第二次乘以0.7，其余的乘以0.9。这个参数可以帮助避免在初始迭代中方差分量超出参数空间，这在NR方法中不经常发生，但可以通过观察可能性的行为来检测
emWeight	一个值的向量(长度等于迭代次数)，用0到1之间的值表示分配给EM信息矩阵的权重。将1 - emWeight值应用于NR/AI信息矩阵，生成联合信息矩阵

1、函数random()详解

指定随机效应名称的公式，即random= ~genotype+year
可以用来拟合异构方差和其他特殊模型

1.1vsr()

vsr是mmer函数中为拟合随机效应建立方差-协方差结构的主要函数
语法：

vsr(..., Gu=NULL, Gti=NULL, Gtc=NULL, reorderGu=TRUE, buildGu=TRUE)

… 方差结构要按照内部克罗内克积中所需的逻辑来指定。例如，如果用户想要定义随机效应“基因型”(g)相对于随机效应“环境”(e)的对角方差结构，则为:$var(g)=G.e @ I.g$
其中：$G.e$是g(基因型)在e(环境)各水平的方差协方差分量的矩阵
          $I.g$是基因型g之间的协方差阵，或基因关系矩阵
          $@$是克罗内克积，这将在mmer函数中指定为$random=$~$vsr(dsr(e),g)$
sommer的一个优点是能够指定非常复杂的结构和尽可能多的克罗内克积。例如：$var(g)=G.e @ G.F @ G.h @ I.g$，R中实现:$random$=~$vsr(e,f,h,g)$

Gu具有第u个随机效应水平的已知方差-协方差值的矩阵(即个体之间的关系矩阵或任何其他已知的协方差矩阵)。如果为NULL，则假设一个单位矩阵。Gu矩阵可以比与之相关的随机效应中的矩阵有更多的能级，但反之则不然。否则，将返回Gu中缺少级别的错误消息

Gti维数为$t*t$($t$等于性状的数量)的矩阵，带有…中指定的随机效应的方差-协方差分量的初始值论点。如果为NULL，程序将提供初始值

Gtc矩阵的维数为$t*t$ ($t$等于性状的数量)，约束为随机效应的方差-协方差分量。参数根据以下规则:
          0：不可估量
          1：估计并约束为正(方差分量)
          2：估计的和无约束的(可以是负的或正的，即协方差分量)
          3：不是估计的，而是固定的(值必须在Gti参数中提供)

• 在多响应场景中，如果没有指定这个参数，默认是构建一个非结构化矩阵(使用unsm()函数)。用于指定约束的有用函数有：Diag ()， unsm()， fixm()。

reorderGu如果Gu矩阵应该根据随机效应的设计矩阵的名称重新排序或使用用户的自定义顺序传递，则为TRUE/FALSE语句。当以定制的方式拟合协方差分量时，这可能很重要。仅限高级用户使用。

buildGu一个TRUE/FALSE语句，指明当参数Gu=NULL时，Gu矩阵是否应该构建在R中。只有当实参Gu的值等于NULL时，才重复。在某些情况下，当关联矩阵很宽时(例如rrBLUP模型)，协方差结构是一个巨大的p × p矩阵，在执行矩阵操作时可以避免。通过设置这个参数为FALSE，它允许跳过形成这个协方差矩阵。

1.2 dsr(x), usr(x), csr(x) and atr(x,levs)

$X$是随机效应的观测向量
可用于指定未知的对角、非结构化和自定义的非结构化和对角协方差结构，由REML估计。

1.3 unsm(x), fixm(x) and diag(x)

可以用来构建简单的矩阵，以在vsr()函数的Gtc参数中指定约束。
$X$指定一个给定的随机效应所拟合的特征的数量。

1.4 gvsr(…,Gu,Guc,Gti,Gtc)

GVSR函数在拟合随机效应中建立一般的方差-协方差结构组合

gvsr(..., Gu=NULL, Guc=NULL, Gti=NULL, Gtc=NULL, form=NULL)

… 用于创建一个基因方差结构的随机效应的名字(属于数据集中的变量)
例如：有两个随机变量mom和progeny
输入$gvsr(mom, progeny)$
得到下面形式的方差结构：
$\begin{pmatrix} sigma2.m & sigma.pm \\ sigma.pm & sigma2.p \end{pmatrix} \quad$
其中不仅估计每个随机效应的方差分量，而且估计两个随机效应之间的协方差分量

Gu在“…”参数中提供的不同随机效应水平的已知方差-协方差值的矩阵列表。如果为NULL，则假设是单位矩阵。例如，具有协方差结构的2个随机效应的模型应为:$gvsr(mom, progeny, Gu=list(Am,Ap))$
其中$Am$和$Ap$分别是母体和子代随机效应的关系矩阵

Gti 是一个$t$阶方阵($t$是性状的数量)，指定"…"中随机效应的方差-协方差分量的初始值点。如果是NULL，R将提供默认初始值

Gtc 是一个$t$阶方阵($t$等于性状的数量)，表示为随机效应的方差-协方差分量。参数根据以下规则:

• 0：不用估计
• 1：估计并约束为正(为方差分量)
• 2：估计且无约束(可为负或正，即协方差分量)
• 3:不是估计的，而是固定的(值必须在Gti参数中提供)

form 在定义的一般协方差结构之上指定一个克罗内克乘积的附加结构

2、特殊的方差结构

2.1 vsr(atr(x,levels),y)

可用于在协变量x的特定水平上指定“y”协变量的异方差
例如：$random=$~$vsr(at(Location,c("A","B")),ID)$
表明：在水平$A$和$B$的位置上，$ID$的方差变量

2.2 vsr(dsr(x),y)

可用于为协变量“x”的所有水平上的“y”协变量，指定为对角协方差结构
例如：$random=$~$vsr(dsr(Location),ID)$
表明：在协变量Location的所有级别上，ID的方差分量

2.3 vsr(usr(x),y)

可用于为协变量“x”的所有级别的“y”协变量指定非结构化协方差结构
例如：$random=$~$vsr(usr(Location),ID)$
表明：在协变量Location的所有级别上为ID拟合方差和协方差分量

2.4 vsr(overlay(…,rlist=NULL,prefix=NULL))

可用于指定在连续随机效应之间的设计矩阵的叠加
例如：$random=$~$vsr(overlay(male,female))$
表明：男性和女性随机效应的关联矩阵，以获得两个效应的单一方差分量。'rlist’参数是一个列表，其中每个元素都是一个数值，它乘以要覆盖的关联矩阵。