加快成对观察计数R中

Question

我有一个数据集，其中的每个条目的测量子集随机丢失：

dat <- matrix(runif(100), nrow=10) 
rownames(dat) <- letters[1:10] 
colnames(dat) <- paste("time", 1:10) 
dat[sample(100, 25)] <- NA 

我很感兴趣，在此数据集计算每一行之间的相关性（即AA ，ab，ac，ad，...）。但是，我想通过在结果相关矩阵中将其值设置为NA来排除少于5个非成对非NA观测值的相关性。

目前，我这样做如下：

cor <- cor(t(dat), use = 'pairwise.complete.obs') 
names <- rownames(dat) 
filter <- sapply(names, function(x1) sapply(names, function(x2) 
    sum(!is.na(dat[x1,]) & !is.na(dat[x2,])) < 5)) 
cor[filter] <- NA 

然而，这种操作的实际数据集包含> 1000项非常缓慢。

是否有方法可以基于矢量化方式中的非NA成对观察数来过滤单元格，而不是在嵌套循环中？

Answer 1

您可以使用矩阵方法计算非NA成对观测值的数量。

让我们使用这个数据生成代码。我使数据更大并增加了更多的NAs。

nr = 1000; 
nc = 900; 
dat = matrix(runif(nr*nc), nrow=nr) 
rownames(dat) = paste(1:nr) 
colnames(dat) = paste("time", 1:nc) 
dat[sample(nr*nc, nr*nc*0.9)] = NA 

然后你过滤代码正在85秒

tic = proc.time() 
names = rownames(dat) 
filter = sapply(names, function(x1) sapply(names, function(x2) 
    sum(!is.na(dat[x1,]) & !is.na(dat[x2,])) < 5)); 
toc = proc.time(); 
show(toc-tic); 
# 85.50 seconds 

我的版本创建与原始数据值1用于非NAS的矩阵。然后使用矩阵乘法计算配对非NAs的数量。它跑了几分之一秒。

tic = proc.time() 
NAmat = matrix(0, nrow = nr, ncol = nc) 
NAmat[ !is.na(dat) ] = 1; 
filter2 = (tcrossprod(NAmat) < 5) 
toc = proc.time(); 
show(toc-tic); 
# 0.09 seconds 

简单的检查，显示结果是相同的：

all(filter == filter2) 
# TRUE