Accélérer beaucoup de GLM dans R

Question

Tout d'abord, désolé pour le long poste. Il est préférable de donner un 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 par paires de variables dans une trame de données. Cela a bien fonctionné à l'époque, mais maintenant un collègue aimerait que je le fasse avec un ensemble de données beaucoup plus grand. Ils ne savent pas combien de variables ils vont avoir au bout du compte, mais ils devinent environ 2 500 à 3 000. Ma fonction ci-dessous est beaucoup trop lente pour cela (4 minutes pour 100 variables). Au bas de ce post, j'ai inclus des timings pour différents nombres de variables et nombre total d'interactions. J'ai les résultats de l'appel de Rprof() sur les 100 variables de ma fonction, donc si quelqu'un veut y jeter un coup d'oeil, 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 regarder directement glm.fit, mais pour autant que je comprenne, pour que cela soit utile, le calcul des matrices de conception et tout ce que je ne comprends pas franchement, doit être le 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.

Toutes les idées pour accélérer le fonctionnement de ce système seraient grandement appréciées. J'ai l'intention d'utiliser la parallélisation pour faire l'analyse à la fin, mais je ne sais pas combien de processeurs j'aurai accès mais je dirais que ce 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, vous trouverez les résultats system.time pour l'augmentation du nombre de variables transmises à travers la fonction. n est le nombre, et Ints, est le nombre d'interactions par paires donné 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 

Un code pour reproduire une trame de données dans le cas où vous voulez regarder les timings ou les résultats Rprof(). S'il vous plaît ne pas exécuter ceci à moins que votre machine soit super rapide, ou votre 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

J'ai donc résolu cela (avec l'aide des listes de diffusion R) et je l'affiche au cas où cela serait utile à tout le monde.

Fondamentalement, où les SNPs ou variables sont indépendantes (non en LD, décorrélée) vous pouvez centrer chaque SNP/variable à cela signifie comme ceci:

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

vous pouvez tester la corrélation entre le phénotype et l'interaction en tant qu'étape de criblage:

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

puis d'étudier complètement en utilisant le glm complet qui semble être corrélé. Le choix du cut-off est, comme toujours, arbitraire. D'autres suggestions devaient utiliser un test de score RAO, qui consiste à ajuster le modèle d'hypothèse nulle en réduisant de moitié le temps de calcul pour cette étape, mais je ne comprends pas vraiment comment cela fonctionne (encore plus de lecture requise.)

En tout cas, c'est parti. Peut-être être utile à quelqu'un un jour.