遇到bug不要慌，都是小場面。

育種資料分析之放飛自我發表於2020-11-17

1. 澄清

今天早上準備跑步，看到外面霧濛濛的，好像是下雨了。

也只是好像下雨了，不一定真下呢，畢竟眼見不一定為實，一定要親身實踐，小馬過河才可以。

所謂不淋雨不知道雨是真實的，不捱揍不知道鐵拳的味道，等到下樓雨滴滴到臉上頭上的那一刻，我就慫了，上去，這天還去跑步就是一個鐵憨憨！

不跑步也不能去睡一個回籠覺，就調bug吧。

2. 事情是這個樣子的

之前培訓GS課程時，有一個HBLUP使用asreml的結果，因為大家大部分都沒有asreml這個軟體，所以僅僅是寫了程式碼，沒有執行。

然後在群裡面收到了報錯的資訊：

我一直都有別人的程式碼報錯不報錯我不在乎，我自己的程式碼報錯我一定解決的態度。

初步判斷：

報錯是沒有報錯，但是方差組分估算為0，肯定是錯誤的：

所以，很有可能是H矩陣的ID定義不一致所致。

我看了一下之前培訓的程式碼，是手動構建H矩陣，然後使用asreml進行ssGBLUP的估算。當時的程式碼如下：

library(asreml)
library(learnasreml)
hinv = write_relation_matrix(Hmat,type="ginv")
hinv = as.matrix(hinv)
str(dd)
attr(hinv,"rowNames") = as.character(dd$ID)
attr(hinv,"INVERSE") = TRUE

3. asreml讀取外部矩陣的規則

    1. 將矩陣變為三元組的逆矩陣形式，這裡用的是我編寫的learnasreml包中的write_relation_matrix函式，將其轉化為廣義逆矩陣，然後變為三元組的稀疏矩陣，對應程式碼library(learnasreml)；hinv = write_relation_matrix(Hmat,type="ginv")
    1. 將三元組轉化為矩陣的形式hinv = as.matrix(hinv)
    1. 用attr函式，將其rowNames變為實際的行名，attr(hinv,"rowNames") = as.character(dd$ID)
    1. 用attr函式，將其INVERSE變為TRUE

4. 如何除錯？

問題出在哪裡？
第三步時，我的程式碼為：attr(hinv,"rowNames") = as.character(dd$ID)，但是Hmat矩陣的行名不一定和dd$ID的ID是一致的，如果不一致，那就是對應關係錯誤，評估方差組分時肯定報錯。

然後我比較一下兩者的異同：

> sum(dd$ID == idh)
[1] 610
> dim(dd)
[1] 2560    4

可以看到，原來共有2560個個體，只有610是對應的，結果不報錯才怪！

5. 更正後結果正常

更正如下：

# attr(hinv,"rowNames") = as.character(dd$ID)
attr(hinv,"rowNames") = as.character(rownames(Hmat))

這樣結果就正常了，執行結果如下：

> mod.as = asreml(phe ~ Sex+Generation, random = ~ vm(ID,hinv), residual = ~ idv(units),workspace="5Gb",data=dd)
Online License checked out Tue Nov 17 07:44:10 2020
Model fitted using the sigma parameterization.
ASReml 4.1.0 Tue Nov 17 07:44:11 2020
Allocating main workspace...done!
          LogLik        Sigma2     DF     wall    cpu
 1     -23152.27           1.0   2554 07:44:47   34.5
 2     -19423.25           1.0   2554 07:44:55    7.1
 3     -15276.23           1.0   2554 07:45:02    7.0
 4     -12315.15           1.0   2554 07:45:08    6.9
 5     -10608.41           1.0   2554 07:45:15    7.0
 6     -10008.30           1.0   2554 07:45:22    7.0
 7      -9857.03           1.0   2554 07:45:29    7.0
 8      -9840.99           1.0   2554 07:45:36    7.0
 9      -9840.70           1.0   2554 07:45:43    7.0
10      -9840.70           1.0   2554 07:45:50    7.0
> summary(mod.as)$varcomp
             component std.error   z.ratio bound %ch
vm(ID, hinv)  53.27475  14.99226  3.553484     P   0
units!units  763.88704  23.60964 32.354877     P   0
units!R        1.00000        NA        NA     F   0
> vpredict(mod.as,h2 ~ V1/(V1+V2))
     Estimate         SE
h2 0.06519487 0.01786531

可以看到asreml執行了1分鐘左右，而sommer執行了10分鐘。

怎麼講？收費的就是香！

6. 調bug心得

  • 先看資料有無問題
  • 再看模型程式碼有無問題
  • 最後看結果有無問題

要有資料敏感性，事出反常必有妖，有時候結果沒有報錯不一定是真的沒有錯誤，要看結果是否正常，不正常肯定有問題。

當然，自己寫的bug，含著淚也要改好。遇到bug不要慌，都是小場面。

相關文章