Don't panic when you encounter a bug. It's all a small scene.

1. Clarification

I was going to run this morning. I saw it was foggy outside. It seemed to rain.

It's just like it's raining. It doesn't have to be true. After all, seeing is not necessarily true. You must practice it yourself and the pony can cross the river.

The so-called don't get caught in the rain, don't know the rain is real, don't get beaten, don't know the taste of iron fist, wait until the moment when the rain drops on my face downstairs, I counselled, go up, and go running this day is an iron Han!

If you don't run, you can't go to sleep. Just adjust the bug.

2. Things are like this

In the previous GS training course, there was a result of using asreml in HBLUP. Because most of us didn't have asreml software, we just wrote code and didn't run it.

Then I received the error message in the group:

I always have other people's code reporting errors. I don't care if I don't report errors. I will solve my own code reporting errors.

Preliminary judgment:

There is no error in error reporting, but the estimation of variance component is 0, which must be wrong:

Therefore, it is likely that the ID definition of H matrix is inconsistent.

I looked at the code of the previous training, which is to manually build the H matrix, and then use asreml to estimate ssGBLUP. The code at that time was as follows:

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 rules for reading external matrix

    1. Change the matrix into the inverse matrix form of triples. Here is the write in the learnasreml package I wrote_ relation_ Matrix function, which is transformed into a generalized inverse matrix, and then into a sparse matrix of triples, corresponding to the code library(learnasreml); hinv = write_relation_matrix(Hmat,type="ginv")
    1. Convert triples to matrix form hinv = as matrix(hinv)
    1. Use the attr function to change its rowNames to the actual row name. Attr (hinv, "rowNames") = as character(dd$ID)
    1. Use the attr function to change its invert to TRUE

4. How to debug?

What's the problem?
In the third step, my code is: attr (hinv, "rownames") = as Character (dd$ID), but the row name of Hmat matrix may not be consistent with the ID of dd$ID. if it is inconsistent, it is the corresponding relationship error, and an error must be reported when evaluating the variance component.

Then I compare the similarities and differences between the two:

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

It can be seen that there are 2560 individuals in the original, and only 610 are corresponding. It's strange that the result doesn't report an error!

5. The corrected result is normal

Correct as follows:

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

In this way, the results are normal, and the operation results are as follows:

> 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

You can see that asreml ran for about a minute, while sommer ran for 10 minutes.

Well? The charge is incense!

6. Experience of de bug ging

  • Let's see if there is any problem with the data
  • Let's see if there is any problem with the model code
  • Finally, see if there is any problem with the result

There must be data sensitivity. If something goes wrong, there must be demons. Sometimes the result does not report an error, it does not necessarily mean there is no error. It depends on whether the result is normal. There must be a problem if it is abnormal.

Of course, the bugs you write should be corrected with tears. Don't panic when you encounter a bug. It's all a small scene.

Posted by newbie5050 on Sat, 07 May 2022 09:59:20 +0300