[R] Genstat into R - Randomisation test

Mike Lawrence Mike.Lawrence at dal.ca
Wed Apr 8 16:31:47 CEST 2009


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




More information about the R-help mailing list