sommer包的预测函数

（1）预测的背景

在混合线性模型中，y是n×1的观测值向量，混合线性模型可以写成：

这里τ是固定效应的t×1的向量；X是n×t的设计矩阵，它将观测值与固定效应联系起来；u是q×1的随机效应向量；Z是n×q的矩阵，它联系观测值与随机效应；e是n×1的残差误差向量。W和β分别表示组合的设计矩阵和效应向量，它假设：

G和R是随机效应和残差的协方差矩阵，函数的参数是γ和φ。协方差矩阵的数据应该写作：

var(y)=σ²(ZGZ’+R)

协方差的参数γ和φ 用有限的最大似然或最大似然估计。

已知D、 γ和φ 的线性组合的BLUP，然后有：

这里

是混合模型公式的解：

这些也可以写作：

这里的τbar是τ最佳线性无偏估计（BLUE），ubar是u的最佳线性无偏预测（BLUP），这些用协方差var(βbar-β)=C^-1。考虑到，为了更加清晰，使用能形成置信区间预测误差方差（Prediction error variance，PEV或var(βbar-β)，而不是通常感兴趣的方差估计者var(βbar)。对于线性组合Dβbar，PEV是DC^-1D^T。因为方差参数是未知的，因此，我们用它们的REML估计替换掉未知方差参数，并使用经验值。

sommer包中的预测函数用混合模型公式和线性组合D矩阵构建C^-1矩阵用于预测。对于均值，它使用线性组合的D矩阵乘以需要固定和随机效应的向量Xτbar+Zubar或Dβbar，这里D矩阵是来自X和/或Z矩阵的固定和/或随机效应的线性组合。

注意，由于这里打不出来βbar，因此用βbar代替。τbar和ubar同理。

（2）预测均值

sommer用predict()函数计算mmer函数中固定和随机效应的均值和标准差。使用yatesoats数据集你和固定和随机效应。

data(DT_yatesoats)
DT <- DT_yatesoats
m3 <- mmer(fixed=Y ~ V + N + V:N,
           random = ~ B + B:MP,
           rcov=~units,
           data = DT, verbose=FALSE)
summary(m3)$varcomp

           VarComp VarCompSE   Zratio Constraint
B.Y-Y     214.4477 168.62790 1.271722   Positive
B:MP.Y-Y  106.0508  67.83280 1.563415   Positive
units.Y-Y 177.0883  37.34293 4.742217   Positive

现在，该模型可以和分类参数一起使用，获得分类的均值。例如，模型包括了固定公式项”V“用于品种，”N“用于指定氮肥处理，”V:N“表示品种和氮肥互作。分类参数可以用于指定需要的均值。下面的例子，氮肥处理的均值：

p0 <- predict.mmer(object=m3, classify = "N")

iteration    LogLik     wall    cpu(sec)   restrained
    1      -0.745921   16:3:8      0           0
    2      -0.745921   16:3:8      0           0
    3      -0.745921   16:3:8      0           0
    4      -0.745921   16:3:8      0           0

p0$pvals

  trait   N predicted.value standard.error
1     Y   0        79.38889       9.006796
2     Y 0.2        98.88889       9.006796
3     Y 0.4       114.22222       9.006796
4     Y 0.6       123.38889       9.006796

结束语

记住，sommer使用直接反演算法（direct inversion，DI），对于大数据集来说可能非常慢。该包主要针p>n的情况（随机效应比观测值多）且模型使用密集的协方差结构。例如，适用于协方差结构密集且重复较少的实验。也适用于随机效应较多的基因组问题。对于重复很多的情况，如200个个体，10次重复，有2000条记录时，asreml的速度要高于sommer。