2015-05-28 7 views
2

У меня есть dataframe со следующими размерами:Быстрое вычисление ANOVA в R

dim(b) 
[1] 974 433685 

Столбцы представляют переменные, которые я хочу запустить ANOVAs на (то есть, я хочу запустить 433,685 ANOVAs). Размер выборки - 974. Последний столбец - это переменная «группа».

Я придумал 3 разных метода, но все они слишком медленны из-за количества тестов.

Во-первых, давайте выдадим небольшую практику набора данных, чтобы играть с:

dat = as.data.frame(matrix(runif(10000*500), ncol = 10000, nrow = 500)) 
dat$group = rep(letters[1:10], 5000) 

Метод 1 (основано на 'sapply'):

system.time(sapply(dat[,-length(dat)], function(x) aov(x~group, data=dat))) 

    user system elapsed 
143.76 0.33 151.79 

Методы 2 (на основе 'mclapply' из «параллельный» пакет):

library(parallel) 
options(mc.cores=3) 
system.time(mclapply(dat[,-length(dat)], function(x) aov(x~group, data=dat))) 

    user system elapsed 
141.76 0.21 142.58 

Методы 3 (на основе 'cbind'-ing the LHS):

formula = as.formula(paste0("cbind(", paste(names(dat)[-length(dat)],collapse=","), ")~group")) 
system.time(aov(formula, data=dat)) 

    user system elapsed 
    10.00 0.22 10.25 

В практике набора данных, метод 3 является явным победителем. Однако, когда я делаю это на моих фактических данных, вычисления на только 10 (из 433685) столбцов, используя метод 3 принимает это долго:

user system elapsed 
119.028 5.430 124.414 

Не знаю, почему он занимает существенно больше времени на моих фактических данных. У меня есть доступ к кластеру Linux с 16 ядрами и 72 ГБ оперативной памяти.

Есть ли способ вычислить это быстрее?

+3

'433,685 ANOVAs'? Для чего вы хотите это сделать? Для вашей проблемы должен быть более эффективный статистический подход. – Roland

+1

Это p-значение, исправленное Bonferonni 1.152911e-07 –

+0

Роланд, я пытаюсь количественно определить «пакетный эффект» (по моей «групповой» переменной), что является вариацией от нежелательного технического артефакта. У меня есть 433K индивидуальных зондов для измерения метилирования ДНК в определенных геномных локусах. Моя мысль заключалась в том, чтобы суммировать индивидуальную F-статистику ANOVA с каждого зонда. Затем я сравнил это число с аналогичными наборами данных, сгенерированными из разных конвейеров предварительной обработки, чтобы найти тот, который удалил больше «пакетного эффекта». –

ответ

0

Для одновременной установки многих общих линейных моделей (таких как ANOVA) с использованием того же design matrix, Bioconductor/R limma package обеспечивает очень быструю функцию lmFit(). Это, как можно подогнать ANOVA модели, используя limma:

library(limma) 

# generate some data 
# (same dimensions as in your question) 
nrows <- 1e4 
ncols <- 5e2 
nlevels <- 10 
dat <- matrix(
    runif(nrows * ncols), 
    nrow = nrows, 
    ncol = ncols 
) 
group <- factor(rep(
    letters[1:nlevels], 
    ncols/nlevels 
)) 

# construct the design matrix 
# (same as implicitly used in your question) 
dmat <- model.matrix(~ group) 
# fit the ANOVA model 
fit <- lmFit(dat, dmat) 

На моем ноутбуке он закончил в 0,4 - 0,45 секунды, по данным тех же размеров, что и данные в вашем вопросе.