Hi Matthew,

First - I fully support Greg Snow proposition.  Sampling is the way to go
here.

But besides that:
1) Try to avoid using data.frames as much as possible (use vectors and
matrixes instead - they are usually faster)
2) Since you are running on a loop, you can try running it in parallel (if
you have more then one core on your computer). I recently wrote how to do
this on Windows (here:
http://www.r-statistics.com/2010/04/parallel-multicore-processing-with-r-on-windows/)
, and there are other ways for doing it on other OS.

But again - sampling is probably going to be the only real solution...

Good luck,
Tal



----------------Contact
Details:-------------------------------------------------------
Contact me: Tal.Galili@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
----------------------------------------------------------------------------------------------




On Tue, Apr 27, 2010 at 1:09 AM, Matthew Finkbeiner <
matthew.finkbeiner@maccs.mq.edu.au> wrote:

> Yes, I suspect that I will end up using a sampling approach, but I'd like
> to use an exact test if it's at all feasible.
>
> Here are two samples of data from 3 subjects:
> Sample  Subj    C1      C2
> 44      1       0.0093  0.0077
> 44      2       0.0089  0.0069
> 44      3       0.051   0.0432
> 44      4       0.014   0.0147
> 44      5       0.0161  0.0117
> 45      1       0.0103  0.0086
> 45      2       0.0099  0.0078
> 45      3       0.0542  0.0458
> 45      4       0.0154  0.0163
> 45      5       0.0175  0.0129
>
>
> and then here is the script I've pieced together from things I've found on
> the web (sorry for not citing the snippets!).  any pointers on how to speed
> it up would be greatly appreciated.
>
> #----------------------------------
> # Utility function
> # that returns binary representation of 1:(2^n) X SubjN
> binary.v <-
> function(n)
> {
>  x <- 1:(2^n)
>  mx <- max(x)
>  digits <- floor(log2(mx))
>  ans <- 0:(digits-1); lx <- length(x)
>  x <- matrix(rep(x,rep(digits, lx)),ncol=lx)
>  x <- (x %/% 2^ans) %% 2
> }
>
> library(plyr)
>
>
> #first some global variables
> TotalSubjects <- 5
> TotalSamples <- 2
> StartSample <- 44
> EndSample <- ((StartSample + TotalSamples)-1)
> maxTs <- NULL
> obsTs <- NULL
>
>
>
>
> #create index array that drives the permuations for all samples
> ind <- binary.v(TotalSubjects)
>
> #transpose ind so that the first 2^N items correspond to S1,
> #the second 2^N correspond to S2 and so on...
> transind <- t(ind)
>
> #get data file that is organized first by sample then by subj (e.g. sample1
> subject1
> # sample1 subject 2 ... sample 1 subject N)
> #sampledatafile <- file.choose()
>
> samples <- read.table(sampledatafile, header=T)
>
> #this is the progress bar
> pb <- txtProgressBar(min = StartSample, max = EndSample, style = 3)
> setTxtProgressBar(pb, 1)
>
> start.t <- proc.time()
>
> #begin loop that analyzes data sample by sample
> for (s in StartSample:EndSample) {
>
>    S <- samples[samples$Sample==s,] #pick up data for current sample
>
>    #reproduce data frame rows once for each permutation to be done
>    expanddata <- S[rep(1:nrow(S), each = 2^TotalSubjects),]
>
>
>    #create new array to hold the flipped (permuted) data
>    permdata = expanddata
>
>    #permute the data
>    permdata[transind==1,3] <- expanddata[transind==1,4] #Cnd1 <- Cnd2
>    permdata[transind==1,4] <- expanddata[transind==1,3] #Cnd2 <- Cnd1
>
>    #create permutation # as a factor in dataframe
>    PermN <- rep(rep(1:2^TotalSubjects, TotalSubjects),2)
>
>    #create Sample# as a factor
>    Sample <- rep(permdata[,1],2) #Sample# is in the 1st Column
>
>    #create subject IDs as a factor
>    Subj <- rep(permdata[,2],2) #Subject ID is in the 2nd Column
>
>    #stack the permutated data
>    StackedPermData <- stack(permdata[,3:4])
>
>    #bind all the factors together
>    StackedPermData <- as.data.frame(cbind(Sample, Subj, PermN,
> StackedPermData))
>
>
>    #sort by perm
>    sortedstack <-
> as.data.frame(StackedPermData[order(StackedPermData$PermN,
>                                StackedPermData$Sample),])
>
>
>    #clear up some memory
>    rm(expanddata, permdata, StackedPermData)
>
>    #pull out data 1 perm at a time
>    res<-ddply(sortedstack, c("Sample", "PermN"), function(.data){
>
>        # Type combinations by Class
>        combs<-t(combn(sort(unique(.data[,5])),2))
>
>        # Applying the t-test for them
>        aaply(combs,1, function(.r){
>        x1<-.data[.data[,5]==.r[1],4] # select first column
>        x2<-.data[.data[,5]==.r[2],4] # select first column
>
>        tvalue <- t.test(x1,x2, paired = T)
>
>        res <- c(tvalue$statistic,tvalue$parameter,tvalue$p.value)
>        names(res) <- c('stat','df','pvalue')
>        res
>        }
>        )
>    }
>    )
>
>     # update progress bar
>     setTxtProgressBar(pb, s)
>
>     #get max T vals
>     maxTs <- c(maxTs, tapply (res$stat, list (res$Sample), max))
>
>     #get observed T vals
>     obsTs <- c(obsTs, res$stat[length(res$stat)])
>
>     #here we need to save res to a binary file
>
>
> }
> #close out the progress bar
> close(pb)
>
> end.t <- proc.time() - start.t
> print(end.t)
>
> #get cutoffs
> #these are the 2-tailed t-vals that maintain experimentwise error at the
> 0.05 level
> lowerT <- quantile(maxTs, .025)
> upperT <- quantile(maxTs, .975)
>
>
>
>
>
>
>
>
>
>
>
>
> On 4/27/2010 6:53 AM, Greg Snow wrote:
>
>> The usual way to speed up permutation testing is to sample from the
>> set of possible permutations rather than looking at all possible
>> ones.
>>
>> If you show some code then we may be able to find some inefficiencies
>> for you, but there is not general solution, poorly written uses of
>> apply will be slower than well written for loops.  In some cases
>> rewriting critical pieces in C or fortran will help quite a bit, but
>> we need to see what you are already doing to know if that will help
>> or not.
>>
>>
> ______________________________________________
> 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.
>

	[[alternative HTML version deleted]]

