R help - Jun 2012 - Help with permutation function from Turner et al. 2010 (Ecology)

Hello,


I'm using R code that includes a residual permutation that was written as a
supplement to the paper:


Turner et al. 2010. A general hypothesis-testing framework for stable isotopes
ratios in ecological studies. Ecology 91:2227-2233.
The supplemental code is available at:
http://www.esapubs.org/archive/ecol/E091/157/suppl-1.htm


When I execute the function, no warnings are given (seems fine). But the
resulting output matrix is simply a matrix of repeating rows (i.e., it appears
as if the function is selecting the same "random" sample every time).
The function is as follows (please keep in mind that many of the internal
functions were written by Turner et al. as well and may not be familiar):


line.1<-c(0,traj.result(traj.size,traj.dir,traj.shape))


permute<-function(y,x2,g,o,p) { 
  p.table<-line.1
  yhat<-predict(x2)
  res<-resid(x2)
  line<-nrow(y)
  
  for(i in 1:p){
    line.rand<-sample(line,replace=FALSE)
    res.temp<-cbind(line.rand,res)
    z<-(order(line.rand))
    res.temp2<-as.matrix(res.temp[z,])
    res.p<-res.temp2[,-1]  # Rows of residuals are now randomized
    y.rand<-yhat+res.p
    
    lm.new<-lm(y.rand~o*g,x=T,model=T) 
    ls<-data.frame(expand.grid(o=levels(o),g=levels(g))) 
    ls[]<-lapply(ls, factor)
    means.r<-predict(lm.new,ls)
    
    traj.r<-arrayspecs(means.r,g,o)
    traj.size.r<-trajsize(traj.r)
    traj.dir.r<-trajorient(traj.r,ang.cor=T)
    traj.shape.r<-trajshape(traj.r)
    result<-traj.result(traj.size.r,traj.dir.r,traj.shape.r)
    p.table<-rbind(p.table,c(i,result))
  }
  p.table
}



But when I try to execute the function the rows in the output matrix repeat the
same numbers:

> test3<-permute(Y,lm.2f.red,group,obs,99)
> test3
        [,1]     [,2]      [,3]      [,4]
p.table    0 3.612322 0.9910104 0.3998703
             1 3.612322 0.9910104 0.3998703
             2 3.612322 0.9910104 0.3998703
             3 3.612322 0.9910104 0.3998703
             4 3.612322 0.9910104 0.3998703
             5 3.612322 0.9910104 0.3998703 
     ......


Can anyone see where there may be a problem with the permutation code? I'm
relatively new to functions, and especially to writing permutation functions, so
I'm having difficulty finding the problem.


Thanks,


Ashley





************************************************
Ashley Stasko
M.Sc. Candidate
Cooperative Freshwater Ecology Unit
Living With Lakes Centre
Laurentian University
Sudbury, ON
Tel: 705-675-1151 ext. 2209
Email: ax_stasko at laurentian.ca

On 09.06.2012 03:57, Ashley Stasko wrote:
> Hello,
>
>
> I'm using R code that includes a residual permutation that was written
as a supplement to the paper:
>
>
> Turner et al. 2010. A general hypothesis-testing framework for stable
isotopes ratios in ecological studies. Ecology 91:2227-2233.
> The supplemental code is available at:
http://www.esapubs.org/archive/ecol/E091/157/suppl-1.htm
>
>
> When I execute the function, no warnings are given (seems fine). But the
resulting output matrix is simply a matrix of repeating rows (i.e., it appears
as if the function is selecting the same "random" sample every time).
The function is as follows (please keep in mind that many of the internal
functions were written by Turner et al. as well and may not be familiar):
>
>
> line.1<-c(0,traj.result(traj.size,traj.dir,traj.shape))
>
>
> permute<-function(y,x2,g,o,p) {
>    p.table<-line.1
>    yhat<-predict(x2)
>    res<-resid(x2)
>    line<-nrow(y)
>
>    for(i in 1:p){
>      line.rand<-sample(line,replace=FALSE)
>      res.temp<-cbind(line.rand,res)
>      z<-(order(line.rand))
>      res.temp2<-as.matrix(res.temp[z,])
>      res.p<-res.temp2[,-1]  # Rows of residuals are now randomized
>      y.rand<-yhat+res.p
>
>      lm.new<-lm(y.rand~o*g,x=T,model=T)
>      ls<-data.frame(expand.grid(o=levels(o),g=levels(g)))
>      ls[]<-lapply(ls, factor)
>      means.r<-predict(lm.new,ls)
>
>      traj.r<-arrayspecs(means.r,g,o)
>      traj.size.r<-trajsize(traj.r)
>      traj.dir.r<-trajorient(traj.r,ang.cor=T)
>      traj.shape.r<-trajshape(traj.r)
>      result<-traj.result(traj.size.r,traj.dir.r,traj.shape.r)
>      p.table<-rbind(p.table,c(i,result))
>    }
>    p.table
> }
>
>
>
> But when I try to execute the function the rows in the output matrix repeat
the same numbers:
>
>
>> test3<-permute(Y,lm.2f.red,group,obs,99)
>> test3
>          [,1]     [,2]      [,3]      [,4]
> p.table    0 3.612322 0.9910104 0.3998703
>               1 3.612322 0.9910104 0.3998703
>               2 3.612322 0.9910104 0.3998703
>               3 3.612322 0.9910104 0.3998703
>               4 3.612322 0.9910104 0.3998703
>               5 3.612322 0.9910104 0.3998703
>       ......
>
>
> Can anyone see where there may be a problem with the permutation code?
I'm relatively new to functions, and especially to writing permutation
functions, so I'm having difficulty finding the problem.


If the results are suspicious (we cannot reproduce easily since you have 
not provided reproducible code as the posting guide asks you), why not 
ask the authors of the code? They hopefully know better what they aimed 
to get than readers of this list.

Best,
Uwe Ligges

>
> Thanks,
>
>
> Ashley
>
>
>
>
>
> ************************************************
> Ashley Stasko
> M.Sc. Candidate
> Cooperative Freshwater Ecology Unit
> Living With Lakes Centre
> Laurentian University
> Sudbury, ON
> Tel: 705-675-1151 ext. 2209
> Email: ax_stasko at laurentian.ca
>
> ______________________________________________
> 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.

Dear Ashley,

I am one of the co-authors for these scripts.  I just downloaded the example
files from the ESA archives and ran in naively on some data and it worked
fine.  My guess is that the problem lies in a portion of the script that you
did not provide.  If I had to be more certain with my guess, I bet your
'y'
data are the predicted values from a linear model instead of the raw data. 
Therefore, when the permutation procedure tries to shuffle residuals, there
are no residuals to shuffle.

If this is not the case, please email me and I can look at the rest of your
adapted script.  My new email address is michael (dot) collyer (at) wku
(dot) edu.  

By the way, one regret that we had after publishing these R scripts in the
ESA archives was to not provide a link for updates.  ESA is going to allow
us to make an amendment.  Look for this change in the near future.  We are
in the process of streamlining the original functions a bit and plan to make
these updates available.

Sincerely,
Michael Collyer
Western Kentucky University, Department of Biology




--
View this message in context:
http://r.789695.n4.nabble.com/Help-with-permutation-function-from-Turner-et-al-2010-Ecology-tp4632865p4632902.html
Sent from the R help mailing list archive at Nabble.com.