Get emmeans in sommer's mmer?

Question

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!

Answer 1

It appears to be possible. Here is one of the package's examples:

library(sommer) # Version 4.1.2
data(DT_cornhybrids)
DT <- DT_cornhybrids
DTi <- DTi_cornhybrids
GT <- GT_cornhybrids
hybrid2 <- DT
A <- GT
K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]
K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]
S <- kronecker(K1, K2, make.dimnames=TRUE)   

ans <- mmer(Yield ~ Location, 
             random = ~ vs(GCA1,Gu=K1) + vs(GCA2,Gu=K2) + vs(SCA,Gu=S),
             rcov=~units,
            data=hybrid2)
summary(ans)

## ...
## Fixed effects:
##   Trait      Effect   Estimate Std.Error    t.value
## 1 Yield (Intercept)  1.379e+02     1.962  7.031e+01
## 2 Yield   Location2  1.776e-14     2.099  8.461e-15
## 3 Yield   Location3  7.835e+00     2.099  3.732e+00
## 4 Yield   Location4 -9.097e+00     2.099 -4.333e+00
## ...

The returned object has elements $Beta and $VarBeta which return the fixed effects and covariances thereof. We can create a reference grid using emmeans::qdrg():

rg <- qdrg(~ Location, data = hybrid2, coef = ans$Beta$Estimate, 
    vcov = ans$VarBeta)
rg
## 'emmGrid' object with variables:
##    Location = 1, 2, 3, 4

emmeans(rg, trt.vs.ctrl1 ~ Location)
## $emmeans
##  Location emmean   SE  df asymp.LCL asymp.UCL
##  1           138 1.96 Inf       134       142
##  2           138 1.96 Inf       134       142
##  3           146 1.96 Inf       142       150
##  4           129 1.96 Inf       125       133

## Confidence level used: 0.95 

## $contrasts
##  contrast estimate  SE  df z.ratio p.value
##  2 - 1        0.00 2.1 Inf  0.000  1.0000 
##  3 - 1        7.84 2.1 Inf  3.732  0.0006 
##  4 - 1       -9.10 2.1 Inf -4.333  <.0001 

## P value adjustment: dunnettx method for 3 tests 

The fact that the EMM for location 1, and its SE, matches the summary() intercept, and that the remaining regression coefficients and SEs match the contrast results, is reassuring.

See the documentation for qdrg for more details.

Answer 2

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.