lem-usp / EvolQG

Tools for Evolutionary Quantitative Genetics
http://cran.r-project.org/web/packages/evolqg/
Other
10 stars 8 forks source link

RandomSkewers #22

Closed alexhubbe closed 11 years ago

alexhubbe commented 11 years ago

Caro,

estava comparando a função RandomSkewers com a minha função de RS (que foi adaptada do Pato) e sistematicamente a minha função gera valores maiores (pouco maiores, mas mesmo assim consistentemente maiores) e não consegui identificar o pq da diferença.

Realizei 1000 comparações entre duas VCV usando os dois scripts e o meu, como já dito gerou valores maiores. Tentei anexar uma figura pra mostrar isso. Nela, os círculos são valores de RS usando a função do GITHUB e os Xs são usando a minha função zzzz

Será que vc pode dar uma olhada na minha função e opinar qual o motivo da diferença? Acho que na prática a diferença é irrisória, mas existe...

A minha função é:

random.skewers <- function (vcv1, vcv2, nsk = 10000) { size <- dim (vcv1) [1] isovec <- normalize (rep (1, times = size)) rvec <- array (rnorm (nsk * size, mean = 0, sd = 1), c(size, nsk)) rvec <- apply (rvec, 2, normalize) dist <- abs (isovec %% rvec) dz1 <- apply (vcv1 %% rvec, 2, normalize) dz2 <- apply (vcv2 %% rvec, 2, normalize) real <- abs (apply (dz1 \ dz2, 2, sum)) ac <- mean (real) stdev <- sd (real) prob <- sum (real < dist) / nsk

mag1=c(mean((apply((vcv1 %*% rvec)^2,2,sum))^(1/2)),sd((apply((vcv1 %*% rvec)^2,2,sum))^(1/2)))
mag2=c(mean((apply((vcv2 %*% rvec)^2,2,sum))^(1/2)),sd((apply((vcv2 %*% rvec)^2,2,sum))^(1/2)))

output <- c(ac, prob, stdev,mag1,mag2)

names(output) <- c("AC","P","SD","magnitude_vcv_1: média","magnitude_vcv_1: sd","magnitude_vcv_2: média","magnitude_vcv_1: sd")

return(output)

}

VALEU!!!!

diogro commented 11 years ago

A única diferença é que a Random Skewers do repositório não toma o valor absoluto das correlações entre os vetores resposta. Vc pode ver uma discussão sobre isso na errada do artigo sobre o SRD: http://link.springer.com/article/10.1007%2Fs11692-012-9192-5?LI=true

sua versão:

real <- abs (apply (dz1 * dz2, 2, sum))

Repo:

real <- (apply (dz1 * dz2, 2, sum))