2017-02-22 19 views
1

Я делаю цикл для вычисления мощности теста Харди-Вайнберга для варианта под названием «snp1» (с командой HWPower из пакета HardyWeinberg). Эта команда требует ввода n (размер выборки) и pA (малая частота аллелей). Я должен делать расчет много раз с большим количеством ns и pAs, потому что они представляют разные образцы населения, поэтому я сделал первые два вручную, но теперь я хочу сделать цикл for для всех остальных. Я начал с простого цикла с первыми двумя, поэтому я могу легко проверить, что результаты в порядке (и, следовательно, цикл работает нормально). Но я столкнулся с проблемой при сравнении результатов обоих вычислений, что заставляет меня думать, что мой код не совсем точен.Несоответствие между результатом вложенного цикла цикла и независимой функцией

install.packages("HardyWeinberg") 
library(HardyWeinberg) 

snp1n=c(661,503) 
snp1pA=c(0.006051,0.174) 
HWpowersnp1<-numeric(2) 

for(i in seq_along(snp1n)) { 
    for(j in seq_along(snp1pA)) { 
    HWpowersnp1[i]<-HWPower(n=snp1n[i],pA=snp1pA[j]) 
    } 
} 
HWpowersnp1 

Это дает мне следующий вектор:

HWpowersnp1 
[1] 0.04109278 0.04253145 

Но когда я рассчитать каждый из них с помощью функции самостоятельно, я получаю:

HWPower(n = 661,pA = 0.006051) 
[1] 0.02107572 
HWPower(n = 503, nA = 175) 
[1] 0.04253145 

Я не знаю, где проблема в том, что вызывает несогласованность. Это странно, потому что его только в первом результате, а не во втором (вторая рассчитанная мощность в порядке, это дает мне тот же результат, но первый не делает).

+0

Вы определяете 'HWpowersnp1 <-numeric (2)' как числовой массив размера два, так что вы должны получить в конце цикла два значения. При втором вызове возвращается скалярное значение, как и должно быть. –

ответ

0

Ваша двойная петля неверна. Если вы запустите это

for(i in seq_along(snp1n)) { 
    for(j in seq_along(snp1pA)) { 
    print(paste(i,j)) 
    HWpowersnp1[i]<-HWPower(n=snp1n[i],pA=snp1pA[j]) 
    } 
} 

Вы увидите, что ваша функция работает 4 раза, а не две, как вы ожидаете. Вы хотите итерации по snp1n и snp1pA одновременно. Поэтому вы должны использовать mapply или Map. Попробуйте это вместо

HWpowersnp1 <- Map(HWPower, n=snp1n, pA=snp1pA) 
+0

Спасибо! Теперь я вижу, что я неправильно понял использование mapply, и именно поэтому я не использовал его в первую очередь. Большое спасибо! Решено. – melunuge92