문제

첫 번째 OFF, 긴 게시물에 대해 유감스럽게 생각합니다. 좋은 답변을 얻기 위해 맥락을주는 것이 더 낫다는 것을 알았습니다 (나는 희망!). 얼마 전 데이터 프레임에서 변수의 모든 쌍의 상호 작용을 얻을 수있는 R 기능을 썼습니다. 이것은 당시에도 잘 작동했지만, 이제 동료가 훨씬 더 큰 데이터 세트로 이것을하고 싶습니다. 그들은 끝에 얼마나 많은 변수가 있을지 모르지만 약 2,500 ~ 3,000 ~ 3,000 명을 추측하고 있습니다. 아래의 내 기능은 너무 느리게 너무 느리게합니다 (100 변수의 경우 4 분). 이 게시물의 맨 아래에는 다양한 변수와 총 상호 작용 수에 대한 몇 가지 타이밍을 포함 시켰습니다. 나는 100 변수의 내 기능을 실행하는 RPROF ()를 호출 한 결과를 가지고 있으므로 누군가가 그것을 봐주고 싶다면 알려줍니다. 나는 그것이 필요하는 것보다 더 오래 슈퍼 오랜 길게 만들고 싶지 않습니다.

내가 알고 싶은 것은이 기능을 빠르게하기 위해 할 수있는 일이있는 경우에 알고 싶습니다. 나는 GLM.FIT에 직접가는 것을 시도했지만, 내가 이해 한 한, 디자인 행렬의 계산과 내가 솔직히 이해하지 못하는 다른 모든 것들의 계산을 유용하기 위해서는 각 모델에 대해 동일해야한다. 아마도 나는 이것에 대해 잘못되었지만 내 분석의 경우가 아닙니다.

이 실행을 더 빨리 만드는 방법에 대한 아이디어가 크게 감사 할 것입니다. 나는 결국 분석을 실행하기 위해 병렬화를 사용하고 있지만, 나는 얼마나 많은 CPU가 액세스 할 수 있는지 모르지만 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)
}
. 아래 표의

은 시스템을 통해 전달되는 변수 수를 증가시키는 데 필요한 시스템입니다. 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)
.

도움이 되었습니까?

해결책

그래서 나는 이것을 해결하고 (r 메일 링리스트에서 도움을 받아 도움을 받아) 누구에게도 유용한 경우를 게시하고 있습니다.

기본적으로 SNP 또는 변수가 독립적 인 경우 (즉, 상관 없음, 상관 없음) 각 SNP / 변수가 의미합니다. 

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

표현형과 상호 작용 간의 상관 관계를 심사 단계로 테스트 할 수 있습니다. 

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

을 누른 다음 전체 GLM을 사용하여 조사하여 상관 관계가있는 것처럼 보입니다.컷오프 선택은 언제든지 임의적입니다.

다른 제안은 RAO 스코어 테스트를 사용하는 것이 었습니다.이 단계는이 단계의 계산 시간을 절반으로하는 NULL 가설 모델만이 포함하지만,이 단계에서 어떻게 작동하는지 이해하지 못합니다 (아직도 더 많은 읽기가 필요합니다.)

어쨌든 당신은 가야합니다.언젠가 누군가에게 사용될 수도 있습니다.

라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top