CODEWITHSUNDEEP
rmixed-modelsasreml
FrosyFeet456
FrosyFeet456

Reputation: 349

Standard error for random effect (BLUP) of a mixed effect model in r

Can we calculate standard errors (SE) for BLUPs from a mixed model in R? I am using a custom package called ASReml and it can calculate SE from the predicted value. However I am not sure where to find the SE for the BLUP values.

So far I have something like this...

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

Use asreml to fit the mixed model.

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.

And then we ask for the predicted values:

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

These std.errors above are for the predicted.value column. Is there a way to get standard errors for the BLUP or random effects column?

The BLUPs for the random effects are:

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

Upvotes: 2

Views: 958

Answers (1)

Kevin Wright
Kevin Wright

Reputation: 2697

Good question. The answer is not obvious. You can get the standard errors for BLUPs using the summary function:

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

Upvotes: 0

Related Questions