У меня есть число c data.frame df с 134946 строками x 1938 столбцами.
99,82% данных составляют NA.
Для каждой пары (отдельных) столбцов "P1" и "P2", мне нужно найти, какие строки имеют значения не NA для обоих, а затем выполнить некоторые операции над этими строками (линейная модель).
Я написал скрипт, который делает это, но, похоже, довольно медленно.
В этом посте , похоже, обсуждается связанная задача, но я не могу сразу понять, может ли она быть адаптирована к моему делу или как она может быть адаптирована.
Заимствование пример из этого поста:
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 минут.
Вы можете себе представить, что происходит, когда мне нужно обработать намного большая (хотя и более разреженная) матрица данных. Мне никогда не удавалось закончить sh вычисления.
Вопрос : как вы думаете, это можно сделать более эффективно?
Дело в том, что я не знаю, какие операции занимают больше времени (подмножество df, и в этом случае я бы удалил дубликаты этих добавляемых матричных данных, в этом случае я бы создал плоский вектор и затем преобразовал его в матрицу в конце? ...).
Спасибо!
РЕДАКТИРОВАТЬ после поста шахты
Как показывает minem, скорость этого вычисления сильно зависела от линейного Параметры регрессии были рассчитаны. Поэтому изменение этой части было самым важным делом.
Мои дальнейшие испытания показали, что: 1) было важно использовать sapply в сочетании с do.call(rbind, а не каким-либо плоским вектором, для хранения данных (Я до сих пор не уверен, почему - я мог бы сделать отдельный пост об этом); 2) в исходной матрице, над которой я работаю, гораздо более разреженной и с гораздо большим отношением nrows / ncolumns, чем в этом примере, используя информацию о векторе x, доступную в начале каждой итерации i, для уменьшение вектора y в начале каждой итерации j увеличивало скорость на несколько порядков, даже по сравнению с оригинальным сценарием 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)
Это, однако, медленнее, чем выше .