Beschleunigung Zählung von paarweise Beobachtungen in R

Question

Ich habe einen Datensatz in die eine Teilmenge von Messungen für jeden Eintrag zufällig fehlen:

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

Ich bin interessiert Korrelationen zwischen jeder Zeile in diesem Datensatz (dh aa Berechnung , ab, ac, ad, ...). Allerdings möchte ich Korrelationen ausschließen, bei denen es weniger als 5 paarweise nicht-NA-Beobachtungen gibt, indem deren Wert auf NA in der resultierenden Korrelationsmatrix gesetzt wird.

Zur Zeit mache ich das wie folgt:

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 

jedoch dieser Vorgang sehr langsam ist als der tatsächliche Datensatz enthält> 1.000 Einträge.

Gibt es eine Möglichkeit, Zellen basierend auf der Anzahl der nicht-NA-paarweisen Beobachtungen in einer vektorisierten Weise zu filtern, anstatt in verschachtelten Schleifen?

Answer 1

Sie können die Anzahl der nicht-NA paarweisen Beobachtungen mithilfe des Matrixansatzes zählen.

Verwenden wir diesen Datengenerierungscode. Ich habe Daten größer gemacht und mehr NAs hinzugefügt.

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 

Dann filtern Sie Code unter 85 Sekunden

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 

Meine Version erstellt eine Matrix mit den Werten 1 für Nicht-NAs in den Originaldaten. Dann berechne ich mit Hilfe der Matrixmultiplikation die Anzahl der paarweisen Nicht-NAs. Es lief in einem Bruchteil einer Sekunde.

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 

Einfache Kontrolle zeigt die Ergebnisse sind die gleichen:

all(filter == filter2) 
# TRUE