1. 澄清

今天早上准备跑步,看到外面雾蒙蒙的,好像是下雨了。

也只是好像下雨了,不一定真下呢,毕竟眼见不一定为实,一定要亲身实践,小马过河才可以。

所谓不淋雨不知道雨是真实的,不挨揍不知道铁拳的味道,等到下楼雨滴滴到脸上头上的那一刻,我就怂了,上去,这天还去跑步就是一个铁憨憨!

不跑步也不能去睡一个回笼觉,就调bug吧。

2. 事情是这个样子的

之前培训GS课程时,有一个HBLUP使用asreml的结果,因为大家大部分都没有asreml这个软件,所以仅仅是写了代码,没有运行。

然后在群里面收到了报错的信息:

遇到bug不要慌,都是小场面。_职场生活

我一直都有别人的代码报错不报错我不在乎,我自己的代码报错我一定解决的态度。

初步判断:

报错是没有报错,但是方差组分估算为0,肯定是错误的:
遇到bug不要慌,都是小场面。_职场生活_02
所以,很有可能是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不要慌,都是小场面。