2

我们可以从 R 中的混合模型计算 BLUP 的标准误差 (SE) 吗？我正在使用一个名为 ASReml 的自定义包，它可以根据预测值计算 SE。但是我不确定在哪里可以找到 BLUP 值的 SE。

到目前为止，我有这样的事情......

df<-data.frame(inbred=c("x1","x2","x3","x4","x5","x1","x2","x3","x4","x5","x1","x2","x3","x4","x5","x1","x2","x3","x4","x5","x1","x2","x3","x4","x5","x1","x2","x3","x4","x5","x1","x2","x3","x4","x5","x1","x2","x3","x4","x5"),
                       trait1=rnorm(40,0,1),
                       block=c(1,1,1,1,1,2,2,2,2,2,1,1,1,1,1,2,2,2,2,2,1,1,1,1,1,2,2,2,2,2,1,1,1,1,1,2,2,2,2,2),
                       rep=c(rep(1,20),rep(2,20)))
df <- transform(df, inbred=factor(inbred), rep=factor(rep), block=factor(block))
head(df)

  inbred      trait1 block rep
1     x1 -1.15668530     1   1
2     x2  0.41492671     1   1
3     x3 -0.08740545     1   1
4     x4  0.37983415     1   1
5     x5  0.27180581     1   1
6     x1  1.22986338     2   1

使用 asreml 来拟合混合模型。

library(asreml)

df.reml < -asreml(trait1 ~  rep + block,
                  random = ~ inbred,
                  data=df)

Model fitted using the gamma parameterization.
ASReml 4.1.0 Fri Mar 12 13:56:29 2021
          LogLik        Sigma2     DF     wall    cpu
 1      -16.4111      0.669903     37 13:56:29    0.0 (1 restrained)
 2      -15.9428      0.692305     37 13:56:29    0.0 (1 restrained)
 3      -15.9183      0.694848     37 13:56:29    0.0 (1 restrained)
 4      -15.9168      0.695017     37 13:56:29    0.0 (1 restrained)
 5      -15.9167      0.695027     37 13:56:29    0.0 (1 restrained)
Warning message:
In asreml(trait1 ~ rep + block, random = ~inbred, data = df) :
  Some components changed by more than 1% on the last iteration.

然后我们要求预测值：

df.pred <- predict(df.reml,classify="inbred")
Model fitted using the gamma parameterization.
ASReml 4.1.0 Fri Mar 12 13:59:52 2021
          LogLik        Sigma2     DF     wall    cpu
 1      -15.9167      0.695028     37 13:59:53    0.0
 2      -15.9167      0.695028     37 13:59:53    0.0
 3      -15.9167      0.695028     37 13:59:53    0.0

df.pred$pvals

Notes:
- The predictions are obtained by averaging across the hypertable
  calculated from model terms constructed solely from factors in
  the averaging and classify sets.
- Use 'average' to move ignored factors into the averaging set.
- The simple averaging set: rep,block

  inbred predicted.value std.error    status
1     x1      -0.2584834 0.1318171 Estimable
2     x2      -0.2584834 0.1318171 Estimable
3     x3      -0.2584836 0.1318171 Estimable
4     x4      -0.2584839 0.1318171 Estimable
5     x5      -0.2584838 0.1318171 Estimable

上面的这些 std.errors 是针对 predict.value 列的。有没有办法获得 BLUP 或随机效应列的标准误差？

随机效应的 BLUP 是：

 df.ran = df.reml$vcoeff$random
 df.ran
 [1] 1.599984e-06 1.599984e-06 1.599984e-06 1.599984e-06 1.599984e-06
4

1 回答 1

0

好问题。答案并不明显。summary您可以使用以下函数获取 BLUP 的标准错误：

summary(df.reml, coef=TRUE)$coef.random
               solution   std.error       z.ratio
inbred_x1  3.231267e-06 0.001054529  3.064180e-03
inbred_x2  3.354557e-06 0.001054529  3.181093e-03
inbred_x3 -9.275731e-08 0.001054529 -8.796086e-05
inbred_x4 -3.410440e-06 0.001054529 -3.234087e-03
inbred_x5 -3.082627e-06 0.001054529 -2.923225e-03
于 2021-03-12T20:09:46.293 回答