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. ~