Блог

Добрый день, я выполняю геномное предсказание (GBLUP) у тетраплоидных видов. Я начал использовать sommer в прошлом году, до изменений. Тогда вы могли использовать функцию mmer с аргументом ZETA, чтобы указать ковариацию случайного члена. Теперь есть спецификация с случайным ~ 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 с доминированием)

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

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

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

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

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

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

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

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

y.hat = Xb Zu

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

y.hat= Xb

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

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

Я надеюсь, что это поможет.

Приветствую, Эдуардо