Ускорение множества GLM в R - PullRequest
Новая
       13

Ускорение множества GLM в R

1 голос
/ 12 марта 2012

Прежде всего, извините за длинный пост. Понял, что лучше дать контекст, чтобы получить хорошие ответы (надеюсь!). Некоторое время назад я написал функцию R, которая будет получать все попарные взаимодействия переменных в кадре данных. В то время это работало нормально, но теперь коллега хотел бы, чтобы я сделал это с гораздо большим набором данных. Они не знают, сколько переменных у них будет в итоге, но они предполагают, что примерно 2500-3000. Моя функция ниже слишком медленная для этого (4 минуты для 100 переменных). В нижней части этого поста я включил некоторые моменты времени для различного количества переменных и общего количества взаимодействий. У меня есть результаты вызова Rprof () для запуска 100 функций моей функции, поэтому, если кто-то захочет взглянуть на него, дайте мне знать. Я не хочу делать супер длиннее, чем нужно.

Что я хотел бы знать, так это то, что я могу сделать что-нибудь, чтобы ускорить эту функцию. Я пытался смотреть прямо на glm.fit, но, насколько я понял, для того, чтобы это было полезно, вычисление матриц проектирования и всех других вещей, которые я, честно говоря, не понимаю, должно быть одинаковым для каждой модели. , что не так для моего анализа, хотя, возможно, я ошибаюсь по этому поводу.

Буду очень признателен за любые идеи о том, как сделать это быстрее. Я планирую использовать распараллеливание для запуска анализа в конце, но я не знаю, сколько процессоров у меня будет, но я бы сказал, что это будет не более 8.

Заранее спасибо, Приветствия
Дэви. 

getInteractions2 = function(data, fSNPcol, ccCol)
{
#fSNPcol is the number of the column that contains the first SNP
#ccCol is the number of the column that contains the outcome variable
  require(lmtest)
  a = data.frame()
  snps =  names(data)[-1:-(fSNPcol-1)]
  names(data)[ccCol] = "PHENOTYPE"
  terms = as.data.frame(t(combn(snps,2)))
  attach(data)

  fit1 = c()
  fit2 = c()
  pval  = c()

  for(i in 1:length(terms$V1))
  {
    fit1 = glm(PHENOTYPE~get(as.character(terms$V1[i]))+get(as.character(terms$V2[i])),family="binomial")
    fit2 = glm(PHENOTYPE~get(as.character(terms$V1[i]))+get(as.character(terms$V2[i]))+I(get(as.character(terms$V1[i]))*get(as.character(terms$V2[i]))),family="binomial")
    a  = lrtest(fit1, fit2)
    pval = c(pval, a[2,"Pr(>Chisq)"])
  }

  detach(data)
  results = cbind(terms,pval) 
  return(results)
}

В таблице ниже приведены результаты system.time для увеличения числа переменных, передаваемых через функцию. n - число, а Ints - количество парных взаимодействий, заданных этим числом переменных. 

      n   Ints     user.self sys.self elapsed
time  10   45      1.20     0.00    1.30
time  15  105      3.40     0.00    3.43
time  20  190      6.62     0.00    6.85
... 
time  90 4005    178.04     0.07  195.84
time  95 4465    199.97     0.13  218.30
time 100 4950    221.15     0.08  242.18

Некоторый код для воспроизведения фрейма данных на случай, если вы захотите посмотреть тайминги или результаты Rprof (). Пожалуйста, не запускайте это, если ваша машина не работает очень быстро или вы не готовы ждать 15-20 минут. 

df = data.frame(paste("sid",1:2000,sep=""),rbinom(2000,1,.5))
gtypes = matrix(nrow=2000, ncol=3000)
gtypes = apply(gtypes,2,function(x){x=sample(0:2, 2000, replace=T);x})
snps = paste("rs", 1000:3999,sep="")
df = cbind(df,gtypes)
names(df) = c("sid", "status", snps)

times = c()
for(i in seq(10,100, by=5)){
  if(i==100){Rprof()}
  time = system.time((pvals = getInteractions2(df[,1:i], 3, 2)))
  print(time)
  times = rbind(times, time)
  if(i==100){Rprof(NULL)}
}
numI = function(n){return(((n^2)-n)/2)}
timings = cbind(seq(10,100,by=5), sapply(seq(10,100,by=5), numI),times)

1 Ответ

0 голосов
/ 13 марта 2012

Итак, я вроде как решил это (с помощью списков рассылки R) и публикую это на тот случай, если это пригодится кому-либо.

В основном, где SNP или переменные независимы (т. Е. Нев LD, не коррелируется) вы можете центрировать каждое SNP / переменную по среднему значению примерно так: 

rs1cent <- rs1-mean(rs1)
rs2cent <- rs2 -mean(rs2)

затем вы можете проверить корреляцию между фенотипом и взаимодействием в качестве шага проверки: 

rs12interaction <- rs1cent*rs2cent
cor(PHENOTYPE, rs12interaction)

, а затем полностью исследуйте, используя полный glm, любой, который кажется коррелированным.Выбор отсечения, как всегда, произвольный.

Другие предложения заключались в использовании критерия оценки RAO, который включает в себя только подгонку модели нулевой гипотезы, что вдвое сокращает время вычисления для этого шага, но я не оченьпонять, как это работает (пока! требуется больше чтения.)

В любом случае, вы идете.Может быть, когда-нибудь кому-нибудь пригодится.

...