【发布时间】：2020-01-21 17:30:13
【问题描述】：
我有一个数字 data.frame df 有 134946 行 x 1938 列。
99.82% 的数据是NA。
对于每一对（不同的）列"P1" 和"P2"，我需要找出哪些行的值都不是NA，然后对这些行进行一些操作（线性模型）。
我写了一个脚本来做这个，但是看起来很慢。
This post 似乎在讨论一项相关任务，但我无法立即看到它是否或如何适应我的情况。
借用那篇文章中的例子：
set.seed(54321)
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
df <- as.data.frame(dat)
df_ps <- names(df)
N_ps <- length(df_ps)
我的脚本是：
tic = proc.time()
out <- do.call(rbind,sapply(1:(N_ps-1), function(i) {
if (i/10 == floor(i/10)) {
cat("\ni = ",i,"\n")
toc = proc.time();
show(toc-tic);
}
do.call(rbind,sapply((i+1):N_ps, function(j) {
w <- which(complete.cases(df[,i],df[,j]))
N <- length(w)
if (N >= 5) {
xw <- df[w,i]
yw <- df[w,j]
if ((diff(range(xw)) != 0) & (diff(range(yw)) != 0)) {
s <- summary(lm(yw~xw))
o <- c(i,j,N,s$adj.r.squared,s$coefficients[2],s$coefficients[4],s$coefficients[8],s$coefficients[1],s$coefficients[3],s$coefficients[7])} else {
o <- c(i,j,N,rep(NA,7))
}
} else {o <- NULL}
return(o)
},simplify=F))
}
,simplify=F))
toc = proc.time();
show(toc-tic);
这在我的机器上大约需要 10 分钟。
您可以想象当我需要处理更大（尽管更稀疏）的数据矩阵时会发生什么。我从来没有完成计算。
问题：您认为这可以更有效地完成吗？
问题是我不知道哪些操作需要更多时间（df 的子集，在这种情况下我会删除重复项？附加矩阵数据，在这种情况下我会创建一个平面向量然后转换它到最后的矩阵？...）。
谢谢！
编辑跟进我的帖子
如 minem 所示，此计算的速度很大程度上取决于计算线性回归参数的方式。因此，更改该部分是最重要的事情。
我自己的进一步试验表明：1) 必须将sapply 与do.call(rbind 结合使用，而不是任何平面向量来存储数据（我仍然不确定为什么 - 我可能会就此单独发表一篇文章); 2) 在我正在处理的原始矩阵上，使用每个i 迭代开始时可用的x 向量信息，比本示例中的矩阵更稀疏且nrows/ncolumns 比率更大在每次 j 迭代开始时减少 y 向量，速度提高了几个数量级，即使与 minem 的原始脚本相比，这已经比我上面的要好得多。
我认为优势来自于先验过滤掉许多行，从而避免了对非常长的向量进行昂贵的 xna & yna 操作。
修改后的脚本如下：
set.seed(54321)
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.90)] = NA
df <- as.data.frame(dat)
df_ps <- names(df)
N_ps <- length(df_ps)
tic = proc.time()
naIds <- lapply(df, function(x) !is.na(x))
dl <- as.list(df)
rl <- sapply(1:(N_ps - 1), function(i) {
if ((i-1)/10 == floor((i-1)/10)) {
cat("\ni = ",i,"\n")
toc = proc.time();
show(toc-tic);
}
x <- dl[[i]]
xna <- which(naIds[[i]])
rl2 <- sapply((i + 1):N_ps, function(j) {
y <- dl[[j]][xna]
yna <- which(naIds[[j]][xna])
w <- xna[yna]
N <- length(w)
if (N >= 5) {
xw <- x[w]
yw <- y[yna]
if ((min(xw) != max(xw)) && (min(yw) != max(yw))) {
# extracts from lm/lm.fit/summary.lm functions
X <- cbind(1L, xw)
m <- .lm.fit(X, yw)
# calculate adj.r.squared
fitted <- yw - m$residuals
rss <- sum(m$residuals^2)
mss <- sum((fitted - mean(fitted))^2)
n <- length(m$residuals)
rdf <- n - m$rank
# rdf <- df.residual
r.squared <- mss/(mss + rss)
adj.r.squared <- 1 - (1 - r.squared) * ((n - 1L)/rdf)
# calculate se & pvals
p1 <- 1L:m$rank
Qr <- m$qr
R <- chol2inv(Qr[p1, p1, drop = FALSE])
resvar <- rss/rdf
se <- sqrt(diag(R) * resvar)
est <- m$coefficients[m$pivot[p1]]
tval <- est/se
pvals <- 2 * pt(abs(tval), rdf, lower.tail = FALSE)
res <- c(m$coefficients[2], se[2], pvals[2],
m$coefficients[1], se[1], pvals[1])
o <- c(i, j, N, adj.r.squared, res)
} else {
o <- c(i,j,N,rep(NA,7))
}
} else {o <- NULL}
return(o)
}, simplify = F)
do.call(rbind, rl2)
}, simplify = F)
out2 <- do.call(rbind, rl)
toc = proc.time();
show(toc - tic)
例如试试nr=100000; nc=100。
我可能应该提到我尝试使用索引，即：
naIds <- lapply(df, function(x) which(!is.na(x)))
然后显然通过交集生成w：
w <- intersect(xna,yna)
N <- length(w)
然而，这比上面的要慢。
【问题讨论】：
-
首先，如果你只有数字数据，我真的建议存储为矩阵，而不是数据框。其次，如果你的数据非常稀疏，为什么不将其存储为稀疏矩阵呢？
-
好的，这两种更改都是可能的。您认为这将有助于加快进程吗？
标签： r performance