1
$\begingroup$

I am using asreml-R to model genotypic effects of crop field trials and I am confused on how to get best linear unbiased estimates for crop varieties of the model.

I've found two different ways how people do this:

The second approach, predict, itself calls update.asreml which leads to some iterative algorithm, it seems to partially refit the model even though it was already fitted to the data before. Even though the result is exactly the same as I get when just reading the coefficients from the data structure of the model using coef, which completes virtually instantaneously. I have shown this with a MVCE below. This confuses me, and leads to my three questions:

  1. Are BLUEs the coefficients/effects of the mixed effects model in all cases, if no, what is the difference?
  2. If yes, why is predict() partially refitting the asreml model, which may take several minutes for ~8000 observations (where a complete fit took ~30 minutes). What needs to be done beyond the point of looking up the coefficients?
  3. So when can I use coef() to get the BLUEs, and when is predict() nessecary?

MCVE

library(asreml)
data(nin89, package = "asreml")

nin89 <- nin89[!is.na(nin89$Rep),]
asr <- asreml(fixed = yield ~ Variety, random =~ Rep, data = nin89)
## ASReml: Wed Sep  5 14:21:56 2018
## 
##      LogLik         S2      DF      wall     cpu
##    -454.8069     50.3285   168  14:21:56     0.0
##    -454.6631     50.1197   168  14:21:56     0.0
##    -454.5323     49.8682   168  14:21:56     0.0
##    -454.4718     49.6374   168  14:21:56     0.0
##    -454.4691     49.5854   168  14:21:56     0.0
##    -454.4691     49.5824   168  14:21:56     0.0
##    -454.4691     49.5824   168  14:21:56     0.0
## 
## Finished on: Wed Sep  5 14:21:56 2018
##  
## LogLikelihood Converged 

# BLUEs Method 1 (coef)
blues1 <- coef(asr)$fixed
blues1 <- blues1 + blues1["(Intercept)",]
blues1 <- blues1[rownames(blues1) != "(Intercept)",,drop = FALSE]
blues1 <- data.frame(Variety = rownames(blues1), predicted.value = blues1[,"effect"], stringsAsFactors = FALSE, row.names = NULL)
blues1$Variety <- gsub("^Variety_","",blues1$Variety)

# BLUEs Method 2 (predict)
blues2 <- predict(asr, classify = "Variety", maxiter = 1)$prediction$pvals
## ASReml: Wed Sep  5 14:21:56 2018
## 
##      LogLik         S2      DF      wall     cpu
##    -454.4691     49.5824   168  14:21:56     0.0
## 
## Finished on: Wed Sep  5 14:21:56 2018

# Compare
blues <- merge(blues1, blues2, by = "Variety", all = TRUE)
blues[is.na(blues$predicted.value.x) | is.na(blues$predicted.value.y)]

plot(predicted.value.y ~ predicted.value.x, data = blues)

dot plot showing each BLUE has same value for both methods

$\endgroup$
2
$\begingroup$

Good question. For trivial predictions with no missing values in your data, you could just use the result of coef indeed but if you want extra things such as the prediction error variances, estimability etc, you'll need to use predict.asreml. For details of the algorithm see Gilmour et al. (2004) An efficient computing strategy for prediction in mixed linear models.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.