rmixed-modelslsmeansemmeans

Get emmeans in sommer's mmer?

Additional keywords: Best linear unbiased Estimator (BLUE), adjusted means, mixed model, fixed effects, linear combination, contrast, R

After fitting a model with mmer() of the sommer package - is it possible to obtain estimated marginal means (emmeans()) / least squares means (LS-means) from an mmer object? Maybe similar to the predict() function with ASReml-R v3?

Actually, I would I want multiple things and maybe it is clearer to ask for them separately:

  1. The emmeans themselves and their
  2. Standard errors (s.e.)
  3. as a column next to the means of each level
  4. variance-covariance matrix of emmeans (see predict(..., vcov=T))
  5. Contrasts between means and their
  6. Standard errors of a difference (s.e.d.)
  7. All pairwise differences between means, preferably with a post hoc test (see emmeans(mod, pairwise ~ effect, adjust="Tukey")
  8. S.e.d. matrix (see predict(..., sed=T))
  9. Minimum, average and maximum s.e.d.
  10. Custom contrasts

So yeah, basically a mix of predict() and emmeans() would be the goal here.

Thanks in advance!

Solution

  • In sommer >= 3.7 the predict function has been implemented to obtain predictions for either fixed or random effects the way asreml does. It takes a model and the classify argument to know which arguments to use for aggregating the hypertable and come up with the right standard errors. For example:

    data(DT_cpdata)
#### create the variance-covariance matrix
A <- A.mat(GT) # additive relationship matrix
#### look at the data and fit the model
head(DT)
mix1 <- mmer(Yield~1,
              random=~vs(id,Gu=A)
                      + Rowf + Colf,
              rcov=~units,
              data=DT)
summary(mix1)

preds <- predict(mix1,classify = "id")
> head(preds$predictions)
    id predicted.value.Yield standard.errors.Yield
1 P003             111.15400              28.16363
2 P004             135.21958              29.81544
3 P005             109.72743              29.68574
4 P006             144.98582              27.99627

preds <- predict(mix1,classify = "Rowf")
> head(preds$predictions)
  Rowf predicted.value.Yield standard.errors.Yield
1    1              81.71645              23.22997
2    2              96.79625              22.92514
3    3             128.89043              22.64216
4    4             132.65795              22.73903

    and so on ... The RtermsToForce and FtermsToForce arguments can be used to force the use of specific fixed or random terms in the predictions. Customized contrasts I guess for the next version.