Acelerando mucho de glms en r

Question

primero, lo siento por el poste largo. Pensé que es mejor dar un contexto para obtener buenas respuestas (¡espero!). Hace algún tiempo, escribí una función R que obtendrá todas las interacciones de las variables en un marco de datos. Esto funcionó bien en ese momento, pero ahora a un colega me gustaría que hiciera esto con un conjunto de datos mucho más grande. No saben cuántas variables van a tener al final, pero están adivinando aproximadamente 2.500 - 3.000. Mi función a continuación es demasiado lenta para esto (4 minutos para 100 variables). En la parte inferior de este post, he incluido algunos tiempos para varios números de variables y números totales de interacciones. Tengo los resultados de llamar a RPROF () en las 100 variables que se ejecutan de mi función, por lo que si alguien quiere echarle un vistazo, hágamelo saber. No quiero hacer un super largo por más tiempo de lo que debe ser.

Qué me gustaría saber es si hay algo que pueda hacer para acelerar esta función. Intenté mirar directamente a GLM.FIT, pero por lo que entendí, para que eso sea útil, el cálculo de las matrices de diseño y todo eso, otras cosas que francamente no entiendo, debe ser lo mismo para cada modelo. , que no es el caso de mi análisis, aunque quizás me equivoque sobre esto.

Cualquier idea sobre cómo hacer que esta carrera sea más rápida sería muy apreciada. Estoy planeando usar la paralelización para ejecutar el análisis al final, pero no sé cuántas CPU voy a tener acceso, pero diría que no será más de 8.

gracias de antemano, Saludos
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)
}

En la tabla a continuación se encuentra el sistema.Tiempo de resultados para aumentar el número de variables que se pasan a través de la función. n es el número, e INTS, es el número de interacciones en pareja dadas por ese número 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

Algún código para reproducir un marco de datos en caso de que desee ver los tiempos o los resultados de RPROF (). Por favor, no ejecute esto a menos que su máquina sea super rápida, o está preparado para esperar unos 15-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

Así que he resuelto esto (con la ayuda de las listas de correo R) y lo estoy publicando en caso de que sea útil para cualquier persona.

Básicamente, donde los SNP o las variables son independientes (es decir, no en LD, no se correlacionan), puede centrar cada SNP / variable en su significado, así:

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

Puede probar la correlación entre el fenotipo y la interacción como una etapa de detección:

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

y luego investigue completamente utilizando el GLM completo que parezca correlacionado.La opción de corte es, como siempre, arbitraria.

Otras sugerencias debían usar una prueba de puntuación RAO, que implica que solo se ajuste a la hipótesis nula. Modelo de la hipótesis. Esto reduce el tiempo de cálculo para este paso, pero realmente no entiendo cómo funciona esto (¡todavía! Más lectura requerida).

de todos modos ahí vayas.Tal vez sea de utilidad para alguien algún día.