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读取外部矩阵的规则
-
- 将矩阵变为三元组的逆矩阵形式,这里用的是我编写的learnasreml包中的
write_relation_matrix
函数,将其转化为广义逆矩阵,然后变为三元组的稀疏矩阵,对应代码library(learnasreml);hinv = write_relation_matrix(Hmat,type="ginv")
- 将矩阵变为三元组的逆矩阵形式,这里用的是我编写的learnasreml包中的
-
- 将三元组转化为矩阵的形式
hinv = as.matrix(hinv)
- 将三元组转化为矩阵的形式
-
- 用
attr
函数,将其rowNames
变为实际的行名,attr(hinv,"rowNames") = as.character(dd$ID)
- 用
-
- 用
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不要慌,都是小场面。