Accélérer beaucoup de GLM dans R

Question

Tout d'abord, désolé pour le long post. J'ai préférable de donner son contexte pour obtenir de bonnes réponses (j'espère!). Il y a quelque temps, j'ai écrit une fonction R qui obtiendra toutes les interactions paires de variables dans un cadre de données. Cela a fonctionné bien à l'époque, mais maintenant un collègue voudrait que je fasse cela avec un jeu de données beaucoup plus grand. Ils ne savent pas combien de variables ils auront à la fin, mais ils devinent environ 2 500 - 3 000. Ma fonction ci-dessous est trop lente pour cela (4 minutes pour 100 variables). Au bas de cet article, j'ai inclus des horaires pour différents nombres de variables et de nombres totaux d'interactions. J'ai les résultats de Calling RProf () sur les 100 variables gérées de ma fonction, alors si quelqu'un veut jeter un coup d'œil à cela, faites le moi savoir. Je ne veux pas faire un super long plus longtemps que nécessaire.

Ce que j'aimerais savoir, c'est s'il y a quelque chose que je peux faire pour accélérer cette fonction. J'ai essayé de chercher directement à Glm.fit, mais aussi loin que j'ai compris, pour que cela soit utile le calcul des matrices de conception et de toutes ces autres choses que je ne comprenne franchement pas, doit être la même pour chaque modèle , ce qui n'est pas le cas pour mon analyse, bien que je me trompe peut-être à ce sujet.

Des idées sur la façon de faire courir plus vite seront grandement appréciées. Je prévois d'utiliser la parallélisation pour exécuter l'analyse à la fin, mais je ne sais pas combien de CPU je vais avoir accès à mais je dirais que cela ne sera pas plus de 8.

Merci d'avance, Cheers
Davy.

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

Dans le tableau ci-dessous est le système System.Time pour augmenter le nombre de variables passées à travers la fonction. n est le nombre et l'INTS est le nombre d'interactions paire-sage données par ce nombre de variables.

      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

Certains code pour reproduire une trame de données au cas où vous souhaitez examiner les horaires ou les résultats RPROF (). S'il vous plaît ne courez pas cela à moins que votre machine soit super rapide, ou si vous êtes prêt à attendre environ 15-20 minutes.

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

Donc, j'ai une sorte de résolution de cela (avec l'aide des listes de diffusion R) et je l'ai publié au cas où il est utile à quiconque.

Fondamentalement, où les SNPS ou les variables sont indépendants (c'est-à-dire non en LD, non corrélé), vous pouvez centrer chaque SNP / variable à son point méchant, comme suit:

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

Vous pouvez ensuite tester une corrélation entre phénotype et interaction comme étape de dépistage:

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

Ensuite, étudiez complètement en utilisant le GLM complet qui semble être corrélé.Le choix de coupure est, comme toujours, arbitraire.

D'autres suggestions devaient utiliser un test de score RAO, qui implique uniquement du modèle d'hypothèse nulle, ce qui a de moitié le temps de calcul de cette étape, mais je ne comprends pas vraiment comment cela fonctionne (encore! Plus de lecture requise.)

Quoi qu'il en soit, vous y allez.Peut-être être utile à quelqu'un un jour.