Можем ли мы вычислить стандартные ошибки (SE) из функции предсказания в r? Я использую пользовательский пакет с именем ASReml, и они могут вычислить SE из прогнозируемого значения, однако я не уверен, где найти SE для значения BLUP. Мои BLUPs происходят от случайных эффектов, поэтому получение стандартных ошибок от случайных эффектов будет таким же, как для 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
inbred trait1 block rep
1 x1 -1.232811824 1 1
2 x2 -0.640402987 1 1
3 x3 1.095322383 1 1
4 x4 -0.005583372 1 1
5 x5 -0.866967153 1 1
6 x1 -0.777962482 2 1
7 x2 0.506640166 2 1
8 x3 -0.685157660 2 1
9 x4 0.332986177 2 1
10 x5 0.473803168 2 1
11 x1 0.620589862 1 1
12 x2 -0.625389932 1 1
13 x3 -1.182449415 1 1
14 x4 -0.329339345 1 1
15 x5 -0.375746779 1 1
16 x1 -0.848761206 2 1
17 x2 0.187828695 2 1
18 x3 2.050570955 2 1
19 x4 -0.136447699 2 1
20 x5 -0.842616216 2 1
21 x1 0.404040616 1 2
22 x2 -0.284649534 1 2
23 x3 1.603720308 1 2
24 x4 -0.123323447 1 2
25 x5 -0.826282508 1 2
26 x1 -0.595282750 2 2
27 x2 0.366978232 2 2
28 x3 -0.739760317 2 2
29 x4 -1.613508072 2 2
30 x5 -0.113346413 2 2
31 x1 0.609426601 1 2
32 x2 -0.430816036 1 2
33 x3 -0.167768992 1 2
34 x4 0.471079912 1 2
35 x5 -0.448943184 1 2
36 x1 0.478253280 2 2
37 x2 0.155989507 2 2
38 x3 -0.082594336 2 2
39 x4 -0.136350248 2 2
40 x5 -0.474254090 2 2
>
df.reml<-asreml(trait1 ~ rep + block,
random = ~ inbred,
data=df)
Model fitted using the gamma parameterization.
ASReml 4.1.0 Thu Feb 27 00:32:42 2020
LogLik Sigma2 DF wall cpu
1 -13.1497 0.561630 37 00:32:42 0.0
2 -12.7446 0.579377 37 00:32:42 0.0 (1 restrained)
3 -12.7013 0.583738 37 00:32:42 0.0 (1 restrained)
4 -12.6989 0.584042 37 00:32:42 0.0 (1 restrained)
5 -12.6987 0.584061 37 00:32:42 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 Thu Feb 27 00:32:42 2020
LogLik Sigma2 DF wall cpu
1 -12.6987 0.584063 37 00:32:42 0.0
2 -12.6987 0.584063 37 00:32:42 0.0
3 -12.6987 0.584063 37 00:32:42 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.
- rep evaluated at average value of 1.500000
- block evaluated at average value of 1.500000
inbred predicted.value std.error status
1 x1 -0.1307322 0.1208373 Estimable
2 x2 -0.1307321 0.1208373 Estimable
3 x3 -0.1307316 0.1208373 Estimable
4 x4 -0.1307322 0.1208373 Estimable
5 x5 -0.1307327 0.1208373 Estimable
$avsed
overall
0.000488976
>
df.ran = df.reml$vcoeff$random
df.ran
[1] 3.236289e-06 3.236289e-06 3.236289e-06 3.236289e-06 3.236289e-06
>
Эти std.errors предназначены для столбца прогнозируемого значения и не связаны со столбцом df.ran, которые являются значениями BLUP из случайные эффекты. Есть ли способ получить стандартные ошибки для столбца BLUP или случайных эффектов, чтобы соответствовать с df.ran? На самом деле я бы включил столбец для стандартной ошибки случайного эффекта, но я не знаю, как она будет выглядеть, поскольку не знаю, как извлечь ее из этого набора данных.
cbind(df.pred$pvals,df.ran)
inbred predicted.value std.error status df.ran
1 x1 -0.1307322 0.1208373 Estimable 3.236289e-06
2 x2 -0.1307321 0.1208373 Estimable 3.236289e-06
3 x3 -0.1307316 0.1208373 Estimable 3.236289e-06
4 x4 -0.1307323 0.1208373 Estimable 3.236289e-06
5 x5 -0.1307327 0.1208373 Estimable 3.236289e-06