Genstat into R - Randomisation test
*Hits head*
Of course, the approach taken by your Genstat code of only shuffling
one variable is sufficient and faster:
n.obs = 100
ID=rnorm(n.obs)
CD=rnorm(n.obs)
obs.cor = cor(ID,CD)
num.permutations = 1e4
perm.cor = rep(NA,num.permutations)
start.time=proc.time()[1]
index = 1:n.obs
for(i in 1:num.permutations){
IDorder = sample(index)
perm.cor[i] = .Internal(cor(ID[IDorder], CD, 4, FALSE))
}
cat('Elapsed time:',start.time-proc.time(1))
sum(perm.cor>obs.cor)/num.permutations
On Wed, Apr 8, 2009 at 11:33 AM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:
"Assuming you have a data frame with columns ID & CD, this should do it" Oops, the code I sent doesn't assume this. It assumes that you have two vectors, ID & CD, as generated in the first 3 lines. On Wed, Apr 8, 2009 at 11:31 AM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:
Looks like that code implements a non-exhaustive variant of the
randomization test, sometimes called a permutation test. Assuming you
have a data frame with columns ID & CD, this should do it:
n.obs = 100
ID=rnorm(n.obs)
CD=rnorm(n.obs)
obs.cor = cor(ID,CD)
num.permutations = 1e4
perm.cor = rep(NA,num.permutations)
start.time=proc.time()[1]
index = 1:n.obs
for(i in 1:num.permutations){
? ? ? ?IDorder = sample(index)
? ? ? ?CDorder = sample(index)
? ? ? ?perm.cor[i] = .Internal(cor(ID[IDorder], CD[CDorder], 4, FALSE))
}
cat('Elapsed time:',start.time-proc.time(1))
sum(perm.cor>obs.cor)/num.permutations
On Wed, Apr 8, 2009 at 10:02 AM, Anne Kempel <kempel at ips.unibe.ch> wrote:
Hello everybody, I have a question. I would like to get a correlation between constitutive and induced plant defence which I messured on 30 plant species. So I have table with Species, Induced defence (ID), and constitutive defence (CD). Since Induced and constitutive defence are not independant (so called spurious correlation) I should do a randomisation test. I have a syntax of my supervisor in Genstat, but I would really like to try this in R. "data from trade-off.IDCD" list variate [nval=1000] slope calc ID1=ID graph ID; CD calc b=corr(ID; CD) calc slope$[1]=b "slope$[1] is the correlation before permutating the data" for i=2...1000 ? ?randomize ID1 ? ?calc b=corr(CD1; ID1) ? ?calc slope$[i]=b endfor hist slope describe slope quantile [proportion=!(0.0005,0.005, 0.025, 0.975, 0.995,0.9995)]slope print slope$[1] corr[p=corr] ID,CD DHISTOGRAM [WINDOW=1; ORIENTATION=vertical; KEY=0; BARCOVERING=1.0] slope; PEN=2 Does anybody have done something similar and has any idea how to make the randomisation part? I would be very grateful for any help!! Thanks in advance, Anne -- Anne Kempel Institute of Plant Sciences University of Bern Altenbergrain 21 CH-3103 Bern kempel at ips.unibe.ch
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
-- Mike Lawrence Graduate Student Department of Psychology Dalhousie University Looking to arrange a meeting? Check my public calendar: http://tinyurl.com/mikes-public-calendar ~ Certainty is folly... I think. ~
-- Mike Lawrence Graduate Student Department of Psychology Dalhousie University Looking to arrange a meeting? Check my public calendar: http://tinyurl.com/mikes-public-calendar ~ Certainty is folly... I think. ~
Mike Lawrence Graduate Student Department of Psychology Dalhousie University Looking to arrange a meeting? Check my public calendar: http://tinyurl.com/mikes-public-calendar ~ Certainty is folly... I think. ~