Получить emmeans в mmer2 Соммера?

Question

Дополнительные ключевые слова: Лучший линейный несмещенный оценщик (СИНИЙ), скорректированное среднее, смешанная модель, фиксированные эффекты, линейная комбинация, контрастность, R

После установки модели с mmer2() из sommer пакет - возможно ли получить расчетное предельное среднее значение (emmeans()) / наименьших квадратов означает (LS-значит) от объекта mmer2? Может быть, похоже на функцию predict() с ASReml-R v3 ?

На самом деле, я бы хотел несколько вещей, и, возможно, было бы яснее спросить их отдельно:

1. Сами эммеи и их
2. Стандартные ошибки (с. Е.)
3. как столбец рядом со средством каждого уровня
4. дисперсионно-ковариационная матрица значений (см. predict(..., vcov=T))
5. Контрасты между средними и их
6. Стандартные погрешности разницы (s.e.d.)
7. Все попарные различия между средними значениями, предпочтительно с помощью теста post hoc (см. emmeans(mod, pairwise ~ effect, adjust="Tukey")
8. S.e.d. матрица (см. predict(..., sed=T))
9. Минимальная, средняя и максимальная s.e.d.
10. Пользовательские контрасты

Так что да, в общем, цель здесь - сочетание predict() и emmeans().

Заранее спасибо!

Answer 1

В sommer> = 3.7 была реализована функция предсказания для получения предсказаний для фиксированных или случайных эффектов, как это делает asreml. Требуется модель и аргумент классификации, чтобы узнать, какие аргументы использовать для агрегирования гипертаблицы, и найти правильные стандартные ошибки. Например:

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

и т. Д. ... Аргументы RtermsToForce и FtermsToForce могут использоваться для принудительного использования определенных фиксированных или случайных членов в предсказаниях. Индивидуальные контрасты, я думаю, для следующей версии.

Answer 2

Это представляется возможным. Вот один из примеров пакета:

data(cornHybrid)
hybrid2 <- cornHybrid$hybrid
A <- cornHybrid$K
K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]
K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]
S <- kronecker(K1, K2, make.dimnames=TRUE)   

ans <- mmer2(Yield ~ Location, 
             random = ~ g(GCA1) + g(GCA2) + g(SCA), 
             rcov = ~ units,
             G=list(GCA1=K1, GCA2=K2, SCA=S),
             data=hybrid2)
summary(ans)
## ...
## Fixed effects:
##
## $`Yield`
##                  Estimate Std. Error       t value
## (Intercept)  1.379120e+02   1.961220  7.031951e+01
## Location2    1.776357e-13   2.099498  8.460865e-14
## Location3    7.835337e+00   2.099498  3.732005e+00
## Location4   -9.097455e+00   2.099498 -4.333158e+00
## ...

Возвращенный объект имеет элементы $beta.hat и $Var.beta.hat, которые возвращают фиксированные эффекты и их ковариации. (Последний находится в структурированной форме, которую необходимо привести к стандарту matrix.) Мы можем создать опорную сетку, используя emmeans::qdrg():

rg <- qdrg(~ Location, data = hybrid2, coef = ans$beta.hat, 
    vcov = as.matrix(ans$Var.beta.hat))
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        137.9120 1.96122 Inf  134.0681  141.7559
##  2        137.9120 1.96122 Inf  134.0681  141.7559
##  3        145.7473 1.96122 Inf  141.9034  149.5913
##  4        128.8146 1.96122 Inf  124.9706  132.6585
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast      estimate       SE  df z.ratio p.value
##  2 - 1     1.776357e-13 2.099498 Inf   0.000  1.0000
##  3 - 1     7.835337e+00 2.099498 Inf   3.732  0.0006
##  4 - 1    -9.097455e+00 2.099498 Inf  -4.333  <.0001
##
## P value adjustment: dunnettx method for 3 tests

Тот факт, что EMM для местоположения 1 и его SE соответствует перехвату summary() и что оставшиеся коэффициенты регрессии и SE соответствуют результатам контрастности, обнадеживает.

См. Документацию для qdrg для получения более подробной информации.