*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.lawre...@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.lawre...@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 <kem...@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
>>> kem...@ips.unibe.ch
>>>
>>> ______________________________________________
>>> R-help@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. ~
______________________________________________
R-help@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.
