Acelerando muitos GLMs em R

Question

Primeiramente, desculpe pela longa postagem.Achei melhor contextualizar para obter boas respostas (espero!).Há algum tempo, escrevi uma função R que obterá todas as interações pares de variáveis ​​em um quadro de dados.Isso funcionou bem na época, mas agora um colega gostaria que eu fizesse isso com um conjunto de dados muito maior.Eles não sabem quantas variáveis ​​terão no final, mas estão adivinhando aproximadamente 2.500 - 3.000.Minha função abaixo é muito lenta para isso (4 minutos para 100 variáveis).No final deste post incluí alguns intervalos para vários números de variáveis ​​e números totais de interações.Eu tenho os resultados da chamada de Rprof() na execução de 100 variáveis ​​da minha função, então se alguém quiser dar uma olhada, me avise.Não quero fazer um superlongo por mais tempo do que o necessário.

O que eu gostaria de saber é se há algo que eu possa fazer para acelerar essa função.Tentei ir diretamente para glm.fit, mas pelo que entendi, para que isso seja útil o cálculo das matrizes de projeto e todas as outras coisas que francamente não entendo, precisam ser iguais para cada modelo , o que não é o caso da minha análise, embora talvez eu esteja errado sobre isso.

Qualquer idéia sobre como tornar isso mais rápido seria muito apreciada.Estou planejando usar paralelização para executar a análise no final, mas não sei a quantas CPUs terei acesso, mas diria que não serão mais do que 8.

Agradecemos antecipadamente
Davi.

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)
}

Na tabela abaixo estão os resultados do system.time para um número crescente de variáveis ​​​​passadas pela função.n é o número, e Ints, é o número de interações entre pares fornecidas por esse número de variáveis.

      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

Algum código para reproduzir um quadro de dados caso você queira ver os tempos ou os resultados do Rprof().Por favor, não execute isso a menos que sua máquina seja super rápida ou esteja preparado para esperar cerca de 15 a 20 minutos.

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)

Answer 1

Então eu meio que resolvi isso (com a ajuda das listas de discussão do R) e estou postando caso seja útil para alguém.

Basicamente, onde os SNPs ou variáveis ​​são independentes (ou seja,Não em LD, não correlacionado), você pode centralizar cada SNP/Variável em sua média, assim:

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

você pode então testar a correlação entre fenótipo e interação como uma etapa de triagem:

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

e então investigue completamente usando o glm completo qualquer um que pareça estar correlacionado.a escolha do ponto de corte é, como sempre, arbitrária.

Outras sugestões foram usar um teste de pontuação RAO, que envolve apenas ajustar o modelo de hipótese nula, reduzindo pela metade o tempo de cálculo para esta etapa, mas eu realmente não entendo como isso funciona (ainda!é necessária mais leitura.)

De qualquer forma, aí está.Talvez seja útil para alguém algum dia.