R的诡异错误 (3)
A pitfall in R (3)
下面的代码有两行，第一行计算Genotype的Minor Allele Frequency，第二行计算Weight matrix （用SKAT paper建议的Beta分布的weight）。
下面的代码有问题吗？
# G is a sample by marker matrix maf <- colMeans(G) / 2 w <- diag(dbeta(maf, 1, 25))
粗粗一看，上面的代码的思路是正确的，不过这里有两个坑。
第一个坑是缺失数据（Missing data），如果Genotype矩阵G有缺失数据，那么maf和w都会有缺失的数据；
第二个坑比较难看出来：当G只有一列的时候，maf也会是一个标量（Scalar），同样dbeta函数的返回值是标量。
你可能想diag函数作用在一个标量上，其返回值是一个1×1的矩阵。
但是，diag()的返回的是一个….不确定的结果，看个例子吧：
diag(3) => 3x3 单位矩阵 diag(3.9) => 还是3x3 单位矩阵 diag(0) => 是0x0 矩阵 diag(-1) => 报错 diag(c(1, 2, 3) => 3x3 矩阵， 对角线是1，2，3
总之，第二个坑告诉我们，如果dbeta返回的是一个长度为1的数值，别指望w一定是1×1的矩阵。
继续吐槽：
diag()这个函数设计的不好用。在R内部，标量就是被实现为长度为1的向量。
在现在的设计下，diag要求使用者同时考虑变量的类型（矩阵，向量，缺失）以及变量的长度和值（向量长度为1且大于等于0，向量长度为1且小于0，向量长度大于1）。
对于一个适应了使用有类型系统的编程语言的程序员，这种要求会增加使用者的脑力负担，同时阅读者在语义上也很容易出现混淆。
我建议的解决方案是：diag(向量)永远输出一个对角线是该向量的矩阵，同时定义一个新函数(可以叫diag_len,类比seq_len；或者叫eye，类比Matlab里的eye函数）专门用来生成单位矩阵。
最后总结一下，我觉得这样写比较
# G is a sample by marker matrix maf <- colMeans(G, na.rm = TRUE) / 2 m <- length(maf) w <- matrix(0, nr = m, nc = m) diag(w) <- dbeta(maf, 1, 25) # 或者：w <- diag(dbeta(maf, 1, 25), nrow = m, ncol = m)