Пакет Sommer: Mmer против MEMMA - PullRequest
Новая
       11

Пакет Sommer: Mmer против MEMMA

0 голосов
/ 09 апреля 2019

Добрый день, я выполняю геномное предсказание (GBLUP) у тетраплоидных видов.Я начал использовать зоммер в прошлом году до изменений.Тогда вы могли использовать функцию mmer с аргументом ZETA , чтобы указать ковариацию случайного члена.Теперь есть спецификация с random ~ vs (...)

Я сравнил это с функцией MEMMA , которая использует ZETA аргумент, и, насколько я понимаю, он используется внутри Mmer.Однако по какой-то причине я получаю разные ответы, например 

    library(devtools)
    install_github("covaruber/sommer")
    library(sommer)
     install.packages("AGHmatrix")
    library(AGHmatrix) 



    data(DT_polyploid)
     # ####=========================================####
     # ####### convert markers to numeric format
     # ####=========================================####
      numo <- atcg1234(data=GT, ploidy=4);

     # ###=========================================####
     # ###### plants with both genotypes and phenotypes
     # ###=========================================####
      common <- intersect(DT$Name,rownames(numo$M))
     #
     # ###=========================================####
     # ### get the markers and phenotypes for such inds
     # ###=========================================####
     marks <- numo$M[common,]; marks[1:5,1:5]
      DT2 <- DT[match(common,DT$Name),];

     # ###=========================================#### 

     G.A <- Gmatrix(marks, method="VanRaden", ploidy=4, missingValue=NA, impute.method = TRUE)
     G.D <- Gmatrix(marks, method="Endelman", ploidy=4, missingValue=NA, impute.method = TRUE)    

     T.pheno <- DT2[,c(1,9)]
     T.pheno$Name <- as.factor(T.pheno$Name)
     T.pheno$Name2 <- T.pheno$Name

    set.seed(1892) 
    rrn <- sample(1:187, 50, replace = F)
    T.pheno$sucrose[rrn] <- NA

    ans.A <-  mmer(sucrose ~ 1,
                   random=~vs(Name, Gu=G.A),
                   rcov = ~units,
                   data=T.pheno)


    cor(ans.A$U$`u:Name`$sucrose[which(is.na(T.pheno$sucrose))], 
        DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")

[1] 0,2205831

Затем применяю тот же анализ с использованием MEMMA 

Z1 <- diag(length(T.pheno$sucrose))
Z2 <- diag(length(T.pheno$sucrose))

ETA.A <- list(list(Z=Z1, K=G.A) )
ETA.AD <- list(list(Z=Z1, K=G.A), list(Z=Z2, K=G.D) )


ans.A <- MEMMA(Y=T.pheno$sucrose, ZETA=ETA.A)

cor(ans.A$fitted.y[which(is.na(T.pheno$sucrose))], 
    DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")

[1] 0,2778689

Почему есть разница?Для добавления матриц доминирования и эпистазиса я нахожу формулировку в MEMMA более простой и менее подверженной ошибкам, и, конечно, я воодушевлен более высокой точностью прогнозирования.Во-вторых, функция MEMMA обеспечивает установленные значения, которые находятся в том же масштабе и, следовательно, сопоставимы с наблюдаемыми значениями.Однако MEMMA довольно медленная ... и все же должно быть объяснение разницы в точности.Ниже приведен код и результаты для mmer и MEMMA при включении доминирования. 

ans.AD <-  mmer(sucrose ~ 1,
                random= ~ vs(Name, Gu=G.A) + vs(Name2, Gu=G.D),
                rcov = ~ units,
                data=T.pheno)

cor(ans.AD$U$`u:Name`$sucrose[which(is.na(T.pheno$sucrose))], 
    DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")


ans.AD <- MEMMA(Y=T.pheno$sucrose, ZETA=ETA.AD)

cor(ans.AD$fitted.y[which(is.na(T.pheno$sucrose))], 
    DT2$sucrose[which(is.na(T.pheno$sucrose))], use = "pairwise")

[1] 0.2357571 (mmer с доминированием)

[1] 0.2785493 (MEMMA с доминированием)

Спасибо, что уделили время

1 Ответ

0 голосов
/ 09 апреля 2019

Причины, по которым функция MEMMA не является частью функции mmer, а доступна только как дополнительная функция:

1) компоненты дисперсии не оцениваются должным образом при переходе к более чем одному случайному эффекту, отличному от ошибки

2) Невозможность использования специальных ковариационных структур в невязках (опять же, в той же задаче, превышающей 1 vc)

3) Помимо самого простого сценария, этот алгоритм ужасно медленен для оценки ВК.

Теперь причина, по которой Ньютон-Рафсон (NR) и Средняя информация (AI) (единственные методы, доступные в функции mmer ()) дают результаты, отличные от MEMMA, должна быть связана с тем, что для более чем один vc оценки REML MEMMA обычно выходят за пределы пространства параметров. Этот алгоритм был первоначально разработан для одного vc рядом с vc для ошибки. Я бы сказал, что NR и AI, скорее всего, дают вам правильные значения VC.

Кроме того, то, что одно значение дает лучшие результаты для MEMMA, не означает, что, скажем, через 200 итераций это будет правдой. Я бы выполнил некоторую перекрестную проверку, чтобы проверить, верно ли это. И даже если это так, NR / AI менее склонны выходить за пределы пространства параметров, давая вам правильные значения.

О слоте установленных значений. МЕММА по-прежнему возвращает подогнанные значения как

y.hat = Xb + Zu

тогда как в результате выполнения функции mmer () слот с установленными значениями имеет вид:

y.hat = Xb

Если вы хотите использовать прогнозы в исходном масштабе, вы должны построить их самостоятельно, построив Xb и Zu, или вы можете использовать доступную функцию предсказания ().

Также, чтобы показать вам, как основание того, какой алгоритм использовать на основе этой точности может быть очень субъективным, вы можете попробовать для того же примера, но теперь укажите в функции mmer () аргумент na.method.Y = "include" и Посмотрите, как ПА поднимается до 0,27. Урок, небольшие изменения в PA не так важны с моей точки зрения. Бумага, показывающая увеличение PA на 5 единиц (0,05), с моей точки зрения, не на том пути, но просто мое скромное мнение.

Надеюсь, это поможет.

Ура, Эдуардо

...